CODEWITHSUNDEEP
rpolygonoverlapshapefileintersect
SophiaL
SophiaL

Reputation: 71

Percentage of overlap between SpatialPolygonsDataFrame

I need to calculate the percentage of overlap between several objects of the type SpatialPolygonsDataFrame in R. I have them as shapefiles and uploaded a selection here. 

## load libraries
library(raster)
library(rgdal)
library(rgeos)
library(maptools)
library(sp)
library(FRK)

setwd("...")
Polar1 <- shapefile("Din_biome_Rock_Ice.shp")
Polar2 <- shapefile("FAO_GEZ_Polar.shp")
Polar3 <- shapefile("KG_Beck_EF.shp")

Take a look: 

> Polar1
class       : SpatialPolygonsDataFrame 
features    : 1 
extent      : -179.9999, 180, -89.89197, 82.3525  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
variables   : 1
names       : dummy 
value       :     0

I could convert them to SpatialPolygons:

T_Polar1 <- as(Polar1, "SpatialPolygons")

As we see... 

> T_Polar1
class       : SpatialPolygons
features    : 1 
extent      : -179.9999, 180, -89.89197, 82.3525  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

Either way I could do the calculation of percentage as easy as that:

intersection_Polar1_Polar2 <- gIntersection(Polar1, Polar2)
round(area(intersection_Polar1_Polar2)/area(Polar1)*100, 2)

But I have over 200 shapefiles to do so. Therefore, I need a more convenient method. My ideal solution is posted here. But what I am struggeling with is to adapt the following line of code from the answer to my files: 

poly <- SpatialPolygons(list(Polygons(list(Polygon(p1)), "a"),Polygons(list(Polygon(p2)), "b"),Polygons(list(Polygon(p3)), "c")),1L:3L)

How can I make a list of my files assigning the variables "a", "b", "c" and so on which are needed for the following calculations? All I need to do is to create an object like poly from the answer with my input SpatialPolygonsDataFrame- or SpatialPolygons-files. Does anybody know how to do this? Thanks a lot!

Upvotes: 0

Views: 1848

Answers (1)

Robert Hijmans
Robert Hijmans

Reputation: 47211

Here an adaptation of that answer

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)
crs(poly) <- "+proj=utm +zone=1, +datum=WGS84"

I understand you to have separate SpatialPolygonDataFrame objects that have a single (multi-) polygon each. You could store these in a list. With the example data:

ss <- list(poly[1,], poly[2,], poly[3,])

And then do something like this:

n <- length(ss)
overlap <- matrix(0, nrow=n, ncol=n)
diag(overlap) <- 1

for (i in 1:n) {
    ss[[i]]$area <- area(ss[[i]])
}


for (i in 1:(n-1)) {
    for (j in (i+1):n) {
        x <- intersect(ss[[i]], ss[[j]])
        if (!is.null(i)) {
            a <- area(x)
            overlap[i,j] <- a / ss[[i]]$area
            overlap[j,i] <- a / ss[[j]]$area
        } 
    }
}

With your data, you could make a list of polygons like this: 

ff <- c("Din_biome_Rock_Ice.shp", "FAO_GEZ_Polar.shp", "KG_Beck_EF.shp")
ss <- lapply(ff, shapefile)

And take it from there.

Upvotes: 2

Related Questions