常用的雷达数据处理Python库中，PyCINRAD不直接计算CAPPI产品,PyCWR和PyART可以计算，但两者出来的效果有所不同，哪个是正确的？哪个更合理一些？这是值得关注的。
下面给出PyCWR和PyART生成3km高度CAPPI产品的代码和效果比对。
所用PyCWR库版本：1.0.8
所用PyCINRAD库版本：2.1.1
1.导入所需的库，并设置全局字体 import pyart from pycwr.io import read_auto from pycwr.qc import apply_dualpol_qc import matplotlib as mpl import matplotlib.pyplot as plt from cinrad.visualize.gpf import _cmap import numpy as np import cartopy.io.shapereader as shpreader import cartopy.crs as ccrs # 设置全局字体为宋体，用于绘图时支持中文标签
'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.data lat_ref = REF0.lat.data # 利用pyCWR生成CAPPI target_h = 3 lon_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() 效果图如下：
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 plt import pyart from pycwr.io import read_auto from cinrad.visualize.gpf import _cmap # 设置全局字体为SimSun，用于绘图时支持中文标签 plt.rcParams[ 'font.family' ] = 'SimSun' plt.rcParams[ 'font.size' ]=10 # 解决负号显示问题 plt.rcParams[ 'axes.unicode_minus' ] = False filename_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' field cappi = pyart.retrieve.create_cappi( radar, fields=[ "reflectivity" ], height=3000 ) # Create RadarMapDisplay objects for both PPI and CAPPI radar_display = pyart.graph.RadarMapDisplay(radar) cappi_display = pyart.graph.RadarMapDisplay(cappi) r_cmap = _cmap( "REF" )[ "cmap" ] # Plotting the PPI and CAPPI for comparison fig, ax = plt.subplots(1, 2, figsize=(15, 5)) # Plot PPI for 'reflectivity' field radar_display.plot_ppi( "reflectivity" , vmin=-10, vmax=60, ax=ax[0],cmap=r_cmap) ax[0].set_title( "PPI Reflectivity" ) # Plot CAPPI for 'reflectivity' field cappi_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 plots plt.show() 目前我们考虑：计算CAPPI产品，PyCWR的算法更合理。
