Detect if raster is within, without, or intersecting a SpatialPolygons object

Question

I have many rasters for which I'd like to check if they are contain fully within a spatial polygon, fully without the spatial polygon, or intersecting with the spatial polygon (this could mean the polygon is fully within the raster, or the polygon and raster overlap). I'm doing this check to save me from time intensive masking when possible.

Here is an example:

                                      # create 3 example rasters
r <- raster()
r[] <- rnorm(n = ncell(r))
e1 <- extent(c(45,55,45,50))
r1 <- crop(r,e1)
e2 <- extent(c(20,25,25,30))
r2 <- crop(r,e2)
e3 <- extent(c(38,55,57,65))
r3 <- crop(r,e3)


                                        #create SpatialPolygons

x <- c(40,60)
y <- c(40,60)
m <- expand.grid(x,y)
m <- m[c(1,2,4,3),]
p1 <- Polygon(m)
p1 <- Polygons(list(p1),1)

x <- c(10,15)
y <- c(10,15)
m <- expand.grid(x,y)
m <- m[c(1,2,4,3),]
p2 <- Polygon(m)
p2 <- Polygons(list(p2),2)

x <- c(30,45)
y <- c(70,80)
m <- expand.grid(x,y)
m <- m[c(1,2,4,3),]
p3 <- Polygon(m)
p3 <- Polygons(list(p3),3)


poly <- SpatialPolygons(list(p1,p2,p3))

plotting these:

I will read in each raster separately and check if it is within, without, or intersecting the SpatialPolygons.

What do you think will be the most efficient way to do this in R? I have thousands of 4mb rasters that I'm planning to mask in parallel and would like this check to speed things up a bit.

Note, there's also this question: https://gis.stackexchange.com/questions/34535/detect-whether-there-is-a-spatial-polygon-in-a-spatial-extent

However, I don't think it gives the detail I'm looking for. For example, all of the rasters are within the extent of the spatial polygons, but not all are within the spatial polygons.

Functions like those in rgeos (gIntersects, gContains), would probably be handy. I'm not sure if these are most efficient, or how I should convert raster (or it's extent), to an sp object.

thanks!

Answer 1

You can also use gRelate for this. It returns a DE-9IM code that describes the relationships between interior, boundary and exterior components of the two geometries.

library(rgeos)
x <- sapply(rlist, function(x) 
  gsub('[^F]', 'T', gRelate(as(extent(x), 'SpatialPolygons'), poly)))

You could then compare the strings to relationships of interest. For example we might define within, disjoint, and overlaps as follows (but note that other some intersections are optional for given relationships - "within" is defined by GEOS as T*F**F***, "disjoint" as FF*FF****, and "overlaps" as T*T***T**):

pat <- c(TFFTFFTTT='within', FFTFFTTTT='disjoint', TTTTTTTTT='overlaps')
pat[x]

##  TFFTFFTTT  FFTFFTTTT  TTTTTTTTT 
##   "within" "disjoint"  "overlaps" 

It seems marginally faster than the gContainsProperly/gIntersects approach, but @Tedward's post is more comprehensible, and more consistent with GEOS definitions (though the power to create specific relationship definitions might be desirable).

---

The elements of the DE-9IM strings represent, in order:

1. Does the interior of geometry A intersect the interior of geometry B?
2. Does the boundary of A intersect the interior of B?
3. Does the exterior of A intersect the interior of B?
4. Does the interior of geometry A intersect the boundary of geometry B?
5. Does the boundary of A intersect the boundary of B?
6. Does the exterior of A intersect the boundary of B?
7. Does the interior of geometry A intersect the exterior of geometry B?
8. Does the boundary of A intersect the exterior of B?
9. Does the exterior of A intersect the exterior of B?