格点插值。
对于创建平滑图形或使用 barbs 或 quiver 绘图时非常有用。当使用 maskoceans 函数时也非常有用。
interp(datain, xin, yin, xout, yout, checkbounds=False, masked=False, order=1)
0 表示 最邻近插值,1 表示双线性插值, 2 表示三次样条插值,需要安装 scipy.ndimage
from mpl_toolkits.basemap import Basemap
from mpl_toolkits.basemap import interp
import matplotlib.pyplot as plt
from osgeo import gdal
import numpy as np
map = Basemap(llcrnrlon=-82., llcrnrlat=28., urcrnrlon=-79., urcrnrlat=29.5,
projection='lcc', lat_1=30., lat_2=60., lat_0=34.83158, lon_0=-98.,
resolution='i')
ds = gdal.Open("../sample_files/wrf.tiff")
lons = ds.GetRasterBand(4).ReadAsArray()
lats = ds.GetRasterBand(5).ReadAsArray()
u10 = ds.GetRasterBand(1).ReadAsArray()
v10 = ds.GetRasterBand(2).ReadAsArray()
x, y = map(lons, lats)
x2 = np.linspace(x[0][0],x[0][-1],x.shape[1]*2)
y2 = np.linspace(y[0][0],y[-1][0],y.shape[0]*2)
x2, y2 = np.meshgrid(x2, y2)
u10_2 = interp(u10, x[0], np.flipud(y[:, 0]), x2, np.flipud(y2),order=1)
v10_2 = interp(v10, x[0], np.flipud(y[:, 0]), x2, np.flipud(y2),order=1)
map.drawmapboundary(fill_color='#9999FF')
map.fillcontinents(color='#cc9955', lake_color='#9999FF', zorder = 0)
map.drawcoastlines(color = '0.15')
map.barbs(x, y, u10, v10,
pivot='middle', barbcolor='#555555')
map.barbs(x2, y2, u10_2, v10_2,
pivot='middle', barbcolor='#FF3333')
plt.show()
译注:以下是关于上述代码的部分解释
2) 新的网格被创建,其元素数是之前的2倍,但拥有相同的边界
3) 新的网格必须是二维数组。上例中使用 meshgrid 方法创建
灰色为原始点,红色为插值点。
如果所选点在 land 则返回 True,如果在 ocean 或 lake 则返回 False。
is_land(xpt, ypt)
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
map = Basemap(projection='aeqd', lon_0 = 10, lat_0 = 50, resolution='h')
x, y = map(0, 0)
print map.is_land(x, y)
#Idea taken from this post at StackOverflow: http://stackoverflow.com/questions/13796315/plot-only-on-continent-in-matplotlib/13811775#13811775
#I've re-written it to work with modern matplotlib versions
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
from matplotlib.path import Path
import numpy as np
map = Basemap(projection='aeqd', lon_0 = 10, lat_0 = 50, resolution='h')
lons = [0., 0., 16., 76.]
lats = [0., 41., 19., 51.]
x, y = map(lons, lats)
locations = np.c_[x, y]
polygons = [Path(p.boundary) for p in map.landpolygons]
result = np.zeros(len(locations), dtype=bool)
for polygon in polygons:
result += np.array(polygon.contains_points(locations))
print result
makegrid 方法会创建随机格点,这些格点在地图坐标系中是等间距的点。用来获取使用地图投影构成的等间距格点的经度和纬度。
makegrid(nx, ny, returnxy=False)
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np
fig=plt.figure(figsize=(9, 3))
map = Basemap(width=12000000,height=8000000,
resolution='l',projection='stere',
lat_ts=50,lat_0=50,lon_0=-107.)
lons, lats, x, y = map.makegrid(30, 30, returnxy=True)
ax = fig.add_subplot(121)
ax.set_title('The regular grid')
map.scatter(x, y, marker='o')
map.drawcoastlines()
ax = fig.add_subplot(122)
ax.set_title('Projection changed')
map = Basemap(width=12000000,height=9000000,projection='aeqd',
lat_0=50.,lon_0=-105.)
x, y = map(lons, lats)
map.scatter(x, y, marker='o')
map.drawcoastlines()
plt.show()
利用数据mask ocean 和 lake。
mpl_toolkits.basemap.maskoceans(lonsin, latsin, datain, inlands=True, resolution=’l’, grid=5)
from mpl_toolkits.basemap import Basemap
from mpl_toolkits.basemap import maskoceans
from mpl_toolkits.basemap import interp
import matplotlib.pyplot as plt
from osgeo import gdal
import numpy as np
map = Basemap(llcrnrlon=-93.7, llcrnrlat=28., urcrnrlon=-66.1, urcrnrlat=39.5,
projection='lcc', lat_1=30., lat_2=60., lat_0=34.83158, lon_0=-98.,
resolution="h")
ds = gdal.Open("../sample_files/wrf.tiff")
lons = ds.GetRasterBand(4).ReadAsArray()
lats = ds.GetRasterBand(5).ReadAsArray()
data = ds.GetRasterBand(3).ReadAsArray()
x, y = map(lons, lats)
plt.figure(0)
mdata = maskoceans(lons, lats, data)
map.drawcoastlines(color = '0.15')
map.contourf(x, y, mdata)
plt.figure(1)
x2 = np.linspace(x[0][0],x[0][-1],x.shape[1]*5)
y2 = np.linspace(y[0][0],y[-1][0],y.shape[0]*5)
x2, y2 = np.meshgrid(x2, y2)
data2 = interp(data, x[0], np.flipud(y[:, 0]), x2, np.flipud(y2),order=1)
lons2, lats2 = map(x2, y2, inverse=True)
mdata = maskoceans(lons2, lats2, data2, resolution = 'c', grid = 10, inlands=True)
map.drawcoastlines(color = '0.15')
map.contourf(x2, y2, mdata)
plt.figure(2)
mdata = maskoceans(lons2, lats2, data2, resolution = 'h', grid = 10, inlands=True)
map.drawcoastlines(color = '0.15')
map.contourf(x2, y2, mdata)
plt.figure(3)
mdata = maskoceans(lons2, lats2, data2, resolution = 'h', grid = 1.25, inlands=True)
map.drawcoastlines(color = '0.15')
map.contourf(x2, y2, mdata)
plt.show()
第一个例子(第22行)直接进行 mask,图中可以看出结果比较粗糙,这是由于输入数据的精度低,而不是 maskoceans 参数设置的原因。
第二个例子显示的图更精细一些(29-36行),从而避免了由于数据导致的大像素块。可查看 interp 部分获取更多细节。但 maskoceans 函数仍使用的 grid 和 resolution 参数的最低配置。
注意 Florida 的 lake 并没有被 mask,这是因为数据粗糙的原因导致。
在此例中,resolution 设置为 'h',lakes已经被 mask了,但 Florida 海岸仍然显示一些像素块(就是有锯齿存在)。
最后,当 grid 参数也设置为最高精度后, Florida 海岸也能显示的很好(译注:其实放大后看,仍存在锯齿)。
绘制指定日期地图中的夜晚地区。
nightshade(date, color=’k’, delta=0.25, alpha=0.5, ax=None, zorder=2)
下例显示了 van der Grinten 投影的阴影区 [注2]
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
from datetime import datetime
map = Basemap(projection='vandg',lon_0=0,resolution='c')
map.drawmapboundary(fill_color="#7777ff")
map.fillcontinents(color="#ddaa66",lake_color="#7777ff")
map.drawcoastlines()
map.nightshade(datetime.now(), delta=0.2)
plt.show()
注1: 原文为 nx and n define the size of the output grid
注2: http://matplotlib.org/basemap/users/vandg.html