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

PS我不是在寻找完整的解决方案.我更愿意自己解决这个问题.相反,我正在寻找建议和提示.我已经使用gnuplot超过10年了,并且在过去几周内只转换到了matplotlib,所以请不要假设我甚至知道关于matplotlib的最简单的事情.

Answer 1

请注意，做相反的事情，即在海洋上放置一个栅格并在大陆上放置一个遮罩，非常简单。只需使用map.fillcontinents(). 因此，该解决方案的基本思想是修改该fillcontinents函数，使其在海洋上放置多边形。

\n\n

步骤是：

\n\n

\n
1. 创建一个覆盖整个地球的大圆形多边形。
\n
2. 为数组中的每个形状创建一个多边形map.coastpolygons。
\n
3. shapely使用及其方法将陆地多边形的形状从圆中切出difference。
\n
4. 添加剩余的多边形，这些多边形具有海洋的形状，位于顶部，高度为zorder。
\n

\n\n

代码：

\n\n

from mpl_toolkits.basemap import Basemap\nimport numpy as np\nfrom scipy import interpolate\nfrom shapely.geometry import Polygon\nfrom descartes.patch import PolygonPatch\n\ndef my_circle_polygon( (x0, y0), r, resolution = 50 ):\n    circle = []\n    for theta in np.linspace(0,2*np.pi, resolution):\n        x = r * np.cos(theta) + x0\n        y = r * np.sin(theta) + y0    \n        circle.append( (x,y) )\n\n    return Polygon( circle[:-1] )\n\ndef filloceans(the_map, color=\'0.8\', ax=None):\n    # get current axes instance (if none specified).\n    if not ax:     \n        ax = the_map._check_ax()\n\n    # creates a circle that covers the world\n    r =  0.5*(map.xmax - map.xmin) # + 50000 # adds a little bit of margin\n    x0 = 0.5*(map.xmax + map.xmin)\n    y0 = 0.5*(map.ymax + map.ymin)\n    oceans = my_circle_polygon( (x0, y0) , r, resolution = 100 )\n\n    # for each coastline polygon, gouge it out of the circle\n    for x,y in the_map.coastpolygons:\n        xa = np.array(x,np.float32)\n        ya = np.array(y,np.float32)\n\n        xy = np.array(zip(xa.tolist(),ya.tolist()))\n        continent = Polygon(xy)\n\n        ##  catches error when difference with lakes \n        try:\n            oceans = oceans.difference(continent)\n        except: \n            patch = PolygonPatch(continent, color="white", zorder =150)\n            ax.add_patch( patch ) \n\n    for ocean in oceans:\n        sea_patch = PolygonPatch(ocean, color="blue", zorder =100)\n        ax.add_patch( sea_patch )\n\n###########  DATA\nx = [3.395833, 31.233333, 28.045556, 45.35   ]\ny = [6.453056, 30.05,    -26.204444, 2.033333]\nz = [0, 90, 180, 270]\n\n# set up orthographic map projection\nmap = Basemap(projection=\'ortho\', lat_0=0, lon_0=20, resolution=\'l\')\n\n## Plot the cities on the map\nmap.plot(x,y,".", latlon=1)\n\n# create a interpolated mesh and set it on the map\ninterpol_func = interpolate.interp2d(x, y, z, kind=\'linear\')\nnewx = np.linspace( min(x), max(x) )\nnewy = np.linspace( min(y), max(y) )\nX,Y = np.meshgrid(newx, newy)\nZ = interpol_func(newx, newy)\nmap.pcolormesh( X, Y, Z, latlon=1, zorder=3)\n\nfilloceans(map, color="blue")\n

\n\n

瞧\xc3\xa0:\n

\n