为雷达 PPI 图像添加地图底图

昨天开始改本科论文，发现了一个问题：如何给自己绘制的雷达 PPI 图添加上地图底图。

画地图底图的话肯定是用 Cartopy 或者是 cnmaps 来画，但是问题是雷达数据存放的形式是极坐标形式，而 Cartopy 或者 cnmaps 绘制地图使用的是经纬度坐标。将雷达数据由极坐标系转换为经纬度坐标系有一个小问题，就是数据量大。haversine 虽然可以根据两个地理位置的距离和角度计算经纬度，但是好像只能计算两个点之间的，一个一个点算的话太慢了。

后来我想到的解决办法是绘制两个图层，然后分别叠加到一起。matplotlib 的 Axes 类有一个 set_alpha 函数可以设置背景透明度，只需要让覆盖在上面的图层背景透明，显示出下面的图，就可以达到图层叠加的目的了。

第一步：绘制雷达数据

绘制雷达数据的库有很多，之前写本科论文时用的 pycwr，但是时隔半年之后竟然读取不了了。所以我转用 xarray 了。由于雷达数据并不像 ERA5 数据那样在 nc 文件里存储的那么整齐，我还多写了一个脚本来进行数据的转换，详细步骤这里略过，最终数据的结构如下。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

<xarray.Dataset>
Dimensions:                    (sweep: 9, azimuth: 360, range: 1000)
Coordinates:
  * sweep                      (sweep) int32 0 1 2 3 4 5 6 7 8
  * azimuth                    (azimuth) int64 0 1 2 3 4 ... 355 356 357 358 359
  * range                      (range) float32 250.0 500.0 ... 2.498e+05 2.5e+05
Data variables: (12/25)
    differential_phase         (sweep, azimuth, range) float32 ...
    time                       (sweep, azimuth) datetime64[ns] ...
    elevation                  (sweep, azimuth) float32 ...
    cross_correlation_ratio    (sweep, azimuth, range) float32 ...
    spectrum_width             (sweep, azimuth, range) float32 ...
    reflectivity               (sweep, azimuth, range) float32 ...
    ...                         ...
    longitude                  float64 ...
    altitude                   float64 ...
    time_coverage_start        |S32 ...
    time_coverage_end          |S32 ...
    time_reference             |S32 ...
    volume_number              int32 ...

第一步就是按照自己熟悉的方式绘制出普通的雷达图像，但是要确保拿到 fig 和 雷达图对应的 ax。

第二步：添加一个 GeoAxes

我们需要根据已有的 Axes 添加一个新的 GeoAxes (GeoAxes 是 Cartopy 中 Axes 的一个子类)。使用 get_position 获得雷达图像 ax 的位置，并使用 add_axes 函数添加一个相同大小的子图到相同的位置。

1
2
3
4
5
6
7
8

# get position
x0 = ax.get_position().x0
x1 = ax.get_position().x1
y0 = ax.get_position().y0
y1 = ax.get_position().y1
proj = ccrs.PlateCarree()
# add another axes
ax2: GeoAxes = fig.add_axes((x0, y0, x1 - x0, y1 - y0), projection=proj)

使用 set_zorder 函数设置图层上下的重叠次序，如果你想显示完整的地图边界，就把绘制地图的图层放置在雷达图层的上方。反之则相反。set_zorder 里数字小的图层位于下方。这里我把地图的图层放在了上方。为了显示出下方的图层，需要用 set_alpha 把地图图层的背景设置为透明的。

1
2
3

ax.set_zorder(0)
ax2.set_zorder(1)
ax2.patch.set_alpha(0)

为了不让显示的经纬度遮盖下面的距离坐标，使用 axis 关闭地图图层的坐标轴。如果要显示经纬度而不显示距离坐标，则关闭雷达图层的坐标轴即可。

1

ax2.axis('off')

第三步：绘制地图

绘制中国的地图需要注意精确性的问题，这里使用 cnmaps 来进行绘制，只需要传入添加好的 ax2 即可。绘制完成后根据雷达的观测范围计算出应该显示的经纬度范围，然后设置坐标轴显示范围即可。

1
2
3

draw_maps(get_adm_maps(level='省'), color='black', ax=ax2)
ax2.set_xlim(109.71952972227395, 115.17847027772608)
ax2.set_ylim(31.83252972227392, 37.29147027772605)

题外话：地图图层的实时修正

虽然我们已经达到了向 PPT 图添加地图底图的效果，但是如果有读者尝试改变 matplotlib 的图片展示窗口或者使用放大镜放大图片的话，就会发现在一些情况下地图的边界线会脱离图片的边缘。如果叠加上真正的地图边缘你就会发现，此时绘制的省界线是偏离了实际位置的，造成这种情况的原因是改变窗口的大小以及使用放大镜放大图片都会导致 Axes 和 GeoAxes 的位置发生变动，如果两个图层位置变动不准确，就会发生这样的情况。

以前的数据和代码找不到了，拿一个最近在画的卫星云图举例子吧

正常的显示，地图海岸线比较符合 (有差别是因为用了 Cartopy 自带的地图数据画的海岸线)

窗口全屏以后，海岸线图层小于卫星云图图层

用放大镜看山东半岛，海岸线图层也变小了

简单的做法是写一个回调函数或者回调类来实时同步图层的位置。

1
2
3
4
5
6
7
8
9
10
11

# 回调类
class EventHandler:
    def __init__(self, ax, ax2):
        self.ax = ax
        self.ax2 = ax2
        
    def __call__(self, *args, **kwargs):
        """
        实现 __call__ 方法让回调类能够像函数一样调用
        """
        self.ax2.set_position(self.ax.get_position())

然后分别设置监听窗口大小改变的事件和坐标轴范围改变的事件就可以了。

1
2
3

fig.canvas.mpl_connect('resize_event', EventHandler(ax, ax2))
ax2.callbacks.connect('xlim_changed', EventHandler(ax, ax2))
ax2.callbacks.connect('ylim_changed', EventHandler(ax, ax2))

注意，这里把回调类和地图图层绑定到了一起。因为我发现绑定到雷达 (卫星) 图层的话是没有效果的