Cartopy 投影比例不一致
Cartopy projection scale not consistent
我正在使用 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() 投影时，您根本看不到世界地图。这是因为世界地图的比例非常大。下面是两个示例的图像：
Plate Carree 投影：
中断的 Goode Homolosine 投影（标准缩放）：
中断的 Goode Homolosine 投影（缩小）：
如果有人能解释为什么会发生这种情况，以及如何解决它以便我可以在不同的投影上绘制相同的数据，那将不胜感激。
编辑：
我还想说明一下，我尝试将 transform=projec 添加到我包含的示例中的第 37 行，即：
ax.add_feature(cfeature.NaturalEarthFeature('physical', 'land', '50m',
edgecolor='black', facecolor="none", transform=projec))
然而这并没有帮助。事实上，添加这个之后，世界地图似乎根本就不再出现了。
编辑：
为了回应 JohanC 的回答，这是我使用该代码时得到的情节：
并缩小：
对您的地块的评论：
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 到轴的投影坐标
的适当转换
- 地图的范围也需要进行适当的转换和设置
- 必须设置长宽比 'equal'，以防止 x 和 y 方向的不等拉伸
关于'gridlines'地块
- 对于位置参考很有用
- latitude/parallels：在这种情况下可以使用 InterruptedGoodeHomolosine
- longitude/meridians: 有问题（不知道怎么解决！！）
这是运行并生成所需地图的修改后的代码。
# 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()
输出图：
我正在使用 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() 投影时，您根本看不到世界地图。这是因为世界地图的比例非常大。下面是两个示例的图像：
Plate Carree 投影：
中断的 Goode Homolosine 投影（标准缩放）：
中断的 Goode Homolosine 投影（缩小）：
如果有人能解释为什么会发生这种情况，以及如何解决它以便我可以在不同的投影上绘制相同的数据，那将不胜感激。
编辑：
我还想说明一下，我尝试将 transform=projec 添加到我包含的示例中的第 37 行，即：
ax.add_feature(cfeature.NaturalEarthFeature('physical', 'land', '50m',
edgecolor='black', facecolor="none", transform=projec))
然而这并没有帮助。事实上，添加这个之后，世界地图似乎根本就不再出现了。
编辑：
为了回应 JohanC 的回答，这是我使用该代码时得到的情节：
并缩小：
对您的地块的评论：
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 到轴的投影坐标 的适当转换
- 地图的范围也需要进行适当的转换和设置
- 必须设置长宽比 'equal'，以防止 x 和 y 方向的不等拉伸
关于'gridlines'地块
- 对于位置参考很有用
- latitude/parallels：在这种情况下可以使用 InterruptedGoodeHomolosine
- longitude/meridians: 有问题（不知道怎么解决！！）
这是运行并生成所需地图的修改后的代码。
# 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()
输出图：