# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3 (ipykernel)
#     language: python
#     name: python3
# ---

# # 空间统计分析之克里金插值
#
# 前面章节主要是讨论一些基础gis数据处理运算和地图可视化方面的内容，从本小节所属章节开始涉及记录一些常用的空间数据分析方法，比如统计计算方法等。
#
# 本小节简单介绍一种经典的空间数据插值方法——克里金插值。
#
# ## 克里金(Kriging)插值的原理

# 克里金插值变种很多，本文主要简单介绍基本克里金插值的原理，内容主要来自[这里](https://xg1990.com/blog/archives/222)。

# ### 0 前言 --从反距离插值(IDW)说起

# + [markdown] tags=[]
# 空间插值问题，就是在已知空间上若干离散点 $(x_i,y_i)$ 的某一属性(如气温，海拔)的观测值$z_i=z(x_i,y_i)$的条件下，估计空间上任意一点$(x,y)$的属性值的问题。
#
# 直观来讲，根据地理学第一定律，  
# ```
# All attribute values on a geographic surface are related to each other, but closer values are more strongly related than are more distant ones.
# ```
#
# 地理属性有空间相关性，相近的事物会更相似。由此人们发明了反距离插值，对于空间上任意一点 $(x,y)$ 的属性 $z=z(x,y)$ ，定义反距离插值公式估计量
# $$\hat{z} = \sum^{n}_{i=1}{\frac{1}{d^\alpha}z_i}$$
#
# 其中$\alpha $通常取1或者2。  
#
# 即，用空间上所有已知点的数据加权求和来估计未知点的值，权重取决于距离的倒数（或者倒数的平方）。那么，距离近的点，权重就大；距离远的点，权重就小。  
# 反距离插值可以有效的基于地理学第一定律估计属性值空间分布，但仍然存在很多问题：  
#     (1)$\alpha $ 的值不确定  
#     (2)用倒数函数来描述空间关联程度不够准确   
#     
# 因此更加准确的克里金插值方法被提出来了。
# -

# ### 1 克里金插值的定义

# 克里金法（Kriging  method）是依据**协方差函数**对**随机过程**/随机场进行空间建模和预测（插值）的回归算法。   
# 相比反距离插值，克里金插值公式更加抽象
# $$\hat{z_o} = \sum^{n}_{i=0}{\lambda_iz_i} $$
# 其中$\hat{z_o}$是点$(x_o,y_o)$处的估计值，即$z_o=z(x_o,y_o)$  
# 这里的$\lambda_i$是权重系数。它同样是用空间上所有已知点的数据加权求和来估计未知点的值。但权重系数并非距离的倒数，而是能够满足点$(x_o,y_o)$处的估计值$\hat{z_o}$与真实值$z_o$的差最小的一套最优系数，即
# $$\min_{\lambda_i} Var(\hat{z_o}-z_o)$$
# 同时满足无偏估计的条件
# $$E(\hat{z_o}-z_o)=0$$
# 无偏指的是估计值和实际值之差的期望等于零，最优指的是估计值和实际值的方差最小。  
# 简单来说克里金（kriging）插值干的事情就是：已经有一些坐标和这些坐标处的真实值，我们可以称之为采样点。然后把这些采样点输入到克里金插值中，我们就可以用来估计其他未知位置处的值。

# ### 2 算法

#   -普通克里金（Ordinary Kriging, OK）  
# **改进算法：**  
#   -泛克里金（Universal Kriging, UK）  
#   -协同克里金（Co-Kriging, CK）  
#   -析取克里金（Disjunctive Kriging, DK）  
# **混合算法：**   
#   -回归克里金（regression-Kriging）  
#   -神经网络克里金（neural Kriging）  
#   -贝叶斯克里金（Bayesian Kriging）  

# ### 3 假设条件

