基于Python和GIS信息的南海气象雷达专题图自动制作

韦凯华， 梁伟杰,2， 韦海宁， 袁宇东,2

(1.广东省气象台，广东 广州 510080；2.中国南海石油联合服务总公司气象服务公司，广东 广州 510080；3.广西水利电力勘测设计研究院有限责任公司，广西 南宁 530023)

0 引言

南海海域是我国重要的战略纵深，地理位置十分重要。随着近年来海洋主权意识的不断增强，在南海进行生产作业和驻守的军民也在增加。受地理环境和气象环境条件影响，南海灾害性天气频发，海洋气象条件复杂，季风和热带云团导致的强对流天气以及路径曲折、强度多变的台风[1-3]，使得在南海航行及进行渔业作业的船舶常受恶劣天气的威胁[4]。天气雷达作为监测近海台风及强对流天气的有效工具，其相关产品对预报员分析判断未来天气动向，以及为公众和行业提供防御决策指导至关重要。但现行的雷达产品图在精细化地理信息(如南海岛礁)的绘制方面的不足限制了其使用价值的进一步挖掘。本文在此前工作的基础上[5-7]利用Python编程语言及其相关类库重新制作南海海域气象雷达可视化专题图，并实现自动化生成，以期进一步推进相关资料应用，充分发挥气象防灾减灾第一道防线的作用。

Python作为一种跨平台的计算机程序设计语言，结合了解释性、编译性、互动性等特点，并且是一种面向对象的高级语言。其设计之初定位自动化脚本(Shell)语言的编写，随着相关类库的不断更新和发展，逐渐被用于大型项目的开发。其提供的便捷的数据分析、拟合等功能，并被应用于金融分析、数据挖掘等多个领域[8]。本研究选NumPy、Pandas进行数据解码，利用强大的绘图类库Matplotlib完成专题图绘制。

1 数据预处理

1.1 气象雷达基数据资料的读取

南海海域以及周边的海南、广东、广西三省沿海共布设有S波段(SA、SC)气象雷达12部，不同型号的雷达其基数据格式不尽相同。中国气象局于2020年度完成了气象雷达数据格式标准化工作，实现标准格式单轨运行。本文选取位于南海中部三沙市永兴岛的西沙SC型雷达作为代表，一方面考虑西沙附近岛礁密集，对于精确地理信息与气象信息结合需求相对较大；另一方面，SC型雷达在上述业务变更前后，基数据格式的改动较大。标准格式雷达基数据的数据结构如图1所示。

图1 标准格式基数据整体结构

其中，N表示第N个仰角；M表示第M个径向；K表示第K个数据类型；I表示径向数据长度，数据读取过程中涉及到上述多个结构字段的循环。本文首先选用NumPy类库中的numpy.dtype函数构建类似结构体的字典类型。并采用NumPy类库的frombuffer函数循环遍历读取雷达基数据的各个数据块并经解码后存入对应numpy.ndarray类型的三维数组(仰角、方位角和径向)中。其关键代码如图2所示。

图2 标准格式基数据读取及解码关键代码

由于Numpy底层采用Fortran和C++ BOOST进行数据处理和运算的工作，本文采用numpy.frombuffer函数可以在确保运行效率的同时充分发挥Python易读性强的特点。数据读取和解码完成后，根据各个径向的方位角信息，以方位角升序的方式对变量三维数组进行排序，形成排序后的反射率因子(ref)和多普勒速度(vel)三维数组数据，用于后续进一步绘图或计算处理。

1.2 GIS信息的预制作

从国家地理信息共享平台(全国地理信息资源目录服务系统)获取1：100万全国地理信息基础数据，由于获取的是分幅数据，故需要通过ArcGIS软件按照所需要的南海区域地理信息范围合并相应的几幅，再进行裁剪，最后得到一组以SouthChinaSea命名的地理信息数据。对于主要岛礁名称的文字标注，需要读取其中的库表关系获取名称和位置信息。为简化步骤，提高读取效率，用ArcGIS软件中的要素转点功能，将.shp文件中原有的面元素和线元素转换为点元素，同时运用软件中的Excel转换工具，将已转换为点元素的.shp文件中包含各元素的属性表转换为.xls文件并另存，每行中就包含了元素所属类别、中英文名称、经纬度坐标等。后续即可通过Python程序读取该.xls文件实现地理信息的标注。

