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 *********************************************** #
# ******************************************************************************************************************* #
参考自
《地理信息系统算法基础》科学出版社
189
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