# + [markdown] tags=[]
# **不同的克里金插值方法的主要差异就是假设条件不同**。本文仅介绍普通克里金插值的假设条件与应用。  
# 普通克里金插值的假设条件为，空间属性z是均一的。对于空间任意一点$(x,y)$，都有同样的期望c与方差$\sigma^2$。  即对任意点$(x,y)$都有$$E[z(x,y)] = E[z] = c$$
#
# $$Var[z(x,y)] = \sigma^2$$
# 换一种说法：任意一点处的值$z(x,y)$，都由区域平均值c和该点的随机偏差$R(x,y)$组成，即
#
# $$z(x,y)=E[z(x,y)] + R(x,y)= c + R(x,y)$$
#
# 其中$R(x,y)$表示点$(x,y)$处的偏差，其方差均为常数
#
# $$Var[R(x,y)] = \sigma^2$$
#
#
#
#
#
# -

# ### 4 无偏约束条件

# 先分析无偏估计条件$E(\hat{z_o}-z_o)=0$，将$\hat{z_o} = \sum^{n}_{i=0}{\lambda_iz_i}$带入则有
#
# $E(\sum^{n}_{i=0}{\lambda_iz_i}- z_o)=0$
#
# 又因为对任意的z都有$E[z] = c$，则
#
# $c \sum^{n}_{i=0}{\lambda_i}- c=0$
#
# 即
#
# $\sum^{n}_{i=0}{\lambda_i} = 1$
#
# 这是$\lambda_i$的约束条件之一。

# ### 5 优化目标/代价函数J

# 再分析估计误差$Var(\hat{z_o}-z_o)$。为方便公式推理，用符号$J$表示，即
#
# $J = Var(\hat{z_o}-z_o)$
#
# 则有

# $$J= Var(\sum^{n}_{i=0}{\lambda_iz_i} - z_o) = Var(\sum^{n}_{i=0}{\lambda_iz_i}) - 2 Cov(\sum^{n}_{i=0}{\lambda_iz_i}, z_o) + Cov(z_o, z_o) = \sum^{n}_{i=0}\sum^{n}_{j=0}{\lambda_i\lambda_jCov( z_i, z_j)} - 2 \sum^{n}_{i=0}{\lambda_iCov(z_i, z_o)} + Cov(z_o, z_o) $$

# 为简化描述，定义符号 $C_{ij} = Cov(z_i,z_j) = Cov(R_i,R_j)$，这里$R_i = z_i - c$，即点$(x_i,y_i)$处的属性值相对于区域平均属性值的偏差。
#
# 则有
#
# $J = \sum^{n}_{i=0}\sum^{n}_{j}{\lambda_i\lambda_jC_{ij}} - 2 \sum^{n}_{i=0}{\lambda_iC_{io}} + C_{oo} $

# ### 6 代价函数的最优解

