Basemap系列教程：绘制子图及小地图

bugsuse

发布于 2020-04-20 13:45:33

4.7K0

发布于 2020-04-20 13:45:33

文章被收录于专栏：气象杂货铺气象杂货铺

使用 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()))

basemask(cs2, axins, map2, r'E:\MATLAB\shp\bou2_4p')

mark_inset(ax, axins, loc1=2, loc2=4, fc = "none", ec = "none")

plt.show()