【万事屋】一个区域斜向经纬网格的补全(上)

以下全文代码均已发布至和鲸社区，复制下面链接，可一键fork跑通：https://www.heywhale.com/mw/project/64f1be79abb16766700babe6?shareby=620ce369c1ae5e001747085e#

这是一个粉丝在万事屋活动中的投稿提问，原文是这样的：

前段时间遇到过一个问题，需要将一个数据进行绘图，数据是一个3列（经度、纬度、变量）的csv文件，是将一个分辨率大约1km的曲线网格（也就是WRF常见的经纬度都有二维的数据类型）展平后得到的结果，原网格的投影信息未知。从绘图的角度来说，当然可以将这些数据当作散点来绘制，但我还是想尝试复原原来的网格，不过在进行了一些简单的尝试后失败了，不知道老师在这方面有没有什么想法呢？

最开始我以为只是一个简单的插值问题，后来经过沟通发现不是，如果用一句话来描述这道题，可能用《一个区域斜向经纬网格的补全》比较恰当。我来用我的理解描述一下这道题：

假如你有一个 csv 表格文件，里面存的是经度、纬度和 DEM 数据。然后这些数据如果用散点图的方式画到地图上大概类似于这样的分布：

可以看出来它大概是遵循一种等距离的网格分布的形式，但又不是我们常见的规则等经纬度的格点坐标，大概可以想象它是把规则等经纬网格进行了一定角度的旋转，但投影未知。现在我们需要在不做格点插值重采样的情况下，将其原坐标补全到一个规则形状的矩阵上，缺失部分可以用 nan 填充，类似这种结果：

这样最后输出的结果的好处是可以作为 NetCDF 文件或矩阵形式的 npy 存储，方便使用常规格点方式画图，同时也不会因插值而引入细微的数据差异。

下面我就来介绍一下我的求解过程，提前预警一下，内容有些干燥，可能需要读者具有一定 numpy 基础的才能比较流畅地阅读。另外整个过程比较繁琐和漫长，需要一定的耐心才能读完。鉴于求解过程的长度，我将分上下两期来呈现完整的过程，另外我所介绍的过程是我做出这道题的第一解，而不是最优解，请理性看待。

题目分析

我们先来大致分析一下这道题，要解决这道题，就需要找出各个点对于行和列的位置关系，把它们分门别类地组织起来。在我们肉眼观察来看，矩阵中带有倾斜特征的行和列一眼就能看出来，但是对于代码来说，它并不知道哪些点应该分到同一行或同一列，需要我们寻找一些方法来建立它们的关联。

一个最直观的方式就是从某个起始点以特定方向进行临近点查找，比如从最左上角的点，沿左下的方向找最近点，再以这个点作为下一次查找的起始点，以此类推，直到没有点满足要求为止。然后再从列首沿右下的方向查找下一列的起始点为止，再以此类推查找下一列的所有点。再以此类推查找所有列和所有行。

先不说这个方法能不能覆盖到所有点，在实际查找过程中可能会遇到的第一个麻烦是断点。如果遇到断点，在查找列的时候就会提前结束，从而无法查到断点之后的点。另外如果我们要补全一个规则网格坐标，也需要把这些断点补上。所以我们解题的上半部分的任务就是先把缺失的断点补上。

代码实现

以下是我们这一阶段所需要用到的包：

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d

读出原始表格浅看一下

df = pd.read_csv('./raw.csv')
df

lon lat RASTERVALU
0 108.001830 25.006561 373
1 108.011730 25.006302 352
2 108.021640 25.006020 360
3 108.031530 25.005741 570
4 108.041435 25.005468 368
... ... ... ...
44310 109.964485 26.991802 228
44311 109.974594 26.991344 195
44312 109.984710 26.990894 292
44313 109.994820 26.990437 355
44314 109.985214 26.999956 226
44315 rows × 3 columns

