How to Calculate the Area of Polygon Overlap in R

How to calculate the area of polygon overlap in R?

EDIT: these days I would use the 'intersect', 'cover', 'erase', 'union' and related functions in the 'raster' package. They do the hard work to keep the top-level object and attributes.

ORIG:
You could use the rgeos package with its gIntersection function. Successive calls between pairs and resulting intersections will get you there. See

library(rgeos)
?gIntersection

You will need to get into the structure of "SpatialPolygons" in the sp package to get the final coordinates. See the vignette("sp").

R: Calculating overlap polygon area

Here is a method slightly different that yours (but only slightly).

Inspection of switzerland@data reveals that, while there are 11 FeatureIDs (representing ethnicitity's), there are only 4 unique named ethnicity's (German, Italian, and French Swiss, and Rhaetoromanians). So the result below is based on the names, not the IDs.

library(rgeos)    # for gIntersection(...), etc.
library(rgdal)    # for readOGR(...)

setwd("<directory to accept your files>")
CH.1903 <- "+proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=600000 +y_0=200000 +ellps=bessel +towgs84=674.374,15.056,405.346,0,0,0,0 +units=m +no_defs"

print(load(url("http://gadm.org/data/rda/CHE_adm1.RData")))
gadm <- spTransform(gadm, CRS(CH.1903))
download.file("http://www.icr.ethz.ch/data/other/greg/GREG.zip","GREG.zip")
unzip("GREG.zip")
greg <- readOGR(dsn=".",layer="GREG")
greg <- spTransform(greg[greg$COW==225,],CRS(CH.1903))

gadm.ids <- gadm@data$ID_1               # Canton Nr.
greg.ids <- unique(greg@data$G1SHORTNAM) # ethnicity
get.area <- Vectorize(function(adm,reg){
  int <- gIntersection(gadm[gadm$ID_1==adm,],greg[greg$G1SHORTNAM==reg,],byid=TRUE)
  if (length(int)==0) return(0)
  gArea(int)
})
result <- outer(gadm.ids,greg.ids,FUN=get.area)
rownames(result) <- gadm.ids
colnames(result)  <- greg.ids
result <- as.data.frame(result)
totals <- rowSums(result)
result <- result/totals
result$totals <- totals/1e6
result$land.area <- sapply(rownames(result),function(p)gArea(gadm[gadm$ID_1==p,])/1e6)
result
#     German Swiss French Swiss Italian Swiss Rhaetoromanians     totals  land.area
# 531  1.000000000   0.00000000    0.00000000    0.000000e+00 1363.27027 1403.22192
# 532  1.000000000   0.00000000    0.00000000    0.000000e+00  244.32279  244.32279
# 533  1.000000000   0.00000000    0.00000000    0.000000e+00  172.40341  172.40341
# 534  1.000000000   0.00000000    0.00000000    0.000000e+00  522.12943  525.73449
# 535  1.000000000   0.00000000    0.00000000    0.000000e+00   70.03116   84.06481
# 536  0.902128658   0.09787134    0.00000000    0.000000e+00 5927.90740 5927.90847
# 537  0.188415744   0.81158426    0.00000000    0.000000e+00 1654.28729 1654.28729

Here we transform both shapefiles to CH.1903 which is a Mercator projection centered on Switzerland with units in meters. Then we identify the Canton Nrs. and the ethnicity's, and use outer(...) to cycle through both lists, calculating the area of intersection in sq.km (sq.m/1e6) using gArea(...) . The final result has one row per Canton, with the percentage of each ethnicity based on land area. $totals is the summation of the intersected areas for each Canton, and $land.area is the total geographic land area in the Canton. So this gives you an idea of the error due to incomplete overlaps between the ethnicity shapfile and the gadm shapefile.

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)

area of polygon intersecting with other polygons within buffer

I've solved my problem through the following code:

data1con<-as(data1, "sf")
geom1<-st_geometry(data1con)
buff1<-st_buffer(geom1,dist=5000)
plot(buff1, border="green")
plot(data1, border="red", add=TRUE)
plot(data2, add=TRUE)

data2con<-as(data2, "sf")
st_is_valid(data2con)

[1]  TRUE  TRUE FALSE  TRUE  TRUE  TRUE FALSE  TRUE FALSE  TRUE FALSE  TRUE  TRUE  TRUE
[15]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
[29]  TRUE  TRUE  TRUE  TRUE

st_is_valid(buff1)

TRUE TRUE TRUE TRUE TRUE TRUE

valid<-st_make_valid(data2con)
intersection<-st_intersection(buff1,valid)
plot(intersection, col="blue", add=TRUE)
area<-st_area(intersection)

picture_update

---

Related Topics