专栏首页气象杂货铺Basemap系列教程:绘制子图及小地图

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

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

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

plt.show()

注意:绘图时一定要设置 vmin 和 vmax 参数,而且两个图的值要设置为一致,否则画出的图可能导致相同区域配色不一致。


注1:http://basemaptutorial.readthedocs.io/en/latest/locator.html#inset-locator

注2:http://stackoverflow.com/questions/3584805/in-matplotlib-what-does-the-argument-mean-in-fig-add-subplot111

注3:https://docs.python.org/2/tutorial/controlflow.html#unpacking-argument-lists

注4:http://matplotlib.org/users/path_tutorial.html

注5:http://basemaptutorial.readthedocs.io/en/latest/locator.html

本文分享自微信公众号 - 气象杂货铺(meteogs),作者:lightning

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2017-03-22

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • Basemap系列教程:绘图

    (1) reduce_C_function 参数用于显示每一个 bin 的最大值,从而代替平均值

    bugsuse
  • Basemap系列教程:3D

    尽管很多人不喜欢3D地图,但是仍可以使用 Basemap 和 matplotlib mplot3d [注1] 工具绘制3D地图。

    bugsuse
  • Basemap工具函数(2)

    对于创建平滑图形或使用 barbs 或 quiver 绘图时非常有用。当使用 maskoceans 函数时也非常有用。

    bugsuse
  • TestNG参数化测试-@DataProvider

    欲将沉醉换悲凉,清歌莫断肠。这混乱的尘世,究竟充斥了多少绝望和悲伤。你想去做一个勇敢的男子,为爱,为信仰,轰轰烈烈的奋斗一场。

    louiezhou001
  • 心理医生妈妈是怎样育儿的?

    用户1756920
  • Shiro集成Spring

    1、Shiro集成Spring,使用maven进行jar包的依赖与管理,pom.xml的配置文件,如下所示:

    别先生
  • SDNLAB群分享(四):利用ODL下发流表创建VxLAN网络

    今天想跟大家分享如何通过ODL控制器下发流表来创建VxLAN网络。ODL作为当前流行的控制器,已经有广泛的应用。基于ODL提供了丰富的北向接口,使得应用对网络有...

    SDNLAB
  • SDNLAB群分享(四):利用ODL下发流表创建VxLAN网络

    今天想跟大家分享如何通过ODL控制器下发流表来创建VxLAN网络。ODL作为当前流行的控制器,已经有广泛的应用。基于ODL提供了丰富的北向接口,使得应用对网络有...

    SDNLAB
  • 数据可视化:用散点图进行数据分析

    导读:散点图的用途有很多,我认为它的核心价值,在于应用相关思维,发现变量之间的关系。

    华章科技
  • 手把手教你估算深度神经网络的最优学习率(附代码&教程)

    来源:机器之心 作者:Pavel Surmenok 学习率(learning rate)是调整深度神经网络最重要的超参数之一,本文作者Pavel Surmen...

    数据派THU

扫码关注云+社区

领取腾讯云代金券