Thunderforest/OpenStreetMap中纬度与经度缩放不一致问题求助

问题分析 #

你当前的问题出在纬度和经度方向的度数/像素比例不匹配：地球是椭球体，相同像素宽度对应的经度差，和相同像素高度对应的纬度差，在非赤道区域是不一样的。你直接复用了经度的度数范围给纬度，导致地理上的圆形区域在地图上被拉伸成椭圆。

根据OSM的说明，赤道处256瓦片的米/像素值，在其他纬度需要乘以该纬度的余弦值——本质是说，纬度方向的单位像素对应的实际距离（度数）是经度方向的1/cos(纬度)倍（因为经度线在两极汇聚，纬度越高，相同经度差的实际距离越短）。

具体修改方案 #

需要调整两个关键部分：

1. 计算纬度方向的地理范围
   把纬度的padding值除以当前发射点纬度的余弦值，这样相同像素高度对应的纬度差会适配实际地理比例：

   import math
   
   # 替换原padding_lat计算代码
   lat_rad = math.radians(launch_latitude)
   padding_lat = map_width_deg / (2 * math.cos(lat_rad))
   

2. 修正matplotlib的aspect设置
   原代码中ax.imshow的aspect='equal'会强制坐标轴比例相等，但实际地图的地理范围已经调整，应该让matplotlib自动适配：

   # 替换原aspect='equal'参数
   ax.imshow(map_img, zorder=0, extent=BBox, aspect='auto')
   

修改后的完整代码（关键部分标注） #

# data
launch_latitude = np.degrees(TOTAL_DATA[0]['TYPE_LATITUDE'][0][0])
launch_longitude = np.degrees(TOTAL_DATA[0]['TYPE_LONGITUDE'][0][0])

# variabels
zoom_level = 14 # 0 <= int <= 20
map_width_pixel = 1100 # max: 2560
save_fig = True
plot_name = 'name'
plot_title = 'NAME'

img_name = f'z{zoom_level}_w{map_width_pixel}_{launch_latitude}N_{launch_longitude}E.png'
img_path = f'../env/images/{img_name}'
tile_width_at_zoom_level = {0:360,1:180,2:90,3:45,4:22.5,5:11.25,6:5.625,7:2.813,8:1.406,9:0.703,10:0.352,11:0.176,12:0.088,13:0.044,14:0.022,15:0.011,16:0.005,17:0.003,18:0.001,19:0.0005,20:0.00025}
map_width_deg = tile_width_at_zoom_level[zoom_level]/256 * map_width_pixel
padding_long = map_width_deg/2

# --- 修改部分开始 ---
import math
lat_rad = math.radians(launch_latitude)
# 纬度方向padding适配纬度余弦值
padding_lat = map_width_deg / (2 * math.cos(lat_rad))
# --- 修改部分结束 ---

BBox = [launch_longitude - padding_long, launch_longitude + padding_long, launch_latitude - padding_lat, launch_latitude + padding_lat]

# download map, if image doesn't allready exist.
if not os.path.isfile(img_path):
    # get URL
    zoom = zoom_level
    style = 'Landscape'
    apikey = '--- API KEY ---'
    lon, lat = launch_longitude, launch_latitude
    width, height = map_width_pixel, map_width_pixel
    URL = f'https://tile.thunderforest.com/static/{style}/{lon},{lat},{zoom}/{width}x{height}.png?apikey={apikey}'

    # save img
    r = requests.get(URL, stream=True)
    r.raw.decode_content = True
    with open(img_path,'wb') as f:
        shutil.copyfileobj(r.raw, f)

# plot
fig, ax = plt.subplots(figsize = (18,7))
map_img = plt.imread(img_path)
attribution = 'Maps © www.thunderforest.com, Data © www.osm.org/copyright'

# --- 修改部分开始 ---
# 取消强制等比例，自动适配地理范围
ax.imshow(map_img, zorder=0, extent = BBox, aspect='auto')
# --- 修改部分结束 ---

ax.scatter(launch_longitude, launch_latitude, zorder=4, marker='^', color='black', label='Launch')

for zorder, raw_data, label in zip([1,2,3], TOTAL_DATA, ['Nominal Flight', 'No Main', 'No Parachute']):
    land_latitude = [np.degrees(flight[-1]) for flight in raw_data['TYPE_LATITUDE']]
    land_longitude = [np.degrees(flight[-1]) for flight in raw_data['TYPE_LONGITUDE']]
    ax.scatter(land_longitude, land_latitude, zorder=zorder, marker='x', label=label)

ax.set_title(f'{plot_title} Dropzone Analysis')
ax.set_xlim(BBox[0],BBox[1])
ax.set_ylim(BBox[2],BBox[3])
ax.set_xlabel('Longitude [WGS84]')
ax.set_ylabel('Latitude [WGS84]')
ax.text(x=1, y=0.01, s=attribution, horizontalalignment='right', verticalalignment='bottom', rotation='vertical', fontsize=8, color='gray', alpha=0.9, transform=ax.transAxes)
ax.legend()
if save_fig:
    plt.savefig(f'../results/{plot_name}_dropzone_analysis.pdf')
plt.show()

原理说明 #

• 地球经度线在所有纬度上的间距相等，但纬度线间距随纬度升高而缩小，相同像素高度对应的纬度差要比经度差大（除以cos(纬度)），这样才能让地理上的圆形区域在地图上保持圆形。
• aspect='auto'会让matplotlib根据地图图片的比例和地理范围自动调整坐标轴比例，避免强制等比例导致的拉伸。

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