# 再定义半方差函数 $r_{ij} = \sigma^2 -C_{ij}$，带入$J$中，有$$J = \sum^{n}_{i=0}\sum^{n}_{j=0}{\lambda_i\lambda_j(\sigma^2 - r_{ij})} - 2 \sum^{n}_{i=0}{\lambda_i(\sigma^2 - r_{io})} + \sigma^2 - r_{oo} =\sum^{n}_{i=0}\sum^{n}_{j=0}{\lambda_i\lambda_j(\sigma^2)} -\sum^{n}_{i=0}\sum^{n}_{j=0}{\lambda_i\lambda_j( r_{ij})}-2\sum^{n}_{i=0}{\lambda_i(\sigma^2)}+2 \sum^{n}_{i=0}{\lambda_i(r_{io})}+\sigma^2 - r_{oo}  $$
# 考虑到$\sum^{n}_{i=0}{\lambda_i} = 1$
# $J = \sigma^2-\sum^{n}_{i=0}\sum^{n}_{j}{\lambda_i\lambda_j(r_{ij})}-2 \sigma^2 +2 \sum^{n}_{i=0}{\lambda_i(r_{io})}+ \sigma^2 - r_{oo}=2 \sum^{n}_{i=0}{\lambda_i(r_{io})} -\sum^{n}_{i=0}\sum^{n}_{j=0}{\lambda_i\lambda_j(r_{ij})} - r_{oo}  $
# 我们的目标是寻找使$J$最小的一组 $\lambda_i$，且$J$是$\lambda_i$的函数，因此直接将$J$对$\lambda_i$求偏导数令其为0即可。即
# $\frac{\partial J}{\partial \lambda_i}= 0;i=1,2,\cdots$
#
# 但是要注意的是，我们要保证求解出来的最优 $\lambda_i$ 满足公式$\sum^{n}_{i=0}{\lambda_i} = 1$，这是一个带约束条件的最优化问题。使用拉格朗日乘数法求解，求解方法为构造一个新的目标函数$$J + \phi(\sum^{n}_{i=0}{\lambda_i}-1)$$
# 其中$\phi$是拉格朗日乘数。求解使这个代价函数最小的参数集${\phi,\lambda_1,\lambda_2,\cdots,\lambda_n}$，则能满足其在$\sum^{n}_{i=0}{\lambda_i} = 1$约束下最小化$J$。即
#
# $$\frac{\partial(J + \phi(\sum^{n}_{i=0}{\lambda_i}-1))}{\partial \lambda_i} = 0;i=1,2,\cdots,n\\ \frac{\partial(J + \phi(\sum^{n}_{i=0}{\lambda_i}-1))}{\partial \phi} = 0  $$
#
# $$ \frac{\partial (2 \sum^{n}_{i=0}{\lambda_i(r_{io})} - \sum^{n}_{i=0}\sum^{n}_{j}{\lambda_i\lambda_j(r_{ij})} - r_{oo})}{\partial \lambda_i}= 0;i=1,2,\cdots,$$$$ \frac{\partial (\partial 2 \sum^{n}_{i=0}{\lambda_i(r_{io})} - \sum^{n}_{i=0}\sum^{n}_{j}{\lambda_i\lambda_j(r_{ij})} - r_{oo})}{\partial \phi} = 0 $$
#
# $$2r_{io} - \sum^{n}_{j=1}{(r_{ij}+r_{ji})\lambda_j}=0;i=1,2,\cdots,$$$$ \sum^{n}_{i=0}{\lambda_i} = 1$$
# 由于$C_{ij}=Cov(z_i,z_j)=C_{ji}$，因此同样地$r_{ij}=r_{ji}$，那么有
#
# $$ r_{io} - \sum^{n}_{j=1}{r_{ij}\lambda_j} = 0;i=1,2,\cdots,$$$$ \sum^{n}_{i=0}{\lambda_i} = 1 $$
# 式子中半方差函数$r_{ij}$十分重要，最后会详细解释其计算与定义
# 在以上计算中我们得到了对于求解权重系数$\lambda_j$的方程组。写成线性方程组的形式就是：
# $$ r_{11}\lambda_1 + r_{12}\lambda_2 + \cdots + r_{1n}\lambda_n + \phi= r_{1o}$$   $$r_{21}\lambda_1 + r_{22}\lambda_2 + \cdots + r_{2n}\lambda_n + \phi= r_{2o}$$ $$ r_{n1}\lambda_1 + r_{n2}\lambda_2 + \cdots + r_{nn}\lambda_n + \phi= r_{no}$$ $$ \lambda_1 + \lambda_2 + \cdots + \lambda_n = 1$$
# 写成矩阵形式即为
# $\begin{bmatrix}r_{11}&r_{12}&\cdots&r_{1n}&1\\ r_{21}&r_{22}&\cdots&r_{2n}&1\\\cdots&\cdots&\cdots&\cdots&\cdots\\r_{n1}&r_{n2}&\cdots&r_{nn}&1\\1&1&\cdots&1&0\end{bmatrix}\begin{bmatrix} \lambda_1\\ \lambda_2\\\cdots\\\lambda_n\\0\end{bmatrix}=\begin{bmatrix} r_{1o}\\ r_{2o}\\\cdots\\r_{no}\\1\end{bmatrix}$
# 对矩阵求逆即可求解。
#
# 唯一未知的就是上文中定义的半方差函数$r_{ij}$，接下来将详细讨论

