IT / 气象/Meteorology · 2022年3月19日 3

如何从中央气象台全国雷达拼图中提取dbz

本文分享的是如何用Python从中央台的全国雷达拼图中提取出纯净的dbz图层(但并不能反解出地理坐标)。这是我很早以前的一段时间出于兴趣研究的一种方法,实现过程中会较多地使用numpy的向量运算技巧,比较适合有比较扎实的numpy基础的朋友进行综合练习,当然,本文也是我自己的一个代码备份。

依赖包

要运行本文的Python程序,需要预先安装以下依赖,使用pip即可安装。

numpy
matplotlib
scipy
webcolors

中央台全国雷达图的特点

中央气象台网站上每6分钟会发布一张全国雷达拼图,该图为PNG格式,也就是说这是一张无损压缩图,这也就为我们对颜色进行精确提取提供了可能性。该图是经过投影转换的,初步分析大概率为兰伯特正形投影,目前我并不知道图片的详细投影参数,因此暂时无法进行反解。全国雷达拼图相对比较干净,底图仅有行政边界以及河流、地名等,这种干净的地图可以大大降低dbz提取的难度。

下图为北京时2021年8月3日05点48分的全国雷达拼图,本文将以此图作为原始输入,剔除其中的底图,从而提取出纯净的dbz。

其最终效果如下:

Here we go~

一些先验信息

在正式开始coding之前,我们需要对图片的一些信息有所了解,首先就是dbz的颜色值,也就是原图右下角图例上显示的13个色块所表示颜色的编码,我们必须要知道这个编码,才能将dbz与底图分离。获取该颜色编码的方式有很多,比如你把图片用Photoshop打开,然后使用取色器直接从色块中取色即可获得它的hex码,或者安装一些小工具(比如MacOS上有一个小软件叫做ColorSlurp)使用取色器取色。由于我已经将色码取出,在这里就可以直接使用:

COLORS = [
    '#AD90F0',
    '#9600B4',
    '#FF00F0',
    '#C00001',
    '#D60100',
    '#FF0200',
    '#FF9000',
    '#E7C000',
    '#FEFF00',
    '#009000',
    '#00D800',
    '#00ECEC',
    '#01A0F6'
]

另外,由于使用RGB进行颜色提取时,图例上色块的颜色会与dbz的颜色混在一起,无法直接区分,因此我们需要使用位置信息将图例块直接剪掉,这就需要知道图例色块的区域范围坐标,获取该坐标的方法也是可以使用Photoshop,直接查看其边界的像素点坐标即可,它们的值如下:

left_x = 955
upper_y = 665
right_x = 973
lower_y = 858

导入所需要的包

import copy
import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import KDTree
import webcolors

将底图与dbz颜色分离

首先我们需要使用matplotlib加载图片文件,由于0-255整数形式的RGB的数值更精确,也较为通用,而plt.imread方法读取的颜色为0-1的浮点数,因此我们需要将其转为int类型。

raw_img_array = plt.imread('./SEVP_AOC_RDCP_SLDAS_EBREF_ACHN_L88_PI_20210711174200001.PNG')
rgb_img_array = (raw_img_array * 255).astype(int)

同时,我们要把颜色的hex码转为RGB整数值。

colors_arrays = np.array([list(webcolors.hex_to_rgb(c)) for c in COLORS])

接下来,我们将会对dbz和底图做分离,我们需要分别存储底图和dbz在两个不同的数组中。而分离过程是通过对原始图片进行过滤和重新赋值的方式完成的,因此初始化阶段我们可以将原始图片数组分别存入这两个数组。存dbz的数组我们取名为data_img_array,存底图的数组我们取名为flaw_img_array,由于后期底图的数组主要用于对缝隙的填补,因此我们使用flaw这个单词。

flaw_img_array = copy.deepcopy(rgb_img_array)
data_img_array = copy.deepcopy(rgb_img_array)

现在我们来看一下要怎么提取出dbz,总体来说我们有两种思路提取出dbz,一类是正向思维,另一类是反向思维,所谓正向思维,就是通过numpy的索引语法,直接筛选出符合dbz颜色的坐标点,再将非dbz的坐标点赋值为空白即 可,实现代码:

for n, colors_array in enumerate(colors_arrays):
    r,g,b = colors_array
    r_bool = np.abs(rgb_img_array[...,0]-r)<=5
    g_bool = np.abs(rgb_img_array[...,1]-g)<=5
    b_bool = np.abs(rgb_img_array[...,2]-b)<=5
    rgb_index = r_bool & g_bool & b_bool
    if n == 0:
        all_dbz_bool = rgb_index
    else:
        all_dbz_bool = all_dbz_bool | rgb_index

data_img_array[~all_dbz_bool] = np.array([255,255,255])

上述代码主要使用numpy布尔矩阵的集合运算(交集、并集运算),对符合条件的RGB通道值进行筛选,获取全图dbz的坐标位置点,并将非dbz坐标点的位置(补集运算)赋值为白色。这里在计算R、G、B的布尔数组的时候,之所以使用减法而非恒等式进行逻辑判断,是因为恒等的精确匹配会出现色彩遗漏(原因可能是读取图片的时候是0-1的小数,转为0-255整数时精度产生偏移),因此放宽对R、G、B的筛选范围,由于做交集时范围的选取产生了新的约束,且原图的色彩丰富度较低,并没有产生负面影响,这种方案可以使用。

得到的结果如下:

上述方法为正向思维,另一种方法是反向思维,即先将图中的dbz剔除掉,得到纯净的底图,然后获取底图带颜色区域的坐标,利用该坐标将底图赋值为白色,从而获取dbz,实现方法如下:

# 剔除dbz像素点,提取出底图部分的坐标
for colors_array in colors_arrays:
    dist = np.sum((rgb_img_array - colors_array) ** 2, axis=2)
    dbz_index = np.where(dist==dist.min())
    flaw_img_array[dbz_index] = np.array([255,255,255])

flaw_index = np.where(flaw_img_array.sum(axis=2)<255*3)

# 将底图部分赋值为空白
data_img_array[flaw_index] = np.array([255,255,255])

上述的实现方法的技巧在于,利用numpy的向量求和运算,得到数组与RGB颜色之间的欧氏距离,取距离最近位置坐标赋值为空白,过滤出底图。计算底图中RGB三个通道之和小于255*3(即非白色)的位置坐标,在原始图片中将该坐标点赋值为白色。

虽然反向思维听起来挺绕的,但实现起来要比正向思维简单,得到的结果和正向方法一样,而它们的右下角都有一个图例,这时候我们就需要使用先验信息里的图例位置将图例剪掉:

data_img_array[upper_y:lower_y+1,left_x:right_x+1] = np.array([255,255,255])

新得到的结果就没有图例了。

这个时候,我们仔细观察会发现,图中会出现由底图留下的空白缝隙(裂痕, flaw),下图中江苏的省界和“南京”二字清晰可辨,这就需要我们用某种方法把它填补上。

填补dbz裂痕

想要填补裂痕,首先要先知道裂痕的坐标。接着上面的思路,如果前面是按照正向方法实现的话,可以通过下面的方法得到裂痕的坐标:

flaw_yx = np.where((~all_dbz_bool) & (rgb_img_array.sum(axis=2)<255*3))

它的原理是先取非dbz坐标点,再与原图中的非空白点取交集,得到的就是缝隙的坐标点。

如果前面是按照反向方法实现的话,可以通过下面的方法得到裂痕的坐标:

flaw_yx = np.where(flaw_img_array.sum(axis=2)<255*3)

为了后面计算的方便起见,我们先来创建一个mask矩阵:

mask = np.full(raw_img_array.shape[:2], False)
mask[flaw_yx] = True

然后利用scipy.spatial.KDTree的方法用最近点的值填补裂缝区域:

flaw_xy = flaw_yx[::-1]
data_xy = np.where(~mask)[::-1]

data_points = np.array(data_xy).T
flaw_points = np.array(flaw_xy).T

data_img_array[mask] = data_img_array[~mask][KDTree(data_points).query(flaw_points)[1]]

由于KDTree函数所需要传入的坐标顺序为x,y而非y,x,因此需要提前做一个顺序的转换。

最后得到的图形:

在江苏地区的效果:

写在最后

这张图虽然可以比较好地提取出dbz,但是由于目前我并不能确切地知道该图片的投影参数,所以没有办法对数据的地理坐标进行反解,曾经对投影参数做过一些探索,但并不顺利,如果有大神知道该地图的精确的投影坐标,可以留言或者私信我,后面我就可以再写一篇反解雷达图的笔记文章了,甚至我也可以直接写一个Python包来简化反解的过程,我的一个小目标是做到将雷达图一键反解为NetCDF文件。