基于pycwr对国产雷达数据进行衰减订正
个人信息
公众号：气python风雨
关注我获取更多学习资料，第一时间收到我的Python学习资料，也可获取我的联系方式沟通合作
温馨提示
由于可视化代码过长隐藏，可点击运行Fork查看
若没有成功加载可视化图，点击运行可以查看
ps:隐藏代码在【代码已被隐藏】所在行，点击所在行，可以看到该行的最右角，会出现个三角形，点击查看即可
前言
项目目标
本项目旨在解决气象c波段雷达衰减订正的问题 本来想写使用pyhail对国产雷达进行冰雹反演的，后发现其算法需要雷达数据有十个仰角才可运行。 于是中途改为写衰减订正了。 本文是基于pycwr的源码qc模块搓出来的，原作者没写那咱们补上吧
项目方法
在以下内容中，将详细介绍
1、基于pycwr的三种气象雷达衰减订正
2、如何对pyart的数据进行查看，数据变量如何自由导出、修改、导入
3、如何使用pyart进行简单可视化
安装与导入库
!pip install pycwr -i https://pypi.mirrors.ustc.edu.cn/simple/from matplotlib import pyplot as plt
import numpy as np
import cartopy.crs as ccrs
from pycwr.io import read_auto
import pyart
from pycwr.qc import correct_attenuation_HB,correct_attenuation,pia_from_kdp
读取c波段数据与转格式
PATH='/home/mw/input/pycwr5461/Z_RADR_I_Z9240_20190703101340_O_DOR_SC_CAP.bin.bz2'
PRD = read_auto(PATH)
radar = PRD.ToPyartRadar()数据信息查看
radar.info()altitude:
data: <ndarray of type: float64 and shape: (1,)>
long_name: Altitude
standard_name: Altitude
units: meters
positive: up
altitude_agl: None
antenna_transition: None
azimuth:
data: <ndarray of type: float64 and shape: (3240,)>
units: degrees
standard_name: beam_azimuth_angle
long_name: azimuth_angle_from_true_north
axis: radial_azimuth_coordinate
comment: Azimuth of antenna relative to true north
elevation:
data: <ndarray of type: float64 and shape: (3240,)>
units: degrees
standard_name: beam_elevation_angle
long_name: elevation_angle_from_horizontal_plane
axis: radial_elevation_coordinate
comment: Elevation of antenna relative to the horizontal plane
fields:
reflectivity:
data: <ndarray of type: float32 and shape: (3240, 1000)>
units: dBZ
standard_name: equivalent_reflectivity_factor
long_name: Reflectivity
coordinates: elevation azimuth range
_FillValue: -9999.0
velocity:
data: <ndarray of type: float32 and shape: (3240, 1000)>
units: meters_per_second
standard_name: radial_velocity_of_scatterers_away_from_instrument
long_name: Mean dopper velocity
coordinates: elevation azimuth range
_FillValue: -9999.0
spectrum_width:
data: <ndarray of type: float32 and shape: (3240, 1000)>
units: meters_per_second
standard_name: doppler_spectrum_width
long_name: Doppler spectrum width
coordinates: elevation azimuth range
_FillValue: -9999.0
fixed_angle:
data: <ndarray of type: float64 and shape: (9,)>
long_name: Target angle for sweep
units: degrees
standard_name: target_fixed_angle
instrument_parameters:
pulse_width:
data: <ndarray of type: float32 and shape: (1,)>
units: seconds
comments: Pulse width
meta_group: instrument_parameters
long_name: Pulse width
radar_beam_width_h:
data: <ndarray of type: float32 and shape: (1,)>
units: degrees
meta_group: radar_parameters
long_name: Antenna beam width H polarization
radar_beam_width_v:
data: <ndarray of type: float32 and shape: (1,)>
units: degrees
meta_group: radar_parameters
long_name: Antenna beam width V polarization
frequency:
data: <ndarray of type: float32 and shape: (1,)>
units: s-1
meta_group: instrument_parameters
long_name: Radiation frequency
nyquist_velocity:
data: <ndarray of type: float64 and shape: (3240,)>
units: meters_per_second
comments: Unambiguous velocity
meta_group: instrument_parameters
long_name: Nyquist velocity
latitude:
data: <ndarray of type: float64 and shape: (1,)>
long_name: Latitude
standard_name: Latitude
units: degrees_north
longitude:
data: <ndarray of type: float64 and shape: (1,)>
long_name: Longitude
standard_name: Longitude
units: degrees_east
nsweeps: 9
ngates: 1000
nrays: 3240
radar_calibration: None
range:
data: <ndarray of type: float64 and shape: (1000,)>
units: meters
standard_name: projection_range_coordinate
long_name: range_to_measurement_volume
axis: radial_range_coordinate
spacing_is_constant: true
comment: Coordinate variable for range. Range to center of each bin.
meters_to_center_of_first_gate: 250
meters_between_gates: 250
scan_rate: None
scan_type: ppi
sweep_end_ray_index:
data: <ndarray of type: int64 and shape: (9,)>
long_name: Index of last ray in sweep, 0-based
units: count
sweep_mode:
data: <ndarray of type: |S20 and shape: (9,)>
units: unitless
standard_name: sweep_mode
long_name: Sweep mode
comment: Options are: "sector", "coplane", "rhi", "vertical_pointing", "idle", "azimuth_surveillance", "elevation_surveillance", "sunscan", "pointing", "manual_ppi", "manual_rhi"
sweep_number:
data: <ndarray of type: int32 and shape: (9,)>
units: count
standard_name: sweep_number
long_name: Sweep number
sweep_start_ray_index:
data: <ndarray of type: int64 and shape: (9,)>
long_name: Index of first ray in sweep, 0-based
units: count
target_scan_rate: None
time:
data: <ndarray of type: float32 and shape: (3240,)>
units: seconds since 2019-07-03T18:06:33Z
standard_name: time
long_name: time_in_seconds_since_volume_start
calendar: gregorian
comment: Coordinate variable for time. Time at the center of each ray, in fractional seconds since the global variable time_coverage_start
metadata:
Conventions: CF/Radial instrument_parameters
version: 1.3
title:
institution:
references:
source:
history:
comment:
instrument_name:
original_container: CINRAD/SAB
site_name: LiaoNing
radar_name: CINRAD/SA/SB/CB/SC以下是提供的雷达数据的主要信息的简洁总结：
• Altitude: 单个值表示雷达高度，单位为米。
• Azimuth: 包含3240个值，表示雷达波束相对于真北方向的角度，单位为度。
• Elevation: 包含3240个值，表示雷达波束相对于水平面的角度，单位为度。
• Fields:
• Reflectivity: 数据形状为(3240, 1000)，单位为dBZ，表示反射率。
• Velocity: 数据形状为(3240, 1000)，单位为m/s，表示平均多普勒速度。
• Spectrum Width: 数据形状为(3240, 1000)，单位为m/s，表示多普勒谱宽。
• Fixed Angle: 包含9个值，表示每个扫掠的目标角度，单位为度。
• Instrument Parameters:
• Pulse Width: 单个值表示脉冲宽度，单位为秒。
• Radar Beam Width H: 单个值表示H极化天线波束宽度，单位为度。
• Radar Beam Width V: 单个值表示V极化天线波束宽度，单位为度。
• Frequency: 单个值表示辐射频率，单位为s⁻¹。
• Nyquist Velocity: 包含3240个值，表示无歧义的速度，单位为m/s。
• Latitude and Longitude: 各包含单个值，分别表示纬度和经度，单位分别为度北和度东。
• Range: 包含1000个值，表示到测量体积的距离，单位为米。
• Scan Type: 表示扫描类型为平面位置指示器（PPI）。
• Sweep Information:
• Number of Sweeps: 9次扫掠。
• Number of Gates per Ray: 每条射线1000个门。
• Number of Rays: 3240条射线。
• Sweep Mode: 扫掠模式信息。
• Sweep Start and End Ray Indices: 每次扫掠的起始和结束射线索引。
• Time: 包含3240个值，表示每条射线中心的时间，单位为自2019-07-03T18:06:33Z以来的秒数。
• Metadata: 遵循CF/Radial标准，版本1.3，提供了雷达站点名称和雷达名称等信息。
sorted(list(radar.fields))['reflectivity', 'spectrum_width', 'velocity']科普
http://kejian2.cmatc.cma.cn/2020/tqradar/study_selftaught2.5.html
在应用新一代天气雷达提供的数据进行预报前，需要了解天气雷达数据覆盖范围，波束展宽，地平线和显示区域的特点和存在的局限性。
新一代天气雷达系统提供较高灵敏度及分辨率的反射率因子、平均径向速度及谱宽三种基数据，从图2.29中可见，新一代天气雷达所探测到的反射率因子数据范围可覆盖0～460km，速度数据范围可覆盖0～230km。基数据通过RPG中的气象水文算法形成的分析产品称为导出产品。一些算法以不同的格式显示基数据，称为无特定算法的导出产品；另一些算法则分析基数据以生成导出产品，称为有特定算法的导出产品。和速度有关的导出产品相应覆盖范围在0～230km，基于反射率因子的算法覆盖范围在0～345km。
在下面的函数中会用到库长这一概念，就扫描半径除以半径的格点数，大伙自个算一算
衰减订正一
correct_attenuation_HB(Ref, a=1.67e-4, b=0.7, gate_length=0.075, thrs=59):
"""
Hitschfeld1954, fillvalue set np.nan
:param Ref: radar ref //dbz [naz, nbins]
:param a: coefficients a
:param b: coefficients b
:param gate_length: bin length //m
:param thrs: thrs //dbz
:return: Zc, pia
corz1,spec_at = correct_attenuation_HB(radar.fields['reflectivity']['data'],
gate_length=0.23)
radar.add_field_like('reflectivity',"specific_attenuation", spec_at)
radar.add_field_like('reflectivity', 'corrected_reflectivity', corz1)/opt/conda/lib/python3.9/site-packages/pycwr/qc/attenuation.py:31: RuntimeWarning: overflow encountered in power
k = a * (10.0 ** ((Ref[:, gate] + ksum) / 10.0)) ** b * 2.0 * gate_length# create the plot
fig = plt.figure(figsize=(15, 5))
ax1 = fig.add_subplot(131)
display = pyart.graph.RadarDisplay(radar)
display.plot(
"reflectivity",
0,
ax=ax1,
vmin=0,
vmax=60.0,
colorbar_label="",
title="Raw Reflectivity",
)
ax2 = fig.add_subplot(132)
display.plot(
"specific_attenuation",
0,
vmin=0,
vmax=1.0,
colorbar_label="",
ax=ax2,
title="Specific Attenuation",
)
ax3 = fig.add_subplot(133)
display = pyart.graph.RadarDisplay(radar)
display.plot(
"corrected_reflectivity",
0,
vmin=0,
vmax=60.0,
colorbar_label="",
ax=ax3,
title="Corrected Reflectivity",
)
plt.suptitle("Attenuation correction using pycwr", fontsize=16)
plt.show()衰减订正二
correct_attenuation(ref, wavelength="C", rscale=0.075):
'''Scientific paper: Ośródka, K., Szturc, J., and Jurczyk, A., 2012.
Chain of data quality algorithms for 3-D single-polarization radar
reflectivity (RADVOL-QC system). Meteorol. Appl. (Early View).
rscale : radar gates length(km)
ref :reflectivity factor (dBZ) fillvalue=np.nan'''
return Z_corr, QI
corz2,spec_at2 = correct_attenuation(radar.fields['reflectivity']['data'],wavelength="C", rscale=0.23)
radar.add_field_like('reflectivity',"specific_attenuation2", spec_at2)
radar.add_field_like('reflectivity', 'corrected_reflectivity2', corz2)# create the plot
fig = plt.figure(figsize=(15, 5))
ax1 = fig.add_subplot(131)
display = pyart.graph.RadarDisplay(radar)
display.plot(
"reflectivity",
0,
ax=ax1,
vmin=0,
vmax=60.0,
colorbar_label="",
title="Raw Reflectivity",
)
ax2 = fig.add_subplot(132)
display.plot(
"specific_attenuation",
0,
vmin=0,
vmax=1.0,
colorbar_label="",
ax=ax2,
title="Specific Attenuation2",
)
ax3 = fig.add_subplot(133)
display = pyart.graph.RadarDisplay(radar)
display.plot(
"corrected_reflectivity",
0,
vmin=0,
vmax=60.0,
colorbar_label="",
ax=ax3,
title="Corrected Reflectivity2",
)
plt.suptitle("Attenuation correction using pycwr", fontsize=16)
plt.show()衰减订正三
pia_from_kdp(kdp, dr, gamma=0.08):
"""Retrieving path integrated attenuation from specific differential \
phase (Kdp).
The default value of gamma is based on :cite:Carey2000.
Parameters
----------
kdp : :class:numpy:numpy.ndarray
array specific differential phase
Range dimension must be the last dimension.
dr : gate length (km)
gamma : float
linear coefficient (default value: 0.08) in the relation between
Kdp phase and specific attenuation (alpha)
Returns
-------
output : :class:numpy:numpy.ndarray
array of same shape as kdp containing the path integrated attenuation
"""
就是基于kdp算偏差
由于手头的c波段雷达数据是单偏振的，
单偏振没kdp
那么只能拿另一个双偏振的s波段试试PATH='/home/mw/input/pycwr5461/Z_RADR_I_Z9898_20190828181529_O_DOR_SAD_CAP_FMT (1).bin.bz2'
PRD = read_auto(PATH)
radar1 = PRD.ToPyartRadar()sorted(list(radar1.fields))['cross_correlation_ratio',
'differential_phase',
'differential_reflectivity',
'horizontal_signal_noise_ratio',
'normalized_coherent_power',
'reflectivity',
'specific_differential_phase',
'spectrum_width',
'total_power',
'velocity']好的，我们看到了我们需要的变量
在这些变量列表中，specific_differential_phase 就是我们所寻找的 KDP 变量
refl_offset = pia_from_kdp(radar1.fields['specific_differential_phase']['data'],0.23)
cal_refl = radar1.fields['reflectivity']['data'] - refl_offset
radar1.add_field_like('reflectivity', 'refl_offset', refl_offset)
radar1.add_field_like('reflectivity', 'corrected_reflectivity', cal_refl)# create the plot
fig = plt.figure(figsize=(15, 5))
ax1 = fig.add_subplot(131)
display = pyart.graph.RadarDisplay(radar1)
display.plot(
"reflectivity",
0,
ax=ax1,
vmin=0,
vmax=60.0,
colorbar_label="",
title="Raw Reflectivity",
)
ax2 = fig.add_subplot(132)
display.plot(
"refl_offset",
0,
vmin=0,
vmax=1.0,
colorbar_label="",
ax=ax2,
title="Specific Attenuation3",
)
ax3 = fig.add_subplot(133)
display = pyart.graph.RadarDisplay(radar1)
display.plot(
"corrected_reflectivity",
0,
vmin=0,
vmax=60.0,
colorbar_label="",
ax=ax3,
title="Corrected Reflectivity3",
)
plt.suptitle("Attenuation correction using pycwr", fontsize=16)
plt.show()小结
前两种方法计算结果相似，后一种没有可相比，但是清除了不少杂波
大家可用更多文件去试试
pyart库在变量导入导出方面可说是独一无二 希望国内雷达库可以学习学习
