两个气象雷达库的简单使用

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



pip install --user --upgrade cinrad

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进行对数据细化处理

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[ns] 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

#抄的官方文档，更多可以自己细化，懒狗记录下
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)

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)

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()
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)

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()
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()
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

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

import cinrad
from cinrad.io import CinradReader, StandardData
f = CinradReader(PATH)
f.available_product(0)

['REF', 'VEL', 'SW', 'azimuth', 'RF']

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])

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

#反射率，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])

/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)

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

/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)

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")

#垂直积分含水量
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

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.
  shading=shading)

# 回波顶
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)

#组合反射率
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)

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[Any, NoneType] = None, norm: Union[Any, NoneType] = None, cmap: Union[Any, NoneType] = None, nlabel: Union[int, NoneType] = None, label: Union[List[str], NoneType] = None, dpi: Union[int, float] = 350, highlight: Union[str, List[str], NoneType] = None, coastline: bool = False, extent: Union[List[Union[int, float]], NoneType] = None, section: Union[xarray.core.dataset.Dataset, NoneType] = 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[str, NoneType] = None)
 |      Call self as a function.
 |  
 |  __init__(self, data: xarray.core.dataset.Dataset, fig: Union[Any, NoneType] = None, norm: Union[Any, NoneType] = None, cmap: Union[Any, NoneType] = None, nlabel: Union[int, NoneType] = None, label: Union[List[str], NoneType] = None, dpi: Union[int, float] = 350, highlight: Union[str, List[str], NoneType] = None, coastline: bool = False, extent: Union[List[Union[int, float]], NoneType] = None, section: Union[xarray.core.dataset.Dataset, NoneType] = 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[int, float] = 0, **kwargs)
 |      Draw grid lines on cartopy axes
 |  
 |  plot_cross_section(self, data: xarray.core.dataset.Dataset, ymax: Union[int, NoneType] = None, linecolor: Union[str, NoneType] = None, interpolate: bool = True)
 |      Plot cross section data below the PPI plot.
 |  
 |  plot_range_rings(self, _range: Union[int, float, list], color: str = 'white', linewidth: Union[int, float] = 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')