0

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

我已经尝试按照本文中的步骤进行操作（如何使用 SpatialPolygon 对 SpatialGrid 进行子集化），但出现以下错误：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' 的边界范围

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

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

4

1 回答 1

0

我现在设法生成了一个 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)

然后，我使用空网格作为 newdata 执行 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)
于 2017-08-01T15:36:05.830 回答