常用的雷达数据处理Python库中，PyCINRAD不直接计算CAPPI产品,PyCWR和PyART可以计算，但两者出来的效果有所不同，哪个是正确的？哪个更合理一些？这是值得关注的。
下面给出PyCWR和PyART生成3km高度CAPPI产品的代码和效果比对。
所用PyCWR库版本：1.0.8
所用PyCINRAD库版本：2.1.1
1.导入所需的库，并设置全局字体
import pyartfrom pycwr.io import read_autofrom pycwr.qc import apply_dualpol_qcimport matplotlib as mplimport matplotlib.pyplot as pltfrom cinrad.visualize.gpf import _cmapimport numpy as npimport cartopy.io.shapereader as shpreaderimport cartopy.crs as ccrs# 设置全局字体为宋体，用于绘图时支持中文标签mpl.rcParams['font.family'] = 'SimSun'mpl.rcParams['font.size']=15# 解决负号显示问题mpl.rcParams['axes.unicode_minus'] = False
2.自定义利用pycwr计算S波段雷达CAPPI的函数
def pycwr_cappi_sband(filename_S, target_height_km=3): try:# 读取基数据 radar=read_auto(filename_S)# 双偏振雷达数据质量控制 qc_radar = apply_dualpol_qc(radar, inplace=False, clear_air_mode="mask",band='S')# 设置产品生成的经纬度范围 lon1d = np.arange(site_lon-0.01*230, site_lon+0.01*230, 0.01) lat1d = np.arange(site_lat-0.01*230, site_lat+0.01*230, 0.01)# 生成CAPPI产品 qc_radar.add_product_CAPPI_lonlat(XLon=lon1d, YLat=lat1d, level_height=target_height_km*1000)# 提取CAPPI产品 ds=qc_radar.product# 提取CAPPI产品经度 lon = ds[f'lon_cappi_{target_height_km*1000}'].values# 提取CAPPI产品纬度 lat = ds[f'lat_cappi_{target_height_km*1000}'].values# 提取 CAPPI 数据值 cappi = ds[f'CAPPI_geo_{target_height_km*1000}_native'].values# 注意这里有转置，因为生成的CAPPI坐标是（纬度，经度），但我需要的是（经度，纬度） ref_data=cappi.T# 生成经纬度网格，方便绘图 lon_grid,lat_grid=np.meshgrid(lon,lat)return lon_grid, lat_grid, ref_data except Exception as e:print(f"PyCWR处理S波段数据失败: {e}")return None, None, None, None
3.自定义利用pyart计算S波段雷达CAPPI的函数
def pyart_cappi_sband(filename_S, target_height_km=3): try:# 读取基数据并转换成pyart适用的数据格式 f=read_auto(filename_S) radar=f.to_pyart_radar()# 基数据质量控制 gatefilter = pyart.filters.GateFilter(radar) gatefilter.exclude_transition()# 生成反射率因子target_height_km高度CAPPI产品 grid = pyart.map.grid_from_radars( (radar,), grid_shape=(1, 500, 500), # 1层高度，500x500的水平网格 grid_limits=( (target_height_km*1000, target_height_km*1000), # 高度限制在 2000 米 (-230000.0, 230000.0), # 纬度范围 (雷达中心向南向北各230公里) (-230000.0, 230000.0) # 经度范围 (雷达中心向西向东各230公里) ), fields=['reflectivity'] # 需要处理的变量 ) # 获取经纬度和数据# 注意：此时的数据已经是笛卡尔网格，形状为 (500, 500) cappi_data_grid = grid.fields['reflectivity']['data'][0] # [0]代表取第0层(2000m) grid_lons = grid.point_longitude['data'][0] grid_lats = grid.point_latitude['data'][0]return grid_lons, grid_lats, cappi_data_grid except Exception as e:print(f"PyART处理S波段数据失败: {e}")return None, None, None, None
4.主程序：读取基数据，用两种方式生成CAPPI产品，并绘图与PPI对比
path=('D:/Z_RADR_I_Z*_O_DOR_SAD_CAP_FMT.bin.bz2')radar = read_auto(path)# 获取雷达站点经纬度site_lon=radar.scan_info.longitude.data.item()site_lat=radar.scan_info.latitude.data.item()# 提取指定仰角序号的变量场数据REF0 = radar.get_sweep_field(0, "dBZ")lon_ref = REF0.lon.datalat_ref = REF0.lat.data# 利用pyCWR生成CAPPItarget_h = 3lon_c,lat_c,cappi_cwr=pycwr_cappi_sband(path, target_height_km=target_h)lon_a,lat_a,cappi_art=pyart_cappi_sband(path, target_height_km=target_h)# 准备绘图参数# 采用pycinrad的画图色标r_cmap = _cmap("REF")["cmap"]# 准备地图文件border_file=r'C:/x.shp'river_file = r'C:/river/river_HN.shp'# 提前读取好shp文件hunan_shp, river_shp = [], []try: hunan_shp = list(shpreader.Reader(border_file).geometries()) river_shp = list(shpreader.Reader(river_file).geometries())except Exception as e:print(f"提示：底图文件未加载。原因: {e}")# 创建 1x3 画布fig, axes = plt.subplots(1, 3, figsize=(24, 8), subplot_kw={'projection': ccrs.PlateCarree()})plt.subplots_adjust(wspace=0.15)# 为每个子图添加底图和经纬度范围限制for ax in axes:if hunan_shp: ax.add_geometries(hunan_shp, ccrs.PlateCarree(), edgecolor='black', facecolor='none', linewidth=0.8)if river_shp: ax.add_geometries(river_shp, ccrs.PlateCarree(), edgecolor='blue', facecolor='none', linewidth=0.5, alpha=0.5)# center_lon,center_lat=112.61, 28.34 ax.set_extent([site_lon-2.3, site_lon+2.3, site_lat-2.3, site_lat+2.3], crs=ccrs.PlateCarree())# 第1幅图：PPI图像axes[0].set_title("PPI", fontsize=15)cf1 = axes[0].pcolormesh(lon_ref, lat_ref, REF0, cmap=r_cmap, vmin=0, vmax=75, transform=ccrs.PlateCarree())cbar1 = plt.colorbar(cf1, ax=axes[0], shrink=0.8, pad=0.02)cbar1.set_label("Reflectivity (dBZ)", fontsize=12)# 第2幅图：pycwr计算出的CAPPI图像axes[1].set_title(f"pycwr {target_h}km CAPPI", fontsize=15)cf1 = axes[1].pcolormesh(lon_c,lat_c,cappi_cwr, cmap=r_cmap, vmin=0, vmax=75, transform=ccrs.PlateCarree())cbar1 = plt.colorbar(cf1, ax=axes[1], shrink=0.8, pad=0.02)cbar1.set_label(f"Reflectivity at {target_h}km (dBZ)", fontsize=12)# 第3幅图：pyart计算出的CAPPI图像axes[2].set_title(f"pyart {target_h}km CAPPI", fontsize=15)cf2 = axes[2].pcolormesh(lon_a,lat_a,cappi_art, cmap=r_cmap, vmin=0, vmax=75, transform=ccrs.PlateCarree())# 设置色标cbar2 = plt.colorbar(cf2, ax=axes[2], shrink=0.8, pad=0.02)cbar2.set_label(f"Reflectivity at {target_h}km (dBZ)", fontsize=12)plt.show()
效果图如下：
可以对比下上面的出图，相比最左边的0.5°仰角的PPI图像，可以看到两个库生成的CAPPI产品的不同。
PyCWR生成的3kmCAPPI产品更合理，PyART生成的3kmCAPPI产品图像中，右上角位置的那团强回波，其0.5°仰角高度已经达到6km了。PyART生成的CAPPI产品截取了≥3km的高度数据。
同时，我们还关注了两点：
1.关注的第1个问题：
软件ROSE出来的CAPPI产品是哪一种形式？
答案：它与PyART产品一致。（我这里用的是ROSE2.1版本，不知道3.0版本有没有改变）
下方左图是0.5°仰角PPI图像，右图是3kmCAPPI图像。右上角强回波的高度在ROSE中显示已经超过5.5km。但依然出现在了3km高度的CAPPI图像上。
2.关注的第2个问题；
PyART本身有一个反演CAPPI的函数：pyart.retrieve.create_cappi，这个函数是否好用呢？
答案：不好用。
请看用同样的基数据生成的PPI图像和该函数生成的CAPPI图像对比。
这个函数仿佛只计算了单个仰角3km高度的CAPPI值，而且位置进行了旋转。
下面是参照原PyART文档里给出的示例编写的代码，大家感兴趣的话可以试验一下自己的基数据文件，现在官方推荐都不用这个函数了，官方推荐用我们上面给出的pyart.map.grid_from_radars函数。
import matplotlib.pyplot as pltimport pyartfrom pycwr.io import read_autofrom cinrad.visualize.gpf import _cmap# 设置全局字体为SimSun，用于绘图时支持中文标签plt.rcParams['font.family'] = 'SimSun'plt.rcParams['font.size']=10# 解决负号显示问题plt.rcParams['axes.unicode_minus'] = Falsefilename_S="D:/Z_RADR_I_Z9731_*_O_DOR_SAD_CAP_FMT.bin.bz2"f=read_auto(filename_S)radar=f.to_pyart_radar()gatefilter = pyart.filters.GateFilter(radar)gatefilter.exclude_transition()# Create CAPPI at 3,000 meters for the 'reflectivity' fieldcappi = pyart.retrieve.create_cappi( radar, fields=["reflectivity"], height=3000)# Create RadarMapDisplay objects for both PPI and CAPPIradar_display = pyart.graph.RadarMapDisplay(radar)cappi_display = pyart.graph.RadarMapDisplay(cappi)r_cmap = _cmap("REF")["cmap"]# Plotting the PPI and CAPPI for comparisonfig, ax = plt.subplots(1, 2, figsize=(15, 5))# Plot PPI for 'reflectivity' fieldradar_display.plot_ppi("reflectivity", vmin=-10, vmax=60, ax=ax[0],cmap=r_cmap)ax[0].set_title("PPI Reflectivity")# Plot CAPPI for 'reflectivity' fieldcappi_display.plot_ppi("reflectivity", vmin=-10, vmax=60, ax=ax[1],cmap=r_cmap)ax[1].set_title("CAPPI Reflectivity at 3000 meters")# Show the plotsplt.show()
目前我们考虑：计算CAPPI产品，PyCWR的算法更合理。
