R'automap'如何创建与 AutoKrige 一起使用的预测网格（例如 meuse.grid）？答案

【问题标题】：R 'automap' how to create prediction grid to use with AutoKrige (e.g. meuse.grid)?R'automap'如何创建与 AutoKrige 一起使用的预测网格（例如 meuse.grid）？
【发布时间】：2018-01-07 09:00:39
【问题描述】：

我在创建预测网格（用于 new_data 参数）以与 automap 包中的 autoKrige 函数一起使用时遇到很多困难。

我已经尝试按照本文 (How to subset SpatialGrid using SpatialPolygon) 中的步骤进行操作，但出现以下错误： x@coords[i, , drop = FALSE] 中的错误： (下标) 逻辑下标太长另外：警告信息： 1：在 min(x) 中：min 没有非缺失参数；返回 Inf 2：在 max(x) 中：max 没有非缺失参数；返回 -Inf

我的（有限）理解是错误与没有非缺失参数有关，因为它是一个空网格。这很好 - 我想要的只是一个由 shapefile 中的多边形约束的空网格。

这是我正在使用的代码：

     shp <-  shapefile("C://path/path/Tobay_Box2.shp")
         shp <-  spTransform (shp,"+proj=utm +ellps=WGS84 +datum=WGS84")
         grid <-        GridTopology(cellcentre.offset=c(731888.0,7457552.0),cellsize=c(2,2),cells.dim=c(122,106))
         grid <- SpatialPixelsDataFrame(grid,
                              data=data.frame(id=1:prod(122,106)),
                              proj4string=CRS("+proj=utm +ellps=WGS84 +    datum=WGS84"))
plot(grid)

[查看保管箱文件夹“Grid.png”]

bound <- shp@polygons
bound <- SpatialPolygons(bound, proj4string=CRS("+proj=utm +ellps=WGS84 +datum=WGS84"))
plot(bound)

[查看保管箱文件夹'Boundary plot.png']

clip_grid <- grid[!is.na(over(grid, bound)),]

到目前为止没有错误或警告。但是后来……

plot(clip_grid)

Error in x@coords[i, , drop = FALSE] : 
  (subscript) logical subscript too long
In addition: Warning messages:
1: In min(x) : no non-missing arguments to min; returning Inf
2: In max(x) : no non-missing arguments to max; returning -Inf

或尝试通过 autokrige 为 new_data 参数传递对象 clip_grid：

PerInkrg <- autoKrige (PerArIn~1, hs1, clip_grid)

Error in predict.gstat(g, newdata = newdata, block = block, nsim = nsim,  : 
  value not allowed for: %s %s newdata empty or only NA's

我在使用非剪裁网格（对象 = 网格）时没有遇到任何问题。

简而言之，我需要这个 [参见 dropbox 文件夹 'Autokrig plot']，但内插曲面约束（剪裁）到 'Torbay_Box2.shp' 的边界范围

附：我试图在此处寻求帮助之前插入我的情节图像和指向我使用过的其他帖子的链接以及指向我的数据的链接，但作为一个新用户，我没有足够的声誉来执行此操作 - 抱歉！

数据和图表可在 Dropbox.com/sh/yqg20z1ibl3h4aa/AACJnHoEuP-S5fTvAXxsnY1za?dl=0 上找到

【问题讨论】：

标签： r gis kriging gstat automap

【解决方案1】：

我现在设法生成了一个 autoKrige [绘图]，它被遮盖到 Torbay_Box2 边界的范围内。但是，我从来没有通过创建像 meuse.grid 这样的预测网格以“传统”方式实现这一点。结果是一样的，所以现在我很高兴，但我最终还是想以传统的方式来做。

我是这样作弊的：

 # Load sample box extent

bx.data <- readOGR (".", "Tobay_Box2")
bx <- spTransform(bx.data,"+proj=utm +ellps=WGS84 +datum=WGS84") #transformsto UTM projection
str(bx)

# Set the boundary extent with that of sample box extent

hs1@bbox <- bx@bbox
#create an empty grid
grd              <- as.data.frame(spsample(hs1, "regular", n=50000))                                                                                                
names(grd)       <- c("X", "Y")
coordinates(grd) <- c("X", "Y")
gridded(grd)     <- TRUE  # Create SpatialPixel object
fullgrid(grd)    <- TRUE  # Create SpatialGrid object

plot(hs1)
plot(grd, pch = ".", add = T)
proj4string(grd) <- proj4string(hs1)

然后，我使用空网格作为新数据执行 IDW 插值，将输出转换为栅格，将其裁剪到 Torbay_Box2 边界，然后将其转换为 SpatialPixelDataFrame，我将其作为 autoKrige 的 new_data 参数传递：

# For PerArIn (% area inhabited) 
#interpolate the grid cells using all points and a power value of 2 

hs1.idw <- gstat::idw(PerArIn ~ 1, hs1, newdata=grd, idp=2.0)



# Convert to raster object then clip to Hollicombe sample box

r       <- raster(hs1.idw)
r.m     <- mask(r, bx)


#Convert and set as prediction grid for Kriging

grd<- rasterToPoints(r.m, spatial=TRUE)
gridded(grd) <- TRUE
grd <- as (grd, "SpatialPixels")

#en voila!
PerInkrg <- autoKrige (PerArIn~1, hs1,grd)

【讨论】：