python - 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()投影时，您根本看不到世界地图。这是因为世界地图的规模很大。下面是两个示例的图像：

板卡里投影：

中断 Goode Homolosine 投影（标准缩放）：

中断的 Goode Homolosine 投影（缩小）：

如果有人能解释为什么会发生这种情况，以及如何解决它，以便我可以在不同的投影上绘制相同的数据，那将不胜感激。

编辑：

我还想指定我尝试transform=projec在我包含的示例中添加到第 37 行，即：

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

然而这并没有帮助。事实上，添加了这个之后，世界地图似乎根本不再出现。

编辑：

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

并缩小：

标签： pythonmatplotlibscipykdecartopy

Answer 1

对你的情节的评论：

Plot1：（参考地图）

• 投影：PlateCarree 投影
• （Zgrid）图像范围覆盖（大约）正方形区域，每边约 40 度
• 图像的左下角在 lat/long: (0,0)

情节2

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

A：该地块占地面积很小，不包括其中任何一个。

• 投影：InterruptedGoodeHomolosine
• 图像数据，Zgrid 被声明为适合网格（mapprojection）坐标（单位：米）
• 该地图在 x 和 y 方向均绘制在几米的小范围内，并且纵横比不相等。

情节3

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

A：绘图覆盖的区域非常大，图像变得太小而无法绘图。

• 投影：InterruptedGoodeHomolosine 投影
• (Zgrid) 图像范围非常小，在此比例下不可见
• 地图绘制范围大，纵横比不相等。

补救措施（对于 Plot2 和 3）

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

关于“网格线”图

• 对位置参考有用
• 纬度/平行：在这种情况下可以使用 InterruptedGoodeHomolosine
• 经度/经线：有问题（不知道如何解决！！）

这是运行并生成所需地图的修改后的代码。

# 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()

输出地图：