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()
最後想說的是，這篇文章算是一個很簡短的介紹，在這個領域只能算是初窺門徑。雷達在氣象領域的應用非常豐富，感興趣的朋友可以深入的學習一下雷達的相關知識，打牢基礎總是沒錯的，這個是南信大官方的課程連結，有需求請自取：學習資源 | 雷達氣象學課程（南信大官方）
獲取方法
示例代碼（py及ipynb）、圖片獲取方式：在「好奇心Log」公眾號後臺留言：pycinrad示例。
除此以外大家也可以去「和鯨社區」看到整個示例。
雷達示例數據暫未直接公開，若有需要可將本篇文章轉至朋友圈截圖發至「好奇心Log」公眾號後臺。
整篇文章是在冬青老師的材料上加以完善，算是比較完整的入門，如果你覺得有幫助的話歡迎轉發，歡迎打賞。打賞代收，最終轉予冬青老師。