python绘图 | 气象雷达入门级讲解&多种雷达图像可视化方法
python绘图 | 气象雷达入门级讲解&多种雷达图像可视化方法
气象雷达简介
气象雷达是专门用于大气探测的雷达。它是一种主动式微波大气遥感设备。 气象雷达是气象观测的重要设备，特别是在突发性、灾害性的监测、预报和警报中具有极为重要的作用，是用于小尺度天气系统（如台风和暴雨云系）的主要探测工具之一。 在国内，我们最常见到和使用的气象雷达，是新一代多普勒天气雷达(CINRAD)。我们在气象局之类建筑楼顶上见到的那些球形建筑，大都属于这一种雷达。这种雷达可以探测反射率因子、多普勒径向速度、谱宽等基本气象要素，从而为短临尺度上的天气预报和预警提供数据支撑。特别是雷达反射率数据，因为其与强对流天气系统直接相关，最常被大家使用。 雷达数据在日常业务科研中的应用非常多，比如雷达数据可以用于数值模式同化中，为数值模式提供一个更加准确的初始场；基于雷达反射率数据的雷达短临预报系统可以预报未来2小时内，雷达探测范围内的强对流天气。例如，眼控科技自主研发的基于深度学习的AI对流临近预报系统就是利用雷达反射率数据，对未来两小时之内强对流天气，进行准确的预报。看了一下，下面的这个预报效果确实很好。
气象雷达的工作原理
要想了解气象雷达的数据结构，首先要知道雷达是如何工作的，为此我特意又把相关知识学习了一遍，感兴趣的小伙伴也可以看一看：学习资源 | 雷达气象学课程（南信大官方） 雷达主要有以下几种扫描形式：
一般来说，业务中雷达的常规扫描方式是VPPI，在扫描时，雷达天线从自体扫模式中最低仰角启动，并以固定仰角从零度方位角（在多普勒天气雷达工作过程中，规定正北方为0°方位角，正东方为90°方位角，天线与水平面平行为0°仰角，与水平面垂直时为90°仰角。）开始沿顺时针方向转动天线，逐方位角完成单层仰角的扫描，在不断发射和接收电磁波的同时记录反射率因子、径向多普勒速度、谱宽等生成基本雷达产品所必需的气象要素。 完成单层仰角扫描之后，雷达天线抬升至对应体扫模式中的下一仰角继续重复上述过程，直到完成所有规定仰角的扫描。 在一次扫描完成之后，我们实际获得的是不同仰角的一整圈的气象要素。如果从三维空间角度来说，就是一个个以雷达站点为定点的，不同倾斜角的圆锥曲面，共同构成了雷达的三维空间观测。下图是其中一个仰角扫描后的示意图。
使用python读取CINRAD雷达数据
一般而言，我们需要使用pyart或wradlib,甚至直接按照雷达数据格式读取二进制文件，然后进行一些必要的处理，才可以顺利读取雷达数据。
幸运的是，对于国内最常见的CINRAD雷达数据，@CyanideCN老师已经开发了一个库PyCINRAD，使得我们可以很简单的对CINRAD数据进行处理和可视化。利用PyCINRAD,我们可以很方便的将雷达基数据读取为xarray.Dataset
#加载所需包
import math
import time
from itkwidgets import view
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import numpy as np
from matplotlib import colors
from scipy.interpolate import griddata
import scipy.misc
import cartopy.crs as ccrs
import cartopy.feature as cfeat
import cinrad
import xarray as xr
from cartopy.mpl.gridliner import LATITUDE_FORMATTER, LONGITUDE_FORMATTER
from cinrad.io import CinradReader, StandardData
from cinrad.visualize import PPI,Section
from mayavi import mlab
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D
import pandas as pd
import plotly.express as px
%matplotlib inline
使用cinrad.io.CinradReader打开二进制雷达基数据文件
f = CinradReader("Z_RADR_I_Z9250_20200612054800_O_DOR_SA_CAP.bin")
查看包含的产品类型
f.available_product(0)
获取雷达扫描仰角
f.get_elevation_angles()
查看数据类型
ele = 0 #选择第1个仰角
radius = 400 #绘制图像的范围大小，单位km
r = f.get_data(ele,radius,"REF") #选择反射率数据
r
通过上面一段代码我们就可以的了解雷达数据的存储结构了。从中可以直观的看出每一层的雷达反射率数据，实际上是由两个维度组成的。一个是方位角(azimuth)，它表示了每一次雷达发出电磁波时，其对准的角度；另一个是距离（distance)，它实际上代表了雷达的每一个观测点与雷达站之间的距离。实际上，我们获取的每一层雷达反射率数据，是由约360条射线，每条射线上每隔1公里的点上的气象要素拼成的一个三维的圆锥面。在这里，PyCinrad库同时计算出了圆锥面上每个点的具体经纬度值和高度值，有了这些值，可以帮助我们更方便的在二维和三维的笛卡尔坐标下进行可视化。
二维可视化
得到一个xarray.Dataset后，自然而然地可以简单地使用cartopy对其进行可视化。
ref = r.REF
proj = ccrs.PlateCarree()#创建投影，选择cartopy的platecarree投影
fig = plt.figure(figsize=(6,6)) #创建页面，可以自己选择大小
ax = fig.subplots(1, 1, subplot_kw={'projection': proj}) #子图
#国界、海岸线、河流、湖泊
ax.add_feature(cfeat.BORDERS.with_scale('50m'), linewidth=0.8, zorder=1)
ax.add_feature(cfeat.COASTLINE.with_scale('50m'), linewidth=0.6, zorder=1)
ax.add_feature(cfeat.RIVERS.with_scale('50m'), zorder=1)
ax.add_feature(cfeat.LAKES.with_scale('50m'), zorder=1)
gl = ax.gridlines(crs=ccrs.PlateCarree(),draw_labels=True,
linewidth=1.2, color='k', alpha=0.5, linestyle='--')
gl.top_labels = False #关闭顶端标签
gl.right_labels = False #关闭右侧标签
gl.xformatter = LONGITUDE_FORMATTER #x轴设为经度格式
gl.yformatter = LATITUDE_FORMATTER #y轴设为纬度格式
color_max=int(ref.max())+1
color_min=int(ref.min())
n_gap=(color_max-color_min)/20
#cmap
levels = np.arange(color_min,color_max+n_gap,n_gap)
plt.contourf( r.longitude, r.latitude,ref,levels=levels,cmap='jet')
#colorbar 左 下 宽 高
l = 0.92
b = 0.23
w = 0.025
h = 0.55
#对应 l,b,w,h；设置colorbar位置；
rect = [l,b,w,h]
cbar_ax = fig.add_axes(rect)
plt.colorbar( cax=cbar_ax)
plt.show()
PPI可视化
也可以使用PyCinrad库中自带的一些可视化的函数，对PPI直接进行可视化,首先定义一个绘制ppi的函数。
tip:使用cinrad.visualize.PPI直接可视化PPI，有时需要运行两次cell才能呈现黑色背景的图片。
def ppi_plot(data):
fig = PPI(data,dpi=75,add_city_names=True)
fig.plot_range_rings(radius, color='white', linewidth=1.0) #绘制圆圈
for i in range(0,radius-30,50):
fig.plot_range_rings(i, color='white', linewidth=1.0) #绘制圆圈
# 设置经纬度
liner = fig.geoax.gridlines(draw_labels=True,linewidth=2, color='gray', alpha=0.5, linestyle='--') #设置经纬度虚线
liner.top_labels = False
liner.right_labels = False
liner.xformatter = LONGITUDE_FORMATTER
liner.yformatter = LATITUDE_FORMATTER
liner.xlabel_style = {'size': 18, 'color': 'red', 'weight': 'bold'} #设置经纬度颜色字号等
liner.ylabel_style = {'size': 18, 'color': 'red', 'weight': 'bold'}
绘制雷达发射率图
ppi_plot(r)
除了雷达反射率以外，PyCinrad还有计算垂直累计液态水(vertically integrated liquid,VIL)，回波顶高(echo tops,ET)，组合反射率(composite reflectivity,CR)等方法
绘制回波顶高
et = cinrad.calc.quick_et([f.get_data(i, 400, 'REF') for i in f.angleindex_r])
ppi_plot(et)
绘制垂直累计液态水
vil = cinrad.calc.quick_vil([f.get_data(i, 400, 'REF') for i in f.angleindex_r])
ppi_plot(vil)
绘制组合反射率
cr = cinrad.calc.quick_cr([f.get_data(i, 400, 'REF') for i in f.angleindex_r])
RHI可视化
由于雷达日常业务数据没有进行RHI扫描，所以要想获得雷达的垂直剖面可以选择cinrad.calc.VCS来计算雷达垂直回波强度后再进行绘图。
rl = [f.get_data(i, 400, 'REF') for i in f.angleindex_r]
vcs = cinrad.calc.VCS(rl)
sec = vcs.get_section(start_cart=(118.5, 32.5), end_cart=(118.5, 34)) # 传入经纬度坐标
fig = Section(sec)
绘制组合图
PyCinrad还支持同时绘制PPI和RHI的组合图
rl = [f.get_data(i, 400, 'REF') for i in f.angleindex_r]
vcs = cinrad.calc.VCS(rl)
sec = vcs.get_section(start_cart=(118.0, 32.5), end_cart=(119.0, 33.5)) # 传入经纬度坐标
fig = PPI(rl[0],dpi=75,add_city_names=True)
fig.settings['is_inline'] = False
fig.plot_range_rings(radius, color='white', linewidth=1.0) #绘制圆圈
for i in range(0,radius-30,50):
fig.plot_range_rings(i, color='white', linewidth=1.0) #绘制圆圈
fig.plot_cross_section(sec) #绘制垂直剖面
liner = fig.geoax.gridlines(draw_labels=True,linewidth=2, color='gray', alpha=0.5, linestyle='--') #设置经纬度虚线
liner.top_labels = False
liner.right_labels = False
liner.xformatter = LONGITUDE_FORMATTER
liner.yformatter = LATITUDE_FORMATTER
liner.xlabel_style = {'size': 18, 'color': 'red', 'weight': 'bold'} #设置经纬度颜色字号等
liner.ylabel_style = {'size': 18, 'color': 'red', 'weight': 'bold'}
三维可视化
雷达数据并非分布在一个曲面上，在经度-维度-高度的笛卡尔坐标系下，一层的PPI数据在三维空间中呈圆锥面分布，因此可以对其进行三维的可视化
matmatplotlib三维静态可视化
fig = plt.figure()
ax = Axes3D(fig) #创建三维绘图空间
X = r.longitude.values.flatten() #读取ppi中经度纬度高度和反射率数值，并转化成一维
Y = r.latitude.values.flatten()
Z = r.height.values.flatten()
value = r.REF.values.flatten()
ax.scatter(X, Y, Z, c=value) #绘制散点图
plotly三维动态可视化
需要在jupyter中才可以进行交互
#取出经度、纬度、高度、反射率
X = r.longitude.values
Y = r.latitude.values
Z = r.height.values
value = r.REF.values
#新建dataframe，保存散点信息
df = pd.DataFrame({"lon":X.flatten(),"lat":Y.flatten(),"height":Z.flatten()/110,"dbz":value.flatten()})
df = df.dropna() #去除nan的点
#生成图像可以拖动
fig = px.scatter_3d(df, x='lon', y='lat', z='height',color='dbz')
fig.show()
最后想说的是，这篇文章算是一个很简短的介绍，在这个领域只能算是初窥门径。雷达在气象领域的应用非常丰富，感兴趣的朋友可以深入的学习一下雷达的相关知识，打牢基础总是没错的，这个是南信大官方的课程链接，有需求请自取：学习资源 | 雷达气象学课程（南信大官方）