如何使用克里格。gstools的扩展漂移以插值2D地理空间数据

Question

如何使用gstools的外部漂移克里格法插值2D地理空间数据。Gstools只有一维外部驱动Kriging插值示例，但我想使用gstool的外部驱动Kriging插值来插值2-D地理空间数据。

这是一个官方的例子：

import numpy as np
from gstools import SRF, Gaussian, krige

# synthetic condtions with a drift
drift_model = Gaussian(dim=1, len_scale=4)
drift = SRF(drift_model, seed=1010)
cond_pos = [0.3, 1.9, 1.1, 3.3, 4.7]
ext_drift = drift(cond_pos)
cond_val = ext_drift * 2 + 1
# resulting grid
gridx = np.linspace(0.0, 15.0, 151)
grid_drift = drift(gridx)
# kriging
model = Gaussian(dim=1, var=2, len_scale=4)
krig = krige.ExtDrift(model, cond_pos, x, ext_drift)
krig(gridx, ext_drift=grid_drift)
ax = krig.plot()
ax.scatter(cond_pos, cond_val, color="k", zorder=10, label="Conditions")
ax.plot(gridx, grid_drift, label="drift")
ax.legend()

Answer 1

我是GSTools的开发人员，当然也可以在2D和3D中执行外部漂移克里格法。

在下面的示例中，我创建了一个在x方向上的人工正弦外漂移。

import numpy as np
import gstools as gs

# conditioning data (x, y, values)
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],
    ]
)

# convert conditioning data
x, y, val = data[:, 0], data[:, 1], data[:, 2]

# grid definition for output field
gridx = np.arange(0.0, 5.5, 0.1)
gridy = np.arange(0.0, 6.5, 0.1)

# external drift at conditioning points
# (given as a sinusodial drift in x direciton)
ext_drift_cond = np.sin(x)

# external drift at the output grid
ext_drift_grid = np.repeat(np.sin(gridx), len(gridy))

# a gaussian covariance model
cov_model = gs.Gaussian(
    dim=2, len_scale=1, anis=0.5, angles=-0.5, var=0.5, nugget=0.1
)

# perform the kriging and plot results
EDK = gs.krige.ExtDrift(
    model=cov_model, 
    cond_pos=(x, y), 
    cond_val=val,
    ext_drift=ext_drift_cond,
)
EDK.structured([gridx, gridy], ext_drift=ext_drift_grid)
EDK.plot()

​

​