一个是idw插值（idw源码），另一个是补缺插值~其实还是一个东西

首先附上idw插值~：

import numpy as np

from scipy.spatial import cKDTree

class tree(object):

def __init__(self, X=None, z=None, leafsize=10):

if not X is None:

self.tree = cKDTree(X, leafsize=leafsize )

if not z is None:

self.z = np.array(z)

def fit(self, X=None, z=None, leafsize=10):

return self.__init__(X, z, leafsize)

def __call__(self, X, k=6, eps=1e-6, p=2, regularize_by=1e-9):

self.distances, self.idx = self.tree.query(X, k, eps=eps, p=p)

self.distances += regularize_by

weights = self.z[self.idx.ravel()].reshape(self.idx.shape)

mw = np.sum(weights/self.distances, axis=1) / np.sum(1./self.distances, axis=1)

return mw

def transform(self, X, k=6, p=2, eps=1e-6, regularize_by=1e-9):

return self.__call__(X, k, eps, p, regularize_by)

然后

正确的补缺idw插值代码：

pmfile= np.genfromtxt(r'D:\Thesis\xiamen\wrwpmn3kripm.csv', delimiter=',',skip_header=True,encoding='gbk')

testpoint2 = np.genfromtxt(r'D:\Thesis\point\2.csv',delimiter=',',skip_header=True)

ds = gdal.Open(r'D:\Thesis\samrast1.tif')

bandg = ds.GetRasterBand(1)

elevationg = bandg.ReadAsArray()

[cols, rows] = elevationg.shape

format = "GTiff"#5

driver = gdal.GetDriverByName(format)#6

pmkrigrid=[]

pmdata=np.array(pmfile)

ii,jj,pm=pmdata[:,3],pmdata[:,5],pmdata[:,9]

pmkridata=np.column_stack((ii,jj,pm)).astype('float')

for i in range(168):

listpm=[]

listpm=pmkridata[i*31:(i+1)*31]

listpm1=[]

for j in range(len(listpm)):

if listpm[j,2]!=999999:

listpm1.append(listpm[j])

listpm1=np.array(listpm1).reshape(len(listpm1),3).astype('float')

idw_tree = test_idw.tree(listpm1[:,0:2], listpm1[:,2])

pre = idw_tree(testpoint2)

# kriging = test_gaussian.SimpleKriging(training_data=listpm1)

# predict = kriging.predict(test_data=testpoint, l=.5, sigma=.2)

#pre=list(predict)

#pre.reverse()#这里面有个顺序问题，我一直没搞懂，到底要不要逆序

pmkrigrid.append(pre)

pe=(np.transpose(np.transpose(pre)[::-1])[::-1])###转置转置

pr=pe.reshape(58,52)[::-1]##转置之后可以把数据恢复到正确的位置

outDataRaster = driver.Create(r'D:\Thesis\xiamen\pm\pmidw1\pm'+str(i)+'.tif', rows, cols, 1, gdal.GDT_Int16)

outDataRaster.SetGeoTransform(ds.GetGeoTransform())

outDataRaster.SetProjection(ds.GetProjection())

outDataRaster.GetRasterBand(1).WriteArray(pr)

outDataRaster.FlushCache()

del outDataRaster

pmkrigrid=np.array(pmkrigrid).reshape(168*3016,)

datapreall=pd.read_csv(r'D:/Thesis/xiamen/pm/datapre/data168pre.csv',encoding='gbk')

datapreall['pmidw']=pd.DataFrame({'pmidw':pmkrigrid})

datapreall.to_csv(r'D:/Thesis/xiamen/pm/datapre/data168preidw.csv',mode='a',index=False,encoding='gbk')

import matplotlib.pyplot as plt

maiac=np.array(datapreall)[:,23]

maiachoy0=maiac[:3016][::-1]

mreshape=maiachoy0.reshape(58,52)

plt.figure(figsize=(5,5)) ##此处是显示

plt.imshow(mreshape,cmap='hsv')

数据的没有！