我使用不同的 basemap 投影有一个奇怪的行为。

我想绘制到世界 map 上的测量网格的形状为 [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”,则生成的插值数据与其原点不同。

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

最佳答案

首先,我鼓励您开始使用 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,您应该始终将坐标转换为 map 的坐标系 - 最终这可能意味着 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

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

HTH

关于python - map 投影和强制插值,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/18186047/

10-12 16:48
查看更多