python使用开源库进行反距离权重插值_反距离加权插值 python

1、反距离权重插值

反距离权重插值算法通用方程如下图；其中在估算中用到的控制点的数目是根据距离来取最近的s个点，在使用计算机编算法时关于最近邻点的获取使用KD-Tree来计算以提高运行效率。

2、KD-Tree

构建kd-tree

拿个例子说明为：

二维样例：{（2,3），（5,4），（9,6），（4,7），（8,1），（7,2）}

构建步骤：

1、初始化分割轴：

发现x轴的方差较大，所以，最开始的分割轴为x轴。

2、确定当前节点：

对{2，5，9，4，8，7}找中位数，发现{5,7}都可以，这里我们选择7，也就是(7,2);

3、划分双支数据：

在x轴维度上，比较和7的大小，进行划分：

左支：{(2,3)，(5,4)，(4,7)}

右支：{(9,6)，(8,1)}

4、更新分割轴：

一共就两个维度，所以，下一个维度是y轴。

5、确定子节点：

左节点：在左支中找到y轴的中位数(5,4)，左支数据更新为{(2,3)}，右支数据更新为{(4,7)}

右节点：在右支中找到y轴的中位数(9,6)，左支数据更新为{(8,1)}，右支数据为null。

6、更新分割轴：

下一个维度为x轴。

7、确定(5,4)的子节点：

左节点：由于只有一个数据，所以，左节点为(2,3)

右节点：由于只有一个数据，所以，右节点为(4,7)

8、确定(9,6)的子节点：

左节点：由于只有一个数据，所以，左节点为(8,1)

右节点：右节点为空。

最终，就可以构建整个的kd-tree了。

示意图如下所示
：

二维空间表示：

二维坐标系下的分割示意图
kd-tree表示：

kd-tree最近邻检索

在构建了完整的kd-tree之后，我们想要使用他来进行高维空间的检索。所以，这里讲解一下比较常用的最近邻检索，其中范围检索也是同样的道理。

最近邻搜索，其实和之前我们曾经学习过的KNN很像。不过，在激光点云章，如果使用常规的KNN算法的话，时间复杂度会空前高涨。因此，为了减少时间消耗，在工程上，一般使用kd-tree进行最近邻检索。

由于kd-tree已经按照维度进行划分了，所以，我们在进行比较的时候，也可以通过维度进行比较，来快速定位到与其最接近的点。由于可能会进行多个最近邻点的检索，所以，算法也可能会发生变化。因此，我将从最简单的一个最近邻开始说起。

一个最近邻
一个最近邻其实很简单，我们只需要定位到对应的分支上，找到最接近的点就可以了。

举个例子：查找(2.1,3.1)的最近邻。

计算当前节点(7,2)的距离，为6.23，并且暂定为(7,2)，根据当前分割轴的维度（2.1 < 7），选取左支。
计算当前节点(5,4)的距离，为3.03，由于3.03 < 6.23，暂定为(5,4)，根据当前分割轴维度（3.1 < 4），选取左支。
计算当前节点(2,3)的距离，为0.14，由于0.14 < 3.03，暂定为(2,3)，根据当前分割轴维度（2.1 > 2），选取右支，而右支为空，回溯上一个节点。
计算(2.1,3.1)与(5,4)的分割轴{y = 4}的距离，如果0.14小于距离值，说明就是最近值。如果大于距离值，说明，还有可能存在值与(2.1,3.1)最近，需要往右支检索。
由于0.14 < 0.9，我们找到了最近邻的值为(2,3)，最近距离为0.14。

3、反距离插值法的python实现

from __future__ import division
import numpy as np
from scipy.spatial import cKDTree as KDTree
import pandas as pd
import gdal
from osgeo import gdal, ogr, osr
import os

# 反距离权重插值
class Invdisttree:

    def __init__( self, X, z, leafsize=10, stat=0 ):
        assert len(X) == len(z), "len(X) %d != len(z) %d" % (len(X), len(z))
        self.tree = KDTree(X, leafsize=leafsize)  # build the tree
        self.z = z
        self.stat = stat
        self.wn = 0
        self.wsum = None

    def __call__( self, q, nnear=6, eps=0.0, p=1, weights=None ):
            # nnear nearest neighbours of each query point --
        q = np.asarray(q)
        qdim = q.ndim
        if qdim == 1:
            q = np.array([q])
        if self.wsum is None:
            self.wsum = np.zeros(nnear)

        self.distances, self.ix = self.tree.query(q, k=nnear, eps=eps )
        interpol = np.zeros( (len(self.distances),) + np.shape(self.z[0]) )

        jinterpol = 0
        for dist, ix in zip(self.distances, self.ix):
            if nnear == 1:
                wz = self.z[ix]
            elif dist[0] < 1e-10:
                wz = self.z[ix[0]]
            else:  # weight z s by 1/dist --
                w = 1 / dist**p
                if weights is not None:
                    w *= weights[ix]  # >= 0
                w /= np.sum(w)
                wz = np.dot( w, self.z[ix] )
                if self.stat:
                    self.wn += 1
                    self.wsum += w
            interpol[jinterpol] = wz
            jinterpol += 1
        return interpol if qdim > 1  else interpol[0]

