matplotlib嵌套南海子图

matplotlib嵌套南海子图

二维可视化的最高境界是可以控制每一个像素的每一个通道。——我瞎说的

简介

目前大多数教程中用matplotlib在中国地图右下角嵌套南海子图主要方案是通过fig.add_axes方法，在指定位置绘制指定宽高的新图轴。这种方案可行但存在两个问题：

1. 大陆图轴和南海图轴在代码逻辑上是独立的：一旦使用会挤压图轴的操作如fig.colorbar，大陆图轴在画布上的位置大小发生改变，但在此之前绘制的南海图轴不会跟随移动。这导致该方案下南海图轴需要在最后绘制，一行一列的图还勉强不算太麻烦，但若是要在一张图上画上个几行几列，必须先for循环绘制所有大陆图轴，再绘制颜色图例、标题等等，最后再用先for循环绘制所有南海图轴，很不优雅。
2. 由于地图投影问题，最终图轴的宽高不一定是创建图轴时所指定的宽高。这导致若不在绘制南海图轴前通过投影关系计算好比例，最终两个图轴的右下两条边界可能不会重合，很不优雅。

我琢磨出一种优雅的方法，可以优雅地解决这两个问题，核心方法为ax.inset_axes。结合cartopy需要matplotlib>=3.6.x。

示例

举个简单例子：

import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import geopandas as gpd
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from matplotlib.font_manager import FontProperties

shp_province = gpd.read_file('E:/map/province.shp')  # 省级边界
shp_10dash = gpd.read_file('E:/map/10dash.shp')  # 十段线
font = FontProperties(fname='C:\Windows\Fonts\simsun.ttc')
platecarree = ccrs.PlateCarree()
extent = (70, 135, 15, 55)  # 大陆图轴经纬度范围
extent_southsea = (107, 122, 0, 22)  # 南海图轴的经纬度范围
height_ratio_southsea = 0.4  # 南海图轴高度与大陆图轴高度的比例
ticks_lon = np.arange(75, 136, 5)  # 经度刻度
ticks_lat = np.arange(20, 51, 5)  # 纬度刻度
figsize = (10, 8)
lon_formatter = LongitudeFormatter(zero_direction_label=True)
lat_formatter = LatitudeFormatter()
# 读取待绘制数据
data = xr.load_dataarray('../data.nc')
lon = data['lon']
lat = data['lat']

目标如下图所示：绘制extent范围的大陆图轴，在其右下角绘制extent_southsea范围的南海图轴，高度占大陆图轴高度的比例为height_ratio_southsea，需要绘制颜色图例。

老方案（fig.add_axes）需要这么实现：

from matplotlib.transforms import Bbox

fig, ax = plt.subplots(figsize=figsize, subplot_kw={'projection': platecarree})
pm = ax.pcolormesh(lon, lat, data, cmap='jet', shading='nearest', vmin=1700, vmax=1880)
ax.add_geometries(shp_province['geometry'], crs=platecarree, facecolor='never', edgecolor='black')
ax.add_geometries(shp_10dash['geometry'], crs=platecarree, edgecolor='black')
ax.set_extent(extent, crs=platecarree)  # 控制绘图经纬度范围
ax.set_xticks(ticks_lon, crs=platecarree)
ax.set_yticks(ticks_lat, crs=platecarree)
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
# 颜色图例
cb = fig.colorbar(pm, ax=ax, location='bottom', pad=0.07)
cb.set_label('colorbar')
# 获取大陆图轴在画布中的位置
ax_pos = ax.get_position()
x1, y0 = ax_pos.x1, ax_pos.y0  # 大陆图轴右下角在在画布中的位置
_, _, width, height = ax_pos.bounds  # 大陆图轴的宽高
# 初步创建南海图轴
ax_southsea = fig.add_axes(ax_pos.bounds, projection=platecarree)
ax_southsea.pcolormesh(lon, lat, data, cmap='jet', shading='nearest', vmin=1700, vmax=1880)
ax_southsea.add_geometries(shp_province['geometry'], crs=platecarree, facecolor='never', edgecolor='black')
ax_southsea.add_geometries(shp_10dash['geometry'], crs=platecarree, edgecolor='black')
ax_southsea.set_extent(extent_southsea, crs=platecarree)  # 控制绘图经纬度范围
_, _, width_southsea, height_southsea = ax_southsea.get_position().bounds  # 当前南海图轴的宽高
# 目标南海图轴的宽高
width_southsea = width_southsea / height_southsea * height * height_ratio_southsea
height_southsea = height * height_ratio_southsea
# 缩放并移动南海图轴到大陆图轴的右下角
ax_southsea.set_position(Bbox([[x1 - width_southsea, y0], [x1, y0 + height_southsea]]))
ax_southsea.set_title('南海诸岛', y=0, pad=3, fontproperties=font)
fig.savefig('matplotlib嵌套南海子图.png', bbox_inches='tight')
plt.close(fig)

