python反距离权重法_先从IDW（反距离权重）插值开始吧

IDW方法是一个很不错，很方便，很快。。。（自行百度）的方法。至少我必须要用到。。。

首先附上idw插值~：

我写的不是很好，如果你们想要更优化的可以再找找其他的版本。

然后

正确的补缺idw插值代码：

pmfile= np.genfromtxt(r'D:Thesisxiamenwrwpmn3kripm.csv', delimiter=',',skip_header=True,encoding='gbk')

testpoint2 = np.genfromtxt(r'D:Thesispoint2.csv',delimiter=',',skip_header=True)

ds = gdal.Open(r'D:Thesissamrast1.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:Thesisxiamenpmpmidw1pm'+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')

文章写的不好，最后打个广告公众号（一个有趣的灵魂W）