超长篇幅！PyCINRAD保姆级教程

文摘 2024-08-11 08:00 北京
前言
近日又有几位读者询问雷达如何处理数据、绘图等等问题
由于pysoer大佬（湖南省益阳市气象局,高级工程师）已经写了非常详细的《PyCINRAD保姆级教程》， 小编就不班门弄斧了，直接私信大佬将文章转载到这里
开始学习吧
简介
PyCINRAD 是一个用于处理和可视化中国新一代天气雷达（CINRAD）数据的开源Python库。它支持多种雷达数据格式的读取，并提供了一些实用的算法和可视化工具。这个库使得用户可以方便地对雷达基数据进行处理和展示，广泛应用于气象学领域。 地址：https://github.com/CyanideCN/PyCINRAD
安装&升级
目前cinrad版本为v1.9.1 由于高版本的cartopy性能反而下降，建议安装0.22版本 numpy 2.0会有冲突，建议安装1.x版本
In [ ]:
pip install numpy==1.26pip install cartopy==0.22  pip install cinrad -U
In [ ]:# 有很多朋友Windows下安装cartopy失败，建议安装anaconda来处理。conda install cartopy -c conda-forge# 同时需要安装 Visual studio 生成工具;# `https://blog.csdn.net/qq_43522889/article/details/127649550`# 参考这个方法一，把右边那些都勾上。windows 10 SDK
如果要以下新功能，则需要到github下载更新（请先安装git）-> 组合反射率拼图:在mergeCR分支
In [ ]:
pip install git+https://github.com/pysoer/PyCINRAD.git@mergeCR -U --force-reinstall --no-deps

支持的文件类型
1. 标准格式基数据StandardData
标准格式基数据: 标准格式文件名中有“FMT”、基数据是O_DOR Z_RADR_I_Z9999_20231025220414_O_DOR_SA_CAP_FMT.bin
2. 标准格式产品StandardPUP
标准格式产品: 标准格式文件名中有“FMT”、产品数据是P_DOR 这里文件名中有HCL，代表是HCL水凝物分类产品 Z_RADR_I_Z9999_20231231160439_P_DOR_SAD_HCL_250_230_5_FMT.bin 下面这种老格式的也是标准格式产品 Z9737_20231126150649Z_CR_00_37
3. 探测中心天气雷达拼图系统v3产品.MocMosaic
ACHN表示是全网拼图，CREF表示是产品类型是CREF组合反射率 Z_RADA_C_BABJ_20231212010615_P_DOR_ACHN_CREF_20231212_010000
4. SWAN
以前的老SWAN3产品.MCR表示拼图组合反射率 Z_OTHE_RADAMCR_20220525085400.bin
5. 老格式的国产基数据.CinradReader
和标准格式的唯一区别是文件名中没有“FMT”（这里指的是大部分新情况） Z_RADR_I_Z9731_20240511070524_O_DOR_SAD_CAP.bin
6. 相控阵PhasedArrayData
标准格式的相控阵雷达基数据(2023以后). 文件名一般有AXPT，就是纳睿雷达的产品。Z_RADR_I_ZBJ02_20210815155836_O_DOR_DXK_CAR
Z_RADR_I_ZA601_20240415183600_O_DOR_AXPT0364_CRA_FMT
Z_RADR_I_ZS999_20231220000000_O_DOR_AXPT0364_CRA_.bz2
7: PUP
NEXRAD Level 3 (NIDS) product files. 已经被淘汰的格式，一言难尽，格式太乱，能用则用吧。
万能读取read_auto
如果你搞不清楚，那请用cinrad.io.read_auto(your_radar_file)来读取。
总结，有以下接口:
cinrad.io.StandardData
cinrad.io.StandardPUP
cinrad.io.MocMosaic
cinrad.io.SWAN(filename,product="CR")
cinrad.io.CinradReader(filename,radar_type="CC")
cinrad.io.PhasedArrayData
cinrad.io.PUP
cinrad.io.read_auto

开始测试
In [2]:
import warningsimport cinradfrom cinrad.visualize import PPIimport numpy as npimport matplotlib.pyplot as pltimport datetime%matplotlib inlinewarnings.filterwarnings("ignore")
basePath = "d:/temp/"

cinrad.__version__
'1.9.1'
标准格式基数据
In [11]:
nFiles = basePath + "/cinrad/bz2/Z_RADR_I_Z9532_20200517124300_O_DOR_SAD_CAP_FMT.bin.bz2"f= cinrad.io.read_auto(nFiles)data = f.get_data(0,230,"REF") # 读取第一层的反射率data

In [7]:
f.available_product(0) #第0个仰角有哪些产品可以读取

['TREF', 'REF', 'SQI', 'ZDR', 'RHO', 'PHI', 'KDP', 'SNRH']
In [40]:
f.available_tilt('REF') #REF产品有哪些仰角可以读取

[0, 2, 4, 5, 6, 7, 8, 9, 10]
In [5]:
f.available_tilt('VEL') #VEL产品有哪些仰角可以读取

[1, 3, 4, 5, 6, 7, 8, 9, 10]
In [3]:
f.el # 雷达的仰角

[0.48339844,
 0.48339844,
 1.4941406,
 1.4941406,
 2.4169922,
 3.2958984,
 4.3066406,
 6.020508,
 9.887695,
 14.589844,
 19.511719]
In [5]:
f.scan_config

