Python CINRAD雷達數據解析與可視化工具

一、CINRAD雷達的基本概述

在介紹Python CINRAD雷達數據解析與可視化工具之前，我們首先需要了解一下CINRAD雷達的基本概念。CINRAD雷達是中國自主研發的一種天氣雷達，被廣泛應用於氣象、公路、航空、水利等領域。CINRAD雷達可以獲取雷達反射率因子、速度以及譜寬等數據，這些數據對於氣象預報和災害預警具有重要的意義。

CINRAD雷達系統主要由雷達機體、天線、發射機和接收機等組成。雷達機體位於地面附近，通過雷達天線向周圍發射電磁波，當電磁波遇到雲層等物體時，會被反射回來並被接收機接收。

二、Python CINRAD雷達數據解析的基本流程

Python CINRAD雷達數據解析的基本流程分為數據讀取、頭文件解析、數據解壓和數據可視化幾個步驟。

1. 數據讀取

我們可以使用Python中的open函數來打開CINRAD雷達數據文件，並用read()函數讀取文件內容。讀取的文件內容為二進位格式的數據，我們需要進行解析。

import numpy as np

filename = 'Z_RADR_I_Z9857_20190809134600_O_DOR_SA_CAP.bin'
with open(filename, 'rb') as f:
    data = f.read()

2. 頭文件解析

在CINRAD雷達數據中，頭文件中包含了一些與數據本身相關的元信息。我們需要使用Python來解析這些信息，例如雷達反射率因子、速度和譜寬等數據的解析度、掃描方式以及探測距離等信息。

header = np.frombuffer(data, dtype='>i4', count=45)
element_number = header[3] #雷達反射率因子、速度和譜寬等數據的元素數
range_number = header[4] #距離元素數
azimuth_number = header[6] #方位角元素數
range_scale = header[7] / 1000.0 #距離解析度
azimuth_scale = header[8] / 1000.0 #方位角解析度
elevation_number = 1
elevation_scale = 1.0

3. 數據解壓

在CINRAD雷達數據中，每個元素都被壓縮為兩個位元組。我們需要使用Python對數據進行解壓縮，並轉換為真實的數據。雷達反射率因子、速度和譜寬等數據都需要進行解壓縮，並乘以各自的比例因子得到真實的數據值。

data_offset = header[11] #數據起始點偏移（byte）
data_block = data[data_offset:]
data_array = np.frombuffer(data_block, dtype='>u2')
data_array = data_array.astype('float32')
data_array[data_array == 0] = np.nan
reflectivity = (data_array[0:element_number*range_number*azimuth_number]
               .reshape((azimuth_number, range_number, element_number))
               .swapaxes(0,1))
reflectivity = reflectivity * header[19] + header[20]
velocity = (data_array[element_number*range_number*azimuth_number:
                       element_number*range_number*azimuth_number*2]
           .reshape((azimuth_number, range_number, element_number))
           .swapaxes(0,1))
velocity = velocity * header[21] + header[22]
spectrum_width = (data_array[element_number*range_number*azimuth_number*2:
                              element_number*range_number*azimuth_number*3]
                  .reshape((azimuth_number, range_number, element_number))
                  .swapaxes(0,1))
spectrum_width = spectrum_width * header[23] + header[24]

4. 數據可視化

最後，我們可以使用Python的數據可視化工具例如matplotlib或者basemap將CINRAD雷達數據進行可視化。我們可以將雷達反射率因子、速度和譜寬等數據繪製成不同的圖像，以便於我們對數據進行分析。

import matplotlib.pyplot as plt

range_array = np.arange(0, range_number) * range_scale
azimuth_array = np.arange(0, azimuth_number) * azimuth_scale
azimuth_mesh, range_mesh = np.meshgrid(azimuth_array, range_array)

#雷達反射率因子
plt.figure()
plt.pcolormesh(azimuth_mesh, range_mesh/1000.0,
               reflectivity.transpose(), cmap='jet', vmin=-32, vmax=94)