用新方案（ax.inset_axes）需要这么实现：

fig, ax = plt.subplots(figsize=figsize, subplot_kw={'projection': platecarree})
pm = ax.pcolormesh(lon, lat, data, cmap='jet', shading='nearest', vmin=1700, vmax=1880)
ax.set_extent(extent, crs=platecarree)  # 控制绘图经纬度范围
ax.add_geometries(shp_province['geometry'], crs=platecarree, facecolor='never', edgecolor='black')
ax.add_geometries(shp_10dash['geometry'], crs=platecarree, edgecolor='black')
ax.set_xticks(ticks_lon, crs=platecarree)
ax.set_yticks(ticks_lat, crs=platecarree)
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
# 创建南海图轴
ax_southsea = ax.inset_axes([0, 0, 1, height_ratio_southsea], anchor=(1, 0), projection=platecarree)
ax_southsea.pcolormesh(lon, lat, data, cmap='jet', shading='nearest', vmin=1700, vmax=1880)
ax_southsea.set_extent(extent_southsea, crs=platecarree)  # 控制绘图经纬度范围
ax_southsea.add_geometries(shp_province['geometry'], crs=platecarree, facecolor='never', edgecolor='black')
ax_southsea.add_geometries(shp_10dash['geometry'], crs=platecarree, edgecolor='black')
ax_southsea.set_title('南海诸岛', y=0, pad=3, fontproperties=font)
# 颜色图例
cb = fig.colorbar(pm, ax=ax, location='bottom', pad=0.07)
cb.set_label('colorbar')
fig.savefig('matplotlib嵌套南海子图.png', bbox_inches='tight')
plt.close(fig)

两个方案绘制结果完全相同。

详解

可以看到新方案的两个明显优势：

1. 代码简洁很多。完全省下了计算南海图轴宽高和位置的工作量。
2. 南海图轴不必在最后绘制。这一点在画多行多列这样的图时非常方便。

这里详细解释一下新方案的核心ax.inset_axes：

ax_southsea = ax.inset_axes([0, 0, 1, 0.4], anchor=(1, 0), projection=platecarree)

最后一个参数projection自然是投影，不多解释。

这里第一个参数[0, 0, 1, 0.4]为新图轴ax_southsea在原图轴ax的相对位置，新创建的ax_southsea从ax左下角(0, 0)到右边边框从下往上4/10处(1, 0.4)，也就是占了大陆图轴ax的下面小半。

为什么这么大？不是应该仅仅占右下角一小块吗？

还记得上文简介中提到的“由于地图投影问题，最终图轴的宽高不一定是创建图轴时所指定的宽高”吗？这里利用了这一特性，在控制高度的时候，把宽度设得大一些，等绘图时宽度会自动缩小以匹配高度。

第二个参数anchor=(1, 0)用相对坐标表示缩放的锚点为ax_southsea的右下角，这保证了自动缩放的时候向右下角缩放。若anchor=(0, 0)，最后南海图轴将与大陆图轴在左下角对齐。

讨论

本文南海图轴用的是默认的四边形外框，有时可能想要左上角缺一块的五边形外框，那么可以把南海图轴背景设为透明，关闭外框，手动创建五边形进行clip（白化）后再叠加显示五边形。可惜这里空白太小，写不下具体实现方式。

该新方案可完全替换老方案，只是依赖的matplotlib版本需要高于3.6.x。