#  ========================================================================================================================================================================================================================================================================

raindata = pd.read_csv("Rain.csv")  # 降雨量数据
maskPath = 'mask.shp'

res = 100 #你自己设置的图像分辨率

# csv中的字段名
longitude = raindata['coox'].tolist()
latitude = raindata['cooy'].tolist()
elevation = raindata['elevation'].tolist()


IDW_path_01 = "IDW_path_01.tif"  # 反距离权重插值得到的TIF原始文件
IDW_path_01_Crop = "IDW_path_01_Crop.tif"   # 反距离权重插值得到的TIF原始文件经掩膜裁剪后的TIF


# ======================================= 输入参数 END ================================================ #


# ---------------------------------------------------------------------------------------------------------------------- #
# --------------------------------------------   反距离权重插值    -------------------------------------------------------- #

longitude = list(map(float,longitude))
latitude = list(map(float,latitude))

elevation = np.array(elevation)
# 生成栅格的x、y
right_add = 300
ud_add = 100
# 由于雨量站的坐标不是在广州地图的边缘因此需要更大的面积  这里在栅格上下增加ud_add 在栅格右侧增加right_add （这里可以根据自身需求更改）
grid_lon = np.linspace(min(longitude)-10, max(longitude)+(10+right_add*100),int(abs(max(longitude) - min(longitude) )/res)+int(right_add*100/res))
grid_lat = np.linspace(min(latitude)-(10+ud_add*100), max(latitude)+(10+ud_add*100),int(abs(max(latitude) - min(latitude))/res)+int(2*ud_add*100/res))
# (min(longitude) - max(longitude))/100

lon_gridmesh, lat_gridmesh = np.meshgrid(grid_lon, grid_lat)


X = np.array([[x,y] for x, y in zip(longitude, latitude)])


invdisttree = Invdisttree(X, elevation, leafsize=10, stat=1)
Data_res = []  # 反距离权重插值后 得到的结果 list存储的每一行插值结果
for i in range(lat_gridmesh.shape[0]-1, 0, -1):
    X_new = np.array([[x, y] for x, y in zip(lon_gridmesh[i], lat_gridmesh[i])])
    interpol = invdisttree(X_new, nnear=12, eps=0.1, p=2)
    Data_res.append(interpol)


# ---------------------保存到tif中------------------------------#
np_data = np.array(Data_res)  # 将反距离权重插值的结果转为array
np_data[np_data == 0] = np.nan

driver = gdal.GetDriverByName("GTiff")

resp = res  # 分辨率 栅格大小
dataset = driver.Create(IDW_path_01, np_data.shape[1], np_data.shape[0], 1, gdal.GDT_Float64)
dataset.SetGeoTransform((min(grid_lon), resp, 0, max(grid_lat), 0, -resp))
#tif影像存放位置

dataset.GetRasterBand(1).WriteArray(np_data)
dataset.GetRasterBand(1).FlushCache()

del dataset

# ******************************************************************************************************************* #
# ***********************************************     加掩膜裁剪      *********************************************** #
# ******************************************************************************************************************* #

from osgeo import gdal

input_raster = IDW_path_01

output_raster = IDW_path_01_Crop


input_raster = gdal.Open(input_raster)

ds = gdal.Warp(output_raster, #输出栅格路径
               input_raster, # 输入栅格路径
               format='GTiff', # 影像保存格式
               cutlineDSName=maskPath, # 输入矢量路径
               cropToCutline=False, #(为True时，结果会与输入矢量大小一致。为False时，结果会与带裁剪的输入栅格大小一致)
               dstSRS='EPSG:4326',  # 获取投影信息, # 参考系：WGS84
               outputType=gdal.GDT_Float64, #数据类型
               dstNodata=0)
ds=None

# ******************************************************************************************************************* #
# ***********************************************     加掩膜裁剪  end *********************************************** #
# ******************************************************************************************************************* #

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
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142

参考自
《地理信息系统算法基础》科学出版社

关于KD-Tree的介绍