Basemap系列教程:绘制子图及小地图
使用 matplotlib 中的 subplots 可以在同一个 figure 中绘制多个地图。有几种方法可以实现这种图形的绘制,而且根据所绘图形的复杂性来选择不同的方法:
- 直接使用 add_subplot 添加 axis
- 使用 pylab.subplots 创建子图
- 使用 subplot2grid
- 创建 inset locators [ 注1 ]
使用 add_subplot
这是大部分情况下一种很好的子图添加方式。
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
fig = plt.figure()
ax = fig.add_subplot(211)
ax.set_title("Hammer projection")
map = Basemap(projection='hammer', lon_0 = 10, lat_0 = 50)
map.drawmapboundary(fill_color='aqua')
map.fillcontinents(color='coral',lake_color='aqua')
map.drawcoastlines()
ax = fig.add_subplot(212)
ax.set_title("Robinson projection")
map = Basemap(projection='robin', lon_0 = 10, lat_0 = 50)
map.drawmapboundary(fill_color='aqua')
map.fillcontinents(color='coral',lake_color='aqua')
map.drawcoastlines()
plt.show()
- 调用 Basemap 构造器之前,先调用 fig.add_subplot 方法 add_subplot 方法中的三个数分别表示:子图的行数,列数,当前是第几个图(从图的 左上方 数起) [ 注2 ]
- 只要创建了 axis,后面绘制地图时就会自动使用(当然也可以通过 ax 参数进行传递)
- 每个子图都可以使用 set_title 方法添加 title
使用 plt.subplots 生成子图
在我看来,使用 add_subplot 方法容易让人困惑。如果 basemap 实例没有使用 ax参数的话,很可能会遇到 bug。因此,可以使用 pyplot.subplots 创建子图:
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
fig, axes = plt.subplots(2, 1)
axes[0].set_title("Hammer projection")
map = Basemap(projection='hammer', lon_0 = 10, lat_0 = 50, ax=axes[0])
map.drawmapboundary(fill_color='aqua')
map.fillcontinents(color='coral',lake_color='aqua')
map.drawcoastlines()
axes[1].set_title("Robinson projection")
map = Basemap(projection='robin', lon_0 = 10, lat_0 = 50, ax=axes[1])
map.drawmapboundary(fill_color='aqua')
map.fillcontinents(color='coral',lake_color='aqua')
map.drawcoastlines()
plt.show()
- 传递给 subplots 方法的参数分别为绘制子图的 行数 和 列数
- subplots 方法会返回 figure 对象及创建的 axes 列表(subplots),仍从左上角开始数起
- 创建 basemap 实例时,必须使用 ax 参数,并将创建的 axes实例 传递给 ax参数( 注 :这样做会使一切操作看起来非常明确)
结果和之前的例子相同。
使用 subplot2grid
当 subplots 数目较多或每个子图的大小不一致时,可以使用 subplots2grid 或 gridspec 。
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
from matplotlib.path import Path
import matplotlib.patches as patches
fig = plt.figure()
ax1 = plt.subplot2grid((2,2), (0,0))
ax2 = plt.subplot2grid((2,2), (1,0))
ax3 = plt.subplot2grid((2,2), (0,1), rowspan=2)
map1 = Basemap(projection='ortho', lon_0 = 0, lat_0 = 40, ax=ax1)
map1.drawmapboundary(fill_color='#9999FF')
map1.fillcontinents(color='#ddaa66',lake_color='#9999FF')
map1.drawcoastlines()
map2 = Basemap(projection='cyl', llcrnrlon=-15,llcrnrlat=30,urcrnrlon=15.,urcrnrlat=50., resolution='i', ax=ax2)
map2.drawmapboundary(fill_color='#9999FF')
map2.fillcontinents(color='#ddaa66',lake_color='#9999FF')
map2.drawcoastlines()
map3 = Basemap(llcrnrlon= -1., llcrnrlat=37.5, urcrnrlon=4.5, urcrnrlat=44.5,
resolution='i', projection='tmerc', lat_0 = 39.5, lon_0 = 3, ax=ax3)
map3.drawmapboundary(fill_color='#9999FF')
map3.fillcontinents(color='#ddaa66',lake_color='#9999FF')
map3.drawcoastlines()
#Drawing the zoom rectangles:
lbx1, lby1 = map1(*map2(map2.xmin, map2.ymin, inverse= True))
ltx1, lty1 = map1(*map2(map2.xmin, map2.ymax, inverse= True))
rtx1, rty1 = map1(*map2(map2.xmax, map2.ymax, inverse= True))
rbx1, rby1 = map1(*map2(map2.xmax, map2.ymin, inverse= True))
verts1 = [
(lbx1, lby1), # left, bottom
(ltx1, lty1), # left, top
(rtx1, rty1), # right, top
(rbx1, rby1), # right, bottom
(lbx1, lby1), # ignored
codes2 = [Path.MOVETO,
Path.LINETO,
Path.LINETO,
Path.LINETO,
Path.CLOSEPOLY,
path = Path(verts1, codes2)
patch = patches.PathPatch(path, facecolor='r', lw=2)
ax1.add_patch(patch)
lbx2, lby2 = map2(*map3(map3.xmin, map3.ymin, inverse= True))
ltx2, lty2 = map2(*map3(map3.xmin, map3.ymax, inverse= True))
rtx2, rty2 = map2(*map3(map3.xmax, map3.ymax, inverse= True))
rbx2, rby2 = map2(*map3(map3.xmax, map3.ymin, inverse= True))
verts2 = [
(lbx2, lby2), # left, bottom
(ltx2, lty2), # left, top
(rtx2, rty2), # right, top
(rbx2, rby2), # right, bottom
(lbx2, lby2), # ignored
codes2 = [Path.MOVETO,
Path.LINETO,
Path.LINETO,
Path.LINETO,
Path.CLOSEPOLY,
path = Path(verts2, codes2)
patch = patches.PathPatch(path, facecolor='r', lw=2)
ax2.add_patch(patch)
plt.show()
- 每一个子图都是用 subplot2grid 创建
三个参数分别是 :
1) 输出矩阵形状,二元素序列,分别为 y 和 x 的大小
2) 子图的位置
3) rowspan 或 colspan, 注 :即每个子图占据多少行多少列,默认只占据一行一列
注 :关于子图绘制的方法会在关于 matplotlib 的相关文章中进行解释。
- 此例中展示了地图的不同放大等级。每一个图都是前一个图的部分放大 1) 每一个图的边界框四角位置都要计算 (1) 使用 basemap 实例的 xmin, xmax, ymin, ymax fields 获取四角位置 (2) 由于这些 fields 都是在投影单元中,因此需要设置 inverse 参数 ( 注 :见 Basemap系列教程:Basemap 使用Basemap实例转换单位 部分) (3) 计算出的拐角经纬度被传递给当前的地图实例,从而获得在当前投影中的坐标,会返回一个序列,可以使用 * 进行解包 [ 注3 ] 2)一旦知道这些点后,可以使用 Path 类创建 polygon [ 注4 ]
嵌入定位器 [ 注5 ]
注 :原文此部分单独成节,因为子图部分包括这部分,因此翻译时将此部分与子图部分合并。
使用嵌入定位器可以在大地图中添加小地图,结果比在主地图中创建子图要好。
嵌入定位器是一个非常酷的类,可以放大一个图的局部,并绘制在这个图上,从而展示某一块区域。 注 :比如用来在地图拐角显示南海地区。
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes
from mpl_toolkits.axes_grid1.inset_locator import mark_inset
import numpy as np
fig = plt.figure()
ax = fig.add_subplot(111)
map = Basemap(projection='cyl',
lat_0=0, lon_0=0)
map.drawmapboundary(fill_color='#7777ff')
map.fillcontinents(color='#ddaa66', lake_color='#7777ff', zorder=0)
map.drawcoastlines()
lons = np.array([-13.7, -10.8, -13.2, -96.8, -7.99, 7.5, -17.3, -3.7])
lats = np.array([9.6, 6.3, 8.5, 32.7, 12.5, 8.9, 14.7, 40.39])
cases = np.array([1971, 7069, 6073, 4, 6, 20, 1, 1])
deaths = np.array([1192, 2964, 1250, 1, 5, 8, 0, 0])
places = np.array(['Guinea', 'Liberia', 'Sierra Leone','United States', 'Mali', 'Nigeria', 'Senegal', 'Spain'])
x, y = map(lons, lats)
map.scatter(x, y, s=cases, c='r', alpha=0.5)
axins = zoomed_inset_axes(ax, 7, loc=1)
axins.set_xlim(-20, 0)
axins.set_ylim(3, 18)
plt.xticks(visible=False)
plt.yticks(visible=False)
map2 = Basemap(llcrnrlon=-20,llcrnrlat=3,urcrnrlon=0,urcrnrlat=18, ax=axins)
map2.drawmapboundary(fill_color='#7777ff')
map2.fillcontinents(color='#ddaa66', lake_color='#7777ff', zorder=0)
map2.drawcoastlines()
map2.drawcountries()
map2.scatter(x, y, s=cases/5., c='r', alpha=0.5)
mark_inset(ax, axins, loc1=2, loc2=4, fc="none", ec="0.5")
plt.show()
- 使用 zoomed_inset_axes 创建嵌入定位器 1) ax 是要添加 嵌入定位器的 axis 2) 7 是放大等级 3) loc 是嵌入定位器的位置(此例中是右上角)
- set_xlim 和 set_ylim 可以改变嵌入定位器所覆盖的范围。因为使用的是 cyl 投影,因此可以直接使用 longitudes 和 latitudes ,而不用进行转换
- xticks 和 yticks 设置为False,从而隐藏新轴的标签(labels),此例中没必要显示
- 使用新的 basemap 实例绘制放大区域
- mark_inset 用于标记放大区域
嵌入南海地区
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes
from mpl_toolkits.axes_grid1.inset_locator import mark_inset
import netCDF4 as nc
import numpy as np
import shapefile
from matplotlib.path import Path
from matplotlib.patches import PathPatch
def basemask(cs, ax, map, shpfile):
sf = shapefile.Reader(shpfile)
vertices = []
codes = []
for shape_rec in sf.shapeRecords():
if shape_rec.record[3] >= 0:
pts = shape_rec.shape.points
prt = list(shape_rec.shape.parts) + [len(pts)]
for i in range(len(prt) - 1):
for j in range(prt[i], prt[i+1]):
vertices.append(map(pts[j][0], pts[j][1]))
codes += [Path.MOVETO]
codes += [Path.LINETO] * (prt[i+1] - prt[i] -2)
codes += [Path.CLOSEPOLY]
clip = Path(vertices, codes)
clip = PathPatch(clip, transform = ax.transData)
for contour in cs.collections:
contour.set_clip_path(clip)
shpfile = u'E:\MATLAB\shp\中国行政区_包含南海九段线'
data = nc.Dataset('E:\python\scripts\sss\pres.mon.ltm.nc')
pres = data.variables['pres'][0,14:35,28:57]
lat = data.variables['lat'][:]
lon = data.variables['lon'][:]
lons, lone = lon[28], lon[57]
lats, late = lat[35], lat[14]
nx = pres.shape[1]
ny = pres.shape[0]
fig = plt.figure()
ax = fig.add_subplot(111)
map = Basemap(llcrnrlon = lon1, llcrnrlat = lat1, urcrnrlon = lon2, urcrnrlat = lat2,
projection = 'cyl', ax = ax)
map.drawparallels(np.arange(lats, late + 1, 5), labels = [1, 0, 0, 0])
map.drawmeridians(np.arange(lons, lone + 1, 10), labels = [0, 0, 0, 1])
x, y = map.makegrid(nx, ny)
# 转换坐标
x, y = map(x, y[::-1])
map.readshapefile(shpfile, 'China')
cs = map.contourf(x, y, pres, vmin = int(pres.min()), vmax = int(pres.max()))
basemask(cs, ax, map, r'E:\MATLAB\shp\bou2_4p')
axins = zoomed_inset_axes(ax, 0.8, loc = 3)
axins.set_xlim(lons + 38, lons + 52)
axins.set_ylim(lons, lons + 20)
map2 = Basemap(llcrnrlon = lons + 38, llcrnrlat = lats, urcrnrlon = lons + 52, urcrnrlat = lats + 20,
ax = axins)
map2.readshapefile(shpfile, 'China')
cs2 = map2.contourf(x, y, pres, vmin = int(pres.min()), vmax = int(pres.max()))