How to mask or crop IDW results My question is: how can I mask or crop the results, using R, of IDW interpolation to only the area containing the original set of data points?

In the example below, 20 random points are used to interpolate a surface using the IDW function in gstat. A convex hull is obtained from the points set and plotted on the interpolated map. But, I would like to crop the map so only areas within the point cloud show the interpolation result.

Thanks for any suggestions!

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