基于 Python（gma） 的 克里金（Kriging）法插值的主要过程
2024-06-12 18:49:49
由于克里金插值的复杂性，本文不再对其原理进行介绍。详情可自行百度。
- 本算法基于 Python 的开源克里金插值包 pykrige。
- 但本算法已对其进行改造，以使其符合 gma 的整体逻辑。
- 本算法已不包含任何 pykrige 的原始代码。
原始代码构建
from gma.algorithm.spmis.interpolate import *class Kriging(IPolate):'''克里金法插值'''# 继承 gma 的标准插值类 IPolate。本处不再做详细说明。def __init__(self, Points, Values, Boundary = None, Resolution = None, InProjection = 'WGS84', VariogramModel = 'Linear',VariogramParameters = None,**kwargs):IPolate.__init__(self, Points, Values, Boundary, Resolution, InProjection)self.eps = eps# 初始化半变异函数及参数self.VariogramModel, self.VParametersList = GetVariogramParameters(VariogramModel, VariogramParameters)self.VariogramFUN = getattr(variogram, self.VariogramModel)if self.VParametersList is None:self.VParametersList = self._INITVariogramModel(**kwargs)# 调整输入数据if self.GType == 'PROJCS':self.Center = (self.Points.min(axis = 0) + self.Points.max(axis = 0)) * 0.5self.AnisotropyScaling = AnisotropyScalingself.AnisotropyAngle = AnisotropyAngleself.DistanceMethod = cdistelse:# 方便后期优化self.Center = np.array([0,0])self.AnisotropyScaling = 1.0self.AnisotropyAngle = 0.0self.DistanceMethod = GreatCircleDistanceself.AdjustPoints = AdjustAnisotropy(self.Points, self.Center, [self.AnisotropyScaling], [self.AnisotropyAngle])self.XYs = AdjustAnisotropy(self.XYs, self.Center,[self.AnisotropyScaling], [self.AnisotropyAngle])def _INITVariogramModel(self, **kwargs):'''初始化参数'''if 'NLags' in kwargs:NLags = kwargs['NLags']initialize.ValueType(NLags, 'pint')else:NLags = 6if 'Weight' in kwargs:Weight = ToNumericArray(kwargs['Weight']).flatten().astype(bool)[0]else:Weight = FalseLags, SEMI = INITVariogramModel(self.Points, self.Values, NLags, self.GType)# 为求解自动参数准备if self.VariogramModel == "Linear":X0 = [np.ptp(SEMI) / np.ptp(Lags), np.min(SEMI)]BNDS = ([0.0, 0.0], [np.inf, np.max(SEMI)])elif self.VariogramModel == "Power":X0 = [np.ptp(SEMI) / np.ptp(Lags), 1.1, np.min(SEMI)]BNDS = ([0.0, 0.001, 0.0], [np.inf, 1.999, np.max(SEMI)])else:X0 = [np.ptp(SEMI), 0.25 * np.max(Lags), np.min(SEMI)]BNDS = ([0.0, 0.0, 0.0], [10.0 * np.max(SEMI), np.max(Lags), np.max(SEMI)])# 最小二乘法求解默认参数def _VariogramResiduals(Params, X, Y, Weight):if Weight:Weight = 1.0 / (1.0 + np.exp(-2.1972 / (0.1 * np.ptp(X)) * (0.7 * np.ptp(X) + np.min(X) - x))) + 1 else:Weight = 1return (self.VariogramFUN(X, *Params) - Y) * WeightRES = least_squares(_VariogramResiduals, X0, bounds = BNDS, loss = "soft_l1",args = (Lags, SEMI, Weight))return RES.xdef _GetKrigingMatrix(self):"""获取克里金矩阵"""LDs = self.DistanceMethod(self.AdjustPoints, self.AdjustPoints)A = -self.VariogramFUN(LDs, *self.VParametersList)A = np.pad(A, (0, 1), constant_values = 1)# 填充主对角线np.fill_diagonal(A, 0.0)return Adef _UKExec(self, A, LDs, SearchRadius):"""泛克里金求解"""Args = LDs.argsort(axis = 1)[:,:SearchRadius]Values = self.Values[Args.T].T# A 的逆矩阵AInv = inv(A)B = -self.VariogramFUN(LDs, *self.VParametersList)B[np.abs(LDs) <= self.eps] = 0.0B = np.pad(B, ((0,0),(0,1)), constant_values = 1)X = np.dot(B, AInv)B = B[np.ogrid[:len(B)], Args.T].TX = X[np.ogrid[:len(X)], Args.T].TX = X / X.sum(axis = 1, keepdims = True)UKResults = np.sum(X * Values, axis = 1), np.sum((X * -B), axis = 1)return UKResultsdef _OKExec(self, A, LDs, SearchRadius):"""普通克里金求解"""Args = LDs.argsort(axis = 1)[:,:SearchRadius]LDs = LDs[np.ogrid[:len(LDs)], Args.T].TB = -self.VariogramFUN(LDs, *self.VParametersList)B[np.abs(LDs) <= self.eps] = 0.0B = np.pad(B, ((0,0),(0,1)), constant_values = 1)OKResults = np.zeros([2, len(LDs)])for i, b in enumerate(B):BSelector = Args[i] ASelector = np.append(BSelector, len(self.AdjustPoints))a = A[ASelector[:, None], ASelector] x = solve(a, b)OKResults[:, i] = x[:SearchRadius].dot(self.Values[BSelector]), -x.dot(b)return OKResultsdef Execute(self, SearchRadius = 12, KMethod = 'Ordinary'):'''克里金插值'''initialize.ValueType(SearchRadius, 'pint')SearchRadius = np.min([SearchRadius, len(self.AdjustPoints)])A = self._GetKrigingMatrix()LDs = self.DistanceMethod(self.XYs, self.AdjustPoints)if KMethod not in ['Universal', 'Ordinary']:raise ValueError("Undefined Kriging method. Please select 'Universal' or 'Ordinary'!")elif KMethod == 'Universal':KResults = self._UKExec(A, LDs, SearchRadius)else:KResults = self._OKExec(A, LDs, SearchRadius)NT = namedtuple('Kriging', ['Data', 'SigmaSQ', 'Transform'])return NT(KResults[0].reshape(self.YLAT, self.XLON),KResults[1].reshape(self.YLAT, self.XLON), self.Transform)
插值应用
在 gma 1.0.13.15 之后的版本可以直接引用。
import gma
import pandas as pdData = pd.read_excel("Interpolate.xlsx")
Points = Data.loc[:, ['经度','纬度']].values
Values = Data.loc[:, ['值']].values# 普通克里金（球面函数模型）插值
KD = Kriging(Points, Values,Resolution = 0.05, VariogramModel = 'Spherical', VariogramParameters = None,InProjection = 'EPSG:4326').Execute(KMethod = 'Ordinary')# 泛克里金类似，这里不做示例gma.rasp.WriteRaster(r'.\gma_OKriging.tif',KD.Data,Projection = 'WGS84',Transform = KD.Transform, DataType = 'Float32')
结果对比
与 ArcGIS Ordinary Kriging 插值结果（重分类后）对比：
与 pykrige 包 Universal Kriging 插值结果（重分类后）对比：
基于 Python（gma） 的 克里金（Kriging）法插值的主要过程相关推荐
- GEE：克里金 Kriging 空间插值（以陕西省2013年生物量为例）
作者:CSDN @ _养乐多_ 本文记录了在Google Earth Engine(GEE)平台上进行 Kriging 插值的介绍和代码案例.本文通过选取的2013年陕西省生物量样本点数据为例,利用 ...
- python 克里金空间插值_Python克里金(Kriging)插值计算及可视化绘制
前面两篇推文我们分别介绍了使用Python和R进行IDW(反距离加权法) 插值的计算及结果的可视化过程,详细内容可见如下: 本期推文,我们将介绍如何使用Python进行克里金(Kriging)插值计算 ...
- 基于主动学习和克里金插值的空气质量推测
基于主动学习和克里金插值的空气质量推测 常慧娟, 於志文, 於志勇, 安琦, 郭斌 西北工业大学计算机学院,陕西 西安 710072 福州大学数学与计算机科学学院,福建 福州 350108 摘要 ...
- 克里金(Kriging)插值的原理与公式推导
转自:http://xg1990.com/blog/archives/222 学过空间插值的人都知道克里金插值,但是它的变种繁多.公式复杂,还有个半方差函数让人不知所云 本文讲简单介绍基本克里金插值的 ...
- matlab克里金插值法,克里金(Kriging)插值的原理与公式推导
学过空间插值的人都知道克里金插值,但是它的变种繁多.公式复杂,还有个半方差函数让人不知所云 本文讲简单介绍基本克里金插值的原理,及其推理过程,全文分为九个部分: 0.引言-从反距离插值说起 1.克里金 ...
- 克里金(kriging)模型的推导详解
Kriging模型理论推导 1.前言 2.条件 3.基础知识 3.1.方差的理解 3.2.概率密度函数 3.3.多元正态分布 4.理论推导 4.1 模型建立 4.2 模型预测 1.前言 简介:Krig ...
- 克里金(Kriging)插值的原理与公式推导_转
转载自 https://xg1990.com/blog/archives/222 如果没有加载出来,请点击查看图片
- [kriging](一)网上下载的kriging克里金的C++程序的初步调试
笔者在网上下载了一份克里金的C++程序,水平有限,正在逐步地调试中. 初步 克里金法现在在许多软件都已经有集成了,据笔者所知: arcgis : 看过有的arcgis培训视频里面简略介绍了里面的插值方 ...
- Kriging 克里金算法Java实现
引入依赖库 import java.util.ArrayList; import java.util.Arrays; import java.util.List; 定义一个类来表示二维坐标点 clas ...
最新文章
- 基于TensorRT优化的Machine Translation
- string生成固定长度的哈希_Redis 选择Hash还是String 存储数据？
- HTML5 未来不可阻挡的力量
- 从一个servlet转发到另一个servlet_javaweb02-创建第一个Servlet
- ubuntu tacacs 服务器安装启动
- Spring+SpringMvc+Mybatis框架集成搭建教程二(依赖配置及框架整合)
- mha数据备份_MHA学习笔记
- C#/.net 中的事件与代理
- 一个非常标准的Java连接Oracle数据库的示例代码
- 【5003】马遍历问题
- 本以为用的MyBatis框架就万无一失了，没想到还是被黑客注入了，我真的无语了！...
- Hibernate 查询缓存
- 一、数字图像处理分析与概述
- mysql冗余_如何合理使用数据库冗余字段的方法
- 如何快速制作App应用软件？国内有哪些比较好用的App制作平台？
- dede {dede:channel currentstyle 中使用~seotitle~
- 【c++中内存拷贝函数(C++ memcpy)详解】
- 【AI】PaddlePaddle实现自动语音识别
- HJ7 取近似值（重点关注）
- [问题已处理]-[nginx]-nginx 报错 could not build server_names_hash