[GIS] Measuring percentage overlap between home ranges and another shape file in R

biologyoverlapping-featuresrsfshapefile

I have tracking data from multiple animals and I have created kernel density estimates of their home ranges. I want to find out the proportion of each home range area that falls within protected areas. I'm using the adehabitatHR package to calculate the home ranges and the wdpar package to download the protected areas from the World Database on Protected Areas. But I'm not sure how to link the two to calculate the proportion of overlap. Here's an example of the code I have so far:

library(amt)
library(wdpar)
library(rgdal)
library(adehabitatHR)

data(deer)
mini_deer <- deer[1:100, ]
mini_deer$id <- ifelse(mini_deer$x_ < 4313000, "A", "B")

mini_deer <-
  mk_track(
    mini_deer,
    .x = x_,
    .y = y_,
    .t = t_,
    id = id,
    crs = CRS(
      "+init=epsg:3035 +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
    )
)

# Keep only the id column and the x and y coordinates
mydata.sp <- mini_deer[, c("id", "x_", "y_")] 
coordinates(mydata.sp) <- c("x_", "y_")

# kde
kernel.ref <- kernelUD(mydata.sp, h = "href")  # href = the reference bandwidth
mydata.kernel.poly <- getverticeshr(kernel.ref, percent = 95, unout = "m2") 
print(mydata.kernel.poly)  # returns the area of each polygon

# get the protected area data
ger_raw_data <- wdpa_fetch("Germany", wait = TRUE)
pa_data <- st_transform(ger_raw_data, 3035)

# scene missing
# proportion of overlap for A and B

Best Answer

In this case it looks like there's zero overlap between your habitats (red hashed) and the protected areas (orange solid). This map produced in QGIS:

enter image description here

but in general once you've got both data sets as sf class objects and as polygons (there's some points in the protected area data) then you do st_intersection of the habitat with the protected area to get the overlapping area. For your data:

> kernel.sf = st_as_sf(mydata.kernel.poly)
> st_intersection(kernel.sf, pa_data)
Simple feature collection with 0 features and 28 fields

it returns 0 features because there's zero intersection. If there was an intersection you could then use st_area to get the area of intersection to work out the percentage of the habitat polygon in the protected area.

Related Question