用 pyart 绘制简单的雷达图像

Python

This article was last updated on <span id="expire-date"></span> days ago, the information described in the article may be outdated.

前言

这是一篇关于使用 Pthon 库 Py-ART 读取雷达数据文件、绘制简单雷达图的教程。

这篇文章是我在去年暑假–也就是 2020年8月–做大创项目的时候写的,想着既然要写自己的博客了那就一块搬过来吧,当然原博客的地址还是要加的

由于中间修改过代码,并不确定代码 100% 正确

原地址:Py-ART 简易中文教程

绘制RHI图

下面是从nc文件中读取数据绘制RHI图的示例

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
#!/usr/bin/python3

#这里解释一下首行注释
#在Unix下,例如Ubuntu,你可以直接通过终端运行文件
#而该行注释则指定了Python解释器的位置
#在Windows下该行注释会自动忽略

import pyart as pa
import matplotlib.pyplot as plt
from sys import exit

if __name__=='__main__':
radar=pa.io.read_cfradial('plot_RHI.nc')

#检查扫描方式

if radar.scan_type!='rhi':
exit('Error: 请使用扫描方式为RHI的数据文件')

#创建一个从雷达数据中绘制图像的对象

display=pa.graph.RadarDisplay(radar)

xlabel='Distance from radar (km)'
ylabel='Height agl (km)'

#figsize指定画布大小,[宽,高],实际绘制出来之后可以用鼠标
#拖动窗口改变大小

fig=plt.figure(figsize=[7,7])
ax=fig.add_subplot(221)

#第一个参数是所要绘制的数据在radar.fields中对应的键名
#colorbar_label不接收''的话会显示默认标题
#set_limits设置横轴和纵轴的范围,ax接受画图的对象

display.plot_rhi('reflectivity',vmin=-20,vmax=80,
title='reflectivity',
axislabels=(xlabel,ylabel),colorbar_label='')
display.set_limits(ylim=[0,16],xlim=[90,140],ax=ax)

ax=fig.add_subplot(222)
display.plot_rhi('differential_reflectivity',vmin=-2,vmax=8,
title='differential_reflectivity',
axislabels=(xlabel,ylabel),colorbar_label='')
display.set_limits(ylim=[0,16],xlim=[90,140],ax=ax)

ax=fig.add_subplot(223)
display.plot_rhi('KDP',vmin=-2,vmax=10,
title='KDP',
axislabels=(xlabel,ylabel),colorbar_label='')
display.set_limits(ylim=[0,16],xlim=[90,140],ax=ax)

ax=fig.add_subplot(224)
display.plot_rhi('cross_correlation_ratio',vmin=0,vmax=1,
title='cross_correlation_ratio',
axislabels=(xlabel,ylabel),colorbar_label='')
display.set_limits(ylim=[0,16],xlim=[90,140],ax=ax)

plt.tight_layout()
plt.show()

绘制出的图像如下图所示

绘制RHI图

绘制PPI图

下面是从nc文件中读取数据绘制PPI图的示例

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
#!/usr/bin/python3

import matplotlib.pyplot as plt
import pyart as pa
from sys import exit

if __name__=='__main__':
radar=pa.io.read_cfradial('plot_PPI.nc')

#检查扫描方式

if radar.scan_type!='ppi':
exit('Error: 请使用扫描方式为PPI的数据文件')

display=pa.graph.RadarDisplay(radar)
fig=plt.figure(figsize=(7,7))

ax=fig.add_subplot(221)
display.plot('reflectivity',0,ax=ax,title='reflectivity',colorbar_label='',
vmin=-10,vmax=60,
axislabels=('','North South distance from radar (km)'))
display.set_limits((-250,250),(-250,250),ax=ax)

ax=fig.add_subplot(222)
display.plot('differential_reflectivity',0,ax=ax,
title='Differential Reflectivity',vmin=-2,vmax=6,colorbar_label='',
axislabels=('',''))
display.set_limits((-250,250),(-250,250),ax=ax)

ax=fig.add_subplot(223)
display.plot('differential_phase',0,ax=ax,
title='Differential Phase',colorbar_label='')
display.set_limits((-250,250),(-250,250),ax=ax)

ax=fig.add_subplot(224)
display.plot('cross_correlation_ratio',0,ax=ax,colorbar_label='',
title='Correlation Coefficient',
axislabels=('East West distance from radar (km)',''))
display.set_limits((-250,250),(-250,250),ax=ax)

plt.show()

绘制出的图像如下图所示

绘制PPI图

提取两个方位角的横截面绘制

下面是在从nc文件所有PPI扫描中提取两个方位角的横截面进行绘制的示例

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
#!/usr/bin/python3

import matplotlib.pyplot as plt
import pyart as pa
from sys import exit

if __name__=='__main__':
radar=pa.io.read_cfradial('plot_cross_section.nc')

#检查扫描方式

if radar.scan_type!='ppi':
exit('Error: 请使用扫描方式为PPI的数据文件')

#从PPI数据中选取方位角数据,第二个参数为天顶角范围