# ### 7 半方差函数

# 上文中对半方差函数的定义为
#
# $r_{ij} = \sigma^2 -C_{ij}$
#
# 其等价形式为
#
# $r_{ij} = \frac{1}{2}E[(z_i-z_j)^2]$
#
# 这也是半方差函数名称的来由，接下来证明这二者是等价的：
#
# 根据上文定义 $R_i = z_i - c$，有$z_i-z_j = R_i - R_j$，则
#
# $$ r_{ij} = \frac{1}{2}E[(R_i-R_j)^2]= \frac{1}{2}E[R_i^2-2R_iR_j+R_j^2]= \frac{1}{2}E[R_i^2]+\frac{1}{2}E[R_j^2]-E[R_iR_j]  $$
#
# 又因为：
#
# $E[R_i^2] =E[R_j^2] = E[(z_i - c)^2] = Var(z_i) = \sigma^2 $
#
# $E[R_iR_j] = E[(z_i - c)(z_j-c)] = Cov(z_i,z_j) = C_{ij}$
#
# 于是有
#
# $$ r_{ij} = \frac{1}{2}E[(z_i-z_j)^2]= \frac{1}{2}E[R_i^2]+\frac{1}{2}E[R_j^2]-E[R_iR_j]= \frac{1}{2}\sigma^2+\frac{1}{2}\sigma^2- C_{ij}=\sigma^2 -C_{ij} $$
#
# $ \sigma^2 -C_{ij} = \frac{1}{2}E[(z_i-z_j)^2]$得证，现在的问题就是如何计算
#
# $$r_{ij} = \frac{1}{2}E[(z_i-z_j)^2]$$
#
# 这时需要用到地理学第一定律，空间上相近的属性相近。$(r_{ij} = \frac{1}{2}(z_i-z_j)^2$表达了属性的相似度；空间的相似度就用距离来表达，定义i与j之间的几何距离
# $d_{ij} = d(z_i,z_j) = d( (x_i,y_i), (x_j,y_j)) = \sqrt{(x_i-x_j)^2 + (y_i - y_j)^2}$
#
# 克里金插值假设$r_{ij}$与$d_{ij}$存在着函数关系，这种函数关系可以是线性、二次函数、指数、对数关系。为了确认这种关系，我们需要首先对观测数据集
#
# $\{z(x_1,y_1),z(x_2,y_2),z(x_3,y_3),\cdots,z(x_{n-1},y_{n-1}),z(x_n,y_n)\}$
#
# 计算任意两个点的 距离$d_{ij}= \sqrt{(x_i-x_j)^2 + (y_i - y_j)^2}$和 半方差 $\sigma^2 -C_{ij} =\frac{1}{2}E[(z_i-z_j)^2]$，这时会得到$n^2$个$(d_{ij}, r_{ij})$的数据对。
#
# 将所有的$d$和$r$绘制成散点图，寻找一个最优的拟合曲线拟合$d$与$r$的关系，得到函数关系式
#
# $r = r(d)$
#
# 那么对于任意两点$(x_i,y_i), (x_j,y_j)$，先计算其距离$d_{ij}$，然后根据得到的函数关系就可以得到这两点的半方差$r_{ij}$

# ### 8 简单克里金（simple kriging）与普通克里金（ordinary kriging）的区别

