地图投影和强制插值答案

【问题标题】：Map projection and forced interpolation地图投影和强制插值
【发布时间】：2013-08-12 11:40:19
【问题描述】：

我在使用不同的底图投影时有一种奇怪的行为。

我想在世界地图上绘制的测量网格的形状是 [181,83]。这意味着我有每个 2°/2° 点的值，范围从 -180° - 180° 经度和 -82° - 82° 纬度。

from mpl_toolkits.basemap import Basemap
import numpy as np
measurementgrid = np.random.random_sample((181,83))
m = Basemap(projection='cyl',llcrnrlon=-180, llcrnrlat=-82, urcrnrlon=180, urcrnrlat=82, resolution='l')
m.drawcoastlines()
m.drawparallels(np.arange(-90,90,30),labels=[1,0,0,0])
m.drawmeridians(np.arange(-180,180,45), labels=[0,0,0,1])
data, x, y = m.transform_scalar(measurementgrid.T, lons=np.arange(-180,182,2), lats=np.arange(-82,84,2), nx = 181, ny = 83, returnxy=True, order=0)
m.imshow(data, origin='lower', interpolation='none')

使用圆柱投影，返回的数据网格等于测量网格，一切都很好。如果我将投影更改为“mill”，则生成的插值数据与其原点不同。

有没有办法按原样绘制测量网格但相对于不断变化的投影？

【问题讨论】：

您的示例代码没有显示您从哪里获得measurementgrid，（它也没有显示import numpy as np，但我设法猜到了！
我可以建议发布您的修复以帮助遇到类似问题的其他人！
我没有解决这个问题，但我根据你的建议编辑了我的帖子。所以问题依然存在！

标签： python matplotlib matplotlib-basemap

【解决方案1】：

首先，我鼓励您开始使用 pcolormesh，而不是 imshow。在尝试在框中绘制网格数据时，Pcolormesh 应该是您的首选可视化工具（就像在将网格数据绘制为等值区域时一样，contourf 也是如此）。

要使用 pcolormesh，您应该通过数据的 x 和 y 角的坐标，所以：

x = np.linspace(-180, 180, 182)
y = np.linspace(-90, 90, 84)
m.pcolormesh(x, y, data)

但是对于 Basemap，您应该始终将坐标转换为地图的坐标系 - 最终这可能意味着 x 和 y 坐标都需要是二维数组，所以我们这样做并转换：

x = np.linspace(-180, 180, 182)
y = np.linspace(-90, 90, 84)
x, y = np.meshgrid(x, y)
converted_x, converted_y = m(x, y)
m.pcolormesh(converted_x, converted_y, data)

这意味着您现在可以继续更改投影，您的数据将绘制在正确的位置。比如我把投影改成“robin”（罗宾逊），得到如下图：

不幸的是，pcolormesh 适用于连续数据块，如果您选择的投影不具有相同的中心经度（又名“lon_0”），那么您将得到不好的结果。比如我把投影改成projection='robin',lon_0=180，得到如下图：

这是因为 Basemap 目前没有处理日期线，而且据我所知，如果不进行重大重写 - 永远不会。

好消息是，这是一个困扰我很长时间的领域，所以我开始编写一个新包来处理这个问题，以及许多其他科学可视化制图的怪癖。结果是一个名为 cartopy 的新包，它为您做了很多工作，因此日期线之类的东西“正常工作”：

import numpy as np
import cartopy.crs as ccrs
import matplotlib.pyplot as plt

x = np.linspace(-180, 180, 182)
y = np.linspace(-90, 90, 84)
measurement_grid = np.random.random_sample((83, 181)) * y[:-1, np.newaxis] ** 2

plt.axes(projection=ccrs.Robinson(central_longitude=180))
plt.pcolormesh(x, y, measurement_grid, transform=ccrs.PlateCarree())
plt.gca().coastlines()
plt.show()

虽然我不是建议您现在应该更改为 cartopy（安装和性能仍在进行中） - 值得知道该软件包存在，并且我预计将来会变得越来越有吸引力遇到这类问题。 http://scitools.org.uk/cartopy/docs/latest

还值得指出的是，网格数据的科学可视化中出现的许多问题都来自于数据、其坐标及其基础坐标系的处理，因此编写了另一个实现数据模型的包将所有这些复杂信息封装到单个对象中，然后可以将其传递给绘图例程以进行简单接口。再次，我鼓励你看看它http://scitools.org.uk/iris/docs/latest。

HTH

【讨论】：