Python3可视化雷达PUP数据产品(CINRAD-PUP)
- Coder: Lai Sheng @ College of Atmospheric Science, Chengdu University of Information Technology.
- E-mail: laish12@lzu.edu.cn
参考metpy.io.Level3File
雷达PUP数据接口的官方样例,编写了这个脚本。
PyCINRAD中的cinrad.io.PUP
坐标映射的BUG已经修复。
请看PUP雷达文件格式说明。
Metpy的官方样例中,将PUP产品文件雷达体扫的方位角az
和极径rng
读出后,把极坐标投影到了平面直角X-Y坐标系中:
# Convert az,range to x,y
xlocs = rng * np.sin(np.deg2rad(az[:, np.newaxis]))
ylocs = rng * np.cos(np.deg2rad(az[:, np.newaxis]))
xlocs
和ylocs
是某点到雷达的经向和纬向距离
这样根据雷达经纬度信息,还有经向和纬向距离,很容易就可以算图上每一点的距离经纬度:
sta_lon=f.lon #站点经纬度
sta_lat=f.lat
a=6.371e3
dy=2*a*np.pi/360.0
lon=np.zeros(xlocs.shape)
lat=np.zeros(xlocs.shape)
for i in range(xlocs.shape[0]):
for j in range(xlocs.shape[1]): #先算纬度,然后再确定经度
lat[i,j]=sta_lat+ylocs[i][j]/dy
dx=np.cos(np.deg2rad(lat[i][j]))*dy
lon[i,j]=sta_lon+xlocs[i][j]/dx
绘图使用的是basemap(1.2.0)
和matplotlib(3.0.3)
,也可以改用cartopy
画等距离圈是直接根据雷达文件读取出的极径rng
确定的
做法大致是用np.where
确定相应等距离圈的下标位置,然后把rng==距离
的所有点的坐标切片出来
最后用plot
将点连成线形成等距离圈
for cir in [50,100,150,200,230]:#画等距离圈
cir_lon=lon[:,np.where(rng==cir)].flatten()
cir_lat = lat[:, np.where(rng == cir)].flatten()
cir_lon,cir_lat=m(cir_lon,cir_lat)
m.plot(cir_lon,cir_lat,'w',linewidth=0.8)
画方位角线也是类似等距离圈的做法,但是画线只用了两个点
一个是雷达站点的经纬度(可以从文件中读取),另一个是径向上最远的一个点
for az_line in np.arange(30,360+30,30):
az_lon=[sta_lon,lon[np.where(np.ceil(az)==az_line),-1]]
az_lat=[sta_lat,lat[np.where(np.ceil(az)==az_line),-1]]
az_lon,az_lat=m(az_lon,az_lat)
m.plot(az_lon, az_lat, 'w', linewidth=0.8)
投影选择等距方位投影(Azimuthal Equidistant Projection)
此时1个纬距和1个经距是相等的,否则等距离圈不是正圆形
(下图是等经度投影cyl,原因自己体会思考吧,发现蛮多人都没发现到这个问题,直接等经纬度投影上边画圆)
我个人不做中尺度方向(被老板狠心抛弃QAQ),这个程序是写给泽儿(我GF可爱多)的。
目前仅做了基本反射率,基本速度还有就是将数据输出成网格的,她有需要的时候我再弄吧~