# 以上介绍的均为普通克里金（ordinary kriging）的公式与推理。
#
# 事实上普通克里金插值还有简化版，即简单克里金（simple kriging）插值。二者的差异就在于如何定义插值形式：
#
# 上文讲到，普通克里金插值形式为
#
# $\hat{z_o} = \sum^{n}_{i=0}{\lambda_iz_i}$
#
# 而简单克里金的形式则为
#
# $\hat{z_o} - c= \sum^{n}_{i=0}{\lambda_i(z_i-c)}$
#
# 这里的符号c在上文介绍过了，是属性值的数学期望，即$E[z] = c$。也就是说，在普通克里金插值中，认为未知点的属性值是已知点的属性值的加权求和；而在简单克里金插值中，假设未知点的属性值相对于平均值的偏差是已知点的属性值相对于平均值的偏差的加权求和，用公式表达即为：
#
# $\hat{R_o} = \sum^{n}_{i=0}{\lambda_iR_i}$
#
# 这里的$R_i$在上文定义过了：$R_i = z_i - c$。
#
# 但是为什么这样的克里金插值称为“简单克里金”呢？由于有假设$E[z] = c$，也就是说$E(R_i + c) = c$，即$E(R_i) = 0$。那么上面的公式$\hat{R_o} = \sum^{n}_{i=0}{\lambda_iR_i}$两边的期望一定相同，那么在求解未知参数$\lambda_i$就不需要有无偏约束条件$\sum^{n}_{i=0}{\lambda_i} = 1$。换句话说，这样的估计公式天生就能满足无偏条件。因此它被称为简单克里金。
#
# 从在上文（第5节优化目标/代价函数J）中可以知道，优化目标的推理和求解过程是通过对属性值相对于期望的偏差量$R_i$进行数学计算而进行的。也就是说这两种克里金插值方法虽然插值形式不一样，求解方法是一样的，重要的区别是简单克里金插值不需要约束条件$\sum^{n}_{i=0}{\lambda_i} = 1$，求解方程组为：
# $$r_{11}\lambda_1 + r_{12}\lambda_2 + \cdots + r_{1n}\lambda_n + \phi= r_{1o}  $$
#
# $$r_{21}\lambda_1 + r_{22}\lambda_2 + \cdots + r_{2n}\lambda_n + \phi= r_{2o}  $$
#
# $$r_{n1}\lambda_1 + r_{n2}\lambda_2 + \cdots + r_{nn}\lambda_n + \phi= r_{no} $$
#
# 还有更重要的一点，简单克里金的插值公式为：
#
# $\hat{z_o} = \sum^{n}_{i=0}{\lambda_i(z_i-c)}+c$
#
# 换句话说，在计算未知点属性值$\hat{z_o}$前，需要知道该地区的属性值期望c。事实上我们在进行插值前很难知道这个地区的真实属性值期望。有些研究者可能会采用对观测数据简单求平均的方法计算期望值c，而考虑到空间采样点位置代表性可能有偏差（比如采样点聚集在某一小片地区，没有代表性），简单平均估计的期望也可能是有偏差的。这是简单克里金方法的局限性。

# ### 9 小结

# 总的来说，进行克里金插值分为这几个步骤：  
# 1、对于**观测数据**，两两计算距离与半方差  
# 2、寻找一个拟合曲线拟合距离与半方差的关系，从而能根据**任意距离计算出相应的半方差**  
# 3、计算出所有已知点之间的半方差$r_{ij}$    
# 4、对于未知点$z_0$，计算它到所有已知点$z_i$的半方差$r_{io}$    
# 5、求解第五节中的方程组，得到最优系数$\lambda_i$    
# 6、使用最优系数对已知点的属性值进行加权求和，得到未知点$z_o$的估计值
#

# # Python Kriging
#
# 在Python里，有两个比较常用的克里金插值包，**pykrige和pykriging**，此外 PySAL中也集成了包括克里金插值在内的很多空间统计方法，后续再补充。
#
# python中Kriging算法: 
#
# ``` OrdinaryKriging```： 2D ordinary kriging with estimated mean  
#   
# ``` UniversalKriging```:2D universal kriging providing drift terms  
#   
# ```OrdinaryKriging3D: ```3D ordinary kriging  
#  
# ```UniversalKriging3D: ```3D universal kriging  
#   
# ```RegressionKriging:``` An implementation of Regression-Kriging 
#   
# ```ClassificationKriging:``` An implementation of Simplicial Indicator Kriging  

# + [markdown] tags=[]
# ## Using pyKrige
# -

# ### Ordinary Kriging Example