2 气象雷达可视化专题图自动制作程序设计

本研究利用Python语言设计程序实现气象雷达可视化专题图自动制作，程序框图及数据处理流程如图3所示。

图3 程序框图及数据处理流程

气象雷达站观测到的雷达回波信息数据上传至国省两级气象信息中心数据库后，再通过推送下发或数据接口访问的方式将基数据提供给气象台等业务单位使用[9]。设计的程序主要分为3大功能模块，即时间调度与实时资料检索匹配模块、雷达数据读取计算处理模块、雷达数据与GIS信息绘图模块。经程序处理后自动制作生成的气象雷达可视化专题图可推送到本地或网络指定位置，供气象部门内及部门外各渠道调用。

2.1 实时资料的检索与匹配

无论通过数据接口访问或是本地推送(或FTP共享)的方式获取气象雷达基数据，要实现实时的专题服务图自动制作，首先需完成最新资料检索下载。以本地(或FTP共享)方式访问为例，读取配置文件settings.ini中的资料路径fileFolder，列出该目录下后缀名为.bz2的文件。标准雷达基数据文件名的形式为Z_RADR_I_Ziiii_YYYYMMDDhhmmss_O_DOR_TT_CAP_FMT.bin.bz2，其中iiii为4位数字的雷达区站号；TT为雷达型号；YYYYMMDDhhmmss为基数据产生时间。本文利用正则表达式提取雷达站号和数据时间信息。通过提取后的雷达站号确定西沙雷达基数据存放的次级目录，进入到目录中后对列表中的所有文件根据文件名时间进行逆序排列，排序后的第一个文件名即为最新的资料。考虑到可能出现的雷达数据传输故障，设定为机器当前时间与最新雷达数据文件名时间的延迟小于15分钟以内时，才进行相应的数据处理工作。具体代码如图4所示。

图4 实时资料检索与匹配关键代码

2.2 气象雷达基数据资料的可视化专题图绘制

在读取雷达基数据资料的基础上进行可视化专题图的绘制。自定义一个带参数的draw_radar_ppi函数，传入的参数为要绘制的专题图数据类型。先利用循环按投影比例绘制雷达距离圈和方位角射线，此后根据传入的参数确定绘制数据类型及其对应的色标(色标预先以数组形式定义在mycolors_变量名.npy文件中，通过自定义的get_cmap函数传递)。以雷达所在坐标经纬度为中心，按方位和径向计算雷达极坐标投影至地理坐标的经纬度位置。本文充分利用NumPy的向量运算的特性，运用Pyproj包中的Geod函数，将极坐标通过WGS84坐标系投影转换为等经纬网上的笛卡尔坐标，得到雷达数据在平面上的lon(经度)与lat(纬度)值，再利用pcolormesh函数按排序处理后的三维数组数据值对应的色标绘制图形。pcolormesh函数可用于均匀/非均匀色标的绘制，通过对BoundaryNorm类型的指定实现，其构造函数为BoundaryNorm(boundaries, ncolors, clip)。各参数含义如下：levels为传入的色标阈值；ncolors为色标的色阶数；clip为控制边界值是否截断进行绘制，选用True表示超出阈值上下边界的均截断为阈值。

再确定子坐标轴，分别以文本形式绘制相关的雷达基本信息和观测信息，以矩形填色方式绘制色标与阈值，并从基数据中提取时间信息，用datetime类完成儒略历、地方时转换，用作保存专题图的文件名使用，完成整幅图的基本绘制。关键代码如图5所示。

图5 专题图绘制关键代码

2.3 GIS信息的绘制

