克里金插值学习笔记

一、类定义：

OrdinaryKriging（x,y,z,variogram_model, variogram_model, variogram_parameters, variogram_function, nlags, weight, anisotropy_scaling, anisotropy_angle, verbose, enable_plotting, enable_statistics, coordinates_type)

参数说明：

x,y,z：坐标x,y，以及数值（z)

variogram_model（方差模型）：7个可选值，包括linear, power, gaussian, spherical, exponential, hole-effect，以及custom

缺省是linear，如选用custom则需提供方差模型的参数以及函数。提醒的是hole-effect模型仅仅在针对一维问题上，存在技术的正确性。

variogram_parameters（方差模型参数）：可选，list或dict类型，由所选模型决定，不提供的话则采用“软”（soft）L1范式最小化方案。

采用dict（字典）形式，参数如下：

linear - {'slope': slope, 'nugget': nugget}

power - {'scale': scale, 'exponent': exponent, 'nugget': nugget}

gaussian - {'sill': s, 'range': r, 'nugget': n} / {'psill': p, 'range': r, 'nugget':n}

spherical - {'sill': s, 'range': r, 'nugget': n} / {'psill': p, 'range': r, 'nugget':n}

exponential - {'sill': s, 'range': r, 'nugget': n} / {'psill': p, 'range': r, 'nugget':n}

hole-effect - {'sill': s, 'range': r, 'nugget': n} / {'psill': p, 'range': r, 'nugget':n}

基台值(Sill):代表变量在空间上的总变异性大小，psill = sill - nugget，所以有两种参数字典方式可选

采用list（列表）形式，参数如下：

linear - [slope, nugget]

power - [scale, exponent, nugget]

gaussian - [sill, range, nugget]

spherical - [sill, range, nugget]

exponential - [sill, range, nugget]

hole-effect - [sill, range, nugget]

列表形式，只能采用sill，不用使用psill值

模型选择 custom（定制），方差参数选项是必须提供，因为方差模型不会自动适应数据，而且参数必须以特定的list（列表）格式提供，按照回调函数（callable function）里使用的顺序提供。

代码不会检查列表中是否提供了合适数量的参数，所以如参数不正确，可能会触发一个难以理解的异常信息。值得注意的是，虽然列表格式期待的参数是全基台值（sill），代码内部使用的是部分基台值（psill）

variogram_function :定制（custom）类型方差模型必须提供此参数，两个参数，一是方差模型列表参数，二是计算方差模型的距离值，前述参数会传递过来。

nlags : int，可选，对于半方差模型的平均箱数（bins），缺省值为6

weight ：bool, 可选，标志半方差的小滞（smaller lags）是否比自动计算方差模型权重更大，常规是权重由逻辑函数计算而来，所以小滞的权重设为1，长滞的权重设为0，逻辑权重的中心设为小滞到大滞距离的70%数值。本参数设为True，意味权重需要考虑，缺省为False。(Kitanidis建议小滞的值在适应方差模型时更重要，所以保留此选项以便于设置权重）

anisotropy_scaling : float, 可选，大面拉伸值，在考虑各向异性时有用，缺省 1（意味没有拉伸），指在旋转数据帧的 y 方向拉伸比例（如果anisotropy_angle !=0时），此参数在coordinate_types值设为'geographic'时无效。

anisotropy_angle : float, 可选，逆时针角度（单位：弧度），在考虑各向异性时，旋转坐标系统的角度值，同上对于coordinate_types值设为'geographic'时无效。

verbose : bool, 可选，设置监测克里金过程文字是否输出，缺省False。

enable_plotting : bool, 可选，结果是否绘图输出，缺省False。

enable_statistics : bool, 可选，缺省False。

coordinates_type : str, 可选，'euclidean' / 'geographic'. 前者平面坐标系，后者球面坐标系，球面坐标系，x为lon，y为lat，lon值为 [0, 360]，lat为 [-90, 90]，缺省：'euclidean'.

 

二、核心方法

execute(self, style, xpoints, ypoints, mask=None, backend='vectorized', n_closest_points = None):

当前确保提供值在结果中属于无误差的，为确切值，未来会允许存在误差

1、输入参数：

style : str，定义如何处理后面的输入点参数，定义为'grid'，则处理xpoints，ypoints为一个矩形网格的x、y坐标值列表，定义为'point‘，则处理为坐标对的列表值。定义为’masked‘，则处理xpoints，ypoints为矩形网格的x、y坐标，且采用掩膜（mask）评估特定的点值

xpoints : 矩阵数组

ypoints : 矩阵数组，如果style 类型为'points'，xpoints，ypoints需有相同的维数

mask : 掩膜矩阵，bool, 维度与xpoints,ypoints相对应，定义哪些网格点数据被排他计算

backend : str, 可选，定义克里计算方法。’vectorized'，向量运算法，速度更快，消耗内存更大，适用于大网格点或大数据量；'loop' 循环计算法，速度慢，消耗内存少。Cpython将使用loop，缺省是'vectorized'.

n_closest_points : int, 可选，对于移窗法计算时，确定计算采用临近点的数量，可加快计算速度（对于大数据量）。但是Kitanidis 提醒，移窗法可能产生不可预期的异常值，如果方差模型没选择好的话

2、返回值：

zvalues : ndarray

sigmasq : ndarray

三、范例

import numpy as np

import pykrige.kriging_tools as kt

from pykrige.ok import OrdinaryKriging



data = np.array([[0.3, 1.2, 0.47],

[1.9, 0.6, 0.56],

[1.1, 3.2, 0.74],

[3.3, 4.4, 1.47],

[4.7, 3.8, 1.74]])



gridx = np.arange(0.0, 5.5, 0.5)

gridy = np.arange(0.0, 5.5, 0.5)



# Create the ordinary kriging object. Required inputs are the X-coordinates of

# the data points, the Y-coordinates of the data points, and the Z-values of the

# data points. If no variogram model is specified, defaults to a linear variogram

# model. If no variogram model parameters are specified, then the code automatically

# calculates the parameters by fitting the variogram model to the binned

# experimental semivariogram. The verbose kwarg controls code talk-back, and

# the enable_plotting kwarg controls the display of the semivariogram.

OK = OrdinaryKriging(data[:, 0], data[:, 1], data[:, 2], variogram_model='linear',

verbose=False, enable_plotting=True)



# Creates the kriged grid and the variance grid. Allows for kriging on a rectangular

# grid of points, on a masked rectangular grid of points, or with arbitrary points.

# (See OrdinaryKriging.__doc__ for more information.)

z, ss = OK.execute('grid', gridx, gridy)



# Writes the kriged grid to an ASCII grid file.

kt.write_asc_grid(gridx, gridy, z, filename="output.asc")