0
$\begingroup$

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 <- rnorm(20)+10
y <- rnorm(20)+10
z <- rnorm(20)+10
#
xyz <- data.frame(x,y,z)
coordinates(xyz) <- ~x+y
xr <- range(x)
yr <- range(y)
b <- round((xr[2]-xr[1])/100,2)
#
grd <- expand.grid(x=seq(from=xr[1],to=xr[2],by=b),
y=seq(from=yr[1],to=yr[2],by=b))
#
coordinates(grd) <- ~ x+y
gridded(grd) <- TRUE
#
idw<-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<-convexhull.xy(x,y)
plot(conv,add=TRUE)
#
$\endgroup$
  • $\begingroup$ I encourage you to use Markdown to format your post in a way that is nicer to read. $\endgroup$ – D.W. Apr 5 '17 at 1:05
0
$\begingroup$

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

| improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ While this link may answer the question, it is better to include the essential parts of the answer here and provide the link for reference. Link-only answers can become invalid if the linked page changes. - From Review $\endgroup$ – oW_ Nov 10 '17 at 0:27
  • $\begingroup$ That makes sense. I just edited the answer. $\endgroup$ – Meo Nov 10 '17 at 16:18
0
$\begingroup$

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)
#
| improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.