当我有一段跨越反子午线的无网格纬度/经度/数据对时，经度从 -180 交换到 +180，如何防止使用 pcolor(mesh) 进行 cartopy绘制充满整个地球的网格单元？我的问题与 here 相同，除了我使用的是 cartopy 而不是 basemap。对链接问题(关于 basemap)的近 5 年前的评论声称有一个 cartopy 解决方案，但尚未发布。
示例代码:
#!/usr/bin/env python3.6
import numpy
import matplotlib.pyplot
import cartopy.crs
lons = numpy.array([[-174.719, -175.297, -175.883],
[-175.164, -175.734, -176.312],
[-175.594, -176.164, -176.734],
[-176.016, -176.578, -177.148],
[-176.43 , -176.984, -177.547],
[-176.836, -177.383, -177.938],
[-177.227, -177.773, -178.312],
[-177.609, -178.148, -178.688],
[-177.984, -178.516, -179.047],
[-178.352, -178.875, -179.398],
[-179.727, 179.766, 179.266],
[ 179.945, 179.445, 178.945],
[ 179.625, 179.133, 178.641],
[ 179.312, 178.828, 178.336],
[ 179.008, 178.523, 178.039],
[ 178.711, 178.234, 177.75 ],
[ 178.414, 177.945, 177.469],
[ 178.133, 177.656, 177.188],
[ 177.844, 177.383, 176.914],
[ 177.57 , 177.109, 176.648]])
lats = numpy.array([[ 67.391, 67.492, 67.586],
[ 67.055, 67.148, 67.25 ],
[ 66.711, 66.812, 66.906],
[ 66.375, 66.469, 66.562],
[ 66.031, 66.125, 66.219],
[ 65.688, 65.781, 65.875],
[ 65.344, 65.438, 65.523],
[ 65. , 65.094, 65.18 ],
[ 64.656, 64.742, 64.836],
[ 64.312, 64.398, 64.484],
[ 62.922, 63. , 63.086],
[ 62.57 , 62.648, 62.734],
[ 62.219, 62.297, 62.383],
[ 61.867, 61.945, 62.023],
[ 61.516, 61.594, 61.672],
[ 61.164, 61.242, 61.32 ],
[ 60.812, 60.891, 60.961],
[ 60.812, 60.891, 60.961],
[ 60.461, 60.531, 60.609],
[ 60.102, 60.18 , 60.25 ]])
data = numpy.array([[ 231.73, 231.56, 231.22],
[ 231.72, 231.72, 231.72],
[ 232.24, 232.73, 233.37],
[ 233.22, 233.69, 234.01],
[ 234.33, 234.94, 235.39],
[ 234.5 , 235.11, 235.71],
[ 235.41, 235.71, 236. ],
[ 235.27, 235.72, 236.31],
[ 234.67, 235.43, 235.73],
[ 235.43, 236.17, 235.88],
[ 236.18, 236.18, 236.18],
[ 236.07, 236.36, 236.79],
[ 235.8 , 236.1 , 235.8 ],
[ 236.84, 236.84, 236.55],
[ 238.27, 238.27, 238.54],
[ 237.72, 237.44, 237.72],
[ 238.42, 238.28, 238.28],
[ 238.57, 238.57, 238.43],
[ 240.17, 240.04, 239.65],
[ 241.21, 241.21, 241.09]])
proj = cartopy.crs.Mollweide()
ax = matplotlib.pyplot.axes(projection=proj)
trans = proj.transform_points(cartopy.crs.Geodetic(), lons, lats)
ax.coastlines()
ax.pcolormesh(trans[:, :, 0], trans[:, :, 1], data, transform=proj)
matplotlib.pyplot.savefig("/tmp/test.png")
预期输出是一张 map ，其中包含一些以北太平洋某处为中心的数据。实际上，我得到了一张非常长的 map ，横跨整个地球的宽度:
我将数据限制为少量点，以便我可以更轻松地将其合并到问题中，但实际上我有一个完整的极地卫星数据轨道，这些数据总是穿过两极，因此总是穿过反子午线。真实轨道的结果可能如下所示:
更改中心经度可以重新定位问题。我可以通过选择远离 map 边缘的中心经度来降低严重性。在此示例中，绘制了与上一张 map 相同的数据，但中心经度为 90°E:
This pull request从 2012 年开始似乎是相关的，所以显然应该有一个相关的功能，但我不知道如何使用它。任何全局 map 投影都会出现这个问题。我正在使用 cartopy 0.15.1。
如何正确绘制此图？
最佳答案
首先，感谢您提供一些数据和一段代码来重现 - 这意味着我可以快速专注于问题本身，而不是重现问题。
cartopy 和 basemap 之间的主要区别在于 cartopy 可以为您处理矢量/栅格转换。完全有可能让 cartopy 以 basemap 的方式运行，用户需要自己转换数据。您提供的示例正是通过将纬度/经度手动转换为目标投影来实现这一点。如果不非常小心，您很快就会发现与您所遇到的类似的逆经络问题。值得庆幸的是，cartopy在数据转换方面非常小心，我鼓励您使用它。
在伪代码中，您的代码执行以下操作:
create a mollweide map
convert your lats/lons to mollweide coordinate system
plot newly converted mollweide data on mollweide map
在实践中，我们希望用 cartopy 改变范例并执行以下操作:
create a mollweide map
plot lat/lon data on mollweide map
通过这样做，我们为 cartopy 提供了正确转换数据所需的上下文。
对代码的主要更改是绘制原始数据(以纬度/经度为单位)，而不是手动转换的坐标:
ax.pcolormesh(lons, lats, data, transform=ccrs.PlateCarree())
在本例中，我使用的是 PlateCarree 投影，而不是大地坐标系，因为我们目前没有实现大地测量 pcolormesh 框(即带有大圆)，并且本质上是生成恒定纬度/经度的框。
使用它，我们最终会生成一个与您问题中的第一张图像非常相似的图，但这并不完全是您想要的。原因是您定义的某些盒子在 PlateCarree 投影空间中的宽度约为 360 度(这是一张平坦的纸，并且对环绕/环绕一无所知)反子午线)。
让我们看一个人为的例子。如果您从测地线的角度思考，您可能会期望以下代码在 map 的两侧生成两个小方框:
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import shapely.geometry as sgeom
box = sgeom.box(minx=170, maxx=-170, miny=40, maxy=60)
proj = ccrs.Mollweide()
ax = plt.axes(projection=proj)
ax.coastlines()
ax.add_geometries([box], ccrs.PlateCarree(), facecolor='coral',
edgecolor='black', alpha=0.5)
plt.show()
唉，这不是我们得到的。如果我们还记得 Plate Carree 投影是二维笛卡尔投影，其中两点之间唯一有效的线是直线，那么这是有道理的 - 它对跨越反子午线一无所知。
(值得注意的是:如果我们将几何投影更改为大地投影，那么我们会在给定点之间绘制大圆并获得所需的框)
因此，为了生成所需的框，我们需要框的坐标具有较小的 x 范围，而不是接近 360 度的范围。值得庆幸的是，cartopy 确实允许我们定义超过 180 度的 PlateCarree 坐标值 - 这是能够定义具有较小 x 范围的 PlateCarree 框的关键。
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import shapely.geometry as sgeom
box = sgeom.box(minx=170, maxx=190, miny=40, maxy=60)
proj = ccrs.Mollweide()
ax = plt.axes(projection=proj)
ax.coastlines()
ax.add_geometries([box], ccrs.PlateCarree(), facecolor='coral',
edgecolor='black', alpha=0.5)
回到你的例子 - 我们有一堆纬度/经度，它们真正定义了大地测量补丁。 Cartopy 还无法 pcolormesh 大地坐标 - 解决方法是 pcolormesh PlateCarree 坐标。尽管大地坐标和 PlateCarree 坐标点可以互换，但它们具有根本不同的拓扑。
在您给出的示例中，可以通过将 360 添加到低于 0 的值来将数据转换为有效的 PlateCarree 拓扑。不幸的是，这不适用于穿过中央子午线的几何图形 - 这会涉及更多一点，并且将是 cartopy IMO 的有用扩展。
最终代码现在如下所示:
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
lons = np.array([[-174.719, -175.297, -175.883],
[-175.164, -175.734, -176.312],
[-175.594, -176.164, -176.734],
[-176.016, -176.578, -177.148],
[-176.43 , -176.984, -177.547],
[-176.836, -177.383, -177.938],
[-177.227, -177.773, -178.312],
[-177.609, -178.148, -178.688],
[-177.984, -178.516, -179.047],
[-178.352, -178.875, -179.398],
[-179.727, 179.766, 179.266],
[ 179.945, 179.445, 178.945],
[ 179.625, 179.133, 178.641],
[ 179.312, 178.828, 178.336],
[ 179.008, 178.523, 178.039],
[ 178.711, 178.234, 177.75 ],
[ 178.414, 177.945, 177.469],
[ 178.133, 177.656, 177.188],
[ 177.844, 177.383, 176.914],
[ 177.57 , 177.109, 176.648]])
lats = np.array([[ 67.391, 67.492, 67.586],
[ 67.055, 67.148, 67.25 ],
[ 66.711, 66.812, 66.906],
[ 66.375, 66.469, 66.562],
[ 66.031, 66.125, 66.219],
[ 65.688, 65.781, 65.875],
[ 65.344, 65.438, 65.523],
[ 65. , 65.094, 65.18 ],
[ 64.656, 64.742, 64.836],
[ 64.312, 64.398, 64.484],
[ 62.922, 63. , 63.086],
[ 62.57 , 62.648, 62.734],
[ 62.219, 62.297, 62.383],
[ 61.867, 61.945, 62.023],
[ 61.516, 61.594, 61.672],
[ 61.164, 61.242, 61.32 ],
[ 60.812, 60.891, 60.961],
[ 60.812, 60.891, 60.961],
[ 60.461, 60.531, 60.609],
[ 60.102, 60.18 , 60.25 ]])
data = np.array([[ 231.73, 231.56, 231.22],
[ 231.72, 231.72, 231.72],
[ 232.24, 232.73, 233.37],
[ 233.22, 233.69, 234.01],
[ 234.33, 234.94, 235.39],
[ 234.5 , 235.11, 235.71],
[ 235.41, 235.71, 236. ],
[ 235.27, 235.72, 236.31],
[ 234.67, 235.43, 235.73],
[ 235.43, 236.17, 235.88],
[ 236.18, 236.18, 236.18],
[ 236.07, 236.36, 236.79],
[ 235.8 , 236.1 , 235.8 ],
[ 236.84, 236.84, 236.55],
[ 238.27, 238.27, 238.54],
[ 237.72, 237.44, 237.72],
[ 238.42, 238.28, 238.28],
[ 238.57, 238.57, 238.43],
[ 240.17, 240.04, 239.65],
[ 241.21, 241.21, 241.09]])
proj = ccrs.PlateCarree(central_longitude=180)
ax = plt.axes(projection=proj)
ax.coastlines('50m')
ax.margins(0.3)
lons[lons < 0] += 360
ax.pcolormesh(lons, lats, data, transform=ccrs.PlateCarree())
plt.show()
如果有兴趣，我鼓励您打开一个 cartopy 功能请求，以添加一个通常将大地测量 pcolormesh 边界转换为板 carree 边界的函数。 cartopy 追踪器可以在 https://github.com/SciTools/cartopy/issues/new 找到。 .
关于python - 防止未网格 pcolor(mesh) 数据出现虚假水平线，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/46527456/