首先附上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')
数据的没有!