Combine and calculate area of geometries that intersect the same polygon in R-CodePudding

I am working with flood data and electoral districts, both are (multi)polygon data. I have different types of flood data: coastal flooding (in blue below) and river flooding (in green below). I have to compute the total surface of each electoral district that is flooded.

I first compute the intersection between each type of flood and electoral districts:

inters_c <- st_intersection(ireED_val, coast_val)
inters_r <- st_intersection(ireED_val, river_val)

This is shown in the picture below; the light blue area is the part of coastal flooding that affects the specific ED considered. Intersection between flood and electoral districts.

An electoral district can be affected by both coastal and river flooding. The two types of flood can overlap, but not necessarily coincide, as shown in this picture. Flood types overlap.

As a result, I cannot calculate the total surface of electoral districts that is flooded by summing the areas of the intersections of each flood type, since this would lead to double counting in those cases where the two flood types overlap.

What I would have to do is to create polygons that are given by the union of the intersections for each electoral district. Initially I though st_union(inters_c, inters_r) would do the trick, but this creates the union of each polygon x with all polygons y, so it is not what I need. It would by great if st_union allowed a by=ED_ID option, but this doesn't seem to be possible. Any suggestion on how I could do this?

Thank you all in advance.

CodePudding user response：

You want to do something like:

# 1. Get intersections of both flood types with base polygon
inters_c <- st_intersection(ireED_val, coast_val)
inters_r <- st_intersection(ireED_val, river_val)

# 2. Create polygon of overlap of floods in the intersecting areas
inters_both <- st_intersection(inters_c, inters_r)

# 3. Compute total area without double-counting
inters_c$AREA   inters_r$AREA - inters_both$AREA

This is pretty direct, so I'm not sure how much more "direct" we can be. However, you could also take the union of the flood polygons then its intersection with the base polygon:

union <- st_union(coast_val, river_val)
inter <- st_intersection(union, ireED_val)
inter$AREA

Or create #2 by getting the intersection of all 3 polygons as such:

inters_both <- st_intersection(st_intersection(ireED_val, coast_val), river_val)

CodePudding user response：

Based on the suggestion by @socialscientist I came up with the following solution.

#Create intersection between flood areas and electoral districts
inters_c <- st_intersection(ireED_val, coast)
inters_r <- st_intersection(ireED_val, river)
inters_cr <- st_intersection(inters_c, inters_r)

#Compute flooding areas of intersection
inters_c$Area_flood_c <- st_area(inters_c)
inters_r$Area_flood_r <- st_area(inters_r)
inters_cr$Area_flood_cr <- st_area(inters_cr)

#Convert intersections to data frames
inters_c.df = data.frame(inters_c)
inters_r.df = data.frame(inters_r)
inters_cr.df = data.frame(inters_cr)

##Aggregate areas of all polygons of a certain flood type intersecting an ED
#Coast
flooded_areas_c = aggregate(
  x = inters_c.df$Area_flood_c,
  by = inters_c.df[c("ED_ID")],
  FUN = sum, na.rm = TRUE
)
#River
flooded_areas_r = aggregate(
  x = inters_r.df$Area_flood_r,
  by = inters_r.df[c("ED_ID")],
  FUN = sum, na.rm = TRUE
)
#Overlap coast-river
flooded_areas_cr = aggregate(
  x = inters_cr.df$Area_flood_cr,
  by = inters_cr.df[c("ED_ID")],
  FUN = sum, na.rm = TRUE
  )

#Merge flood areas to ED based on ID
ireED_val <-  merge(ireED_val, flooded_areas_c, by = "ED_ID", all.x= TRUE)
ireED_val <-  merge(ireED_val, flooded_areas_r, by = "ED_ID", all.x= TRUE)
ireED_val <-  merge(ireED_val, flooded_areas_cr, by = "ED_ID", all.x= TRUE)

#Replace NA from merge with 0 area
ireED_val$Area_flood_c[is.na(ireED_val$Area_flood_c)] <- 0
ireED_val$Area_flood_r[is.na(ireED_val$Area_flood_r)] <- 0
ireED_val$Area_flood_cr[is.na(ireED_val$Area_flood_cr)] <- 0

#Compute total flooded area removing double counting
ireED_val$Flooded_area <- ireED_val$Area_flood_c   ireED_val$Area_flood_r - ireED_val$Area_flood_cr

The code is not optimized yet. I will update it when I implement the iteration over the many flood combinations (types, scenarios, likelihoods of event) that I have. But, although not elegant, at least it is working for now.