python - 使用双三次插值的彩色 matplotlib map

Question

我知道 matplotlib 和 scipy 可以进行双三次插值: http://matplotlib.org/examples/pylab_examples/image_interp.html http://docs.scipy.org/doc/scipy/reference/tutorial/interpolate.html http://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp2d.html

我还知道可以用 matplotlib 绘制世界地图: http://matplotlib.org/basemap/users/geography.html http://matplotlib.org/basemap/users/examples.html http://matplotlib.org/basemap/api/basemap_api.html

但是我可以根据 4 个数据点进行双三次插值并只为陆地着色吗？

例如将这些用于 4 个数据点(经度和纬度)和颜色:

Lagos: 6.453056, 3.395833; red HSV 0 100 100 (or z = 0)
Cairo: 30.05, 31.233333; green HSV 90 100 100 (or z = 90)
Johannesburg: -26.204444, 28.045556; cyan HSV 180 100 100 (or z = 180)
Mogadishu: 2.033333, 45.35; purple HSV 270 100 100 (or z = 270)

我在想，一定可以在纬度和经度范围内进行双三次插值，然后在该层的顶部添加海洋、湖泊和河流吗？我可以用 drawmapboundary 来做到这一点。实际上有一个选项 maskoceans 为此: http://matplotlib.org/basemap/api/basemap_api.html#mpl_toolkits.basemap.maskoceans

我可以像这样插入数据:

xnew, ynew = np.mgrid[-1:1:70j, -1:1:70j]
tck = interpolate.bisplrep(x, y, z, s=0)
znew = interpolate.bisplev(xnew[:,0], ynew[0,:], tck)

或者使用 scipy.interpolate.interp2d: http://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp2d.html

这里解释了如何转换为 map 投影坐标: http://matplotlib.org/basemap/users/mapcoords.html

但我需要弄清楚如何为计算出的表面而不是单个点执行此操作。实际上有一个使用外部数据的地形图示例，我应该可以复制它: http://matplotlib.org/basemap/users/examples.html

附言我不是在寻找一个完整的解决方案。我更愿意自己解决这个问题。相反，我正在寻找建议和提示。我已经使用 gnuplot 超过 10 年，并且在过去几周内才切换到 matplotlib，所以请不要以为我知道关于 matplotlib 的最简单的事情。

Answer 1

我认为这就是您正在寻找的(大致)。请注意，在绘制 pcolor 并传入 hsv 颜色图(文档:cmap parameter for pcolormesh 和 available colormaps)之前，重要的事情是屏蔽数据数组。

我将绘制 map 的代码与示例保持得非常接近，因此应该很容易理解。出于同样的原因，我保留了您的插值代码。请注意，插值是线性的而不是三次插值 - kx=ky=1 - 因为你没有提供足够的点数来进行三次插值(你至少需要 16 - scipy 会提示少说“m 必须 >= (kx+1)(ky+1)”，尽管 documentation 中未提及该约束)。

我还扩展了您的网格的范围，并始终保持 x 和 y 的纬度/经度。

代码

from mpl_toolkits.basemap import Basemap,maskoceans
import matplotlib.pyplot as plt
import numpy as np
from scipy import interpolate
# set up orthographic map projection with
# perspective of satellite looking down at 0N, 20W (Africa in main focus)
# use low resolution coastlines.
map = Basemap(projection='ortho',lat_0=0,lon_0=20,resolution='l')

# draw coastlines, country boundaries
map.drawcoastlines(linewidth=0.25)
map.drawcountries(linewidth=0.25)
# Optionally (commented line below) give the map a fill colour - e.g. a blue sea
#map.drawmapboundary(fill_color='aqua')
# draw lat/lon grid lines every 30 degrees.
map.drawmeridians(np.arange(0,360,30))
map.drawparallels(np.arange(-90,90,30))


data = {'Lagos': (6.453056, 3.395833,0),
        'Cairo': (30.05, 31.233333,90),
        'Johannesburg': (-26.204444, 28.045556,180),
        'Mogadishu': (2.033333, 45.35, 270)}

x,y,z = zip(*data.values())

xnew, ynew = np.mgrid[-30:60:0.1, -50:50:0.1]

tck = interpolate.bisplrep(x, y, z, s=0,kx=1,ky=1)
znew = interpolate.bisplev(xnew[:,0], ynew[0,:], tck)
znew = maskoceans(xnew, ynew, znew)

col_plot = map.pcolormesh(xnew, ynew, znew, latlon=True, cmap='hsv')
plt.show()

输出