两个气象雷达库的简单使用
前置任务：安装与升级!pip install pycwr -i <a href="https://pypi.mirrors.ustc.edu.cn/simple" target="_blank">https://pypi.mirrors.ustc.edu.cn/simple</a>
pip install --user --upgrade cinrad
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)
#提取变量，可以尝试更改数字查看，改变的是仰角
dbz=PRD.fields["dBZ"].values
#是data array格式，可以学习xarray进行对数据细化处理
Warning: ecCodes 2.21.0 or higher is recommended. You are running version 2.14.1
<xarray.Dataset>
Dimensions: (time: 361, range: 920)
Coordinates:
azimuth (time) float64 319.5 320.5 321.4 322.4 ... 317.2 318.2 319.2
elevation(time) float64 1.527 1.527 1.527 1.527 ... 1.527 1.527 1.527
x (time, range) float64 -162.3 -324.5 ... -1.498e+05 -1.499e+05
y (time, range) float64 190.1 380.1 570.2 ... 1.738e+05 1.74e+05
z (time, range) float64 145.5 152.1 158.8 ... 9.363e+03 9.377e+03
lat (time, range) float64 32.19 32.19 32.2 32.2 ... 33.74 33.74 33.74
lon (time, range) float64 118.7 118.7 118.7 ... 117.1 117.1 117.1
* range (range) float64 250.0 500.0 750.0 ... 2.295e+05 2.298e+05 2.3e+05
* time (time) datetime64 2020-06-12T05:49:35.324000 ... 2020-06-1...
Data variables:
dBZ (time, range) float32 nan nan nan -5.5 -5.5 ... nan nan nan nan
V (time, range) float64 nan nan nan nan 4.5 ... nan nan nan nan nan
W (time, range) float64 nan nan nan nan 8.5 ... nan nan nan nan nan
1.1.1 无地图版PPI DBZ
#抄的官方文档，更多可以自己细化，懒狗记录下
fig, ax = plt.subplots()
graph = Graph(PRD)
graph.plot_ppi(ax, 0, "dBZ", cmap="CN_ref") ## 0代表第一层, dBZ代表反射率产品
graph.add_rings(ax, )
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.
1.1.2 有地图版PPI DBZ
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()
plt.show()
/opt/conda/lib/python3.7/site-packages/cartopy/mpl/geoaxes.py:1598: 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.
shading=shading)
1.1.3 RHI
In :
fig, ax = plt.subplots()
graph = Graph(PRD)
graph.plot_rhi(ax, 0, field_name="dBZ", cmap="CN_ref", clabel="Radar Reflectivity")
ax.set_ylim() #设置rhi的高度范围 （units：km）
ax.set_xlabel("distance from radar (km)", fontsize=14)
ax.set_ylabel("Height (km)", fontsize=14)
plt.tight_layout()
plt.show()/opt/conda/lib/python3.7/site-packages/pycwr/draw/RadarPlot.py:103: 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.
norm=norm, shading='auto', **kwargs)
说明一下，普通的业务用双偏振雷达是不开RHI模式的，所以画成这鸟样
不过这个数据是单偏振格式的，双偏振的数据会多几个变量
什么，你说RHI是什么
当然是垂直扫描
1.1.4 剖面图
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()
ax.set_xlim()
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()
plt.show()
1.1.5 CAPPI
CAPPI是等高平面位置显示产品
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()
plt.show()
<xarray.Dataset>
Dimensions: (x_cappi_3000: 301, y_cappi_3000: 301)
Coordinates:
* x_cappi_3000(x_cappi_3000) int64 -150000 -149000 -148000 ... 149000 150000
* y_cappi_3000(y_cappi_3000) int64 -150000 -149000 -148000 ... 149000 150000
Data variables:
CAPPI_3000 (x_cappi_3000, y_cappi_3000) float64 nan nan nan ... nan nan
1.1.6 组合反射率（CR)
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()
plt.show()
/opt/conda/lib/python3.7/site-packages/pycwr/core/NRadar.py:166: RuntimeWarning: All-NaN slice encountered
radar_height, GridX.astype(np.float64), GridY.astype(np.float64), -999.)
<xarray.Dataset>
Dimensions: (x_cappi_3000: 301, y_cappi_3000: 301, x_cr: 301, y_cr: 301)
Coordinates:
* x_cappi_3000(x_cappi_3000) int64 -150000 -149000 -148000 ... 149000 150000
* y_cappi_3000(y_cappi_3000) int64 -150000 -149000 -148000 ... 149000 150000
* x_cr (x_cr) int64 -150000 -149000 -148000 ... 148000 149000 150000
* y_cr (y_cr) int64 -150000 -149000 -148000 ... 148000 149000 150000
Data variables:
CAPPI_3000 (x_cappi_3000, y_cappi_3000) float64 nan nan nan ... nan nan
CR (x_cr, y_cr) float64 nan nan nan nan ... 26.3 22.94 19.25
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函数参数依次是仰角，半径，变量.。这时候就有不懂的小伙伴问了，仰角有哪些？当然是自己看拉。
In :
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')
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()/opt/conda/lib/python3.7/site-packages/cartopy/mpl/geoaxes.py:1598: 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.
shading=shading)
In :
#速度产品
VEL = f.get_data(3, 230, 'VEL')
fig = cinrad.visualize.PPI(VEL, dpi=50,add_city_names=True)
fig.plot_range_rings()
/opt/conda/lib/python3.7/site-packages/cartopy/mpl/geoaxes.py:1598: 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.
shading=shading)
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)
/home/mw/.local/lib/python3.7/site-packages/cinrad/io/level2.py:430: UserWarning: Requested data range exceed max range in this tilt
warnings.warn("Requested data range exceed max range in this tilt")
1.2.4 其他可计算的变量
In :#垂直积分含水量
vil = cinrad.calc.quick_vil(rl)
print(vil)
<xarray.Dataset>
Dimensions: (azimuth: 361, distance: 230)
Coordinates:
* azimuth (azimuth) float64 0.0 0.01745 0.03491 ... 6.248 6.266 6.283
* distance (distance) float64 1.0 2.0 3.0 4.0 ... 227.0 228.0 229.0 230.0
Data variables:
VIL (azimuth, distance) float64 nan nan nan nan ... nan nan nan nan
longitude(azimuth, distance) float64 118.7 118.7 118.7 ... 118.7 118.7
latitude (azimuth, distance) float64 32.2 32.21 32.22 ... 34.25 34.26
Attributes:
elevation: 0
range: 230
scan_time: 2020-06-12 05:48:04
site_code: Z9250
site_name: 南京
site_longitude: 118.69777777777777
site_latitude: 32.19083333333333
tangential_reso:1.0
nyquist_vel: 8.43
task: VCP21
In :fig = cinrad.visualize.PPI(vil, dpi=50)
/opt/conda/lib/python3.7/site-packages/cartopy/mpl/geoaxes.py:1598: 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.
s
In :
# 回波顶
et = cinrad.calc.quick_et(rl)
fig = cinrad.visualize.PPI(et, dpi=50)
/opt/conda/lib/python3.7/site-packages/cartopy/mpl/geoaxes.py:1598: 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.
shading=shading)
In :
#组合反射率
cr = cinrad.calc.quick_cr(rl)
fig = cinrad.visualize.PPI(cr, dpi=50)
/home/mw/.local/lib/python3.7/site-packages/cinrad/calc.py:110: RuntimeWarning: All-NaN axis encountered
cr = np.nanmax(r_data, axis=0)
1.2.5 透明图层设置
help(cinrad.visualize.PPI)
Help on class PPI in module cinrad.visualize.ppi:
class PPI(builtins.object)
|PPI(data: xarray.core.dataset.Dataset, fig: Union = None, norm: Union = None, cmap: Union = None, nlabel: Union = None, label: Union, NoneType] = None, dpi: Union = 350, highlight: Union, NoneType] = None, coastline: bool = False, extent: Union], NoneType] = None, section: Union = None, style: str = 'black', add_city_names: bool = False, plot_labels: bool = True, **kwargs)
|
|Create a figure plotting plan position indicator
|
|By default, norm, cmap, and colorbar labels will be determined by the
|data type.
|
|Args:
| data (xarray.Dataset): The data to be plotted.
|
| fig (matplotlib.figure.Figure): The figure to plot on. Optional.
|
| norm (matplotlib.colors.Normalize): Customized norm data. Optional.
|
| cmap (matplotlib.colors.Colormap): Customized colormap. Optional.
|
| nlabel (int): Number of labels on the colorbar. Optional.
|
| dpi (int): DPI of the figure. Optional.
|
| highlight (str, list(str)): Areas to be highlighted. Optional.
|
| coastline (bool): Plot coastline on the figure if set to True. Default False.
|
| extent (list(float)): The extent of figure. Optional.
|
| add_city_names (bool): Label city names on the figure if set to True. Default True.
|
| plot_labels (bool): Text scan information on the side of the plot. Default True.
|
|Methods defined here:
|
|__call__(self, fpath: Union = None)
| Call self as a function.
|
|__init__(self, data: xarray.core.dataset.Dataset, fig: Union = None, norm: Union = None, cmap: Union = None, nlabel: Union = None, label: Union, NoneType] = None, dpi: Union = 350, highlight: Union, NoneType] = None, coastline: bool = False, extent: Union], NoneType] = None, section: Union = None, style: str = 'black', add_city_names: bool = False, plot_labels: bool = True, **kwargs)
| Initialize self.See help(type(self)) for accurate signature.
|
|gridlines(self, draw_labels: bool = True, linewidth: Union = 0, **kwargs)
| Draw grid lines on cartopy axes
|
|plot_cross_section(self, data: xarray.core.dataset.Dataset, ymax: Union = None, linecolor: Union = None, interpolate: bool = True)
| Plot cross section data below the PPI plot.
|
|plot_range_rings(self, _range: Union, color: str = 'white', linewidth: Union = 0.5, **kwargs)
| Plot range rings on PPI plot.
|
|storm_track_info(self, filepath: str)
| Add storm tracks from Nexrad Level III (PUP) STI product file
|
|----------------------------------------------------------------------
|Data descriptors defined here:
|
|__dict__
| dictionary for instance variables (if defined)
|
|__weakref__
| list of weak references to the object (if defined)
|
|----------------------------------------------------------------------
|Data and other attributes defined here:
|
|data_crs = <cartopy.crs.PlateCarree object>
#绘制透明图层
fig = cinrad.visualize.PPI(cr, dpi=50,style='transparent')
小结
1.pycwr自定义度高，绘图相对好看
2.cinrad入门快，代码简便
完整代码与文件在此
文章来源于微信公众号：气ython风雨
页:
[1]