使用双三次插值为Matplotlib着色图像

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

但是，我能否基于四个数据点进行双三次插值，并且只给陆地上上色呢？

比如说，使用这四个数据点（经度和纬度）和颜色：

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

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

但我需要弄清楚如何针对计算出的表面而不是单个点来做这个。实际上，有一个使用外部数据的地形图示例，我应该能够复制这个：http://matplotlib.org/basemap/users/examples.html

附言：我并不在寻找一个完整的解决方案。我更希望自己解决这个问题。其实我在使用 gnuplot 已经超过十年了，最近几周才转到 matplotlib，所以请不要认为我对 matplotlib 的任何简单知识都了解。

Answer 1

注意，做相反的事情，即在海洋上放一个栅格，然后在大陆上覆盖一个遮罩，其实非常简单。只需使用 map.fillcontinents() 就可以了。所以这个解决方案的基本思路是修改 fillcontinents 函数，让它在海洋上覆盖多边形。

具体步骤如下：

1. 创建一个像大圆一样的多边形，覆盖整个地球。
2. 为 map.coastpolygons 数组中的每个形状创建一个多边形。
3. 使用 shapely 和它的 difference 方法，从圆形中切掉陆地的形状。
4. 将剩下的多边形（形状是海洋的）放在最上面，设置一个较高的 zorder。

代码如下：

from mpl_toolkits.basemap import Basemap
import numpy as np
from scipy import interpolate
from shapely.geometry import Polygon
from descartes.patch import PolygonPatch

def my_circle_polygon( (x0, y0), r, resolution = 50 ):
    circle = []
    for theta in np.linspace(0,2*np.pi, resolution):
        x = r * np.cos(theta) + x0
        y = r * np.sin(theta) + y0    
        circle.append( (x,y) )

    return Polygon( circle[:-1] )

def filloceans(the_map, color='0.8', ax=None):
    # get current axes instance (if none specified).
    if not ax:     
        ax = the_map._check_ax()

    # creates a circle that covers the world
    r =  0.5*(map.xmax - map.xmin) # + 50000 # adds a little bit of margin
    x0 = 0.5*(map.xmax + map.xmin)
    y0 = 0.5*(map.ymax + map.ymin)
    oceans = my_circle_polygon( (x0, y0) , r, resolution = 100 )

    # for each coastline polygon, gouge it out of the circle
    for x,y in the_map.coastpolygons:
        xa = np.array(x,np.float32)
        ya = np.array(y,np.float32)

        xy = np.array(zip(xa.tolist(),ya.tolist()))
        continent = Polygon(xy)

        ##  catches error when difference with lakes 
        try:
            oceans = oceans.difference(continent)
        except: 
            patch = PolygonPatch(continent, color="white", zorder =150)
            ax.add_patch( patch ) 

    for ocean in oceans:
        sea_patch = PolygonPatch(ocean, color="blue", zorder =100)
        ax.add_patch( sea_patch )

###########  DATA
x = [3.395833, 31.233333, 28.045556, 45.35   ]
y = [6.453056, 30.05,    -26.204444, 2.033333]
z = [0, 90, 180, 270]

# set up orthographic map projection
map = Basemap(projection='ortho', lat_0=0, lon_0=20, resolution='l')

## Plot the cities on the map
map.plot(x,y,".", latlon=1)

# create a interpolated mesh and set it on the map
interpol_func = interpolate.interp2d(x, y, z, kind='linear')
newx = np.linspace( min(x), max(x) )
newy = np.linspace( min(y), max(y) )
X,Y = np.meshgrid(newx, newy)
Z = interpol_func(newx, newy)
map.pcolormesh( X, Y, Z, latlon=1, zorder=3)

filloceans(map, color="blue")

完成了：

Answer 2

我觉得这就是你想要的（大致上）。注意，关键的地方是在你绘制 pcolor 之前要对数据数组进行遮罩处理，并且要使用 hsv 这个颜色映射（文档链接：pcolormesh 的 cmap 参数 和 可用的颜色映射）。

我把绘制地图的代码保持得很接近示例，这样应该容易理解。我也保留了你的插值代码，原因是一样的。注意，这里的插值是线性的，而不是立方的 - kx=ky=1 - 因为你提供的点不够做立方插值（至少需要16个点 - 如果少于这个，scipy会报错，提示 "m must be >= (kx+1)(ky+1)"，虽然这个限制在文档中没有提到）。

我还扩展了你的网格范围，并在整个过程中保留了经纬度作为 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()

输出