#下面代码给出了使用普通克里金进行插值的一个简单例子，其他类型的克里金插值类似。
from pykrige.ok import OrdinaryKriging
import numpy as np
from matplotlib import pyplot as plt

# 已知采样点的数据，是坐标（x，y）和坐标对应的值
# 矩阵中第一列是x,第二列是y,第三列是坐标对应的值
data = np.array(
    [
        [0.1, 0.1, 0.9],
        [0.2, 0.1, 0.8],
        [0.1, 0.3, 0.9],
        [0.5, 0.4, 0.5],
        [0.3, 0.3, 0.7],
    ])

# 网格
x_range = 0.6
y_range = 0.6
range_step = 0.1
gridx = np.arange(0.0, x_range, range_step) #三个参数的意思：范围0.0 - 0.6 ，每隔0.1划分一个网格
gridy = np.arange(0.0, y_range, range_step)

ok3d = OrdinaryKriging(
    data[:, 0],
    data[:, 1], 
    data[:, 2], 
    variogram_model="linear") # 模型
# variogram_model是变差函数模型，pykrige提供 linear, power, gaussian, spherical, exponential, hole-effect几种variogram_model可供选择，默认的为linear模型。
# 使用不同的variogram_model，预测效果是不一样的，应该针对自己的任务选择合适的variogram_model。

k3d1, ss3d = ok3d.execute("grid", gridx, gridy) # k3d1是结果，给出了每个网格点处对应的值
print(np.round(k3d1,2))

# 绘图
fig, (ax1) = plt.subplots(1)
ax1.imshow(k3d1, origin="lower")
ax1.set_title("ordinary kriging")
plt.tight_layout()
plt.show()
#黄色是值比较大的区域

# ### Universal Kriging Example

from pykrige  import UniversalKriging
UK = UniversalKriging(
    data[:, 0],
    data[:, 1],
    data[:, 2],
    variogram_model="linear",
     drift_terms=["regional_linear"],
)
z, ss = UK.execute("grid", gridx, gridy)
plt.imshow(z)
plt.show()


# ### Three-Dimensional Kriging Example

# +
from pykrige import OrdinaryKriging3D
from pykrige  import UniversalKriging3D
data = np.array(
    [
        [0.1, 0.1, 0.3, 0.9],
        [0.2, 0.1, 0.4, 0.8],
        [0.1, 0.3, 0.1, 0.9],
        [0.5, 0.4, 0.4, 0.5],
        [0.3, 0.3, 0.2, 0.7],
    ]
)
gridx = np.arange(0.0, 0.6, 0.05)
gridy = np.arange(0.0, 0.6, 0.01)
gridz = np.arange(0.0, 0.6, 0.1)
ok3d = OrdinaryKriging3D(
    data[:, 0], data[:, 1], data[:, 2], data[:, 3], variogram_model="linear"
)
k3d1, ss3d = ok3d.execute("grid", gridx, gridy, gridz)

uk3d = UniversalKriging3D(
    data[:, 0],
    data[:, 1],
    data[:, 2],
    data[:, 3],
    variogram_model="linear",
    drift_terms=["regional_linear"],
)
k3d2, ss3d = uk3d.execute("grid", gridx, gridy, gridz)

# -

zg, yg, xg = np.meshgrid(gridz, gridy, gridx, indexing="ij")
uk3d = UniversalKriging3D(
    data[:, 0],
    data[:, 1],
    data[:, 2],
    data[:, 3],
    variogram_model="linear",
    drift_terms=["specified"],
    specified_drift=[data[:, 0], data[:, 1], data[:, 2]],
)
k3d3, ss3d = uk3d.execute(
    "grid", gridx, gridy, gridz, specified_drift_arrays=[xg, yg, zg]
)

func = lambda x, y, z: x
uk3d = UniversalKriging3D(
    data[:, 0],
    data[:, 1],
    data[:, 2],
    data[:, 3],
    variogram_model="linear",
    drift_terms=["functional"],
    functional_drift=[func],
)
k3d4, ss3d = uk3d.execute("grid", gridx, gridy, gridz)

