利用Metpy绘制雷达CAPPI图像

国内各位大神都开发了多个雷达数据解析绘图的模块，PyCinrad、Pycwr以及Pyart应该是很多同学使用频繁的模块，其中本人使用Pycinrad次数比较多，而且PyCinrad模块功能很强大，除了能够读取各种雷达基数据外，还能读取PUP和SWAN数据，但是目前只支持径向类型数据，具体读取程序见开发者的Github  https://github.com/CyanideCN/PyCINRAD/blob/master/README_zh.md

但是本人想做的分析是PUP的cappi产品文件，目前没能找到合适的模块对产品进行分析，很多大神建议我还是利用pycinrad对基数据进行处理并计算cappi，但是我实在是不愿意看文献找算法，无奈之下就搜寻了pycinrad的源文件，由于github中给出了读取PUP的语句是:

import cinrad.io import PUP

而最终在PyCINRAD-master/cinrad/io/level3.py中找到了PUP的类。

 根据这个PUP类，我有了个偷懒的法子，我决定利用Metpy模块对cappi产品文件进行读取。

首先分享metpy的两个指导方法手册：

https://github.com/Unidata/MetPy

MetPy — MetPy 1.3

用class PUP 和 metpy level3的例子（NEXRAD Level 3 File — MetPy 1.3）相结合，程序如下，文件读取为天擎下载的雷达单站PUP产品-等高面反射率bin文件：

from metpy.io.nexrad import Level3File
import numpy as np
from metpy.plots import add_metpy_logo, add_timestamp, colortables
import matplotlib.pyplot as plt
import matplotlib.colors as colors

path = "*********"
name = "Z_RADR_I_Z****_2022**********_P_DOR_SA_CAR_10_230_NUL.***.bin"

filepath = path + name
f = Level3File(filepath)
ctable = (('NWSStormClearReflectivity', -5, 20))
fig, ax = plt.subplots(1, 1, figsize=(15, 8))
with open(filepath, "rb") as buf:
    buf.seek(12)
    code = np.frombuffer(buf.read(2), ">i2")[0]
    cds = str(code).zfill(3)
    code = "Z*" + cds  # 站点号
    product_code = f.prod_desc.prod_code   # 产品编号
    data_block = f.sym_block[0][0]
    data = f.map_data(data_block['data'])
    
    az = np.array(data_block['start_az'] + [data_block['end_az'][-1]])
    rng = np.linspace(0, f.max_range, data.shape[-1] + 1)
    
    xlocs = rng * np.sin(np.deg2rad(az[:, np.newaxis]))
    ylocs = rng * np.cos(np.deg2rad(az[:, np.newaxis]))
    
   
    ax.pcolormesh(xlocs, ylocs, data, norm=norm, cmap=cmap)
    ax.set_aspect('equal', 'datalim')
    ax.set_xlim(-100, 100)
    ax.set_ylim(-100, 100)
    #add_timestamp(ax, f.metadata['prod_time'], y=0.02, high_contrast=True)
    
    plt.show()

程序比较简单，而且我还没有来得及坐标系转换，也没有搞colorbar，先写这么多，等有空再继续。

图片出来的效果：

 

 和我用PUP看的CAPPI形状还是挺像的。（由于资料尽可能不外传，仅截取部分）

 另外，稍微提一下metpy的其他用法，程序中 data_block = f.sym_block[0][0] 这个语句中，吸引我注意的是sym_block这一部分，在metpy中还有对这个the wide array of NEXRAD Level 3 (NIDS) 文件有其他的字节解码，具体：Level3File — MetPy 1.3，可以自己试着读取。