This article was last updated on <span id="expire-date"></span> days ago, the information described in the article may be outdated.
昨天开始改本科论文，发现了一个问题：如何给自己绘制的雷达 PPI 图添加上地图底图。
画地图底图的话肯定是用 Cartopy 或者是 cnmaps 来画，但是问题是雷达数据存放的形式是极坐标形式，而 Cartopy 或者 cnmaps 绘制地图使用的是经纬度坐标。将雷达数据由极坐标系转换为经纬度坐标系有一个小问题，就是数据量大。haversine 虽然可以根据两个地理位置的距离和角度计算经纬度，但是好像只能计算两个点之间的，一个一个点算的话太慢了。
后来我想到的解决办法是绘制两个图层，然后分别叠加到一起。matplotlib 的 Axes 类有一个 set_alpha 函数可以设置背景透明度，只需要让覆盖在上面的图层背景透明，显示出下面的图，就可以达到图层叠加的目的了。
第一步：绘制雷达数据
绘制雷达数据的库有很多，之前写本科论文时用的 pycwr，但是时隔半年之后竟然读取不了了。所以我转用 xarray 了。由于雷达数据并不像 ERA5 数据那样在 nc 文件里存储的那么整齐，我还多写了一个脚本来进行数据的转换，详细步骤这里略过，最终数据的结构如下。
1 | <xarray.Dataset> |
第一步就是按照自己熟悉的方式绘制出普通的雷达图像，但是要确保拿到 fig 和 雷达图对应的 ax。
第二步：添加一个 GeoAxes
我们需要根据已有的 Axes 添加一个新的 GeoAxes (GeoAxes 是 Cartopy 中 Axes 的一个子类)。使用 get_position 获得雷达图像 ax 的位置，并使用 add_axes 函数添加一个相同大小的子图到相同的位置。
1 | # get position |
使用 set_zorder 函数设置图层上下的重叠次序，如果你想显示完整的地图边界，就把绘制地图的图层放置在雷达图层的上方。反之则相反。set_zorder 里数字小的图层位于下方。这里我把地图的图层放在了上方。为了显示出下方的图层，需要用 set_alpha 把地图图层的背景设置为透明的。
1 | ax.set_zorder(0) |
为了不让显示的经纬度遮盖下面的距离坐标，使用 axis 关闭地图图层的坐标轴。如果要显示经纬度而不显示距离坐标，则关闭雷达图层的坐标轴即可。
1 | ax2.axis('off') |
第三步：绘制地图
绘制中国的地图需要注意精确性的问题，这里使用 cnmaps 来进行绘制，只需要传入添加好的 ax2 即可。绘制完成后根据雷达的观测范围计算出应该显示的经纬度范围，然后设置坐标轴显示范围即可。
1 | draw_maps(get_adm_maps(level='省'), color='black', ax=ax2) |
Author: Syize
Permalink: https://blog.syize.cn/2023/01/13/plot-map-on-radar/
本博客所有文章除特别声明外，均采用 CC BY-NC-SA 4.0 许可协议。转载请注明来自 Syizeのblog！
Comments