0

我正在尝试使用gstat()and创建反距离加权栅格raster::interpolate()。我在将列名传递给 gstat 函数的公式参数时遇到问题。

对列名进行硬编码就可以了：

gs <- gstat(formula=v1~1, locations = data)
r <- raster(shape, res=1000, crs = crs(data))
idw <- raster::interpolate(r, gs)
idwr <- mask(r, gs)
plot(idwr)

在此处输入图像描述

将其包装到一个函数中，以便我可以遍历多个列名会引发错误：

apply_gstat <- function(col_name, data = data, shape = shape) {
    gs <- gstat(formula=col_name~1, locations = data)
    r <- raster(shape, res=1000, crs = crs(data))
    ## interpolate() throws an error because of issue with gstat
    idw <- raster::interpolate(r, gs)
    idwr <- mask(r, gs)
    plot(idwr)
}

col_names <- c("v1", "v2", "v3")
lapply(col_names, function(x) {
    gstat_apply(col_name = x, data = data, shape = shape)
}


Error in predict.gstat(model, blockvals, debug.level = debug.level, ...) : 
  too many spatial dimensions: 18
In addition: Warning message:
In predict.gstat(model, blockvals, debug.level = debug.level, ...) :
  NAs introduced by coercion
4

1 回答 1

1

如果您使用它可能会起作用as.formula(paste(col_name, '~1')) 这是一个独立的示例：

library(raster)
library(gstat)
data(meuse)
r <- raster(system.file("external/test.grd", package="raster"))

gfun <- function(v) {
    mg <- gstat(formula = as.formula(paste(v, "~1")), locations= ~x+y, data=meuse, nmax=7, set=list(idp = .5))
    z <- interpolate(r, mg)
    mask(z, r)
}

vars <- c('cadmium', 'copper', 'lead', 'zinc')
x <- sapply(vars, gfun)
s <- stack(x)
s
#class       : RasterStack 
#dimensions  : 115, 80, 9200, 4  (nrow, ncol, ncell, nlayers)
#resolution  : 40, 40  (x, y)
#extent      : 178400, 181600, 329400, 334000  (xmin, xmax, ymin, ymax)
#coord. ref. : +init=epsg:28992 +towgs84=565.237,50.0087,465.658,-0.406857,0.350733,-1.87035,4.0812 +proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs 
#names       :     cadmium,      copper,        lead,        zinc 
#min values  :   0.3131223,  17.3430502,  42.3218862, 123.9369493 
#max values  :    11.91185,    97.28761,   470.70461,  1387.50615
于 2018-10-06T08:25:05.347 回答