[ScanConfig(process_mode=1, wave_form=0, PRF1=322.0, PRF2=322.0, dealias_mode=1, azimuth=0.0, elev=0.48339844, start_angle=0.0, end_angle=0.0, angular_reso=1.0, scan_spd=11.34, log_reso=250, dop_reso=250, max_range1=460000, max_range2=460000, start_range=0, sample1=28, sample2=28, phase_mode=0, atmos_loss=0.011, nyquist_spd=8.37801, moments_mask=69286, moments_size_mask=1152, misc_filter_mask=255, SQI_thres=0.4, SIG_thres=5.0, CSR_thres=60.0, LOG_thres=3.0, CPA_thres=0.0, PMI_thres=0.45, DPLOG_thres=5.0, res_thres=void(b'\x00\x00\x00\x00'), dBT_mask=1, dBZ_mask=1, vel_mask=1, sw_mask=1, DP_mask=32, res_mask=void(b'\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00'), scan_sync=0, direction=1, ground_clutter_classifier_type=3, ground_clutter_filter_type=1, ground_clutter_filter_notch_width=3, ground_clutter_filter_window=1, res4=void(b'\x01\x03\x03\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00')),
 ScanConfig(process_mode=1, wave_form=1, PRF1=1014.0, PRF2=1014.0, dealias_mode=1, azimuth=0.0, elev=0.48339844, start_angle=0.0, end_angle=0.0, angular_reso=1.0, scan_spd=15.84, log_reso=250, dop_reso=250, max_range1=300000, max_range2=300000, start_range=0, sample1=64, sample2=64, phase_mode=2, atmos_loss=0.011, nyquist_spd=26.382925, moments_mask=67108888, moments_size_mask=0, misc_filter_mask=63, SQI_thres=0.4, SIG_thres=5.0, CSR_thres=60.0, LOG_thres=3.0, CPA_thres=0.0, PMI_thres=0.45, DPLOG_thres=5.0, res_thres=void(b'\x00\x00\x00\x00'), dBT_mask=1, dBZ_mask=1, vel_mask=1, sw_mask=1, DP_mask=32, res_mask=void(b'\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00'), scan_sync=0, direction=0, ground_clutter_classifier_type=3, ground_clutter_filter_type=1, ground_clutter_filter_notch_width=3, ground_clutter_filter_window=1, res4=void(b'\x01\x03\x03\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00')),
 ScanConfig(process_mode=1, wave_form=0, PRF1=322.0, PRF2=322.0, dealias_mode=1, azimuth=0.0, elev=1.4941406, start_angle=0.0, end_angle=0.0, angular_reso=1.0, scan_spd=11.34, log_reso=250, dop_reso=250, max_range1=460000, max_range2=460000, start_range=0, sample1=28, sample2=28, phase_mode=0, atmos_loss=0.011, nyquist_spd=8.37801, moments_mask=69254, moments_size_mask=1152, misc_filter_mask=127, SQI_thres=0.4, SIG_thres=5.0, CSR_thres=60.0, LOG_thres=3.0, CPA_thres=0.0, PMI_thres=0.45, DPLOG_thres=5.0, res_thres=void(b'\x00\x00\x00\x00'), dBT_mask=1, dBZ_mask=1, vel_mask=1, sw_mask=1, DP_mask=32, res_mask=void(b'\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00'), scan_sync=0, direction=0, ground_clutter_classifier_type=3, ground_clutter_filter_type=1, ground_clutter_filter_notch_width=3, ground_clutter_filter_window=1, res4=void(b'\x01\x03\x03\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00')),
 ScanConfig(process_mode=1, wave_form=1, PRF1=1014.0, PRF2=1014.0, dealias_mode=1, azimuth=0.0, elev=1.4941406, start_angle=0.0, end_angle=0.0, angular_reso=1.0, scan_spd=15.84, log_reso=250, dop_reso=250, max_range1=300000, max_range2=300000, start_range=0, sample1=64, sample2=64, phase_mode=2, atmos_loss=0.011, nyquist_spd=26.382925, moments_mask=67108888, moments_size_mask=0, misc_filter_mask=63, SQI_thres=0.4, SIG_thres=5.0, CSR_thres=60.0, LOG_thres=3.0, CPA_thres=0.0, PMI_thres=0.45, DPLOG_thres=5.0, res_thres=void(b'\x00\x00\x00\x00'), dBT_mask=1, dBZ_mask=1, vel_mask=1, sw_mask=1, DP_mask=32, res_mask=void(b'\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00'), scan_sync=0, direction=0, ground_clutter_classifier_type=3, ground_clutter_filter_type=1, ground_clutter_filter_notch_width=3, ground_clutter_filter_window=1, res4=void(b'\x01\x03\x03\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00')),
 ScanConfig(process_mode=1, wave_form=4, PRF1=1014.0, PRF2=446.0, dealias_mode=1, azimuth=0.0, elev=2.4169922, start_angle=0.0, end_angle=0.0, angular_reso=1.0, scan_spd=11.18, log_reso=250, dop_reso=250, max_range1=146000, max_range2=330000, start_range=0, sample1=72, sample2=8, phase_mode=0, atmos_loss=0.011, nyquist_spd=26.382925, moments_mask=69278, moments_size_mask=1152, misc_filter_mask=63, SQI_thres=0.4, SIG_thres=5.0, CSR_thres=60.0, LOG_thres=3.0, CPA_thres=0.0, PMI_thres=0.45, DPLOG_thres=5.0, res_thres=void(b'\x00\x00\x00\x00'), dBT_mask=1, dBZ_mask=1, vel_mask=1, sw_mask=1, DP_mask=32, res_mask=void(b'\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00'), scan_sync=0, direction=0, ground_clutter_classifier_type=3, ground_clutter_filter_type=1, ground_clutter_filter_notch_width=3, ground_clutter_filter_window=1, res4=void(b'\x03\x03\x03\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00')),
 ScanConfig(process_mode=1, wave_form=4, PRF1=1014.0, PRF2=446.0, dealias_mode=1, azimuth=0.0, elev=3.2958984, start_angle=0.0, end_angle=0.0, angular_reso=1.0, scan_spd=11.18, log_reso=250, dop_reso=250, max_range1=146000, max_range2=330000, start_range=0, sample1=72, sample2=8, phase_mode=0, atmos_loss=0.011, nyquist_spd=26.382925, moments_mask=69278, moments_size_mask=1152, misc_filter_mask=63, SQI_thres=0.4, SIG_thres=5.0, CSR_thres=60.0, LOG_thres=3.0, CPA_thres=0.0, PMI_thres=0.45, DPLOG_thres=5.0, res_thres=void(b'\x00\x00\x00\x00'), dBT_mask=1, dBZ_mask=1, vel_mask=1, sw_mask=1, DP_mask=32, res_mask=void(b'\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00'), scan_sync=0, direction=0, ground_clutter_classifier_type=3, ground_clutter_filter_type=1, ground_clutter_filter_notch_width=3, ground_clutter_filter_window=1, res4=void(b'\x03\x03\x03\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00')),
 ScanConfig(process_mode=1, wave_form=4, PRF1=1014.0, PRF2=446.0, dealias_mode=1, azimuth=0.0, elev=4.3066406, start_angle=0.0, end_angle=0.0, angular_reso=1.0, scan_spd=11.18, log_reso=250, dop_reso=250, max_range1=146000, max_range2=330000, start_range=0, sample1=72, sample2=8, phase_mode=0, atmos_loss=0.011, nyquist_spd=26.382925, moments_mask=69278, moments_size_mask=1152, misc_filter_mask=63, SQI_thres=0.4, SIG_thres=5.0, CSR_thres=60.0, LOG_thres=3.0, CPA_thres=0.0, PMI_thres=0.45, DPLOG_thres=5.0, res_thres=void(b'\x00\x00\x00\x00'), dBT_mask=1, dBZ_mask=1, vel_mask=1, sw_mask=1, DP_mask=32, res_mask=void(b'\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00'), scan_sync=0, direction=0, ground_clutter_classifier_type=3, ground_clutter_filter_type=1, ground_clutter_filter_notch_width=3, ground_clutter_filter_window=1, res4=void(b'\x03\x03\x03\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00')),
 ScanConfig(process_mode=1, wave_form=4, PRF1=1014.0, PRF2=644.0, dealias_mode=1, azimuth=0.0, elev=6.020508, start_angle=0.0, end_angle=0.0, angular_reso=1.0, scan_spd=11.18, log_reso=250, dop_reso=250, max_range1=146000, max_range2=228000, start_range=0, sample1=78, sample2=8, phase_mode=0, atmos_loss=0.011, nyquist_spd=26.382925, moments_mask=69278, moments_size_mask=1152, misc_filter_mask=63, SQI_thres=0.4, SIG_thres=5.0, CSR_thres=60.0, LOG_thres=3.0, CPA_thres=0.0, PMI_thres=0.45, DPLOG_thres=5.0, res_thres=void(b'\x00\x00\x00\x00'), dBT_mask=1, dBZ_mask=1, vel_mask=1, sw_mask=1, DP_mask=32, res_mask=void(b'\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00'), scan_sync=0, direction=0, ground_clutter_classifier_type=3, ground_clutter_filter_type=1, ground_clutter_filter_notch_width=3, ground_clutter_filter_window=1, res4=void(b'\x01\x03\x03\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00')),
 ScanConfig(process_mode=1, wave_form=2, PRF1=1181.0, PRF2=1181.0, dealias_mode=1, azimuth=0.0, elev=9.887695, start_angle=0.0, end_angle=0.0, angular_reso=1.0, scan_spd=18.0, log_reso=250, dop_reso=250, max_range1=124000, max_range2=124000, start_range=0, sample1=65, sample2=65, phase_mode=0, atmos_loss=0.011, nyquist_spd=30.728043, moments_mask=69278, moments_size_mask=1152, misc_filter_mask=63, SQI_thres=0.4, SIG_thres=5.0, CSR_thres=60.0, LOG_thres=3.0, CPA_thres=0.0, PMI_thres=0.45, DPLOG_thres=5.0, res_thres=void(b'\x00\x00\x00\x00'), dBT_mask=1, dBZ_mask=1, vel_mask=1, sw_mask=1, DP_mask=32, res_mask=void(b'\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00'), scan_sync=0, direction=0, ground_clutter_classifier_type=3, ground_clutter_filter_type=1, ground_clutter_filter_notch_width=3, ground_clutter_filter_window=1, res4=void(b'\x03\x03\x03\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00')),
 ScanConfig(process_mode=1, wave_form=2, PRF1=1181.0, PRF2=1181.0, dealias_mode=1, azimuth=0.0, elev=14.589844, start_angle=0.0, end_angle=0.0, angular_reso=1.0, scan_spd=18.0, log_reso=250, dop_reso=250, max_range1=124000, max_range2=124000, start_range=0, sample1=65, sample2=65, phase_mode=0, atmos_loss=0.011, nyquist_spd=30.728043, moments_mask=69278, moments_size_mask=1152, misc_filter_mask=63, SQI_thres=0.4, SIG_thres=5.0, CSR_thres=60.0, LOG_thres=3.0, CPA_thres=0.0, PMI_thres=0.45, DPLOG_thres=5.0, res_thres=void(b'\x00\x00\x00\x00'), dBT_mask=1, dBZ_mask=1, vel_mask=1, sw_mask=1, DP_mask=32, res_mask=void(b'\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00'), scan_sync=0, direction=0, ground_clutter_classifier_type=3, ground_clutter_filter_type=1, ground_clutter_filter_notch_width=3, ground_clutter_filter_window=1, res4=void(b'\x03\x03\x03\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00')),
 ScanConfig(process_mode=1, wave_form=2, PRF1=1181.0, PRF2=1181.0, dealias_mode=1, azimuth=0.0, elev=19.511719, start_angle=0.0, end_angle=0.0, angular_reso=1.0, scan_spd=18.0, log_reso=250, dop_reso=250, max_range1=124000, max_range2=124000, start_range=0, sample1=65, sample2=65, phase_mode=0, atmos_loss=0.011, nyquist_spd=30.728043, moments_mask=69278, moments_size_mask=1152, misc_filter_mask=63, SQI_thres=0.4, SIG_thres=5.0, CSR_thres=60.0, LOG_thres=3.0, CPA_thres=0.0, PMI_thres=0.45, DPLOG_thres=5.0, res_thres=void(b'\x00\x00\x00\x00'), dBT_mask=1, dBZ_mask=1, vel_mask=1, sw_mask=1, DP_mask=32, res_mask=void(b'\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00'), scan_sync=0, direction=0, ground_clutter_classifier_type=3, ground_clutter_filter_type=1, ground_clutter_filter_notch_width=3, ground_clutter_filter_window=1, res4=void(b'\x03\x03\x03\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00'))]