后面为了处理方便和性能方面的考虑，我们要抛弃 pandas，纯粹用 numpy 来进行矩阵运算。

# 将表格按经纬度增序排列
df = df.sort_values(['lon', 'lat'])

# 将经纬度分别读取为 ndarray 矩阵
lons = df['lon'].values
lats = df['lat'].values

画一张坐标全景分布图浅看一下

plt.figure(figsize=(15,15))
plt.scatter(lons, lats, s=1, c='grey')
plt.savefig('./first-scatter.png', bbox_inches='tight')

考验眼神的时候到了，你能看到缺失的三个点在哪吗？我们现在的任务就是把这三个缺失的点的具体坐标找出来，并在表格中用 nan 填回。但是由于坐标并不是以规则的网格点排列，所以无法直接用索引坐标获取，需要一些技巧。首先我们要对网格点的一些基本信息有所认知，至少我们需要知道每个网格点紧邻点与之的距离。

distance = sorted(((lons - lons[500]) ** 2 + (lats - lats[500]) ** 2) ** 0.5)
distance[:10]
[0.0,
 0.009054971010443344,
 0.009066301395831013,
 0.010073836607766624,
 0.010078563836182674,
 0.013505179598955783,
 0.01351177027631989,
 0.013588866766582372,
 0.013594153154937977,
 0.018111275521064726]

上面是表格中第500个点与所有表格点之间距离的升序前10位，第1个0是本尊，第2至第5个是上下左右紧邻的点与之的距离，第6到9个是斜方向上4个点的距离。为了防止后续在查找点的过程中出现意外情况，我们需要取一个介于紧邻点和斜向点之间的距离选一个数字作为过滤阈值，用于保证后续的搜索是按照预期的方向进行。

DIST_LIMIT = 0.011  # 紧邻点距离不超过该数，而斜向点的距离必超过该数

我们先来找到最左上角的坐标。

# 查找最左上角的点
distance = ((lons - lons.min()) ** 2 + (lats - lats.max()) ** 2) ** 0.5
first_idx = np.where(distance==distance.min())
first_idx
(array([184]),)

画图验证一下位置是否正确

# 验证左上角的点是否正确
plt.figure(figsize=(15,15))
plt.scatter(lons, lats, s=0.5, c='grey')
plt.scatter(lons[first_idx], lats[first_idx], s=5, c='r')
plt.savefig('./find-first-point.png', bbox_inches='tight')

位置正确，下面需要编写一个查找左下角限定范围内最近点索引的函数

def find_nearest_of_left_lower(lons, lats, lon, lat):
    filter_idx = np.where((lons < lon)&(lats < lat))
    lons1 = lons[filter_idx]
    lats1 = lats[filter_idx]
    
    distance = ((lons1 - lon) ** 2 + (lats1 - lat) ** 2) ** 0.5
    sorted_dist = np.sort(distance)
    sorted_dist = sorted(sorted_dist[(sorted_dist<=DIST_LIMIT)&(sorted_dist>0)])
    idx = None
    if len(sorted_dist) > 0:
        idx_ = np.where(distance==sorted_dist[0])
        lon_ = lons1[idx_]
        lat_ = lats1[idx_]
        idx = np.where((lons==lon_)&(lats==lat_))

        return idx

上述函数向左下方向查找距离自己最近且在 DIST_LIMIT 范围内的点，并返回其索引值。

接下来，我们在选定的左上角的第一个点上，沿左下方向循环迭代查找沿途点的索引。

idx = first_idx
first_col_lons = []
first_col_lats = []
while idx is not None:
    first_col_lons.append(lons[idx])
    first_col_lats.append(lats[idx])
    idx = find_nearest_of_left_lower(lons, lats, lons[idx], lats[idx])

