（geopandas）如何在墨卡托投影中正确输出经度/纬度比例？答案

【问题标题】：(geopandas) How to output longitude/latitude scale correctly in Mercator projetion?（geopandas）如何在墨卡托投影中正确输出经度/纬度比例？
【发布时间】：2021-04-28 12:19:38
【问题描述】：

在使用 geopandas (Python3) 创建的墨卡托地图中，坐标轴旁边的经度和纬度比例不正确。

我用这样的地理熊猫制作了一张墨卡托投影地图。

import geopandas as gpd
import matplotlib.pyplot as plt

world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
world = world[(world.name != "Antarctica") & (world.name != "Fr. S. Antarctic Lands")]
world = world.to_crs("EPSG:3395")

ax = world.plot()
ax.set_title("Mercator")

plt.show()

地图本身输出正确，但坐标轴旁边的经度/纬度比例不正确，如下图所示。

墨卡托地图

经度的范围应该是 -180 到 180，纬度的范围应该是 -90 到 90。但是，图像中的比例不是这样。

我怎样才能得到正确的输出比例？

【问题讨论】：

请直接在帖子中添加您的图片，而不是作为链接
看起来您的绘图中的坐标以米为单位，地球的周长约为 4x10^7 米。事实上 EPSG:3395 有米的单位，见epsg.io/3395

标签： python python-3.x matplotlib geopandas

【解决方案1】：

EPSG:3395 是一个以米为单位的墨卡托（圆柱）投影。您需要以度为单位的圆柱投影（称为 Plate Carree 或等距圆柱投影）。使用 EPSG:4326 - 例如参见 https://geopandas.readthedocs.io/en/latest/docs/user_guide/projections.html

更多关于 EPSG:4326 的信息在这里：https://epsg.io/4326

更多关于 EPSG:3395 的信息在这里：https://epsg.io/3395

讨论墨卡托与 Plate Carree（又名 Platte Carre）https://idvux.wordpress.com/2007/06/06/mercator-vs-well-not-mercator-platte-carre/

【讨论】：

【解决方案2】：

要绘制地理标线（子午线和纬线），并用度数标记它们，您可以使用cartopy gridlines() 函数。这是修改后的代码和输出图。

cartopy.gridlines document

import geopandas as gpd
import matplotlib.pyplot as plt
import numpy as np
import cartopy.crs as ccrs

myproj = ccrs.Mercator()
pcproj = ccrs.PlateCarree()

fig = plt.figure(figsize=(16, 7))
extent =[-179,179,  -78, 85]  #lonmin, lonmax, latmin, latmax
ax = plt.axes(projection=myproj )
ax.set_extent(extent, crs=pcproj)

lon_grid = np.arange(-180, 180, 30)
lat_grid = np.arange(-75, 83, 15)
gl = ax.gridlines(draw_labels=True,
                  xlocs=lon_grid, ylocs=lat_grid,
                  x_inline=False, y_inline=False,
                  color='k', linestyle='dotted')

world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
antarctica = world[world.name == "Antarctica"].to_crs("EPSG:3395")  #just in case

world = world[(world.name != "Antarctica") & (world.name != "Fr. S. Antarctic Lands")]
world = world.to_crs("EPSG:3395")

ax = world.plot(ax=ax, edgecolor='k', lw=0.3)

# if you prefer to plot `Antarctica`
if True:
    antarctica.plot(ax=ax, color='lightgray')

ax.set_title("Mercator")

plt.show()

【讨论】：