In [7]:
f.scantime

datetime.datetime(2023, 4, 2, 7, 14, 19, tzinfo=datetime.timezone.utc)
In [6]:
vel0 = f.get_data(1,230,"VEL") # 速度只有1,3,4....10vel0

相关系数CC
In [12]:
cc = f.get_data(0,230,"RHO") # KDP/ZDR/RHOfig = cinrad.visualize.PPI(cc,style="black",add_city_names=True)

# 将图片保存
# fig("d:/")
# 或者 fig("d:/abc.png")
# 或者 imgName = fig("d:/")
[<matplotlib.lines.Line2D object at 0x000001E90328D9D0>, <matplotlib.lines.Line2D object at 0x000001E90328DA00>, <matplotlib.lines.Line2D object at 0x000001E90328DB20>, <matplotlib.lines.Line2D object at 0x000001E90328DC10>, <matplotlib.lines.Line2D object at 0x000001E90328DD00>, <matplotlib.lines.Line2D object at 0x000001E90328DDF0>, <matplotlib.lines.Line2D object at 0x000001E90328DEE0>, <matplotlib.lines.Line2D object at 0x000001E90328DFD0>, <matplotlib.lines.Line2D object at 0x000001E903356100>, <matplotlib.lines.Line2D object at 0x000001E9033561F0>, <matplotlib.lines.Line2D object at 0x000001E9033562E0>, <matplotlib.lines.Line2D object at 0x000001E9033563D0>, <matplotlib.lines.Line2D object at 0x000001E9033564C0>, <matplotlib.lines.Line2D object at 0x000001E9033565B0>, <matplotlib.lines.Line2D object at 0x000001E9033566A0>, <matplotlib.lines.Line2D object at 0x000001E903356790>, <matplotlib.lines.Line2D object at 0x000001E903356880>, <matplotlib.lines.Line2D object at 0x000001E903356970>, <matplotlib.lines.Line2D object at 0x000001E903356A60>, <matplotlib.lines.Line2D object at 0x000001E903356B50>, <matplotlib.lines.Line2D object at 0x000001E903356C40>, <matplotlib.lines.Line2D object at 0x000001E903356D30>, <matplotlib.lines.Line2D object at 0x000001E903356E20>, <matplotlib.lines.Line2D object at 0x000001E903356F10>, <matplotlib.lines.Line2D object at 0x000001E90335E040>, <matplotlib.lines.Line2D object at 0x000001E90335E130>, <matplotlib.lines.Line2D object at 0x000001E90335E220>, <matplotlib.lines.Line2D object at 0x000001E90335E310>, <matplotlib.lines.Line2D object at 0x000001E90335E400>, <matplotlib.lines.Line2D object at 0x000001E90335E4F0>, <matplotlib.lines.Line2D object at 0x000001E90335E5E0>, <matplotlib.lines.Line2D object at 0x000001E90335E6D0>, <matplotlib.lines.Line2D object at 0x000001E90335E7C0>, <matplotlib.lines.Line2D object at 0x000001E90335E8B0>, <matplotlib.lines.Line2D object at 0x000001E90335E9A0>, <matplotlib.lines.Line2D object at 0x000001E90335EA90>, <matplotlib.lines.Line2D object at 0x000001E90335EB80>, <matplotlib.lines.Line2D object at 0x000001E90335EC70>, <matplotlib.lines.Line2D object at 0x000001E90335ED60>, <matplotlib.lines.Line2D object at 0x000001E90335EE50>, <matplotlib.lines.Line2D object at 0x000001E90335EF40>, <matplotlib.lines.Line2D object at 0x000001E903369070>, <matplotlib.lines.Line2D object at 0x000001E903369160>, <matplotlib.lines.Line2D object at 0x000001E903369250>, <matplotlib.lines.Line2D object at 0x000001E903369340>, <matplotlib.lines.Line2D object at 0x000001E903369430>, <matplotlib.lines.Line2D object at 0x000001E903369520>, <matplotlib.lines.Line2D object at 0x000001E903369610>, <matplotlib.lines.Line2D object at 0x000001E903369700>, <matplotlib.lines.Line2D object at 0x000001E9033697F0>, <matplotlib.lines.Line2D object at 0x000001E9033698E0>, <matplotlib.lines.Line2D object at 0x000001E9033699D0>, <matplotlib.lines.Line2D object at 0x000001E903369AC0>, <matplotlib.lines.Line2D object at 0x000001E903369BB0>, <matplotlib.lines.Line2D object at 0x000001E903369CA0>, <matplotlib.lines.Line2D object at 0x000001E903369D90>, <matplotlib.lines.Line2D object at 0x000001E903369EE0>, <matplotlib.lines.Line2D object at 0x000001E903369FD0>, <matplotlib.lines.Line2D object at 0x000001E903372100>, <matplotlib.lines.Line2D object at 0x000001E9033721F0>]
计算CR
In [4]:
r_list = list(f.iter_tilt(230, 'REF'))cr = cinrad.calc.quick_cr(r_list)cr

计算ET
In [5]:
et = cinrad.calc.quick_et(r_list)et

In [ ]:
# 计算VIL、VILD# vil = cinrad.calc.quick_vil(r_list)# vild = cinrad.calc.quick_vild(r_list)

