Spatial subsetting SpatialPolygonsDataFrames completely within selecting area using R

Question

I am trying to spatially subset in R:

• the forest stand data (many small polygons) by the
• watershed data (one polygon)

My datasets are both in SpatialPolygonDataFrame format. I don't want to intersect or clip the data, because this creates small polygons on the border. Instead, I would like to select polygons completely within or within, such Select Layer By Location tool in ArcGIS. The output, unlike in ArcGIS, should be a new shapefile.

Here is how it is looks like: forest stands = green, watershed: blue border, selected stands completely within: yellow.

However, I cannot find, how to make this in R?

Using sp::over(x = watershed, y = forest) the output is a a data.frame, not a spatial feature. The tools sf::st_within(x, y, sparse = TRUE, prepared = TRUE) seems promissing, but they operates on sf objects, not on spatialPolygonDataFrame.

Is there any quick spatial subset in R solution that I have overlooked?

Answer 1

**** Update (2019-10-15), I added a spatial.select function to the spatialEco package that emulates a spatial select in ArcGIS/Pro. The predicate argument controls the type of intersection (intersect", "contains", "covers", "touches", "proximity"). If predicate = "proximity" then the distance argument needs to be specified as well. This will return features within a specified distance. This is currently in the development version but will be on CRAN at a future date. You can install the development version using the remotes package: remotes::install_github("jeffreyevans/spatialEco") ****

Take a look at some of the identity functions in rgeos such as gTouches, gIntersects, gContains, gRelate, gWithin, ect...

At their most basic, these functions can return a Boolean that will allow you to perform the equivalent of a spatial select. To control output, take a close look at the function(s) arguments and experiment a bit with outputs from various functions to ensure that your results are as expected. The specific function you choose depends on the application. This Lin.ear th.inking blog gives some guidance on the different rgeos functions and insight to the Dimensionally-Extended 9 Intersection Model (DE-9IM) topology model that GEOS uses.

You do not want a "new spatial feature object" resulting from a given function, per se as this would represent an intersection of some sort. ArcGIS does not return a new feature either, just the subset query of spatially selected features within an existing feature class. You can use results from rgeos or the sp::over function to subset a new feature class based on an index query of the Boolean.

Here is a simple worked example.

Add packages and make example data

library(sp)
library(rgeos)
  p1 = readWKT("POLYGON((0 0,1 0,1 1,0 1,0 0))")
  p2 = readWKT("POLYGON((0.5 1,0 2,1 2,0.5 1))")
  p3 = readWKT("POLYGON((0.5 0.5,0 1.5,1 1.5,0.5 0.5))")
  p2 <- rbind(p2,p3, makeUniqueIDs = TRUE)

plot(p1, xlim=c(0,1), ylim=c(0,2))
  plot(p2,add=TRUE)

Here are the Boolean results from gOverlaps, note the byid = TRUE argument.

gOverlaps(p1, p2, byid=TRUE)

Now, we can wrap which in a row index query using gOverlaps to create the subset object and plot the result.

p2.over <- p2[which(gOverlaps(p1, p2, byid=TRUE)),]

plot(p1, xlim=c(0,1), ylim=c(0,2))
  plot(p2.over, add=TRUE)