Create polygons encompassing specific raster values

Question

I have 4 raster layers contained in a raster stack and want to generate polygons that encompass all cells of a specified value.

The raster package can be use to generate an example data set.

library(raster)
filename <- system.file("external/test.grd", package="raster")
r <- raster(filename)

Like the raster below, my real data are akin to maps of animal habitat and have a patchy distribution of 'good' and 'bad' areas.

To more closely mirror my real data, we can then add a bit of variation to three other rasters and make a stack.

s <- stack(r, r+250, r-250, r+100)

Working with the stack s is it possible to create polygons that surround all cells less than 300 in all stack layers?

As an extension, my end goal is to then calculate the area (or percent) overlap between the resulting polygons.

Any suggestions (specific or general) would be greatly appreciated.

coletl coletl · Accepted Answer · 2017-07-12T10:54:13

Since you're working with a raster stack, all your cells should have the same area. In that case, I don't think you need to use polygons at all. (Note that I adjusted your example data a bit.)

library(raster)
filename <- system.file("external/test.grd", package="raster")
r <- raster(filename)
s <- stack(r, r + 50, r - 50, r + 100)

# Create a new raster stack with results of a logical test
s2 <- s < 300

# Create a raster indicating which cells of the new stack
# have values that are all TRUE
r2 <- sum(s2) == length(unstack(s2))
# Multiply by the area of a single cell
r3 <- r2 * area(r2)[1]

# Sum the area for all raster values
sum(values(r3), na.rm = TRUE)
## 124800

If you'd like to use polygons and your rasters aren't too large, converting the stack to a SpatialPolygonsDataFrame should be pretty quick. Here's an analogous method that yields the same result:

# Create a new raster stack with results of a logical test
s2 <- s < 300

# Convert to sp object
spdf <- as(s2, "SpatialPolygonsDataFrame")

# Index to the rows/features where the values in s2 were all TRUE
spdf2 <- spdf[which(rowSums(spdf@data) == length(unstack(s))), ]

rgeos::gArea(spdf2)
## 124800

Create polygons encompassing specific raster values

1 Answers