python：克里金插值

python：克里金插值
最近写代码遇到了使用样本数据做克里金插值的事情。于是将Excel保存的【x坐标，y坐标，样本值】数据结合tif数据做了克里金插值，并将代码记录下来。

克里金插值结果：

输入数据格式：

import os, sys
import pandas as pd
from sklearn.gaussian_process import GaussianProcessRegressor
import numpy as np
from numpy import inf
from osgeo import gdal
from tqdm import tqdm
import pykrige.kriging_tools as kt
from pykrige.ok import OrdinaryKriging
from numpy import genfromtxt# 读写tif的类
class GRID:# 读图像文件def read_img(self, filename):dataset = gdal.Open(filename)  # 打开文件im_width = dataset.RasterXSize  # 栅格矩阵的列数im_height = dataset.RasterYSize  # 栅格矩阵的行数im_geotrans = dataset.GetGeoTransform()  # 仿射矩阵im_proj = dataset.GetProjection()  # 地图投影信息im_data = dataset.ReadAsArray(0, 0, im_width, im_height)  # 将数据写成数组，对应栅格矩阵im_data = np.array(im_data)del dataset  # 关闭对象，文件datasetreturn im_proj, im_geotrans, im_data, im_width, im_height# 写文件，写成tiffdef write_img(self, filename, im_proj, im_geotrans, im_data):# 判断栅格数据的数据类型if 'int8' in im_data.dtype.name:datatype = gdal.GDT_Byteelif 'int16' in im_data.dtype.name:datatype = gdal.GDT_UInt16else:datatype = gdal.GDT_Float32# 判读数组维数if len(im_data.shape) == 3:im_bands, im_height, im_width = im_data.shapeelse:im_bands, (im_height, im_width) = 1, im_data.shape# 创建文件driver = gdal.GetDriverByName("GTiff")  # 数据类型必须有，因为要计算需要多大内存空间dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)dataset.SetGeoTransform(im_geotrans)  # 写入仿射变换参数dataset.SetProjection(im_proj)  # 写入投影if im_bands == 1:dataset.GetRasterBand(1).WriteArray(im_data)  # 写入数组数据else:for i in range(im_bands):dataset.GetRasterBand(i + 1).WriteArray(im_data[i])del dataset#读取每个tiff图像的属性信息，和上面函数相似，懒得改了，混着用。
def Readxy(RasterFile):ds = gdal.Open(RasterFile,gdal.GA_ReadOnly)if ds is None:print ('Cannot open ',RasterFile)sys.exit(1)cols = ds.RasterXSizerows = ds.RasterYSizeband = ds.GetRasterBand(1)# data = band.ReadAsArray(0,0,cols,rows)noDataValue = band.GetNoDataValue()projection=ds.GetProjection()geotransform = ds.GetGeoTransform()return rows,cols,geotransform,projection,noDataValue# 克里金插值
def KrigingInterpolate(InputExcelPathAndName, X_coordinateField, Y_coordinateField, Sample_valueField,tifDataPath, tifName, resultDataSetPath, outputKrigingTifName):Data = pd.read_excel(InputExcelPathAndName)Points = Data.loc[:, [X_coordinateField, Y_coordinateField]].valuesValues = Data.loc[:, [Sample_valueField]].valuesPoints1 = np.array(Points)Values1 = np.array(Values)# 数据去空值Points1 = np.nan_to_num(Points1)Points1[Points1 == inf] = np.finfo(np.float16).maxValues1 = np.nan_to_num(Values1)Values1[Values1 == inf] = np.finfo(np.float16).max# 读取遥感影像数据run = GRID()# 这一个没有参与运算，主要为了读取它的行列数、投影信息、坐标系和noData值rows, cols, geotransform, projection, noDataValue = Readxy(tifDataPath + tifName)print(rows, cols, geotransform, projection, noDataValue)nXsize = colsnYsize = rows# **********************************//dataset = gdal.Open(tifDataPath + tifName)  # 打开tifadfGeoTransform = dataset.GetGeoTransform()  # 读取地理信息# 左上角地理坐标# print(adfGeoTransform[0])# print(adfGeoTransform[3])# 右下角地理坐标px = adfGeoTransform[0] + nXsize * adfGeoTransform[1] + nYsize * adfGeoTransform[2]py = adfGeoTransform[3] + nXsize * adfGeoTransform[4] + nYsize * adfGeoTransform[5]OK = OrdinaryKriging(Points1[:, 0],Points1[:, 1],Values1[:, 0],variogram_model="linear",verbose=False,enable_plotting=False,)gridx = np.arange(adfGeoTransform[0], px, adfGeoTransform[1])gridy = np.arange(adfGeoTransform[3], py, adfGeoTransform[5])print(type(gridy), gridx)z, ss = OK.execute("grid", gridx, gridy)run.write_img(resultDataSetPath + '//' + outputKrigingTifName + '.tif',projection, geotransform, z)def removeFile(removeTifFileFolder, removeTifName):filePathAndName = removeTifFileFolder + '//' + removeTifNameos.remove(filePathAndName)if __name__ == "__main__":# 样本点数据InputExcelPathAndName = "E:\\RemoteSensing\\2018Sample_ALL.xlsx"# 保存x坐标的表头X_coordinateField = 'x_coor'# 保存y坐标的表头Y_coordinateField = 'y_coor'# 插值点的值Sample_valueField = 'TN'# 参考tif（提供插值的栅格数和行列号）# 参考tif所在的文件夹tifDataPath = 'E:\\RemoteSensing\\DataSet\\'# 参考tif的名字tifName = 'NDVI_Resample150m.tif'# 克里金插值结果保存的文件夹resultDataSetPath = 'E:\\RemoteSensing\\resultDataSet'# 克里金插值结果保存的名字outputKrigingTifName = 'KrigingInterpolate'# ************************** 运行克里金插值函数 ************************ #KrigingInterpolate(InputExcelPathAndName, X_coordinateField, Y_coordinateField, Sample_valueField,tifDataPath, tifName, resultDataSetPath, outputKrigingTifName)# 删除不需要的tifremoveTifFileFolder = 'E:\\RemoteSensing\\resultDataSet'removeTifName = 'regressionTN.tif'removeFile(removeTifFileFolder, removeTifName)

