两个库读取雷达数据

pycwr和pycinrad

前置任务：安装与升级

!pip install pycwr -i https://pypi.mirrors.ustc.edu.cn/simple

1.1 pycwr

项目地址：

https://pycwr.readthedocs.io/en/latest/draw.html

导入库和看变量

from pycwr.io import read_auto
import numpy as np
import matplotlib.pyplot as plt
from pycwr.draw.RadarPlot import Graph
from pycwr.draw.RadarPlot import GraphMap
import cartopy.crs as ccrs
PATH='/home/mw/input/data5692/Z_RADR_I_Z9250_20200612054800_O_DOR_SA_CAP.bin'
#读取
PRD = read_auto(PATH)
#查看变量
print(PRD.fields[1])
#提取变量，可以尝试更改数字查看，改变的是仰角
dbz=PRD.fields[1]["dBZ"].values

#是data array格式，可以学习xarray训练营进行对数据细化处理

1.1.1 无地图版PPI DBZ

findfont: Font family ['Arial'] not found. Falling back to DejaVu Sans.

#抄的官方文档，更多可以自己细化，懒狗记录下
fig, ax = plt.subplots()
graph = Graph(PRD)
graph.plot_ppi(ax, 0, "dBZ", cmap="CN_ref") ## 0代表第一层, dBZ代表反射率产品
graph.add_rings(ax, [0, 50, 100, 150, 200, 250, 300])
ax.set_title("PPI Plot", fontsize=16)
ax.set_xlabel("Distance From Radar In East (km)", fontsize=14)
ax.set_ylabel("Distance From Radar In North (km)", fontsize=14)
plt.show()
/opt/conda/lib/python3.7/site-packages/pycwr/draw/RadarPlot.py:57: UserWarning: The input coordinates to pcolormesh are interpreted as cell centers, but are not monotonically increasing or decreasing. This may lead to incorrectly calculated cell edges, in which case, please supply explicit cell edges to pcolormesh.
  zorder=0, norm=norm, shading='auto', **kwargs)
findfont: Font family ['Arial'] not found. Falling back to DejaVu Sans.

1.1.2 有地图版PPI DBZ

plt.show()

from pycwr.io import read_auto
import matplotlib.pyplot as plt
from pycwr.draw.RadarPlot import GraphMap
import cartopy.crs as ccrs
 
ax = plt.axes(projection=ccrs.PlateCarree())
graph = GraphMap(PRD, ccrs.PlateCarree())
graph.plot_ppi_map(ax, 0, "dBZ", cmap="CN_ref") ## 0代表第一层, dBZ代表反射率产品，cmap
ax.set_title("PPI Plot with Map", fontsize=16)
plt.tight_layout()

1.1.3 RHI

plt.show()

fig, ax = plt.subplots()
graph = Graph(PRD)
graph.plot_rhi(ax, 0, field_name="dBZ", cmap="CN_ref", clabel="Radar Reflectivity")
ax.set_ylim([0, 6]) #设置rhi的高度范围 （units：km）
ax.set_xlabel("distance from radar (km)", fontsize=14)
ax.set_ylabel("Height (km)", fontsize=14)
plt.tight_layout()

说明一下，普通的业务用双偏振雷达是不开RHI模式的，所以画成这鸟样 不过这个数据是单偏振格式的，双偏振的数据会多几个变量 什么，你说RHI是什么 当然是垂直扫描

1.1.4 剖面图

plt.show()

from pycwr.io import read_auto
import matplotlib.pyplot as plt
from pycwr.draw.RadarPlot import Graph
 
fig, ax = plt.subplots()
graph = Graph(PRD)
graph.plot_vcs(ax, (0,0), (150, 0), "dBZ", cmap="copy_pyart_NWSRef")
ax.set_ylim([0, 10])
ax.set_xlim([0, 230])
ax.set_ylabel("Height (km)", fontsize=14)
ax.set_xlabel("Distance From Section Start (Uints:km)", fontsize=14)
ax.set_title("VCS Plot", fontsize=16)
plt.tight_layout()

1.1.5 CAPPI

CAPPI是等高平面位置显示产品

plt.show()

