2

我正在尝试使用 R 包中的krige函数在gstatR 中插入一些空间海洋深度数据。我发现超过约 1000 个点，该函数开始花费不合理的时间来完成（即，几小时到几天到从未完成）。这是正常的还是我做错了什么？我特别担心，因为我的最终目标是对一个非常大的数据集（> 30,000 个数据点）进行时空克里金法，我担心考虑到这些运行时间它是不可行的。

我正在运行 gstat-1.1-3 和 R-3.3.2。下面是我正在运行的代码：

library(sp); library(raster); library(gstat)
v.utm # SpatialPointsDataFrame with >30,000 points

# Remove points with identical positons
zd = zerodist(v.utm)
nzd = v.utm[-zd[,1],] # Layer with no identical positions

# Make a raster layer covering point layer
resolution=1e4
e = extent(as.matrix(v.utm@coords))+resolution
r = raster(e,resolution=resolution) 
proj4string(r) = proj4string(v.utm)

# r is a 181x157 raster

# Fit variogram
fv = fit.variogram(variogram(AVGDEPTH~1, nzd),model=vgm(6000,"Exp",1,5e5,1))

# Krige on random sample of 500 points - works fine
size=500
ss=nzd[sample.int(nrow(nzd),size),]
depth.krig = krige(AVGDEPTH~1,ss,as(r,"SpatialPixelsDataFrame"),
               model=depth.fit)

# Krige on random sample of 5000 points - never seems to end
size=5000
ss=nzd[sample.int(nrow(nzd),size),]
depth.krig = krige(AVGDEPTH~1,ss,as(r,"SpatialPixelsDataFrame"),
               model=depth.fit)
4

2 回答 2

1

对于大型数据集，一种比克里金法更快的替代方法是marmap包中的 griddify。我花了一段时间才找到这个，但效果很好。它使用双线性插值，虽然它是为测深地图设计的，但它适用于任何 xyz 数据。

于 2020-07-03T07:02:57.150 回答
1

Choleski 分解（或类似分解）的复杂度为 O(n^3)，这意味着如果将点数乘以 10，则所需时间将增加 1000 倍。解决此问题有两种方法，至少就 gstat 而言：

  1. 安装 BLAS 的优化版本（例如 OpenBLAS 或 MKL）——这不能解决 O(n^3) 问题，但可以在 n 可用内核数的情况下最大程度地加速 n
  2. 通过选择局部邻域（参数 maxdist 和/或 nmax）避免分解完整的协方差矩阵
于 2016-12-10T13:32:25.930 回答