I want to calculate the percentage area of habitat suitability of a species that overlaps with protected area polygons. I do not know the R language very well, but here is what I have so far.
These are the attributes of the area of habitat suitability derived from a maxent prediction:
class : RasterLayer
dimensions : 6480, 8520, 55209600 (nrow, ncol, ncell)
resolution : 0.008333333, 0.008333333 (x, y)
extent : -103, -32, -36, 18 (xmin, xmax, ymin, ymax)
crs : +proj=longlat +ellps=WGS84
of the protected areas:
Simple feature collection with 5667 features and 2 fields (with 8 geometries empty)
geometry type: GEOMETRY
dimension: XY
bbox: xmin: -118.6344 ymin: -59.85538 xmax: -25.29094 ymax: 32.48333
CRS: +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
Does someone know a way to calculate the percentage area of habitat suitability that overlaps with protected area polygons?
Sorry, I really do not know so much about how to work with these data. I hope I gave all the relevant information.
I would appreciate any input.
To answer your first question, you should be able to use zonal statistics to calculate the area of potential habitat found in protected areas using the spatialEco package:
zonal.stats(x, y, stats = c("min", "mean", "max"))
#x = Polygon object of class SpatialPolygonsDataFrame
#y = rasterLayer object of class raster
https://www.rdocumentation.org/packages/spatialEco/versions/1.3-0/topics/zonal.stats
Here is a reproducible example from the spatialEco package that first calculates the percentage of pixels in each polygon >= a threshold value and second calculates the sum of pixels in each polygon >= the threshold value used to reclassify the input raster. You might be interested in both avenues for your work.
library(spatialEco)
library(raster)
library(sp)
# here the fxn will calculate the percentage of cells >= 0.5
# percent x >= p function
pct <- function(x, p=0.50, na.rm = FALSE) {
if ( length(x[x >= p]) < 1 ) return(0)
if ( length(x[x >= p]) == length(x) ) return(1)
else return( length(x[x >= p]) / length(x) )
}
# create some example data
p <- raster(nrow=10, ncol=10)
p[] <- runif(ncell(p)) * 10
p <- rasterToPolygons(p, fun=function(x){x > 9})
r <- raster(nrow=100, ncol=100)
r[] <- runif(ncell(r))
plot(r)
plot(p, add=TRUE, lwd=4)
# run zonal statistics using pct functions
z.pct <- zonal.stats(x=p, y=r, stats = "pct")
z.pct
#Alternatively, reclassify the raster based on a threshold
r.c<-reclassify(r, c(-Inf, 0.5, 0, 0.5, Inf, 1)) #all values >0.5 reclassified to 1
plot(r.c)
plot(p, add=TRUE, lwd=4) #add poly to the plot
# run zonal stats and calculate sum of cells in each poly
z.sum <- zonal.stats(x=p, y=r.c, stats = "sum")
z.sum