xsect=pa.util.cross_section_ppi(radar,[45,90])

display=pa.graph.RadarDisplay(xsect)
fig=plt.figure(figsize=(7,7))

#说明一下set_limits的作用,如果不设置,纵轴会因为20km以上
#有过多空白而影响图片效果
#这里的第二个参数不再是绘制PPI图时代表仰角层的含义了,由于
#官方介绍并不清晰,这里猜测其代表方位角,起始方位角对应0

ax=fig.add_subplot(211)
display.plot('reflectivity',0,vmin=-10,vmax=60)
display.set_limits((0,250),(0,20),ax=ax)

ax=fig.add_subplot(212)
display.plot('reflectivity',1,vmin=-10,vmax=60)
display.set_limits((0,250),(0,20),ax=ax)

#这里的tight_layout是一个神奇的函数,他会自动调整子图参数
#使之填充整个图像区域。带来的好处之一就是可以让你避免绘制
#出的子图之间标签出现重叠现象

plt.tight_layout()
plt.show()

绘制出的图像如下图所示

绘制方位角横截面图

绘制两个RHI图

下面是在同一张图中绘制两个RHI的示例
前面已经展示过如何绘制RHI图,再次说一遍的原因在于下面会介绍几个新的函数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
#!/usr/bin/python3

import matplotlib.pyplot as plt
import pyart as pa
import netCDF4 as net
from sys import exit

if __name__=='__main__':
radar=pa.io.read_cfradial('CfRadial_nuistCDP_20150712_074636.nc')

#检查扫描方式

if radar.scan_type!='rhi':
exit('Error: 请使用扫描方式为RHI的数据文件')

display=pa.graph.RadarDisplay(radar)

fig=plt.figure(figsize=[10,4])

ax=fig.add_subplot(1,2,1)
display.plot_rhi('reflectivity',0,vmin=-20,vmax=60,title_flag=0,
colorbar_label='')
display.set_limits(xlim=[90,140],ylim=[0,20],ax=ax)

ax=fig.add_subplot(1,2,2)
display.plot_rhi('velocity',0,vmin=-17,vmax=17,title_flag=0,
colorbar_label='')
display.set_limits(xlim=[90,140],ylim=[0,20],ax=ax)

#读取设备名称
instrument_name=radar.metadata['instrument_name']

#读取时间,radar.time['units']中包含两个信息,单位和起始时间
#第一个参数接收数值形式的时间值

time_start=net.num2date(radar.time['data'][0],
radar.time['units'])

#strftime接收格式字符串,将时间格式化表示

time_text=' '+time_start.strftime('%Y-%m-%d%H:%M:%SZ')

#获取角度

azimuth=radar.fixed_angle['data'][0]

title='RHI '+instrument_name+time_text+'Azimuth %.2f'%(azimuth)

plt.suptitle(title)
plt.tight_layout()
plt.show()

绘制出的图像如下图所示

绘制两个RHI子图

绘制包含等值线的图

下面是从nc文件中读取RHI数据并绘制包含等值线的示例

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
#!/usr/bin/python3

import matplotlib.pyplot as plt
import pyart as pa
import numpy as np
import scipy.ndimage as scn
from sys import exit

if __name__=='__main__':
radar=pa.io.read_cfradial('CfRadial_nuistCDP_20150712_074636.nc')

if radar.scan_type!='rhi':
exit('Error: 请使用扫描方式为RHI的数据文件')

#sweep表示绘制哪一个扫描层

sweep=0
display=pa.graph.RadarDisplay(radar)
fig=plt.figure(figsize=[20,5])
ax=fig.add_subplot(111)

#alpha用来控制图像透明度,数值越低图像越透明
#linewidth用来控制网格线的宽窄
#antialiased参数控制是否消除摩尔纹

display.plot('reflectivity',sweep=sweep,vmin=-10,vmax=60,
cmap='jet',fig=fig,ax=ax,
colorbar_label='Reflectivity (dB)',alpha=0.75,
linewidth=0.001,
antialiased=True)
display.set_limits(xlim=[90,140],ylim=[0,16])

#通过读取不同类型的数据,你还可以绘制不同的数据的等值线

data=radar.get_field(sweep,'differential_reflectivity')

#edges参数控制是否加上边界的数据

x,y,z=radar.get_gate_x_y_z(sweep,edges=False)

x/=1000
y/=1000
z/=1000

#此处的sigma参数的选取会影响等值线的分布

data=scn.gaussian_filter(data,sigma=3)

#sign函数用来保证结果符号的准确

R=np.sqrt(x**2+y**2)*np.sign(y)

#设置等值线的范围,应注意,levels的选取应根据data中数据的范围选取

plt.clabel(contours,levels,fmt='%r',inline=True,fontsize=7)
plt.show()

绘制出的图像如下图所示

包含等值线的RHI图

Author: Syize

Permalink: https://blog.syize.cn/2021/09/05/pyart-introduction/

本博客所有文章除特别声明外,均采用 CC BY-NC-SA 4.0 许可协议。转载请注明来自 Syizeのblog

Comments