Create polygons based on raster values

Question

I have 4 raster layers contained in a raster stack within program R. The layers represent occurrence surfaces for two species in two seasons (i.e. species one summer and winter and species 2 in summer and winter). Values for all raster layers range from 0 to 10 and represent a relative ranking of habitat, worst to best, respectively. For each layer I want to make polygon/s that encompass raster values greater than a specified number, for example 8.

Within a given layer there is a patchy distribution of 'good' habitat (i.e. > 8) which will results in multiple polygons for each layer.

This process is similar to creating animal home ranges from kernel densities and I was surprised to not see previous posts on the topic, although I may be searching the wrong terms...

My hope is to then calculate the overlap of seasonal polygons among species (i.e. overlap of species 1 and 2 in summer), likley using intersect from the raster package, although I welcome any suggestions.

Answer 1

Why vectorize the data when you can get at this, much more efficiently, in a raster data structure. It is always easier, and less computationally expensive, to overlay cells rather than vector topology. This also removes the spatial uncertainty of converting an n-resolution raster to a vector and then intersecting the results. You could always convert the final raster overlay results to polygons.

First, lets create some data emulating what you have described. I am using "raster.gaussian.smooth" to create a correlated "patchy" landscape and then stretching back to a 0-10 scale thus, the extra steps.

library(raster)
library(spatialEco)

s = 3   # Gaussian sigma parameter for raster.gaussian.smooth function
w = 21  # Matrix (window) size for raster.gaussian.smooth function

s1.summer <- raster(nrows=500, ncols=500)
  s1.summer[] <- round(runif(ncell(s1.summer), 0,10),0) 
  s1.summer <- raster.gaussian.smooth(s1.summer, sigma = s, n = w)
  s1.summer <- raster.transformation(s1.summer, "stretch", 0, 10)
s2.summer <- raster(nrows=500, ncols=500)
  s2.summer[] <- round(runif(ncell(s2.summer), 0,10),0) 
  s2.summer <- raster.gaussian.smooth(s2.summer, sigma = s, n = w)
  s2.summer <- raster.transformation(s2.summer, "stretch", 0, 10)

s1.winter <- raster(nrows=500, ncols=500)
  s1.winter[] <- round(runif(ncell(s1.winter), 0,10),0) 
  s1.winter <- raster.gaussian.smooth(s1.winter, sigma = s, n = w)
  s1.winter <- raster.transformation(s1.winter, "stretch", 0, 10)
s2.winter <- raster(nrows=500, ncols=500)
  s2.winter[] <- round(runif(ncell(s2.winter), 0,10),0) 
  s2.winter <- raster.gaussian.smooth(s2.winter, sigma = s, n = w)
  s2.winter <- raster.transformation(s2.winter, "stretch", 0, 10)

par(mfrow=c(2,2))
  plot(s1.summer)
  plot(s2.summer)
  plot(s1.winter)
  plot(s2.winter)

We can now reclassify each raster to [0,1] based on the occupancy threshold and then sum the results thus, illustrating seasonal overlap.

occ <- 5 # occupancy threshold 

s1.summer <- setValues(s1.summer, ifelse(getValues(s1.summer) >= occ, 1, 0) )
s2.summer <- setValues(s2.summer, ifelse(getValues(s2.summer) >= occ, 1, 0) )
s1.winter <- setValues(s1.winter, ifelse(getValues(s1.winter) >= occ, 1, 0) )
s2.winter <- setValues(s2.winter, ifelse(getValues(s2.winter) >= occ, 1, 0) )

summer.overlap <- s1.summer + s2.summer
winter.overlap <- s1.winter + s2.winter

par(mfrow=c(2,1))
 plot(summer.overlap, main="summer species overlap")
 plot(winter.overlap, main="winter species overlap")

The resulting summer and winter overlap rasters could then be vectorized based on the summed values or, by reclassifying to [0,1]. However, even in treating the overlap as a binomial process, I would not expect a single polygon.