(array([177]),)
(array([170]),)
(array([164]),)
(array([156]),)
(array([149]),)
(array([142]),)
(array([135]),)
(array([128]),)
(array([121]),)
(array([114]),)
(array([107]),)
(array([100]),)
(array([93]),)
(array([86]),)
(array([79]),)
(array([72]),)
(array([65]),)
(array([58]),)
(array([51]),)
(array([44]),)
(array([37]),)
(array([30]),)
(array([24]),)
(array([18]),)
(array([12]),)
(array([6]),)
(array([0]),)
None

可以看到输出了一长串的索引坐标，我们画个图来验证一下结果是否符合预期

# 验证是否成功识别出最左侧一列的点
plt.figure(figsize=(15,15))
plt.scatter(lons, lats, s=0.5, c='grey')
plt.scatter(first_col_lons, first_col_lats, s=5, c='r')
plt.savefig('./find-first-col.png', bbox_inches='tight')

可以看到它已经识别出了最左侧第一列的点，我们还需要查找到下一列的第一个点，于是我们需要编写一个查找右下方向最近点的函数，参考上面的函数改成过滤右下侧格点即可。

def find_nearest_of_right_lower(lons, lats, lon, lat):
    filter_idx = np.where((lons > lon)&(lats < lat))
    lons1 = lons[filter_idx]
    lats1 = lats[filter_idx]

    distance = ((lons1 - lon) ** 2 + (lats1 - lat) ** 2) ** 0.5
    sorted_dist = sorted(distance[(distance<=DIST_LIMIT)&(distance>0)])
    idx = None
    if len(sorted_dist) > 0:
        idx_ = np.where(distance==sorted_dist[0])
        lon_ = lons1[idx_]
        lat_ = lats1[idx_]
        idx = np.where((lons==lon_)&(lats==lat_))

        return idx

我们可以从第一列第一行的点开始，查找以它为起点的右下侧最近点：

next_idx = find_nearest_of_right_lower(lons, lats, lons[first_idx], lats[first_idx])
next_idx
(array([413]),)

再画张图来验证一下：

next_col_first_lons = lons[next_idx]
next_col_first_lats = lats[next_idx]
# 验证是否成功识别出左侧第二列的第一个点（蓝色）
plt.figure(figsize=(15,15))
plt.scatter(lons, lats, s=0.5, c='grey')
plt.scatter(first_col_lons, first_col_lats, s=5, c='r')
plt.scatter(next_col_first_lons, next_col_first_lats, s=10, c='b')
plt.savefig('./find-second-col.png', bbox_inches='tight')

可以看到，已经正确地定位了第二列的第一个点的位置（蓝色）。然后我们需要再编写一个函数，把上面的功能再组装一下，自动查找每一列的所有点并给出下一列初始点的坐标。

def seek_col_elements(idx, lons, lats):
    first_idx = idx
    col_lons = []
    col_lats = []
    while idx is not None:
        col_lons.append(lons[idx][0])
        col_lats.append(lats[idx][0])
        idx = find_nearest_of_left_lower(lons, lats, lons[idx], lats[idx])

    next_idx = find_nearest_of_right_lower(lons, lats, lons[first_idx], lats[first_idx])
    return col_lons, col_lats, next_idx

剩下的就可以交给循环了，循环查找所有列的数据，为了后续得到缺失点的信息，我们需要在循环过程中获取列的序号以及每列识别出来的元素个数。

col_lons, col_lats, next_idx = seek_col_elements(first_idx, lons, lats)
columns = [{
    'col_n': 0,
    'lons': col_lons, 
    'lats': col_lats,
    'length': len(col_lats)
}]
already_set = set([(float(c[0]), float(c[1])) for c in list(zip(col_lons, col_lats))])
col_n = 1
while next_idx is not None:
    col_lons, col_lats, next_idx_= seek_col_elements(next_idx, lons, lats)
    already_set = set([(float(c[0]), float(c[1])) for c in list(zip(col_lons, col_lats))])
    columns.append({
        'col_n': col_n,
        'lons': col_lons, 
        'lats': col_lats,
        'length': len(col_lats)
    })
    next_idx = next_idx_
    col_n += 1

