Syize の blog
home tags archives categories about links

用 pyart 绘制简单的雷达图像

2021-09-05 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！

updated at 2024-11-19

# Python# pyart
Previous Post
用 pyart 绘制雷达图像并进行反射率修正

Comments