How to mask or crop IDW results

Here I am answering my own question:

To obtain an IDW surface interpolation restricted to the area containing the original set of xyz data points:

• Create a rectangular grid based on the xy point ranges using expand.grid {base}
• Create a convex hull using the function convexhull.xy{spatstat}
• Cut an irregular grid using the function inout{splancs} on the rectangular grid masked by the convex hull
• Perform the interpolation using the function idw{gstat} with the new grid.
• Be sure to detach spatstat, because it has a different function for IDW with the same name

Code:

###
# IDW in R
#
# Load gstat:
library(gstat)
# 
# An example XYZ data set:
set.seed(1234)
x <- rnorm(20)*1000+3000
y <- rnorm(20)*1000+3000
z <- rnorm(20)+10
xyz <- data.frame(x,y,z)
# 
# Calculate ranges and spacing:
xr <- round(range(x),1);xr
yr <- round(range(y),1);yr
b <- round((xr[2]-xr[1])/60,2);b
#
# Create a grid:
grd1 <- expand.grid(x=seq(from=xr[1],to=xr[2],by=b),
    y=seq(from=yr[1],to=yr[2],by=b))
plot(grd1)
#
# Load spatstat and splancs
library(spatstat)
library(splancs)
#
# Extract convex hull 
w2 <- convexhull.xy(x,y)
#
# detach spatstat because it has a different "idw" function:
detach("package:spatstat", unload=TRUE)
#
# Create a new cropped grid:
grd2 <- grd1[inout(grd1,w2$bdry[[1]]),]
plot(grd2)
#
# Convert grids and dataset to sp objects:
coordinates(grd1) <- ~x+y
gridded(grd1)   <- TRUE
coordinates(grd2) <- ~x+y
gridded(grd2)   <- TRUE
coordinates(xyz) <- ~x+y
#
# Run the interpolations:
idw1<-idw(formula=z~1, locations=xyz, newdata=grd1)
idw2<-idw(formula=z~1, locations=xyz, newdata=grd2)
#
# Plot the results:
image(idw1, col=terrain.colors(16))
contour(idw1, add=T)
points(x,y, col=4, pch=19)
image(idw2, col=terrain.colors(16))
contour(idw2, add=T)
points(x,y, col=4, pch=19)
#

Here is what I would do:

idw_output <- as.data.frame(idw)[, 1:3]

# Adding the IDW raster
# The first step to incorporate the IDW results into the interactive map is to turn it into a spatial object
idw_raster <- rasterFromXYZ(idw_output)

# We only want to show the interpolation results within the area of which we have data points
# To do this, we will clip the raster to data region by creating convex hull, buffer, and mask (clip) raster
# Perform convex hull and buffer operations on projected UTM data
xyz_hull <- gConvexHull(xyz)
xyz_buff <- gBuffer(xyz_hull, width = 25000) # arbitrary 1km buffer

# Mask IDW raster to confine to regions that we have information for
idw_raster_crop <- mask(idw_raster, xyz)

spplot(idw_raster_crop)

# Create contours of the IDW
idw_contour <- rasterToContour(idw_raster_crop, nlevels = 15)

When you use your real data, you can use mapview, ie mapview(idw_raster_crop) to see the interactive map.

Cite: https://rstudio-pubs-static.s3.amazonaws.com/212793_f130ecc723da4ed98680d6d3d5c4aff9.html