使用双三次插值的彩色 matplotlib 地图
Posted
技术标签:
【中文标题】使用双三次插值的彩色 matplotlib 地图【英文标题】：color matplotlib map using bicubic interpolation 【发布时间】：2014-08-21 03:35:06 【问题描述】：我知道 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 的最简单的事情。
【问题讨论】：
我自己也想弄清楚同样的事情。运气好吗？ 这对我来说将是一个圣诞节项目。所以如果你能等到 2015 年…… 您真的想要一个表面，而不是一组（可能非常密集的）点？您是指 interp2d 输出的函数类型吗？ 是的，在这个阶段有一组密集的点就可以了。我似乎找不到时间来做这件事。我应该删除这个问题。为此，我将为自己设定一个截止日期为 4 月 1 日。 看看这个例子：webrian.ch/2012/11/heat-maps-with-python.html 【参考方案1】：请注意，做相反的事情，即在海面上放置一个栅格并在大陆上放置一个掩膜，这很容易。只需使用map.fillcontinents()。所以这个解决方案的基本思想是修改fillcontinents函数，让它在海洋上放置多边形。
步骤如下：
-
创建一个覆盖整个地球的大型圆形多边形。
为
map.coastpolygons 数组中的每个形状创建一个多边形。
使用shapely 及其difference 方法将陆地多边形的形状从圆中剪掉。
在顶部添加具有海洋形状的剩余多边形，并带有高 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")
瞧：
【讨论】：
此解决方案并非没有问题。首先它很慢，然后地图边缘的小陆地没有正确渲染，最后没有测试鲁棒性，它不适用于其他投影。我会尝试改进这些方面并更新答案，但请随时提出建议。 什么是descartes？我好像没有。
拉各斯在您的地图上似乎是青色而不是红色。
好吧，插值栅格可能与您的不同。我的重点是掩盖海洋。至于双三次插值，我认为 4 个点不足以让您这样做。您将需要更大的集合或线性插值。这能回答你的评论吗？【参考方案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()
输出
【讨论】：
@tommy.carstensen 这是否给了您想要的东西，或者这里缺少什么？如果有什么不适合您的，请随时联系我。 我会给你赏金。如果过期了，我会再发一次。我只想根据您的答案提供一个示例，其中使用双三次插值提供更多数据点。您的示例因 16 个或更多数据点而失败。抱歉耽搁了。我会在这个周末完成它。感谢您的出色回答，我已经投了赞成票。 酷 - 没问题。如果你有 16 个数据点（或者可以告诉我从哪里得到它们），我很高兴用那个数据集来看看它。我无法弄清楚这些数据的来源，但我不是非洲城市的专家。即使您授予此赏金，我也会重新考虑 - 这是一个有趣的问题，我不希望您浪费更多的代表。以上是关于使用双三次插值的彩色 matplotlib 地图的主要内容，如果未能解决你的问题，请参考以下文章