Percentage of Overlap Between Polygons

How to compute all pairwise interaction between polygons and the the percentage coverage in R with sf?

Without a real example to benchmark, I'm not sure it's faster than your solution. But it's simpler and easier to understand (at least for my brain).

sf::st_intersection() is vectorized. So it will find & return all the intersections of the first & second argument for you. In this case, the two arguments are the same set of polygons.

sf::st_intersection(s.f.2, s.f.2) %>% 
dplyr::mutate(
area = sf::st_area(.),
proportion = area / area.1
) %>%
tibble::as_tibble() %>%
dplyr::select(
id_1 = id,
id_2 = id.1,
proportion,
) %>%
# tidyr::complete(id_1, id_2, fill = list(proportion = 0))
tidyr::pivot_wider(
names_from = id_1,
values_from = proportion,
values_fill = 0
)

Output:

# A tibble: 3 x 4
id_2 `1` `2` `3`
<dbl> <dbl> <dbl> <dbl>
1 1 1 0.316 0
2 2 0.270 1 0
3 3 0 0 1

Things to consider:

  • Keep the areas as proportions, instead of percentages. It's usually better for calculations later.
  • Stay long, and don't pivot. It's usually better for calculations later, because you can join on id_1 and id_2.
  • If you pivot wide, you might want it as a Matrix, instead of a data.frame.

Get percentage of overlap between two multipolygons in R

Let's create a function to do it by year.

library(sf)

## storing years in a variable
years <- unique(district.df$year)

## creating an "auxiliary" geometry
flood.df$geom_2 <- st_geometry(flood.df)

## function to calculate this "proportion of intersection"
## inputs: y = year; dt1 = dataset1 (district); dt2 = dataset2(flood)
prop_intersect <- function(y, dt1, dt2) {
## filtering year
out1 <- dt1[dt1$year == y,]
out2 <- dt2[dt2$year == y,]

## calculating the areas for the first dataset
out1 <- transform(out1, dist_area = as.numeric(st_area(geometry)))

## joining the datasets (polygons that intersect will be joined).
## It would be nice to have id variables for both datasets
output <- st_join(
x = out1,
y = out2,
join = st_intersects
)

## calculating area of intersection
output <- transform(output,
inter_area = mapply(function(x, y) {
as.numeric(sf::st_area(
sf::st_intersection(x, y)
))}, x = geometry, y = geom_2))

## calculating proportion of intersected area
output <- transform(output, prop_inter = inter_area/dist_area)

return(output)
}

So, now we have a function to do what you request. It is very likely that there exists a better (more efficient and "clean") way to do it. However, this is what I could come up with. Also, since this is not a reprex, it is hard for me to know whether this code works or not.

That being said, now we can iterate on "year" as follows

final_df <- lapply(years, prop_intersect, dt1 = district.df, dt2 = flood.df)

final_df <- do.call("rbind", final_df)

Percentage of overlap between SpatialPolygonsDataFrame

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.

Calculating the percent overlap of two polygons in JavaScript

To compute the overlapping percentage

  1. Compute the intersection of the two polygons

    Intersection = intersect(Precinct, Block)
  2. Divide the area of Intersection by the area of the parent polygon of interest.

    Overlap = area(Intersection) / area(Parent)
  3. It is a little unclear what you mean by the percent overlap. The parent polygon could be one of several possibilities

    a) area(Intersection) / area(Precinct)

    b) area(Intersection) / area(Block)

    c) area(Intersection) / area(Precinct union Block)

As for a javascript library, this one seems to have what you need Intersection.js

There's also the JSTS Topology Suite which can do geospatial processing in JavaScript. See Node.js examples here.

Percentage overlap of spatial polygons for a sensitivity analysis of convex hull

Here is a solution to finding the interior without any errors using spatstat
and the underlying polyclip package.

library(spatstat)

# Data from OP
set.seed(11)
dt <- data.frame(x = rnorm(1e3, 10, 3) + sample(-5:5, 1e3, replace = TRUE))
dt$y <- (rnorm(1e3, 3, 4) + sample(-10:10, 1e3, replace = TRUE)) + dt$x
dt <- rbind(dt, data.frame(x = -dt$x, y = dt$y))

# Converted to spatstat classes (`ppp` not strictly necessary -- just a habit)
X <- as.ppp(dt, W = owin(c(-25,25),c(-15,40)))
p1 <- owin(poly = dt[rev(chull(dt)),])

# Plot of data and convex hull
plot(X, main = "")
plot(p1, add = TRUE, border = "green")

# Convex hulls of sampled points in spatstat format
polys <- lapply(1:100, function(i) {
tmp <- dt[sample(rownames(dt), 1e2),]
owin(poly = tmp[rev(chull(tmp)),])
})

# Plot of convex hulls
for(i in seq_along(polys)){
plot(polys[[i]], add = TRUE, border = "red")
}

# Intersection of all convex hulls plotted in transparent blue
interior <- do.call(intersect.owin, polys)
plot(interior, add = TRUE, col = rgb(0,0,1,0.1))

Sample Image

It is not clear to me what you want to do from here, but at least this approach
avoids the errors of polygon clipping.

To do the grid based solution in spatstat I would convert the windows to
binary image masks and then work from there:

Wmask <- as.im(Window(X), dimyx = c(200, 200))
masks <- lapply(polys, as.im.owin, xy = Wmask, na.replace = 0)
maskmean <- Reduce("+", masks)/100
plot(maskmean)

Sample Image

The speed depends on the resolution you choose, but I would guess it is much
faster than the current suggestion using sp/raster (which can probably
be improved a lot using the same logic as here, so that would be another
option to stick to raster).



Related Topics



Leave a reply



Submit