Calculating Intersections of shape files

I have merged two shape files: A ZIP code SF and a SF of the voting districts in Germany. Many zip codes are in one voting district and some are on the border of two (or more) voting districts. I need to find out what proportion (of the total ZIP code area) is in a voting district.

I tried it with sf_intersect:

overlap_areas <- st_intersection( shape_wahl,shape_PLZ) %>%
  mutate(overlap_area2 = st_area(.))

but the problem is that this code calculates the entire intersection points and not just those that are not 100% in a voting district (the area where grey and black-outlined areas are intersected)

I plotted the "intersections" and as it shows the whole map is red

pi <- st_intersection(shape_wahl,shape_PLZ)

ggplot() +
  geom_sf(data = shape_PLZ, aes(fill = category), color = "lightgrey") +
  geom_sf(data = shape_wahl, aes(), color = "black", fill = NA)+
  geom_sf(data = pi, aes(), fill = "red")

in the end I would like to have something like:

I'd be very grateful for some tips :)

You can add areas to your ZIP object, after going through st_intersecation(), you can do the same to get new areas, from there you can get ratios.

Generic example, in this case based on example dataset from sf package, might look something like this:

library(sf)
#> Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE
library(dplyr, warn.conflicts = FALSE)
library(ggplot2)

# example dataset, "ZIP polygons"
sf_1 <- 
  read_sf(system.file("shape/nc.shp", package="sf"))[, "NAME"] |>
  mutate(area_orig = st_area(geometry))
print(sf_1, n = 3)
#> Simple feature collection with 100 features and 2 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
#> Geodetic CRS:  NAD27
#> # A tibble: 100 × 3
#>   NAME                                                        geometry area_orig
#> * <chr>                                             <MULTIPOLYGON [°]>     [m^2]
#> 1 Ashe      (((-81.47276 36.23436, -81.54084 36.27251, -81.56198 36.2…    1.14e9
#> 2 Alleghany (((-81.23989 36.36536, -81.24069 36.37942, -81.26284 36.4…    6.11e8
#> 3 Surry     (((-80.45634 36.24256, -80.47639 36.25473, -80.53688 36.2…    1.42e9
#> # ℹ 97 more rows

# generate grid for "voting districts"
sf_2 <- 
  st_sf(geometry = st_make_grid(sf_1, n = c(9,3))) |> 
  mutate(id = row_number())
print(sf_2, n = 3)
#> Simple feature collection with 27 features and 1 field
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
#> Geodetic CRS:  NAD27
#> First 3 features:
#>                         geometry id
#> 1 POLYGON ((-84.32385 33.8819...  1
#> 2 POLYGON ((-83.33864 33.8819...  2
#> 3 POLYGON ((-82.35344 33.8819...  3

sf_3 <- 
  st_intersection(sf_2, sf_1) |>
  mutate(percentage = as.numeric(st_area(geometry) / area_orig)) 
#> Warning: attribute variables are assumed to be spatially constant throughout
#> all geometries
print(sf_3, n = 3)
#> Simple feature collection with 194 features and 4 fields
#> Geometry type: GEOMETRY
#> Dimension:     XY
#> Bounding box:  xmin: -84.32385 ymin: 33.88252 xmax: -75.45698 ymax: 36.58965
#> Geodetic CRS:  NAD27
#> First 3 features:
#>      id      NAME        area_orig                       geometry percentage
#> 21   21      Ashe 1137107793 [m^2] POLYGON ((-81.45289 36.2395...  0.8582486
#> 22   22      Ashe 1137107793 [m^2] POLYGON ((-81.36823 36.5740...  0.1417514
#> 22.1 22 Alleghany  610916077 [m^2] POLYGON ((-81.17667 36.4154...  1.0000000

ggplot() +
  geom_sf(data = sf_3, aes(fill = percentage)) +
  geom_sf(data = sf_1, color = "red", fill = NA) +
  geom_sf(data = sf_2, color = "gray80", linewidth = 1, fill = NA) +
  scale_fill_viridis_b() +
  theme_minimal()

^{Created on 2024-06-04 with reprex v2.1.0}