使用双三次插值的彩色 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

但是我可以根据 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

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

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

附：我不是在寻找一个完整的解决方案。我更愿意自己解决这个问题。相反，我正在寻找建议和提示。我使用 gnuplot 已经 10 多年了，并且在过去几周内才切换到 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 颜色图之前，关键的事情是屏蔽数据数组（文档：cmap parameter for pcolormesh 和available colormaps）。

我将用于绘制地图的代码保留在与示例非常接近的地方，因此应该很容易理解。出于同样的原因，我保留了您的插值代码。请注意，插值是线性的而不是三次的 - kx=ky=1 - 因为你没有给出足够的点来进行三次插值（你至少需要 16 - scipy 会抱怨说 "m must be >= (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()

输出