kri插值实例，批量气象插值（文件对应关系）

这几天懵了。。。

懵的不懂逻辑了，好吧废话不多说，这次解决的问题其实也比较基础，但却是非常常用和实用，对于入门简直神器。。。通常我们遇到的数据，不会整理的十分友好，需要我们对数据进行进一步处理，才能应用，特别是。。。如果数据之间排列跟预期的差别很大的时候。。。那就。。。虽然我说的也不是很清楚，但是明白的自然明白，就是这么佛系的自己：

贴代码咯

from __future__ import division

import numpy as np

from scipy import stats

from scipy import spatial

import test_gaussian

import matplotlib.pyplot as plt

import os

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

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

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

bandg = ds.GetRasterBand(1)

elevationg = bandg.ReadAsArray()

[cols, rows] = elevationg.shape

format = "GTiff"#5

driver = gdal.GetDriverByName(format)#6

time=train[:,13]

timeuni=np.unique(time)

for i in range(int(len(timeuni))):

train_hour=[]

x_hour=[]

y_hour=[]

for j in range(int(len(train))):

if timeuni[i] == train[j,13]:##对应序号或者对应属性，例如x=1，2 y=1，1，1，2，2，2 依次对应

# train_hour=list(train_hour)

# x_hour=list(x_hour)

# y_hour=list(y_hour)

i1=train[j,9]

i2=train[j,2]

i3=train[j,3]

train_hour.append(i1)

x_hour.append(i2)

y_hour.append(i3)

train_hour=np.array(train_hour)

x_hour=np.array(x_hour)

y_hour=np.array(y_hour)

conqx=np.column_stack((x_hour,y_hour,train_hour))

kriging = test_gaussian.SimpleKriging(training_data=conqx)

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

pre=list(predict)

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

pe=np.array(pre)

pr=pe.reshape(58,52)

outDataRaster = driver.Create("D:/ren0606/windraster/"+str(timeuni[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

以上是导出

plt.figure(figsize=(10,10))

im=plt.imshow(pr,cmap='YlGn',interpolation = 'sinc', vmin = 10, vmax = 0)

plt.colorbar()

plt.show

以上是显示

批量效果贼好