python：克里金插值相关推荐

Python进行克里金插值可视化
处理数据并进行克里金插值 # 处理数据后进行克里金插值 import pandas as pd pm25_region=pd.read_csv(r"C:\Users\LHW\Desktop\ ...
【Matlab 克里金】克里金插值
前言早在去年就开始要学 matlab,直到昨天,有明确的需要,和可视化的效果,才定下心来学: 从克里金插值上手,把该会的都捋了一遍. 学习路径克里金插值,有代码,有dacefit函数包,但是第一次 ...
基于主动学习和克里金插值的空气质量推测
基于主动学习和克里金插值的空气质量推测常慧娟, 於志文, 於志勇, 安琦, 郭斌西北工业大学计算机学院,陕西西安 710072 福州大学数学与计算机科学学院,福建福州 350108 摘要 ...
gstat | 空间插值（四）——克里金插值之协同克里金和交叉验证
本篇是介绍克里金插值的第三篇推文,也是最后一篇. 因为前面两篇使用的数据中已知点的样本太少,本篇使用gstat工具包说明文档中的数据集.该数据集来自sp工具包. 实际上,gstat工具包的方法对sp和 ...
gstat | 空间插值（三）——克里金插值之泛克里金和简单克里金
本篇接着上篇继续介绍克里金插值.首先加载相关工具包和上篇使用的示例数据: library(gstat) library(sf) library(tidyverse) library(readxl) l ...
gstat | 空间插值（二）——克里金插值之普通克里金
说明:昨天的推文误把可吸入颗粒物当作PM2.5,实应该为PM10,这里修正后重发. 从本篇开始计划分三篇介绍克里金插值.与反距离权重插值不同,克里金插值是无偏估计,其中也涉及到模型估计.本篇先对普通克 ...
ArcGIS中使用协同克里金插值（co-kriging interplotation ）对气象数据插值
ArcGIS中如何使用协同克里金插值(co-kriging interplotation )对气象数据插值 ANUSPLIN气象站点数据插值局限性百度搜索ArcGIS 克里金插值搭建梯子搜索Arc ...
克里金插值---MATLAB程序
最近在研究克里金插值,拜读了@lanainluv的笔记,备受启发, 在这里做一些补充,并分享自己的代码,希望对各位有所帮助,有误的地方请批评指正!(有帮助的话点个赞吧~) @lanainluv的原文链 ...
克里金插值详细步骤_openlayers4 入门开发系列之前端动态渲染克里金插值 kriging 篇（附源码下载）...
前言 openlayers4 官网的 api 文档介绍地址 openlayers4 api,里面详细的介绍 openlayers4 各个类的介绍,还有就是在线例子:openlayers4 官网在线例子 ...

python：克里金插值

python：克里金插值相关推荐

最新文章

热门文章