如何从中央气象台全国雷达拼图中提取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文件。