Cartopy投影比例不一致

Question

我正在使用cartopy显示叠加在世界地图上的KDE。起初，我使用的是ccrs.PlateCarree投影，没有任何问题，但当我尝试使用另一个投影时，它似乎爆炸了投影的规模。为了便于参考，我在下面提供了一个可以在自己的机器上测试的示例（只需注释掉两行projec即可在投影之间切换）

from scipy.stats import gaussian_kde
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature

projec = ccrs.PlateCarree()
#projec = ccrs.InterruptedGoodeHomolosine()


fig = plt.figure(figsize=(12, 12))



ax = fig.add_subplot(projection=projec)

np.random.seed(1)

discrete_points = np.random.randint(0,10,size=(2,400))

kde = gaussian_kde(discrete_points)
x, y = discrete_points
# https://www.oreilly.com/library/view/python-data-science/9781491912126/ch04.html
resolution = 1
x_step = int((max(x)-min(x))/resolution)
y_step = int((max(y)-min(y))/resolution)
xgrid = np.linspace(min(x), max(x), x_step+1)
ygrid = np.linspace(min(y), max(y), y_step+1)
Xgrid, Ygrid = np.meshgrid(xgrid, ygrid)
Z = kde.evaluate(np.vstack([Xgrid.ravel(), Ygrid.ravel()]))
Zgrid = Z.reshape(Xgrid.shape)

ext = [min(x)*5, max(x)*5, min(y)*5, max(y)*5]
earth = plt.cm.gist_earth_r

ax.add_feature(cfeature.NaturalEarthFeature('physical', 'land', '50m', 
                                                edgecolor='black', facecolor="none"))

ax.imshow(Zgrid,
    origin='lower', aspect='auto',
    extent=ext,
    alpha=0.8,
    cmap=earth, transform=projec)

ax.axis('on')
ax.get_xaxis().set_visible(True)
ax.get_yaxis().set_visible(True)

ax.set_xlim(-30, 90)
ax.set_ylim(-60, 60)

plt.show()

您会注意到，当使用ccrs.PlateCarree()投影时，KDE被很好地放置在非洲上方，但是当使用ccrs.InterruptedGoodeHomolosine()投影时，您根本看不到世界地图。这是因为世界地图的比例尺很大。下面是两个示例的图像：

板块卡里投影：

中断的古德同色线投影（标准缩放）：

中断的古德同色线投影（缩小）：

如果有人能解释为什么会发生这种情况，以及如何修复它，这样我就可以在不同的预测上绘制相同的数据，那将不胜感激

编辑：

我还想具体说明，在我所包含的示例中，我尝试将transform=projec添加到第37行，即：

ax.add_feature(cfeature.NaturalEarthFeature('physical', 'land', '50m', 
                                                edgecolor='black', facecolor="none", transform=projec))

然而，这没有帮助。事实上，加上这一点，世界地图似乎不再出现

编辑：

作为对约翰回答的回应，这是我使用该代码时得到的图：

并缩小：

Answer 1

对您的绘图的评论：

绘图1：（参考地图）

• 投影：PlateCarree投影
• （Zgrid）图像范围覆盖（约）正方形区域，每侧约40度
• 图像的左下角位于横向/纵向：（0,0）

Plot2

问：为什么地图上没有显示地形特征

答：该地块占地面积很小，不包括任何一块

• 投影：中断的去氧甘草碱
• 图像数据Zgrid声明为符合栅格（地图投影）坐标（单位：米）
• 地图在x和y方向上绘制在几米的范围内，纵横比不相等

Plot3

问：为什么地图上看不到Zgrid图像

答：绘图覆盖的区域非常大，以至于图像变得太小而无法绘图

• 投影：中断的古德霍洛辛投影
• （Zgrid）图像范围非常小，在此比例下不可见
• 地图在较大范围内打印，且纵横比不相等

补救措施（第2和第3部分）

• Zgrid需要从lat/long到轴的投影坐标进行适当的转换
• 地图的范围也需要进行适当的变换和设置
• 纵横比必须设置为“相等”，以防止x和y方向的拉伸不相等

关于“网格线”绘图

• 用于位置参考
• 纬度/平行线：在这种情况下，适用于中断的古德霍洛辛
• 经度/子午线：有问题（不知道如何修复！！）

下面是运行并生成所需映射的修改代码

# proposed code

from scipy.stats import gaussian_kde
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature

fig = plt.figure(figsize=(7, 12))
ax = plt.axes(projection=ccrs.InterruptedGoodeHomolosine())

np.random.seed(1)
discrete_points = np.random.randint(0,10,size=(2,400))

kde = gaussian_kde(discrete_points)
x, y = discrete_points
# https://www.oreilly.com/library/view/python-data-science/9781491912126/ch04.html
resolution = 1
x_step = int((max(x)-min(x))/resolution)
y_step = int((max(y)-min(y))/resolution)
xgrid = np.linspace(min(x), max(x), x_step+1)
ygrid = np.linspace(min(y), max(y), y_step+1)
Xgrid, Ygrid = np.meshgrid(xgrid, ygrid)
Z = kde.evaluate(np.vstack([Xgrid.ravel(), Ygrid.ravel()]))
Zgrid = Z.reshape(Xgrid.shape)

ext = [min(x)*5, max(x)*5, min(y)*5, max(y)*5]
earth = plt.cm.gist_earth_r


ocean110 = cfeature.NaturalEarthFeature('physical', 'ocean', \
        scale='110m', edgecolor='none', facecolor=cfeature.COLORS['water'])
ax.add_feature(ocean110, zorder=-5)

land110 = cfeature.NaturalEarthFeature('physical', 'land', '110m', \
        edgecolor='black', facecolor="silver")
ax.add_feature(land110, zorder=5)

# extents used by both Zgrid and axes
ext = [min(x)*5, max(x)*5, min(y)*5, max(y)*5]

# plot the image's data array
# note the options: `extent` and `transform`
ax.imshow(Zgrid,
    origin='lower', aspect='auto',
    extent=ext,  #set image's extent
    alpha=0.75,
    cmap=earth, transform=ccrs.PlateCarree(),
    zorder=10)

# set the plot's extent with proper coord transformation
ax.set_extent(ext, ccrs.PlateCarree())

ax.coastlines()
#ax.add_feature(cfeature.BORDERS) #uncomment if you need

ax.gridlines(linestyle=':', linewidth=1, draw_labels=True, dms=True, zorder=30, color='k')
ax.set_aspect('equal')  #make sure the aspect ratio is 1

plt.show()

输出映射：