from pycwr.draw.RadarPlot import plot_xy
x1d = np.arange(-150000, 150001, 1000) ##x方向1km等间距， -150km～150km范围
y1d = np.arange(-150000, 150001, 1000) ##y方向1km等间距， -150km～150km范围
PRD.add_product_CAPPI_xy(XRange=x1d, YRange=y1d, level_height=3000) ##level height units:meters
print(PRD.product)
grid_x, grid_y = np.meshgrid(x1d, y1d, indexing="ij")
fig, ax = plt.subplots()
plot_xy(ax, grid_x, grid_y, PRD.product.CAPPI_3000) ##画图显示
ax.set_xlabel("Distance From Radar In East (km)", fontsize=14)
ax.set_ylabel("Distance From Radar In North (km)", fontsize=14)
plt.tight_layout()

1.1.6 组合反射率（CR)

plt.show()

x1d = np.arange(-150000, 150001, 1000) ##x方向1km等间距， -150km～150km范围
y1d = np.arange(-150000, 150001, 1000) ##y方向1km等间距， -150km～150km范围
PRD.add_product_CR_xy(XRange=x1d, YRange=y1d)
print(PRD.product)
grid_x, grid_y = np.meshgrid(x1d, y1d, indexing="ij")
fig, ax = plt.subplots()
plot_xy(ax, grid_x, grid_y, PRD.product.CR) ##画图显示
ax.set_xlabel("Distance From Radar In East (km)", fontsize=14)
ax.set_ylabel("Distance From Radar In North (km)", fontsize=14)
plt.tight_layout()

1.2 pycinrad

项目地址为https://gitee.com/CyanideCN/PyCINRAD/

1.2.1 看变量

import cinrad
from cinrad.io import CinradReader, StandardData
f = CinradReader(PATH)
f.available_product(0)
['REF', 'VEL', 'SW', 'azimuth', 'RF']

f.get_data函数参数依次是仰角，半径，变量.。这时候就有不懂的小伙伴问了，仰角有哪些？当然是自己看拉。

f.get_elevation_angles()

array([ 0.5657959 , 0.51635742, 1.44470215, 1.52709961, 2.31811523, 3.24645996, 4.26818848, 5.9765625 , 9.83825684, 14.51843262, 19.46777344])

按照python的从零开始的特点，填入自己需要的仰角顺位即可

REF = f.get_data(2, 230, 'REF')
REF

1.2.2 PPI

#反射率，add_city_names加城市名字，fig.plot_range_rings加环
fig = cinrad.visualize.PPI(REF, dpi=50,add_city_names=True)
fig.plot_range_rings([50, 100, 150, 200, 230])

fig.plot_range_rings([50, 100, 150, 200, 230])

#速度产品
VEL = f.get_data(3, 230, 'VEL')
fig = cinrad.visualize.PPI(VEL, dpi=50,add_city_names=True)

1.2.3 剖面图(不是rhi)

from cinrad.visualize import Section
rl = list(f.iter_tilt(230, 'REF'))
vcs = cinrad.calc.VCS(rl)
sec = vcs.get_section(start_cart=(118, 32.5), end_cart=(118, 34)) # 传入经纬度坐标
fig = Section(sec)

1.2.4 其他可计算的变量

#垂直积分含水量
vil = cinrad.calc.quick_vil(rl)
print(vil)

回波顶

et = cinrad.calc.quick_et(rl)
fig = cinrad.visualize.PPI(et, dpi=50)

组合反射率

cr = cinrad.calc.quick_cr(rl)
fig = cinrad.visualize.PPI(cr, dpi=50)

1.2.5 透明图层设置

#绘制透明图层
fig = cinrad.visualize.PPI(cr, dpi=50,style='transparent')

小结

1. pycwr自定义度高，绘图相对好看

2. cinrad入门快，代码简便

更多具体图像访问以下链接

https://www.heywhale.com/mw/project/6468caf98cc5d51b4c7bbc1a?shareby=620dbd752e0e510017d2f245#

我分享了一个项目给你《两个气象雷达库的简单使用》，快来看看吧

因为换了新编辑器，不知道怎么搞图床，懒得一个个贴图