I am busy developing a module that implements autoKrige (automap library) on soil sample data with a 20x20m cell-size. Once the kriging is complete I would like to crop/mask the kriging result with the field's boundary.
The problem with cropping/masking is that (because of the 20x20 cells) the border of the field gives a "step" effect. What I am looking for is a smooth boundary (cutting through the cells).
Below is the code to generate the two scenarios explained above:
library(sp)
library(rgeos)
#create polygon
r1 <- cbind(c(641777, 642290, 642276, 641794), c(7036885, 7036743, 7036154, 7036146))
r2 <- cbind(c(642320, 642478, 642494, 642314), c(7036732, 7036699, 7036296, 7036290))
sr1 <- Polygons(list(Polygon(r1)),"r1")
sr2 <- Polygons(list(Polygon(r2)),"r2")
boundary.sp <- SpatialPolygons(list(sr1,sr2))
boundary.sp@proj4string <- CRS('+proj=utm +zone=35 +south +datum=WGS84 +units=m +no_defs')
#create bounding box grid
bbox <- bbox(boundary.sp)
boundary.grid <- expand.grid(x = seq(from = bbox[1], to = bbox[3], by = 20), y = seq(from = bbox[2], to = bbox[4], by = 20))
coordinates(boundary.grid) <- ~x + y
gridded(boundary.grid) <- TRUE
boundary.grid@proj4string <- boundary.sp@proj4string
#create SpatialPixels grid
boundary.grid.stepped <- boundary.grid[!is.na(over(boundary.grid, boundary.sp)),]
plot(boundary.grid.stepped)
#cut grid with polygon to create SpatialPolygons grid
boundary.poly.grid <- as.SpatialPolygons.GridTopology(getGridTopology(boundary.grid), proj4string = CRS('+proj=utm +zone=35 +south +datum=WGS84 +units=m +no_defs'))
boundary.grid.smooth <- gIntersection(boundary.poly.grid, boundary.sp, byid=TRUE)
plot(boundary.grid.smooth)
Currently the above grid (boundary.grid.stepped) is passed as the 'new_data' parameter when calling autoKrige.
Which is a better approach and how do I go about implementing it:
1) Prepare the target grid beforehand and use it as 'new_data' or,
2) Do the kriging onto a bounding box grid and cut/crop/mask afterwards?
I'd say do the gridding directly on the target grid, via feeding it to new_data. The kriging result depends only on the input data, so I suspect there is little difference between cropping/calculating and calculating/cropping.