python反距离权重法_先从IDW(反距离权重)插值开始吧
IDW方法是一个很不错,很方便,很快。。。(自行百度)的方法。至少我必须要用到。。。首先附上idw插值~:我写的不是很好,如果你们想要更优化的可以再找找其他的版本。然后正确的补缺idw插值代码:pmfile= np.genfromtxt(r'D:Thesisxiamenwrwpmn3kripm.csv', delimiter=',',skip_header=True,encoding='gbk')

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)
更多推荐
所有评论(0)