Reputation: 41
I'm new to R and it's my first time posting here, thought to let you guys know I'm still a newbie in the R language with the hope it will simplify your answers. I'm currently doing my masters dissertation which includes a lot of work with shapefiles containing distributions of different species in South Africa. My research includes 4 families of animals; 1 prey family (consisting of 29 species) and 3 predatory families (each containing even more species). For each species I have a shapefile with a polygon of their distribution, so I have a LOT of shapefiles. My problem now is, I have to calculate the percentage of overlap each predator has with each prey species (100% overlap is when the distribution/polygon of the predator is within the distribution/polygon of the prey, so percentage should be calculated relative to the polygon of the prey species).
I know I could probably do this one by one for each species separately, but this would literally take me weeks and I don't have time for that. Even my promotor at uni hasn't done anything like this before and is trying to figure it out himself. He originally gave me this code:
my_rangemaps <- list.files(path = "imagine_rangemaps", pattern = ".shp", full.names = TRUE)
my_rangemaps
rangemap_matrix <- pairwiseRangemaps(my_rangemaps, projection = 3035,
Ncpu = 1, nchunks = 1, filename = "rangemap_matrix.csv")
But the result was not what we expected as it doesn't calculate percentages relative to whatever polygon. Does anyone know a way to get percentages of overlap using a relative easy code without too much fuzz? Maybe resulting in the form of a matrix containing all species/shapefiles/polygons?
Thank you all in advance!
Upvotes: 4
Views: 2942
Reputation: 54287
For the fun of it, here's an example. Some folks here or on gis.stackexchange.com might have better ways in petto:
library(raster)
library(sp)
## example data:
p1 <- structure(c(0, 0, 0.4, 0.4, 0, 0.6, 0.6, 0), .Dim = c(4L, 2L), .Dimnames = list(NULL, c("x", "y")))
p2 <- structure(c(0.2, 0.2, 0.6, 0.6, 0, 0.4, 0.4, 0), .Dim = c(4L, 2L), .Dimnames = list(NULL, c("x", "y")))
p3 <- structure(c(0, 0, 0.8, 0.8, 0, 0.8, 0.8, 0), .Dim = c(4L, 2L), .Dimnames = list(NULL, c("x", "y")))
poly <- SpatialPolygons(list(Polygons(list(Polygon(p1)), "a"),Polygons(list(Polygon(p2)), "b"),Polygons(list(Polygon(p3)), "c")),1L:3L)
plot(poly)
## areas for the original shapes:
areas_poly <- vector(length = length(poly))
for (x in seq_along(poly)) areas_poly[x]<-area(poly[x])
## areas for the overlapping regions:
idx <- combn(seq_along(poly),2)
areas_intersect <- sapply(1:ncol(idx), function(x) {
area(intersect(poly[idx[1,x]], poly[idx[2,x]]))
})
## get overlaps in percentages:
overlap_perc <-
round(do.call(cbind, lapply(seq_len(ncol(idx)), function(x)
rbind(
areas_intersect[x] / areas_poly[idx[1,x]] * 100,
areas_intersect[x] / areas_poly[idx[2,x]] * 100
)
)), 2)
## into matrix form:
m <- matrix(100, ncol=length(poly), nrow=length(poly))
m[rbind(t(idx),t(idx)[,2:1])] <- as.vector(t(overlap_perc))
m
# [,1] [,2] [,3]
# [1,] 100.0 33.33 100
# [2,] 50.0 100.00 100
# [3,] 37.5 25.00 100
Upvotes: 6