fig, (ax1, ax2, ax3, ax4) = plt.subplots(4)
ax1.imshow(k3d1[:, :, 0], origin="lower")
ax1.set_title("ordinary kriging")
ax2.imshow(k3d2[:, :, 0], origin="lower")
ax2.set_title("regional lin. drift")
ax3.imshow(k3d3[:, :, 0], origin="lower")
ax3.set_title("specified drift")
ax4.imshow(k3d4[:, :, 0], origin="lower")
ax4.set_title("functional drift")
plt.tight_layout()
plt.show()

# ### Geometric example

#A small example script showing the usage of the 'geographic' coordinates type for ordinary kriging on a sphere.
from pykrige  import OrdinaryKriging
import numpy as np
from matplotlib import pyplot as plt
# Make this example reproducible:
np.random.seed(89239413)
N = 7
lon = 360.0 * np.random.random(N)
lat = 180.0 / np.pi * np.arcsin(2 * np.random.random(N) - 1)
z = 3.5 * np.random.rand(N) + 2.0

# Generate a regular grid with 60° longitude and 30° latitude steps:
grid_lon = np.linspace(0.0, 360.0, 7)
grid_lat = np.linspace(-90.0, 90.0, 7)


# Create ordinary kriging object:
OK = OrdinaryKriging(
    lon,
    lat,
    z,
    variogram_model="linear",
    verbose=False,
    enable_plotting=False,
    coordinates_type="geographic",
)
OK

# Execute on grid:
z1, ss1 = OK.execute("grid", grid_lon, grid_lat)
z1

# Create ordinary kriging object ignoring curvature:
OK = OrdinaryKriging(
    lon, lat, z, variogram_model="linear", verbose=False, enable_plotting=False
)
# Execute on grid:
z2, ss2 = OK.execute("grid", grid_lon, grid_lat)
z2

# Print data at equator (last longitude index will show periodicity):
print("Original data:")
print("Longitude:", lon.astype(int))
print("Latitude: ", lat.astype(int))
print("z:        ", np.array_str(z, precision=2))
print("\nKrige at 60° latitude:\n======================")
print("Longitude:", grid_lon)
print("Value:    ", np.array_str(z1[5, :], precision=2))
print("Sigma²:   ", np.array_str(ss1[5, :], precision=2))
print("\nIgnoring curvature:\n=====================")
print("Value:    ", np.array_str(z2[5, :], precision=2))
print("Sigma²:   ", np.array_str(ss2[5, :], precision=2))

fig, (ax1, ax2) = plt.subplots(1, 2)
ax1.imshow(z1, extent=[0, 360, -90, 90], origin="lower")
ax1.set_title("geo-coordinates")
ax2.imshow(z2, extent=[0, 360, -90, 90], origin="lower")
ax2.set_title("non geo-coordinates")
plt.show()

# ## Using pyKriging

# pyKriging 旨在简化创建代理模型的过程。以下示例演示了如何创建抽样计划,评估这些位置的测试函数,创建和训练克里金模型并添加填充点以减少模型的均方误差 (MSE)。

import pyKriging  
from pyKriging.krige import kriging  
from pyKriging.samplingplan import samplingplan

# +
# The Kriging model starts by defining a sampling plan, we use an optimal Latin Hypercube (最优的拉丁超立方体)here
sp = samplingplan(2) 

X = sp.optimallhc(20)
X
# -

# Next, we define the problem we would like to solve
testfun = pyKriging.testfunctions().branin  
y = testfun(X)

# Now that we have our initial data, we can create an instance of a Kriging model
k = kriging(X, y, testfunction=testfun, name='simple')  
k.train()
# Now, five infill points are added. Note that the model is re-trained after each point is added
numiter = 5 
for i in range(numiter):
    #print "Infill iteration {0} of {1}....".format(i + 1, numiter)
    newpoints = k.infill(1)
    for point in newpoints:
        k.addPoint(point, testfun(point)[0])
    k.train()
# And plot the results
k.plot()  

StarcoderdataJupyterScriptsDedupFiltered