防止虚假的pcolor(mesh)数据出现伪水平线 [英] preventing spurious horizontal lines for ungridded pcolor(mesh) data
问题描述
当我有一排跨过子午线的无纬纬度/经度/数据对，从而使经度从-180变为+180时，如何防止pcolor(mesh)的矩形绘制填充整个地球的网格单元?我的问题与此处相同，除了我使用的是cartopy而不是basemap.对链接的问题(大约与basemap有关)已有近5年的评论，声称有一个cartopy解决方案，但尚未发布.
When I have a stretch of ungridded lat/lon/data pairs that cross the antimeridian, such that longitudes swap from -180 to +180, how can I prevent cartopy with pcolor(mesh) from drawing grid cells filling the entire globe? My problem is identical to the one here, except I'm using cartopy rather than basemap. A nearly 5 year old comment to the linked question (which is about basemap) claims there is a cartopy solution but such has not been posted.
示例代码:
#!/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")
预期输出将是一张地图，其中有一些数据以北太平洋某处为中心.实际上，我得到了一张拉长的地图，横跨整个地球的宽度:
Expected output would be a map with a bit of data centred somewhere in the North Pacific Ocean. In reality, I get a very elongated map spanning the width of the entire Earth:
我将数据限制在少数几个点上，以便我可以更轻松地将其合并到问题中，但实际上，我有一个完整的极地卫星数据轨道，该轨道始终穿过两个极点，因此始终穿过antimeridian.真实轨道的结果可能看起来像这样:
I've limited the data to a small number of points such that I can more easily incorporate it in the question, but in reality I have a full orbit of polar satellite data that always crosses both poles, and therefore always crosses the antimeridian. The result for a real orbit may look like this:
更改中心经度可以解决此问题.通过选择远离地图边缘的中心经度，可以降低严重程度.在此示例中，绘制了与上一幅地图相同的数据，但中心经度为90°E:
Changing the central longitude relocates the problem. I can reduce the severity by choosing the central longitude away from where I cross the map edge. In this example, the same data as in the previous map are plotted but with a central longitude of 90°E:
2012年的此拉取请求似乎是相关的，因此显然应该是一个相关的功能，但我不知道如何使用它.任何全局地图投影都会出现问题.我正在使用Cartopy 0.15.1.
This pull request from 2012 appears related, so apparently there is supposed to be a related feature, but I have no clue how to use it. The problem appears with any global map projection. I'm using cartopy 0.15.1.
如何正确绘制?
推荐答案
首先，感谢您提供一些数据和一段可重现的代码-这意味着我可以快速专注于问题本身，而不是重现问题
Firstly, thanks for providing some data and a piece of code to reproduce - it meant that I could quickly focus on the issue itself, and not on reproducing the problem.
cartopy和底图之间的主要区别在于，cartopy可以为您处理矢量/栅格转换.完全有可能让底牌以底图的方式进行操作，这取决于用户自己转换数据.您提供的示例正是通过将纬度/经度手动转换为目标投影来完成此操作的.无需过多照顾，您将快速找到诸如遇到的问题之类的子午线问题.值得庆幸的是，cartopy在数据转换方面已经非常注意了，我鼓励您使用它.
The major difference between cartopy and basemap is that cartopy can handle vector/raster transformations for you. It is entirely possible to get cartopy to operate in basemap's fashion where it is beholden on the user to transform their data themselves. The example you have provided is doing precisely this by transforming the lats/lons manually into the target projection. Without a great deal of care, you are going to quickly find antimeridian issues such as the one you have encountered. Thankfully, cartopy has taken a great deal of care with regards to data transformation, and I encourage you to make use of it.
在伪代码中，您的代码可以:
In pseudo code, your code does:
create a mollweide map
convert your lats/lons to mollweide coordinate system
plot newly converted mollweide data on mollweide map
在实践中，我们想用cartopy改变范例并这样做:
In practice, we want to change the paradigm with cartopy and do:
create a mollweide map
plot lat/lon data on mollweide map
这样做，我们为cartopy提供了必要的上下文，以正确地转换您的数据.
By doing so, we are giving cartopy the necessary context to transform your data correctly.
代码的主要变化是绘制原始数据(以经度/经度为单位)，而不是您手工转换的坐标:
The major change to your code is to plot the original data (in lats/lons), not the coordinates you transformed by hand:
ax.pcolormesh(lons, lats, data, transform=ccrs.PlateCarree())
在这种情况下，我使用的是PlateCarree投影，而不是大地坐标系，因为我们目前尚未实现大地测量的pcolormesh框(即带有大圆的框)，并且实际上是在制作经纬度恒定的框.
In this instance I've used the PlateCarree projection not the Geodetic coordinate system as we don't currently implement geodetic pcolormesh boxes (i.e. with great circles) and are essentially producing boxes of constant lat/lon.
使用此方法，我们最终生成了一个与您问题中的第一张图片非常相似的图，这并非您想要的.原因是您要定义的某些盒子在PlateCarree投影空间中的宽度为〜360度(这是一张平坦的纸，不知道任何环绕效果/antimeridian).
Using this, we end up producing a plot very similar to the first image in your question, which isn't exactly what you wanted. The reason for this is that some of the boxes you are defining have a width of ~360 degrees in the PlateCarree projected space (which is a flat piece of paper, and knows nothing of wrap-arounds/the antimeridian).
让我们看一个人为的例子.如果您以测地学的角度来考虑，您可能会 期望以下代码在地图的任一侧产生两个小方框:
Let's take a look at a contrived example. If you think in terms of geodesics, you might expect the following code to produce two small boxes on either side of the 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()
A，那不是我们得到的.如果我们还记得Plate Carree投影是2D笛卡尔投影，则这是有道理的，其中两个点之间的唯一有效线是直线-它对绕过子午线一无所知.
Alas, that is not what we get. This makes sense if we remember that the Plate Carree projection is a 2d cartesian projection where the only valid line between two points is a straight line - it knows nothing about wrapping across the antimeridian.
(值得注意的是:如果将几何投影更改为大地测量，则在给定点之间绘制大圆并获得所需的框)
(it is worth noting: if we were to change the geometry projection to geodetic then we draw great circles between the given points and get the desired boxes)
因此，要生成所需的盒子，我们需要盒子的坐标具有较小的x范围，而不是接近360度的范围.幸运的是，cartopy确实使我们能够定义超过180度的PlateCarree坐标值-这是能够定义x范围较小的PlateCarree框的关键.
So to produce the desired boxes we need the box's coordinates to have a small x range, not one nearing 360 degrees. Thankfully cartopy does allow us to define PlateCarree coordinate values beyond 180 degrees - this is the key to being able to define a PlateCarree box with a small x range.
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坐标 points 是可以互换的，但它们的拓扑结构却根本不同.
So going back to your example - we have a bunch of lat/lons, which are really defining geodetic patches. Cartopy can't pcolormesh geodetic coordinates yet - the workaround is to pcolormesh PlateCarree coordinates. Despite geodetic coordinates and PlateCarree coordinate points being interchangeable, they have fundamentally different topology.
在此示例中，您可以通过将360加到0以下的值来将数据转换为有效的PlateCarree拓扑.不幸的是，这不适用于跨越中央子午线的几何图形-这将涉及更多点，这将是对Cartopy IMO的有用扩展.
In the example you have given it is possible to convert your data to valid PlateCarree topology by adding 360 to the values below 0. Unfortunately this wouldn't work for geometries that cross the central meridian - that would be a little more involved, and would be a useful extension to 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边界转换为carrees边界的函数.可以在> https://github.com/SciTools/cartopy/issues/new中找到Cartopy跟踪器.
If interested, I encourage you to open a cartopy feature request to add a function that generally converts geodetic pcolormesh bounds into plate carree ones. The cartopy tracker can be found at https://github.com/SciTools/cartopy/issues/new.
这篇关于防止虚假的pcolor(mesh)数据出现伪水平线的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
