防止未网格 pcolor(mesh) 数据出现虚假水平线
ger*_*rit 5 python matplotlib cartopy
当我有一段跨越反子午线的未网格纬度/经度/数据对时，经度从 -180 交换到 +180，如何防止 cartopy 绘制pcolor(mesh)填充整个地球的网格单元？我的问题与此处的问题相同，只是我使用的是cartopy而不是basemap. 对链接问题（关于basemap）的一条近 5 年的评论声称有一个cartopy解决方案，但尚未发布。
示例代码：
\n\n#!/usr/bin/env python3.6\n\nimport numpy\nimport matplotlib.pyplot\nimport cartopy.crs\n\nlons = numpy.array([[-174.719, -175.297, -175.883],\n [-175.164, -175.734, -176.312],\n [-175.594, -176.164, -176.734],\n [-176.016, -176.578, -177.148],\n [-176.43 , -176.984, -177.547],\n [-176.836, -177.383, -177.938],\n [-177.227, -177.773, -178.312],\n [-177.609, -178.148, -178.688],\n [-177.984, -178.516, -179.047],\n [-178.352, -178.875, -179.398],\n [-179.727, 179.766, 179.266],\n [ 179.945, 179.445, 178.945],\n [ 179.625, 179.133, 178.641],\n [ 179.312, 178.828, 178.336],\n [ 179.008, 178.523, 178.039],\n [ 178.711, 178.234, 177.75 ],\n [ 178.414, 177.945, 177.469],\n [ 178.133, 177.656, 177.188],\n [ 177.844, 177.383, 176.914],\n [ 177.57 , 177.109, 176.648]])\n\nlats = numpy.array([[ 67.391, 67.492, 67.586],\n [ 67.055, 67.148, 67.25 ],\n [ 66.711, 66.812, 66.906],\n [ 66.375, 66.469, 66.562],\n [ 66.031, 66.125, 66.219],\n [ 65.688, 65.781, 65.875],\n [ 65.344, 65.438, 65.523],\n [ 65. , 65.094, 65.18 ],\n [ 64.656, 64.742, 64.836],\n [ 64.312, 64.398, 64.484],\n [ 62.922, 63. , 63.086],\n [ 62.57 , 62.648, 62.734],\n [ 62.219, 62.297, 62.383],\n [ 61.867, 61.945, 62.023],\n [ 61.516, 61.594, 61.672],\n [ 61.164, 61.242, 61.32 ],\n [ 60.812, 60.891, 60.961],\n [ 60.812, 60.891, 60.961],\n [ 60.461, 60.531, 60.609],\n [ 60.102, 60.18 , 60.25 ]])\n\ndata = numpy.array([[ 231.73, 231.56, 231.22],\n [ 231.72, 231.72, 231.72],\n [ 232.24, 232.73, 233.37],\n [ 233.22, 233.69, 234.01],\n [ 234.33, 234.94, 235.39],\n [ 234.5 , 235.11, 235.71],\n [ 235.41, 235.71, 236. ],\n [ 235.27, 235.72, 236.31],\n [ 234.67, 235.43, 235.73],\n [ 235.43, 236.17, 235.88],\n [ 236.18, 236.18, 236.18],\n [ 236.07, 236.36, 236.79],\n [ 235.8 , 236.1 , 235.8 ],\n [ 236.84, 236.84, 236.55],\n [ 238.27, 238.27, 238.54],\n [ 237.72, 237.44, 237.72], \n [ 238.42, 238.28, 238.28],\n [ 238.57, 238.57, 238.43],\n [ 240.17, 240.04, 239.65],\n [ 241.21, 241.21, 241.09]])\n\nproj = cartopy.crs.Mollweide() \nax = matplotlib.pyplot.axes(projection=proj)\ntrans = proj.transform_points(cartopy.crs.Geodetic(), lons, lats)\nax.coastlines()\nax.pcolormesh(trans[:, :, 0], trans[:, :, 1], data, transform=proj)\n\nmatplotlib.pyplot.savefig("/tmp/test.png")\n预期输出将是一张地图，其中包含一些以北太平洋某处为中心的数据。实际上，我得到了一张非常长的地图，横跨整个地球的宽度：
\n\n
我将数据限制为少量点，以便我可以更轻松地将其合并到问题中，但实际上我有一个完整的极地卫星数据轨道，这些数据总是穿过两极，因此总是穿过反子午线。真实轨道的结果可能如下所示：
\n\n
改变中心经度可以重新定位问题。我可以通过选择远离地图边缘的中心经度来降低严重性。在此示例中，绘制了与上一个地图中相同的数据，但中心经度为 90\xc2\xb0E：
\n\n
这个 2012 年的拉取请求看起来是相关的，所以显然应该有一个相关的功能，但我不知道如何使用它。任何全球地图投影都会出现这个问题。我正在使用 cartopy 0.15.1。
\n\n我怎样才能正确地绘制这个？
\n首先，感谢您提供一些数据和一段代码来重现 - 这意味着我可以快速专注于问题本身，而不是重现问题。
cartopy 和底图之间的主要区别在于 cartopy 可以为您处理矢量/栅格转换。完全有可能让 cartopy 以底图的方式运行，用户需要自己转换数据。您提供的示例正是通过将纬度/经度手动转换为目标投影来实现这一点。如果不非常小心，您很快就会发现与您所遇到的类似的逆经络问题。值得庆幸的是，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 度。
让我们看一个人为的例子。如果您从测地线的角度思考，您可能会期望以下代码在地图的两侧生成两个小框：
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 的一个有用的扩展。
最终代码现在如下所示：
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找到。
编辑：请注意，对此有两个后续问题：
这将总共显示给定问题的解决方法**
原答案：
我不确定以下内容是否有很大帮助，但让我们尝试一下。如果您能够在 -180 和 180 之间发生偏移的位置分割数据，您可以绘制两个不同的 pcolor 图，并以这种方式防止 pcolorshape 绕地球旋转一次。
使用测试数据，这相当容易；我刚刚更改了其中一个数据点的符号。但也可以省略一行数据左右。然后分别绘制数据即可得到所需的结果。请注意，需要进行标准化以确保两个图的颜色映射是同步的。
norm = plt.Normalize(data.min(), data.max())
ax.pcolormesh(trans[:10, :, 0], trans[:10, :, 1], data[:10,:], transform=proj, norm=norm)
ax.pcolormesh(trans[10:, :, 0], trans[10:, :, 1], data[10:,:], transform=proj, norm=norm)
……
完整代码：
import numpy
import matplotlib.pyplot as plt
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.999, 179.766, 179.266], #<- changed sign here
[ 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]])
print lons.shape, lats.shape, data.shape
proj = cartopy.crs.Mollweide()
ax = plt.axes(projection=proj)
trans = proj.transform_points(cartopy.crs.Geodetic(), lons, lats)
ax.coastlines()
norm = plt.Normalize(data.min(), data.max())
ax.pcolormesh(trans[:10, :, 0], trans[:10, :, 1], data[:10,:], transform=proj, norm=norm)
ax.pcolormesh(trans[10:, :, 0], trans[10:, :, 1], data[10:,:], transform=proj, norm=norm)
plt.show()
| 归档时间： |
|
| 查看次数： |
2701 次 |
| 最近记录： |