plt.xlabel('Azimuth (degree)')
plt.ylabel('Range (km)')
plt.title('Reflectivity')
plt.colorbar()

#速度
plt.figure()
plt.pcolormesh(azimuth_mesh, range_mesh/1000.0,
               velocity.transpose(), cmap='seismic', vmin=-30, vmax=30)
plt.xlabel('Azimuth (degree)')
plt.ylabel('Range (km)')
plt.title('Velocity')
plt.colorbar()

#譜寬
plt.figure()
plt.pcolormesh(azimuth_mesh, range_mesh/1000.0,
               spectrum_width.transpose(), cmap='plasma', vmin=0, vmax=8)
plt.xlabel('Azimuth (degree)')
plt.ylabel('Range (km)')
plt.title('Spectrum Width')
plt.colorbar()

plt.show()

三、Python CINRAD雷達數據解析與可視化實例

通過以上的步驟，我們可以使用Python CINRAD雷達數據解析與可視化工具來對CINRAD雷達數據進行處理。下面我們以解析通州站雷達2019年8月9日13時46分的數據為例，進行數據可視化。

import numpy as np
import matplotlib.pyplot as plt

filename = 'Z_RADR_I_Z9857_20190809134600_O_DOR_SA_CAP.bin'
with open(filename, 'rb') as f:
    data = f.read()

header = np.frombuffer(data, dtype='>i4', count=45)
element_number = header[3] #雷達反射率因子、速度和譜寬等數據的元素數
range_number = header[4] #距離元素數
azimuth_number = header[6] #方位角元素數
range_scale = header[7] / 1000.0 #距離解析度
azimuth_scale = header[8] / 1000.0 #方位角解析度
elevation_number = 1
elevation_scale = 1.0

data_offset = header[11] #數據起始點偏移（byte）
data_block = data[data_offset:]
data_array = np.frombuffer(data_block, dtype='>u2')
data_array = data_array.astype('float32')
data_array[data_array == 0] = np.nan

reflectivity = (data_array[0:element_number*range_number*azimuth_number]
               .reshape((azimuth_number, range_number, element_number))
               .swapaxes(0,1))
reflectivity = reflectivity * header[19] + header[20]

velocity = (data_array[element_number*range_number*azimuth_number:
                       element_number*range_number*azimuth_number*2]
           .reshape((azimuth_number, range_number, element_number))
           .swapaxes(0,1))
velocity = velocity * header[21] + header[22]

spectrum_width = (data_array[element_number*range_number*azimuth_number*2:
                              element_number*range_number*azimuth_number*3]
                  .reshape((azimuth_number, range_number, element_number))
                  .swapaxes(0,1))
spectrum_width = spectrum_width * header[23] + header[24]

range_array = np.arange(0, range_number) * range_scale
azimuth_array = np.arange(0, azimuth_number) * azimuth_scale
azimuth_mesh, range_mesh = np.meshgrid(azimuth_array, range_array)

#雷達反射率因子
plt.figure()
plt.pcolormesh(azimuth_mesh, range_mesh/1000.0,
               reflectivity.transpose(), cmap='jet', vmin=-32, vmax=94)
plt.xlabel('Azimuth (degree)')
plt.ylabel('Range (km)')
plt.title('Reflectivity')
plt.colorbar()

#速度
plt.figure()
plt.pcolormesh(azimuth_mesh, range_mesh/1000.0,
               velocity.transpose(), cmap='seismic', vmin=-30, vmax=30)
plt.xlabel('Azimuth (degree)')
plt.ylabel('Range (km)')
plt.title('Velocity')
plt.colorbar()

#譜寬
plt.figure()
plt.pcolormesh(azimuth_mesh, range_mesh/1000.0,
               spectrum_width.transpose(), cmap='plasma', vmin=0, vmax=8)
plt.xlabel('Azimuth (degree)')
plt.ylabel('Range (km)')
plt.title('Spectrum Width')
plt.colorbar()

plt.show()