把所有能查找出来的行列重新组织存入 columns 列表以后，我们可以再画出一张图来验证一下结果是否符合预期

# 验证是否能把每一列都正确地识别出来，做到顺序对，不串色
plt.figure(figsize=(15,15))
plt.scatter(lons, lats, s=0.5, c='grey', alpha=0.5)
colors = ['r', 'g', 'b']
for n, column in enumerate(columns):
    color = colors[n % len(colors)]
    col_lons = column['lons']
    col_lats = column['lats']
    plt.scatter(col_lons, col_lats, s=1, c=color)
plt.savefig('./rgb-grid.png', bbox_inches='tight')

我们按照红、绿、蓝三个颜色轮流涂色的方式来区分识别出来的每一列互相之间的关系是否有错位（原底色用灰色占位）。可以看到识别的结果中，每一列的颜色都是单一颜色，没有串色的情况发生，说明每一列的识别都是正确的，但是右上和右下部分是识别不出来的，不过没关系，这个后面我们有办法解决。现在的任务是找到缺失的三个点的位置，从图中可以比较明显地看出来有3列不正常，“彩带”的尽头，就是断点，我们要利用这个现象把断点揪出来。

可以通过比长短的方式识别出断裂列，即找出那些长度比自己前后列长度都要短的列。

for n, col in enumerate(columns):
    if n == 0 or n == len(columns)-1:
        continue
    col_n = col['col_n']
    previous_n = columns[n-1]['length']
    current_n = columns[n]['length']
    next_n = columns[n+1]['length']
    if current_n < previous_n and current_n < next_n:
        print(n,col_n)
        
82 82
90 90
120 120

我们现在找出了缺失点在第82、第90和第120列。现在我们要拿到这3列的坐标点，用1维插值算法插出缺失点的坐标，先写个函数来实现这个外插1个点的功能：

def get_next_one_coord(lons, lats):
    lats = np.array(lats)
    deltay_mean = (lats[1:] - lats[:-1]).mean()  # 计算纬度平均间隔

    interp = interp1d(lats, lons, fill_value='extrapolate')
    
    y1 = deltay_mean+lats[-1]
    x1 = interp(y1)

    return float(x1), y1
    
missing_lons = []
missing_lats = []

for i in [82, 90, 120]:
    missing_lon, missing_lat = get_next_one_coord(columns[i]['lons'], columns[i]['lats'])
    missing_lons.append(missing_lon)
    missing_lats.append(missing_lat)

missing_df = pd.DataFrame({'lon':missing_lons, 'lat':missing_lats, 'RASTERVALU':np.nan})
missing_df


lon lat RASTERVALU
0 108.826256 26.691993 NaN
1 108.887980 26.254511 NaN
2 109.151330 25.429875 NaN

通过插值计算，我们推算出了3个缺失点的经纬度坐标，并将其生成了一个和原表格同构的 DataFrame 对象，但是值列表是用 nan 替代的，后续将用于二者的拼接，现在我们先来画图验证一下这三个缺失点的位置：

# 验证左上角的点是否正确
plt.figure(figsize=(15,15))
plt.scatter(lons, lats, s=0.5, c='grey')
plt.scatter(missing_lons, missing_lats, s=5, c='r')
plt.savefig('./show-missing-point.png', bbox_inches='tight')

你一开始的时候找对了吗？

好了，后面我们需要做的，就是把这个增补的 DataFrame 和原始 DataFrame 进行一个合并、保存，供后续使用：

df = pd.concat([missing_df, df])
df.to_csv('./complete.csv', index=False)

最后，我们把补全后的坐标点再画出一张图来，看看效果。

可以看到这是一个把缺失点填充上以后的散点图，中间已经没有任何空缺了，我们的第一阶段任务也算完成了。

下一期我将介绍如何识别出完整原始矩阵，以及对其进行扩展、去重、最后保存成规则矩阵的过程。