In [7]:
# 画个图试试
fig = cinrad.visualize.PPI(et,style="black")

# 将图片保存
# fig("d:/")
# 或者 fig("d:/abc.png")
# 或者 imgName = fig("d:/")
计算HCL
CINRAD使用的论文中的算法（详见github底部）
In [13]:
nFiles = basePath + "/cinrad/bz2/Z_RADR_I_Z9735_20240511082558_O_DOR_SAD_CAP_FMT.bin.bz2"f = cinrad.io.read_auto(nFiles)ref = f.get_data(0, 230, "REF")  # 这里全部使用的是第一个仰角的zdr = f.get_data(0, 230, "ZDR")rho = f.get_data(0, 230, "RHO")kdp = f.get_data(0, 230, "KDP")cHCL = cinrad.calc.hydro_class(ref, zdr, rho, kdp, band="S")  # band手动输入S/C/XcHCLfig = cinrad.visualize.PPI(cHCL, add_city_names=True, dpi=300, style="black")

导出为pyart
In [4]:
from cinrad.io.export import standard_data_to_pyartnFiles = basePath + "/cinrad/bz2/Z_RADR_I_Z9731_20240510070522_O_DOR_SAD_CAP_FMT.bin.bz2"f = cinrad.io.read_auto(nFiles)print(f.available_product(0))sradar = standard_data_to_pyart(f)sradar

## You are using the Python ARM Radar Toolkit (Py-ART), an open source
## library for working with weather radar data. Py-ART is partly
## supported by the U.S. Department of Energy as part of the Atmospheric
## Radiation Measurement (ARM) Climate Research Facility, an Office of
## Science user facility.
##
## If you use this software to prepare a publication, please cite:
##
##     JJ Helmus and SM Collis, JORS 2016, doi: 10.5334/jors.119

['TREF', 'REF', 'ZDR', 'RHO', 'PHI', 'KDP', 'SNRH']
<pyart.core.radar.Radar at 0x249d5a0fdc0>
非标准格式基数据
SB/CC/CA等少部分数据最好是用CinradReader并手动输入参数radar_type
In [3]:
nFiles = basePath + "/cinrad/bz2/Z_RADR_I_Z9636_20200625140000_O_DOR_SB_CAP.bz2"f = cinrad.io.read_auto(nFiles)# f = cinrad.io.CinradReader(nFiles, radar_type="SB")print(type(f).__name__)data = f.get_data(0, 230, "REF")data

CinradReader
In [4]:
f.get_elevation_angles()

array([ 0.49987793,  0.49987793,  1.44470215,  1.44470215,  2.39501953,
        3.34533691,  4.2956543 ,  5.99853516,  9.89868164, 14.59533691,
       19.49523926])
In [5]:
f.angleindex_r

array([0, 2, 4, 5, 6, 7, 8, 9])
相控阵雷达基数据
请更新到最新版本 ≥ 1.9.1
AXPT0364
In [6]:
nFiles = basePath + "/cinrad/bz2/Z_RADR_I_ZS409_20231220000000_O_DOR_AXPT0364_CRA_.bz2"f = cinrad.io.read_auto(nFiles)print(type(f).__name__)data = f.get_data(0, 230, "RHO")data

StandardData
In [ ]:
fig = cinrad.visualize.PPI(data, style="white") # 画图

In [8]:
f.available_product(0)

['TREF', 'REF', 'VEL', 'SW', 'ZDR', 'PHI', 'KDP', 'RHO', 'SNRH', 'SQI']
In [11]:
str(f.el)

'[0.9, 1.8, 2.7, 3.6, 4.5, 5.4, 6.3, 7.2, 8.1, 9.0, 9.9, 10.8, 11.7, 12.6, 13.5, 14.4, 15.3, 16.2, 17.1, 18.0, 18.9, 19.8, 20.7, 21.6, 22.5, 23.4, 24.3, 25.2, 26.1, 27.0, 27.9, 28.8, 29.7, 30.6, 31.5, 32.4, 33.3, 34.2, 35.1, 36.0]'
In [12]:
str(f.available_tilt('REF')) #REF产品有哪些仰角可以读取

'[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39]'
DXK
In [14]:
nFiles = basePath + "/cinrad/bz2/Z_RADR_I_ZBJ02_20210815155836_O_DOR_DXK_CAR.bin"f = cinrad.io.read_auto(nFiles)print(type(f).__name__)data = f.get_data(0, 230, "REF")data

PhasedArrayData
In [15]:
f.available_product(0)
['TREF', 'REF', 'VEL', 'SW', 'ZDR', 'PHI', 'KDP', 'RHO']
In [16]:
f
ig = cinrad.visualize.PPI(data, style="black")

In [17]:
# dbz<=0 的会在画图的时候有一个圈，修改为nan即可
import numpy as npdata["REF"].values = np.ma.masked_less(data["REF"].values, 0)fig = cinrad.visualize.PPI(data, style="black")

AXPT0364_FMT
AXPT貌似都是标准格式的
In [20]:
nFiles = basePath + "/cinrad/bz2/Z_RADR_I_ZA601_20240415183600_O_DOR_AXPT0364_CRA_FMT.bin.bz2"f = cinrad.io.read_auto(nFiles)print(type(f).__name__)data = f.get_data(0, 230, "REF")data

StandardData
In [21]:
fig = cinrad.visualize.PPI(data, style="black")

标准格式产品ROSE2.0
已测试支持的产品：HCL、HI、V、MAX、HSR、WER、STI、CC、KDP、PDP、R、ZDR、CR、M、OHP、CAR、DOHP、DSTP、STP、VIL、SRM、THP、TVS、ET、VWP、DTHP、UAM、SWP、CAPPI
已测试不支持的产品：ML、SS、CS
BR产品
In [6]:
nFiles = basePath + "/cinrad/bz2/Z_RADR_I_Z9094_20230402071419_P_DOR_CAD_R_300_400_5_FMT.bin"f = cinrad.io.read_auto(nFiles)dt = f.get_data()dt

CR产品
In [51]:
nFiles = basePath + "/cinrad/bz2/Z_RADR_I_Z9730_20240510000142_P_DOR_SB_CR_250X250_230_NUL_FMT.bin"f = cinrad.io.read_auto(nFiles)dt = f.get_data()dt

MAX产品
这个产品实际上读取出来只能看到俯视最大数据投影区的数据，因为只读取了1遍。
In [33]:
nFiles = basePath + "/cinrad/bz2/Z_RADR_I_ZN701_20230805000000_P_DOR_YLD2-D_CR_150X150_230_NUL_FMT.bin.bz2"f = cinrad.io.read_auto(nFiles)data = f.get_data()data

WER产品
In [34]:
nFiles = basePath + "/cinrad/bz2/Z_RADR_I_Z9731_20240519000105_P_DOR_SAD_WER_250_50X50_NUL_FMT.bin"f= cinrad.io.read_auto(nFiles)data = f.get_data()data

NEXRAD Level 3产品
In [52]:
# 用Metpy读取nFiles = basePath + "/cinrad/bz2/KCYS_SDUS25_N3UCYS_202405160000"from metpy.io.nexrad import Level3Filef = Level3File(nFiles)print(f)