运用Matplotlib类库中的Basemap库指定数据绘图的投影方式以及经纬度范围，并通过readshapefile函数依次读入国省市及县界的.shp地理信息，本设计中，使用的.shp地理信息基础数据包括了南海所有岛礁的信息以及南海九段线，符合地理信息规范化制图的规定。drawmapboundary函数使用指定颜色绘制底图，此处指定为海洋颜色，fillcontinents函数将在底图上用指定颜色填充绘制陆地地形。此外，南海岛礁名称的标注，通过Pandas类库中的read_excel函数，按指定格式读取第2.2节中预先制作好的SouthChinaSea.xls中按行列存储的地理信息，提取其中的名称、经纬度值，最后通过循环的方式逐条通过text函数将岛礁名称按前述绘图投影方式的经纬度以标注文字的形式绘制在底图上。

3 专题图效果展示与应用

原程序制作的专题图如图6(a)所示。设计的Python程序自动运行后制作的南海气象雷达专题图如图6(b)、图6(c)所示。其中图6(b)为绘制的雷达反射率因子(ref)资料。相较于原有VB(Visual Basic)程序绘制的专题图(图6(a))，可以发现本文设计程序制图的若干优点。其一是增加了精确的地理信息，三沙市范围内西沙群岛、中沙群岛的岛礁位置及名称均清晰标注在图内，同时在海南岛范围内将地理信息精确绘制至县级(原来仅为市级)，并增加标注沿海重要港口及雷达探测距离圈范围，帮助终端用户能清楚获悉恶劣天气的影响区域、范围和过去及未来一段时间的移动趋势，特别是在台风临近登陆时，上述功能可以辅助确定精确到县级的可能登陆地段，以便更有针对性地采取防御措施。其二是改进了雷达数据处理的逻辑，原有的VB专题图制作程序设计直接采取径向数据方位角取整的方法处理和绘制图像，造成了个别径向的缺失，如图6(a)中150°-180°方位角的部分径向缺失影响了雷达回波图像的完整性。本文的程序由于在绘制图像前先进行了排序，并考虑到可能缺失径向的插补，避免了绘制缺失径向问题的出现。其三是高精度图像绘制，由于雷达进行了技术升级与数据格式标准化统一工作，分辨率从原来的1.2 km提升了4倍至300 m，相应程序也进行了修改，使得绘制的图像更精细。同时本文的程序设计可以通过改变图像dpi的方式调整输出图像的分辨率(默认dpi为100)。其四是在原有基础上新增了径向速度(vel)资料的专题图，如图6(c)所示。

图6 各程序制作的专题图

所选个例为2020年第22号台风“环高”过程，通过制作的专题图，预报员可以分析计算得到专题图代表的雷达探测0.5°仰角上，该台风底层中心附近的最大风速为51 m/s(15级上限、16级下限)，考虑到台风中心距离雷达约150 km，根据测高公式算得回波离地高度约2.5 km，故结合预报经验公式折算至地面估计的风应在14到15级(46-48 m/s)左右。本台预报员最终制作发布的台风路径和强度在未来18 h内均维持在45-50 m/s，同时继续发布了台风紧急警报，指出“受冷空气和台风‘环高’的共同影响，南海中西部海面：旋转风10到11级阵风13级，其中台风中心经过的附近海面14级阵风15级”，预报结论与观测实况较为一致，应用专题图起到了较好的服务效果。

此外程序自动生成的专题图可以按配置文件settings.ini中的目标路径输出到指定目录下或FTP服务器，除可以生成网站及App端常用的.jpg或.png格式外，还可以生成矢量图.eps及.svg格式文件，为研究应用印刷出版提供支持。

4 总结

本文基于Python语言，结合精细的GIS信息实现了南海气象雷达专题服务图的自动制作，支持雷达基数据自动检索下载、解码读取、图像投影绘制等功能。程序流程简明清晰，从数据检索到成图分发耗时在5 s以内，满足业务时效需求。制作的专题服务图供气象部门内预报分析网站及对内/外App等多渠道应用，特别是为近海台风、强对流等天气的气象防灾减灾工作提供了有力的资料支撑。添加的精确地理信息对于天气影响区域的研判准确性较此前大为提升，且可以通过修改相关预制GIS文件的方式增加或突出显示，在后续也可以按需增加近海航线、渔场等信息[10]，为更具体的应用场景服务。此外，程序设计是基于标准格式的雷达基数据处理，可扩展应用到南海海域及周边任意一部业务雷达进行专题服务图的制作，具有一定的推广应用价值。