R 'automap' how to create prediction grid to use with AutoKrige (e.g. meuse.grid)?

Question

I'm having a lot of difficulty creating a prediction grid (for the new_data argument) to use with the autoKrige function in the automap package.

I've already tried following the steps in this post (How to subset SpatialGrid using SpatialPolygon) but get the following error : 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

My (limited) understanding is the error relates to there being no non-missing arguments because it is an empty grid. This is fine - all I want is an empty grid constrained by a polygon from a shapefile.

Here is the code I'm working with:

     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)

[see dropbox folder 'Grid.png']

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

[see dropbox folder 'Boundary plot.png']

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

No errors or warnings up to this point. But then...

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

or attempting to pass the object clip_grid through autokrige for the new_data argument:

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

I've had no issues using the non-clipped grid (object = grid).

In a nutshell, I require this [see dropbox folder 'Autokrig plot'] but with the interpolated surfaced constrained (clipped) to the boundary extent of 'Torbay_Box2.shp'

P.S. I attempted to insert images of my plots and links to other posts I've used before asking for help here and a link to my data but as a new user I don't have enough reputation to do this - sorry!

Data and plots can be found on Dropbox.com/sh/yqg20z1ibl3h4aa/AACJnHoEuP-S5fTvAXxsnY1za?dl=0

score 0 · Answer 1 · answered Aug 01 '17 at 15:36

I've now managed to produce an autoKrige [plot] which is masked to the extent of the Torbay_Box2 boundary. However, I never achieved this in the 'conventional' way by creating a prediction grid like meuse.grid. The result is the same so for now I'm happy but I would still like to do it the conventional way eventually.

Here's how I cheated it:

 # 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)

I then performed an IDW interpolation using the empty grid as the newdata, converted the output to raster, clipped this to the Torbay_Box2 boundary and then converted this to a SpatialPixelDataFrame which I passed through as the new_data argument for autoKrige:

# 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)

R 'automap' how to create prediction grid to use with AutoKrige (e.g. meuse.grid)?

1 Answers1