Python PIL经纬度转地图像素点绘制图标问题求助

问题背景 #

将经纬度坐标转换为像素坐标，使用PIL库在地图图片上粘贴自定义图标，采用地图中心点相对绘制法，但图标与预期标记偏差过大。

现有条件 #

• 地图图片：宽度2800，高度2600
• 缩放级别：16.045383428216734
• 地图中心点：[149.7985724, -35.4420469]
• 待绘制点：[(149.79889, -35.44060), (149.79772, -35.44349)]

当前问题 #

生成的图片中橙色框图标与预期圆形标记偏差过大（已知滚动条会导致少量偏移，但当前偏差远超预期）。

原始代码 #

import math
from math import *
from math import pi, atan, sinh, log, tan, cos

import pandas as pd
import geopandas as gpd

from PIL import Image, ImageOps, ImageChops, ImageDraw

def getDistance(y,x,lat,lng):
    p1 = (float(lat), float(lng))
    p2 = (float(y),float(x))
    distance = round(geodesic(p1, p2).meters,0)
    return distance

mapboxZoom = 16.045383428216734
snagDpi = 4
mapWidth = 700
mapHeight = 650
centroid = [149.7985724, -35.4420469]

schDfDict = {'lng': {1: 149.79889, 2: 149.7977248}, 'lat': {1: -35.4406, 2: -35.4434938}}
schDf = pd.DataFrame(schDfDict)
schDf = gpd.GeoDataFrame(schDf, geometry=[Point(xy) for xy in zip(schDf['lng'].astype(float), schDf['lat'].astype(float))])

im = Image.open(io.BytesIO(base64.b64decode(img)))

metersPerPx = 156543.03392 * cos(centroid[1] * pi / 180) / 2.0 ** mapboxZoom
mapWidthKm = metersPerPx * mapWidth
mapHeightKm = metersPerPx * mapHeight

## multiply distances by snagDpi
snagMetersPerPx = metersPerPx / snagDpi

centroidPx = ((mapWidth * snagDpi) / 2, (mapHeight * snagDpi) / 2)

foreground = Image.open("C:\\Users\\ronan\\Downloads\\testIcon.png")

for index, row in schDf.reset_index().iterrows():

    prop = schDf.iloc[[index]]
    propX = prop.geometry.x.values[0]
    propY = prop.geometry.y.values[0]

    angle = degrees(atan((propX - centroid[0]) / (propY - centroid[1])))

    mDist = getDistance(centroid[1],centroid[0],propY,propX)

    pxDist = mDist / snagMetersPerPx

    ## radian
    alpha = angle * pi / 180

    xx = centroidPx[0] + (pxDist * cos(alpha))
    yy = centroidPx[1] + (pxDist * sin(alpha))

    im.paste(foreground, (round(xx),round(yy)), foreground)

问题排查与修复方案 #

1. 角度计算逻辑错误 #

当前用atan((propX - centroid[0]) / (propY - centroid[1]))计算角度，无法区分四个象限的方向，且未处理分母为0的极端情况，导致方向判断错误。

修复： 替换为atan2函数，自动处理象限问题：

# 替换原角度计算代码
dx = propX - centroid[0]
dy = propY - centroid[1]
angle = degrees(atan2(dy, dx))

2. 像素坐标系方向不匹配 #

PIL的图片坐标系Y轴向下递增，而地理纬度Y轴向上递增（纬度越大越靠北），导致Y方向偏移完全反向。

修复： 计算Y轴像素坐标时取反偏移量：

yy = centroidPx[1] - (pxDist * sin(alpha))

3. 地图尺寸变量与实际不符 #

代码中mapWidth=700、mapHeight=650与实际地图尺寸2800x2600不匹配，即使通过snagDpi缩放，也会引入不必要的计算误差。

修复： 直接使用实际地图尺寸：

mapWidth = 2800
mapHeight = 2600
snagDpi = 1  # 无需额外缩放

4. 缺失关键依赖导入 #

geodesic函数未导入，代码运行时会报错，同时Point类也未明确导入。

修复： 添加必要导入：

from geopy.distance import geodesic
from shapely.geometry import Point
import io
import base64

5. 图标锚点未对齐中心 #

im.paste默认以图标左上角为定位点，但预期标记通常是图标中心对齐目标位置，导致视觉偏移。

修复： 计算图标中心点偏移：

# 获取图标尺寸
icon_w, icon_h = foreground.size
# 调整坐标为图标中心位置
paste_x = round(xx - icon_w / 2)
paste_y = round(yy - icon_h / 2)

修复后的完整代码 #

import math
from math import pi, degrees, atan2, cos
from geopy.distance import geodesic
import pandas as pd
import geopandas as gpd
from shapely.geometry import Point
import io
import base64
from PIL import Image

def getDistance(y,x,lat,lng):
    p1 = (float(lat), float(lng))
    p2 = (float(y),float(x))
    distance = round(geodesic(p1, p2).meters, 0)
    return distance

mapboxZoom = 16.045383428216734
mapWidth = 2800
mapHeight = 2600
centroid = [149.7985724, -35.4420469]

schDfDict = {'lng': {1: 149.79889, 2: 149.7977248}, 'lat': {1: -35.4406, 2: -35.4434938}}
schDf = pd.DataFrame(schDfDict)
schDf = gpd.GeoDataFrame(schDf, geometry=[Point(xy) for xy in zip(schDf['lng'].astype(float), schDf['lat'].astype(float))])

im = Image.open(io.BytesIO(base64.b64decode(img)))

metersPerPx = 156543.03392 * cos(centroid[1] * pi / 180) / (2.0 ** mapboxZoom)
centroidPx = (mapWidth / 2, mapHeight / 2)

foreground = Image.open("C:\\Users\\ronan\\Downloads\\testIcon.png")
icon_w, icon_h = foreground.size

for _, row in schDf.iterrows():
    propX = row.geometry.x
    propY = row.geometry.y

    dx = propX - centroid[0]
    dy = propY - centroid[1]
    angle = degrees(atan2(dy, dx))
    alpha = math.radians(angle)

    mDist = getDistance(centroid[1], centroid[0], propY, propX)
    pxDist = mDist / metersPerPx

    xx = centroidPx[0] + (pxDist * cos(alpha))
    yy = centroidPx[1] - (pxDist * sin(alpha))

    paste_x = round(xx - icon_w / 2)
    paste_y = round(yy - icon_h / 2)

    im.paste(foreground, (paste_x, paste_y), foreground)

im.save("result_map.png")

内容的提问来源于stack exchange，提问作者ronan