d:/temp//cinrad/bz2/KCYS_SDUS25_N3UCYS_202405160000: Base Velocity Data Array
MsgHdr(code=99, date=19860, time=247, msg_len=18257, src_id=335, dest_id=0, num_blks=3)
ProdDesc(divider=-1, lat=41152, lon=-104806, height=6192, prod_code=99, op_mode=2, vcp=215, seq_num=6741, vol_num=52, vol_date=19860, vol_start_time=28, prod_gen_date=19860, prod_gen_time=247, dep1=0, dep2=0, el_num=6, dep3=31, thr1=-635, thr2=5, thr3=254, thr4=0, thr5=0, thr6=0, thr7=0, thr8=0, thr9=0, thr10=0, thr11=0, thr12=0, thr13=0, thr14=0, thr15=0, thr16=0, dep4=-92, dep5=96, dep6=0, dep7=6912, dep8=1, dep9=6, dep10=27294, version=0, spot_blank=0, sym_off=60, graph_off=0, tab_off=0)
[-635, 5, 254, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
[0, 0, 31, -92, 96, 0, 6912, 1, 6, 27294]
{'msg_time': datetime.datetime(2024, 5, 16, 0, 4, 7), 'vol_time': datetime.datetime(2024, 5, 16, 0, 0, 28), 'prod_time': datetime.datetime(2024, 5, 16, 0, 4, 7), 'el_angle': 3.1, 'min': -92, 'max': 96, 'delta_time': 216, 'supplemental_scan': 'Non-supplemental scan', 'compression': 1, 'uncompressed_size': 420510}
CYS
In [53]:
# 调用read_auto，会跳转到PUP接口，但是支持的产品类型比较少。nFiles = basePath + "/cinrad/bz2/KCYS_SDUS25_N3UCYS_202405160000"f = cinrad.io.read_auto(nFiles)print(f)

---------------------------------------------------------------------------
RadarDecodeError                          Traceback (most recent call last)
Cell In[53], line 3
      1 # 调用read_auto，会跳转到PUP接口，但是支持的产品类型比较少。
      2 nFiles = basePath + "/cinrad/bz2/KCYS_SDUS25_N3UCYS_202405160000"
----> 3 f = cinrad.io.read_auto(nFiles)
      4 print(f)

File d:\software\anaconda3\envs\fix\lib\site-packages\cinrad\io\__init__.py:58, in read_auto(filename, **kwargs)
     56     return SWAN(filename, **kwargs)
     57 elif flag[:4] == b"NXUS" or flag[:4] == b"SDUS" or flag[:4] == b"NOUS":
---> 58     return PUP(filename)
     59 sc_flag = flag[100:106]
     60 cc_flag = flag[116:122]

File d:\software\anaconda3\envs\fix\lib\site-packages\cinrad\io\level3.py:58, in PUP.__init__(self, file)
     56 self._update_radar_info()
     57 product_code = f.prod_desc.prod_code
---> 58 self.dtype = self._det_product_type(product_code)
     59 self.radial_flag = self._is_radial(product_code)
     60 data_block = f.sym_block[0][0]

File d:\software\anaconda3\envs\fix\lib\site-packages\cinrad\io\level3.py:189, in PUP._det_product_type(spec)
    187     return "OHP"
    188 else:
--> 189     raise RadarDecodeError("Unsupported product type {}".format(spec))

RadarDecodeError: Unsupported product type 99
HCL产品
In [101]:
nFiles = basePath + "/cinrad/bz2/Z_RADR_I_Z9735_20240511082558_P_DOR_SAD_HCL_250_230_5_FMT.bin"f = cinrad.io.read_auto(nFiles)hcl = f.get_data()print(hcl)fig = cinrad.visualize.PPI(hcl, add_city_names=True, dpi=300, style="black")

<xarray.Dataset> Size: 11MB
Dimensions:    (azimuth: 360, distance: 920)
Coordinates:
  * azimuth    (azimuth) float64 3kB 0.0 0.0175 0.035 ... 6.248 6.266 6.283
  * distance   (distance) float64 7kB 0.25 0.5 0.75 1.0 ... 229.5 229.8 230.0
Data variables:
    HCL        (azimuth, distance) float64 3MB nan nan nan 4.0 ... 6.0 6.0 6.0
    longitude  (azimuth, distance) float64 3MB 113.0 113.0 113.0 ... 113.0 113.0
    latitude   (azimuth, distance) float64 3MB 25.74 25.74 25.74 ... 27.8 27.81
    height     (azimuth, distance) float64 3MB 0.3672 0.3694 ... 5.475 5.484
Attributes:
    elevation:        0.5
    range:            230.0
    scan_time:        2024-05-11 08:25:58
    site_code:        Z9735
    site_name:        郴州
    site_longitude:   112.975
    site_latitude:    25.734167
    tangential_reso:  0.25
    task:             VCP21D
STI产品
这是目前为止，唯一一个返回JSON数据的接口。
In [ ]:
# rose stiimport jsonimport numpy as np
# 用于json序列化numpy数据，确保前端可以正常解析class NumpyEncoder(json.JSONEncoder):    def default(self, obj):        if isinstance(obj, np.float32):            return float(obj)        elif isinstance(obj, np.int32):            return float(obj)        return json.JSONEncoder.default(self, obj)

f = cinrad.io.read_auto(basePath+'cinrad/bz2/Z9737_20230424171836Z_STI_00_58')data = f.get_data()js = json.dumps(data, cls = NumpyEncoder)js

In [ ]:
{    "data": [        {            "id": "217",            "current_position": [                111.7,                28.0            ],            "current_speed": 21.4,            "current_direction": 260.35,            "forecast_position": [                [                    111.9,                    28.1                ],                [                    112.0,                    28.1                ]            ],            "history_position": [                [                    111.6,                    28.0                ],                [                    111.5,                    28.0                ]            ],            "max_ref": 51.0,            "max_ref_height": 1587.0,            "vil": 13.26,            "top_height": 5585.0        }    ],    "attrs": {        "scan_time": "2023-04-24 17:18:36",        "site_code": "Z9737",        "site_name": "益阳",        "site_longitude": 112.5,        "site_latitude": 28.4,        "task": "VCP21",        "sti_count": 64.0    }}

UAM产品
In [13]:
nFiles = basePath + "/cinrad/bz2/Z9737_20231128135254Z_UAM_00_hail"f = cinrad.io.read_auto(nFiles)dt = f.get_data()dt

VEL速度产品
In [99]:
nFiles = basePath + "/cinrad/bz2/Z_RADR_I_Z9735_20240510000045_P_DOR_SAD_V_500_230_5_FMT.bin"f = cinrad.io.read_auto(nFiles)dt = f.get_data()print(dt)dt = dt.rename({'Vc': 'VEL'})fig = PPI(dt, add_city_names=True, dpi=300, style="black")

<xarray.Dataset> Size: 5MB
Dimensions:    (azimuth: 357, distance: 460)
Coordinates:
  * azimuth    (azimuth) float64 3kB 0.002444 0.02009 0.03774 ... 6.268 0.002444
  * distance   (distance) float64 4kB 0.5 1.0 1.5 2.0 ... 229.0 229.5 230.0
Data variables:
    Vc         (azimuth, distance) float64 1MB nan nan nan nan ... nan nan nan
    longitude  (azimuth, distance) float64 1MB 113.0 113.0 113.0 ... 113.0 113.0
    latitude   (azimuth, distance) float64 1MB 25.74 25.74 25.75 ... 27.8 27.81
    height     (azimuth, distance) float64 1MB 0.3694 0.3738 ... 5.466 5.484
Attributes:
    elevation:        0.5
    range:            230.0
    scan_time:        2024-05-10 00:00:45
    site_code:        Z9735
    site_name:        郴州
    site_longitude:   112.975
    site_latitude:    25.734167
    tangential_reso:  0.5
    task:             VCP21D
ZDR产品
In [15]:
nFiles = basePath + "/cinrad/bz2/Z_RADR_I_Z9743_20240510000107_P_DOR_SAD_ZDR_250_230_5_FMT.bin"f = cinrad.io.read_auto(nFiles)dt = f.get_data()dt

CC产品
In [97]:
nFiles = basePath + "/cinrad/bz2/Z_RADR_I_Z9739_20240510000327_P_DOR_SAD_CC_250_230_15_FMT.bin"f = cinrad.io.read_auto(nFiles)dt = f.get_data()print(dt)fig = PPI(dt, add_city_names=True, dpi=300, style="black")

<xarray.Dataset> Size: 11MB
Dimensions:    (azimuth: 366, distance: 920)
Coordinates:
  * azimuth    (azimuth) float64 3kB 0.005061 0.02228 0.03949 ... 6.271 0.005061
  * distance   (distance) float64 7kB 0.25 0.5 0.75 1.0 ... 229.5 229.8 230.0
Data variables:
    RHO        (azimuth, distance) float64 3MB nan nan nan nan ... nan nan nan
    longitude  (azimuth, distance) float64 3MB 111.4 111.4 111.4 ... 111.5 111.5
    latitude   (azimuth, distance) float64 3MB 27.21 27.21 27.21 ... 29.28 29.28
    height     (azimuth, distance) float64 3MB 0.2975 0.3041 ... 9.41 9.423
Attributes:
    elevation:        1.5
    range:            230.0
    scan_time:        2024-05-10 00:03:27
    site_code:        Z9739
    site_name:        邵阳
    site_longitude:   111.44583
    site_latitude:    27.207222
    tangential_reso:  0.25
    task:             VCP21D
CAPPI产品
cinrad目前不支持计算CAPPI，这里指的是ROSE生成的产品
In [2]:
nFiles = basePath + "/cinrad/bz2/Z_RADR_I_Z9094_20230402073754_P_DOR_CAD_CAR_150_200_NUL_FMT.bin"f = cinrad.io.read_auto(nFiles)print(type(f).__name__)dt = f.get_data()dt# CAPPI是可以多层的，因此画图前需要将height维度降维# dt0 = dt.sel(height=dt.coords["height"][0])

StandardPUP
HI产品
In [63]:
nFiles = basePath + "/cinrad/bz2/Z_RADR_I_Z9094_20230402070825_P_DOR_CAD_HI_NUL_200_NUL_FMT.bin"f = cinrad.io.read_auto(nFiles)data = f.get_data()data

OHP产品
In [64]:
nFiles = basePath + "/cinrad/bz2/Z_RADR_I_Z9094_20230402072013_P_DOR_CAD_OHP_150_200_NUL_FMT.bin"f = cinrad.io.read_auto(nFiles)data = f.get_data()data

VWP产品
In [25]:
# 读取VWP
import warningsimport cinrad import numpy as np
warnings.filterwarnings("ignore")nFiles = "d:/temp/cinrad/bz2/Z_RADR_I_Z9731_20240126000025_P_DOR_SAD_VWP_NUL_NUL_NUL_FMT.bin"f = cinrad.io.read_auto(nFiles)vwp = f.get_data()# 有部分雷达的产品输出时，风速有小于0的情况，需要过滤掉# vwp.wind_speed.values = np.ma.masked_less(vwp.wind_speed.values, 0.0)# vwp.wind_speed.values = np.round(vwp.wind_speed.values, 2)vwp

VWP 画图
In [26]:
import numpy as npimport matplotlib.pyplot as pltimport xarrayimport datetime
height = np.round(np.array(vwp.height) / 1000, 1)times = np.array(vwp.times)times = [datetime.datetime.fromtimestamp(time, datetime.timezone.utc) for time in times]times = [time.strftime("%H:%M") for time in times]wind_direction = vwp.wind_directionwind_speed = vwp.wind_speedrms = vwp.rmsu = -wind_speed * np.sin(np.radians(wind_direction))v = -wind_speed * np.cos(np.radians(wind_direction))

def _get_vwp_color(rms: xarray.DataArray) -> list:    """    风羽的颜色是由RMS值决定的。    """    data = rms.data    color_map = [        (0, "#00FF00"),        (2, "#FFFF00"),        (4, "#FF0000"),        (6, "#00EFFF"),        (8, "#FF7BFF"),        (10, "#FFFFFF"),    ]
    color = []    for value in data:        cr = color_map[0][1]        for i in range(len(color_map)):            if value > color_map[i][0]:                cr = color_map[i][1]        color.append(cr)    return color

fig, ax = plt.subplots(1, 1, figsize=(12, 15))ax.set_xlabel("Time(UTC)")plt.style.use("dark_background")ax.set_ylabel("Height (km)")nums = np.arange(1, 31, dtype=int)ax.set_yticks(nums, labels=height)plt.ylim((0.5, 30))plt.grid(True, which="both", axis="y", linestyle="--")for i in range(len(times)):    x = [times[i] for _ in range(len(height))]    colors = _get_vwp_color(rms[i])    ax.barbs(        x,        nums,        u[i],        v[i],        rounding=False,        barb_increments=dict(half=2, full=4, flag=20),        sizes=dict(emptybarb=0.01, spacing=0.23, height=0.5, width=0.25),        color=colors,    )plt.tight_layout()

SWAN3产品
三维拼图 MOSAIC
In [5]:
import numpy as np
nFiles = basePath + "/cinrad/bz2/Z_OTHE_RADAMOSAIC_20240604020600.bin.bz2"f = cinrad.io.read_auto(nFiles)# f = cinrad.io.SWAN(nFiles, product="3DREF")  # 用这种接口加上参数的话，可以直接PPI画图data = f.get_data(level=5)# data[data<0] = np.nandata

CR
In [4]:
nFiles = basePath + "/cinrad/bz2/Z_OTHE_RADAMCR_20220525085400.bin.bz2"f = cinrad.io.read_auto(nFiles)# f = cinrad.io.SWAN(nFiles, product="CR")  # 用这种接口加上参数的话，可以直接PPI画图data = f.get_data()data

TREC
In [5]:
nFiles = basePath + "/cinrad/bz2/Z_TREC_20240520233600.BIN.BZ2"f = cinrad.io.read_auto(nFiles)data = f.get_data()data

ET/VIL/VILD/QPE
这几个数据都是扩大了10倍存储的，读取后记得除以10
In [ ]:
nFiles = basePath + "/cinrad/bz2/Z_OTHE_RADAMTOP_20240604014800.BIN.BZ2"f = cinrad.io.read_auto(nFiles)data = f.get_data()data = data / 10datafig = PPI(data, add_city_names=True, dpi=300, style="black")

探测中心天气雷达拼图系统v3产品
拼图
In [6]:
# MocMosaicnFiles = (    basePath    + "cinrad/bz2/Z_RADA_C_BABJ_20240424000701_P_DOR_ACHN_CAP_20240424_000000_23.bin")f = cinrad.io.read_auto(nFiles)print("CAPPI高度:" + str(f.header["height"][0])) # CAPPI才有意义dt = f.get_data()print(dt)

CAPPI高度:16000
<xarray.Dataset> Size: 208MB
Dimensions:    (latitude: 4200, longitude: 6200)
Coordinates:
  * latitude   (latitude) float64 34kB 12.2 12.21 12.22 ... 54.18 54.19 54.2
  * longitude  (longitude) float64 50kB 73.0 73.01 73.02 ... 135.0 135.0 135.0
Data variables:
    CAP        (latitude, longitude) float64 208MB nan nan nan ... nan nan nan
Attributes:
    scan_time:        2024-04-24 00:00:00
    site_code:        MOC
    site_name:        MOC
    tangential_reso:  nan
    range:            nan
    elevation:        0
单站
In [8]:
# MocMosaicnFiles = (    basePath    + "cinrad/bz2/Z_RADA_C_BABJ_20240520001131_P_DOR_Z9734_CAP_20240520_000532.bin")f = cinrad.io.read_auto(nFiles)dt = f.get_data()print(dt)# 此时发现height维度，如果要画图的话，需要选择一个高度层if "height" in dt.coords:    dt0 = dt.sel(height=dt.coords["height"][18])    dt0 = dt0.rename({"CAP": "CR"})  # 重命名，因为PPI并不支持CAP    fig = PPI(dt0, style="black")

<xarray.Dataset> Size: 71MB
Dimensions:    (height: 24, latitude: 575, longitude: 645)
Coordinates:
  * height     (height) int32 96B 500 1000 1500 2000 ... 13000 14000 15000 16000
  * latitude   (latitude) float64 5kB 24.06 24.07 24.08 ... 29.79 29.8 29.81
  * longitude  (longitude) float64 5kB 109.4 109.4 109.4 ... 115.8 115.8 115.8
Data variables:
    CAP        (height, latitude, longitude) float64 71MB nan nan ... nan nan
Attributes:
    scan_time:        2024-05-20 00:05:00
    site_code:        Z9734
    site_name:        衡阳g_Z9734
    site_longitude:   112.5822
    site_latitude:    26.935
    tangential_reso:  0.1
    range:            320
    elevation:        0
    task:             VCP21
雷达拼图
基本反射率拼图（老方法）
其实所有单仰角的数据都可以拼
In [8]:
mergePath = "d:/temp/cinrad/cr/"file1 = mergePath + "Z_RADR_I_Z9730_20240526075538_O_DOR_SB_CAP_FMT.bin.bz2"file2 = mergePath + "Z_RADR_I_Z9731_20240526075702_O_DOR_SAD_CAP_FMT.bin.bz2"f1 = cinrad.io.read_auto(file1)  # 读取文件br1 = f1.get_data(0, 230, "REF")f2 = cinrad.io.read_auto(file2)br2 = f2.get_data(0, 230, "REF")
gm_br = cinrad.calc.GridMapper([br1, br2], max_dist=0.005)grid_br = gm_br(step=0.005)print(grid_br)fig = cinrad.visualize.PPI(grid_br, style="black", add_city_names=True)

<xarray.Dataset>
Dimensions:    (latitude: 997, longitude: 966)
Coordinates:
  * latitude   (latitude) float64 26.39 26.39 26.4 26.4 ... 31.36 31.36 31.37
  * longitude  (longitude) float64 110.7 110.7 110.7 110.7 ... 115.5 115.5 115.5
Data variables:
    REF        (latitude, longitude) float64 nan nan nan nan ... nan nan nan nan
Attributes:
    elevation:        0
    range:            230
    scan_time:        2024-05-26 08:00:20
    site_code:        RADMAP
    site_name:        RADMAP
    tangential_reso:  nan
    task:             VCP21
基本反射率拼图（新方法）
此方法需要升级到github上的最新版本,理论上ZDR、KDP这些也可以拼图
In [ ]:
pip install git+https://github.com/pysoer/PyCINRAD.git@mergeCR -U --force-reinstall --no-deps

In [9]:
mergePath = "d:/temp/cinrad/cr/"file1 = mergePath + "Z_RADR_I_Z9730_20240526075538_O_DOR_SB_CAP_FMT.bin.bz2"file2 = mergePath + "Z_RADR_I_Z9731_20240526075702_O_DOR_SAD_CAP_FMT.bin.bz2"f1 = cinrad.io.read_auto(file1)  # 读取文件br1 = f1.get_data(0, 230, "REF")f2 = cinrad.io.read_auto(file2)br2 = f2.get_data(0, 230, "REF")
br1_cr = cinrad.calc.polar_to_xy(br1)br2_cr = cinrad.calc.polar_to_xy(br2)
gm_br = cinrad.calc.GridMapper([br1_cr, br2_cr], max_dist=0.005)grid_br = gm_br(step=0.005)print(grid_br)fig = cinrad.visualize.PPI(grid_br, style="black", add_city_names=True)

<xarray.Dataset>
Dimensions:    (latitude: 997, longitude: 966)
Coordinates:
  * latitude   (latitude) float64 26.39 26.39 26.4 26.4 ... 31.36 31.36 31.37
  * longitude  (longitude) float64 110.7 110.7 110.7 110.7 ... 115.5 115.5 115.5
Data variables:
    CR         (latitude, longitude) float64 nan nan nan nan ... nan nan nan nan
Attributes:
    elevation:        0
    range:            230
    scan_time:        2024-05-26 08:00:20
    site_code:        RADMAP
    site_name:        RADMAP
    tangential_reso:  nan
    task:             VCP21
组合反射率拼图（新方法）
此方法需要升级到github上的最新版本，其实适用于所有经过calc方法的数据
In [7]:
mergePath = "d:/temp/cinrad/cr/"file1 = mergePath + "Z_RADR_I_Z9730_20240526075538_O_DOR_SB_CAP_FMT.bin.bz2"file2 = mergePath + "Z_RADR_I_Z9731_20240526075702_O_DOR_SAD_CAP_FMT.bin.bz2"f1 = cinrad.io.read_auto(file1)  # 读取文件list1 = list(f1.iter_tilt(230, "REF"))f2 = cinrad.io.read_auto(file2)list2 = list(f2.iter_tilt(230, "REF"))
cr1 = cinrad.calc.quick_cr(list1)cr2 = cinrad.calc.quick_cr(list2)
gm_cr = cinrad.calc.GridMapper([cr1, cr2], max_dist=0.005)grid_cr = gm_cr(step=0.005)print(grid_cr)fig = cinrad.visualize.PPI(grid_cr, style="black", add_city_names=True)

<xarray.Dataset>
Dimensions:    (latitude: 997, longitude: 966)
Coordinates:
  * latitude   (latitude) float64 26.39 26.39 26.4 26.4 ... 31.36 31.36 31.37
  * longitude  (longitude) float64 110.7 110.7 110.7 110.7 ... 115.5 115.5 115.5
Data variables:
    CR         (latitude, longitude) float64 nan nan nan nan ... nan nan nan nan
Attributes:
    elevation:        0
    range:            230
    scan_time:        2024-05-26 08:00:20
    site_code:        RADMAP
    site_name:        RADMAP
    tangential_reso:  nan
    task:             VCP21
画图
cinrad.visualize.PPI
|参数|功能|
|:-:|:-:|
|cmap|色阶|
|norm|色阶范围|
|nlabel|色阶条标注个数|
|label|色阶条标注|
|highlight|地区边界高亮|
|dpi|分辨率|
|extent|绘图的经纬度范围 e.g. extent=[90, 91, 29, 30]|
|section|在ppi图中绘制的剖面的数据，为xarray.Dataset类型|
|style|背景颜色，可设置为黑色black或者白色white或者透明transparent|
|add_city_names|标注城市名|
PPI支持以下xarray.Dataset类型数据：
"REF":"Base Reflectivity"
"VEL":"Base Velocity"
"CR":"Composite Ref."
"ET":"Echo Tops"
"VIL":"V Integrated Liquid"
"ZDR":"Differential Ref."
"PHI":"Differential Phase"
"RHO":"Correlation Coe."
"TREF":"Total Reflectivity"
"KDP":"Spec. Diff. Phase"
"VILD":"VIL Density"
"OHP":"One-Hour Precip."
"HCL":"Hydrometeor Class"
"VELSZ":"Velocity SZ Recovery"
有些数据读取出来并不是以上名称
例如产品V读取出来可能是Vc,此时重命名即可data=data.rename({"Vc":"VEL"})
画图基本操作
In [10]:
# 读取数据，这里以BR为例nFiles = basePath + "/cinrad/bz2/Z_RADR_I_Z9737_20231025220414_O_DOR_SA_CAP_FMT.bin.bz2"f = cinrad.io.read_auto(nFiles)
# 想把时间改成北京时？# f.scantime = f.scantime + datetime.timedelta(hours=8)data = f.get_data(0, 230, "REF")
# 开始画图fig = PPI(    data, dpi=300, style="black", add_city_names=True)  # 背景颜色支持：white /black / transparent
# 附加操作fig.plot_range_rings([10, 20, 30, 50, 150], color="white", linewidth=1)  # 用这个来画圈fig.gridlines(draw_labels=True, linewidth=1, color="white")  # 用这个来画经纬度网格线
# 要在地图上画一个红色大圆点？import cartopy.crs as ccrs
fig.geoax.scatter(    x=112.34, y=29.22, s=500, c="r", marker=".", transform=ccrs.PlateCarree())


# 将图片保存
# fig("d:/")
# 或者 fig("d:/abc.png")
# 或者 imgName = fig("d:/")
<matplotlib.collections.PathCollection at 0x1e96d6d7ee0>
剖面图
In [9]:
nFiles = basePath + "/cinrad/bz2/Z_RADR_I_Z9737_20231025220414_O_DOR_SA_CAP_FMT.bin.bz2"f = cinrad.io.read_auto(nFiles)rl = list(f.iter_tilt(230, "REF"))vcs = cinrad.calc.VCS(rl)sec = vcs.get_section(    start_cart=(111.97, 28.05), end_cart=(112.59, 28.14))  # 传入经纬度坐标# sec = vcs.get_section(start_polar=(0, 30), end_polar=(90, 30)) # 传入极坐标fig = cinrad.visualize.Section(sec)plt.scatter(    0.5, 3, s=100, c="r", marker="o")  # 在x轴的0.5位置（用两端坐标经纬度换算一下），3公里高度位置画个点plt.axline((0, 2), (1, 2), linewidth=2, color="r", marker="o")
  # 在2公里高度画根线
<matplotlib.lines.AxLine at 0x1463fe65c40>
在PPI图下方添加剖面图
In [10]:
# VCS画图测试nFiles = basePath + "/cinrad/bz2/Z_RADR_I_Z9737_20231025220414_O_DOR_SA_CAP_FMT.bin.bz2"f = cinrad.io.read_auto(nFiles)rl = list(f.iter_tilt(230, "REF"))cr = cinrad.calc.quick_cr(rl)fig = PPI(cr, dpi=300, style="black")fig.settings["is_inline"] = False  # 强行跳过notebook模式检查vcs = cinrad.calc.VCS(rl)sec = vcs.get_section(start_cart=(111.97, 28.05), end_cart=(112.59, 28.14))  # 传入经纬度坐标# sec = vcs.get_section(start_polar=(113, 250), end_polar=(114, 28)) # 传入极坐标fig.plot_cross_section(sec, linecolor="red")

自定义色标
Cinrad自带了很多的标准色标，一般不需要自定义；
如果你需要使用其他色标，则可以使用下面的方法自定义一个
In [10]:
nFiles = basePath + "/cinrad/bz2/Z_RADR_I_Z9735_20240511082558_P_DOR_SAD_HCL_250_230_5_FMT.bin"f = cinrad.io.read_auto(nFiles)data = f.get_data()import matplotlib.colors as cmxfrom cinrad.visualize.gpf import _cmap
cmapFile = (    basePath + "HCL"  # 色标文件参照/data/colormap/目录下的格式，也放到这个目录下)cmap = _cmap(cmapFile)["cmap"]norm = cmx.Normalize(0, 10)  # 取值区间label = ["小雨", "大雨", "冰雹", "大雨滴", "晴空回波", "地物", "干雪", "湿雪", "冰晶", "霰", "未知",""]fig = PPI(    data,    add_city_names=True,    dpi=300,    style="black",    cmap=cmap,    norm=norm,    label=label,)# fig("d:/temp/")# HCL.map文件内容如下："""*NAME:HCL*AUTHOR:CMA*TYPE:LISTED*UNIT:None*UNDER:0/251/144*OVER:231/0/2550 0/251/1441 0/187/02 255/0/03 208/208/964 156/156/1565 118/118/1186 0/255/2557 0/144/2558 255/176/1769 210/132/13210 231/0/255"""

'\n*NAME:HCL\n*AUTHOR:CMA\n*TYPE:LISTED\n*UNIT:None\n*UNDER:0/251/144\n*OVER:231/0/255\n0 0/251/144\n1 0/187/0\n2 255/0/0\n3 208/208/96\n4 156/156/156\n5 118/118/118\n6 0/255/255\n7 0/144/255\n8 255/176/176\n9 210/132/132\n10 231/0/255\n'
番外篇
基数据编辑
修改标准格式基数据，替换里面的站名站号等基本信息
In [ ]:
import structimport bz2



def prepare_file(file):    if hasattr(file, "read"):        return file    f = open(file, "rb")    magic = f.read(3)    f.close()    if magic.startswith(b"BZh"):        return bz2.BZ2File(file, "rb")    return open(file, "rb")

def modify_fmt_site_config(    file: str,    newfile: str,    site_code: str = None,    site_name: str = None,    latitude: float = None,    longitude: float = None,    atenna_height: float = None,    ground_height: float = None,    radar_type: int = None,):    """    修改标准格式天气雷达数据的站点信息@pysoer
    :param file: 原始基数据文件名    :param newfile: 修改后的基数据文件名    :param site_code: 站点代码    :param site_name: 站点名称    :param latitude: 纬度    :param longitude: 经度    :param atenna_height: 天线高度    :param ground_height: 地面高度    :param radar_type: 雷达类型        1–SA 2–SB        3–SC 4–SAD        33–CA 34–CB        35–CC 36–CCJ        37–CD 38–CAD        65–XA 66–XAD    :return:    """
    # 打开原始文件，输出数据    f = cinrad.io.StandardData(file)    print(f.get_data(0, 230, "REF"))
    # 打开原始文件，读取所有字节    ff = prepare_file(file)    content = ff.read()    ff.close()    print(content[:100])    # 修改指定位置数据    if site_code:        site_code_encoded = struct.pack("8s", site_code.encode("utf-8"))        content = content[:32] + site_code_encoded + content[40:]    if site_name:        site_name_encoded = struct.pack(            "32s", site_name.encode("Unicode_escape", errors="ignore")        )        content = content[:40] + site_name_encoded + content[72:]    if latitude:        latitude_encoded = struct.pack("f", latitude)        content = content[:72] + latitude_encoded + content[76:]    if longitude:        longitude_encoded = struct.pack("f", longitude)        content = content[:76] + longitude_encoded + content[80:]    if atenna_height:        atenna_height_encoded = struct.pack("f", atenna_height)        content = content[:80] + atenna_height_encoded + content[84:]    if ground_height:        ground_height_encoded = struct.pack("f", ground_height)        content = content[:84] + ground_height_encoded + content[88:]    if radar_type:        radar_type_encoded = struct.pack("h", radar_type)        content = content[:104] + radar_type_encoded + content[106:]    print(content[:100])    # 写入新文件    if file.lower().endswith("bz2"):        with bz2.open(newfile, "wb") as f:            f.write(content)    else:        with open(newfile, "wb") as f0:            f0.write(content)
    # 输出修改后的数据    f1 = cinrad.io.StandardData(newfile)    dt1 = f1.get_data(0, 230, "REF")    print(dt1)

modify_fmt_site_config(    file="D:/temp/bz2/Z_RADR_I_Z9737_20231025220414_O_DOR_SA_CAP_FMT.bin.bz2",    newfile="D:/temp/bz2/Z_RADR_I_Z9737_20231025220414_O_DOR_SA_CAP_FMT_EDIT.bin.bz2",    site_code="T9999",    site_name="hello",    latitude=32.1,    longitude=118.9,    atenna_height=100,    ground_height=100,    radar_type=1,)
http://mp.weixin.qq.com/s?__biz=MzkzMDYxNDYzOQ==&mid=2247486581&idx=1&sn=e86cf8549020e61bcd0e4a2030e45a0d
气python风雨
主要发一些涉及大气科学的Python文章与个人学习备忘录