Sum/aggregation of polygon's values within other polygons using dyplr

Question

Let's say you have two sf objects made of polygons. In my case, one is administrative units (ADM) and the other is buildings (SHAPES). I would like to calculate the sum of the surface area of SHAPES by administrative units. Administrative units have an identifier which is ADM$ID

Here is a reproducible example with two administrative units and some shapes within those.

structure(list(geometry = structure(list(structure(list(structure(c(6.14080991669186, 
6.13947697270692, 6.16645297360113, 6.16780165721562, 6.14080991669186, 
43.3518254218006, 43.3317355383674, 43.3309488103689, 43.351153846062, 
43.3518254218006), .Dim = c(5L, 2L))), class = c("XY", "POLYGON", 
"sfg")), structure(list(structure(c(6.16890637059613, 6.16740543717258, 
6.19429519052069, 6.19526756384816, 6.16890637059613, 43.3510636393805, 
43.3309787504817, 43.3300719009479, 43.3504615102039, 43.3510636393805
), .Dim = c(5L, 2L))), class = c("XY", "POLYGON", "sfg"))), n_empty = 0L, crs = structure(list(
    epsg = 4326L, proj4string = "+proj=longlat +datum=WGS84 +no_defs"), class = "crs"), class = c("sfc_POLYGON", 
"sfc"), precision = 0, bbox = structure(c(xmin = 6.13947697270692, 
ymin = 43.3300719009479, xmax = 6.19526756384816, ymax = 43.3518254218006
), class = "bbox")), ID = 1:2), row.names = 1:2, sf_column = "geometry", agr = structure(c(ID = NA_integer_), class = "factor", .Label = c("constant", 
"aggregate", "identity")), class = c("sf", "data.frame"))

structure(list(geometry = structure(list(structure(list(structure(c(6.14610631720318, 
6.1444416896701, 6.14698673458045, 6.14872736274246, 6.14610631720318, 
43.3486637269545, 43.3472082709604, 43.3459195725539, 43.3473149345871, 
43.3486637269545), .Dim = c(5L, 2L))), class = c("XY", "POLYGON", 
"sfg")), structure(list(structure(c(6.14685762308362, 6.14559823469476, 
6.14883926645652, 6.15082025971966, 6.14685762308362, 43.3421058187015, 
43.3408121210182, 43.3392139218913, 43.3406600795448, 43.3421058187015
), .Dim = c(5L, 2L))), class = c("XY", "POLYGON", "sfg")), structure(list(
    structure(c(6.15692525266947, 6.15575738447906, 6.15884632257924, 
    6.16050796769147, 6.15692525266947, 43.3436054388965, 43.3425405435208, 
    43.3410621584899, 43.3424596372303, 43.3436054388965), .Dim = c(5L, 
    2L))), class = c("XY", "POLYGON", "sfg")), structure(list(
    structure(c(6.17287816794718, 6.17142829167592, 6.17667469587113, 
    6.17837798466806, 6.17287816794718, 43.3368337767714, 43.3349675471364, 
    43.3338302197931, 43.3359781405233, 43.3368337767714), .Dim = c(5L, 
    2L))), class = c("XY", "POLYGON", "sfg"))), n_empty = 0L, crs = structure(list(
    epsg = 4326L, proj4string = "+proj=longlat +datum=WGS84 +no_defs"), class = "crs"), class = c("sfc_POLYGON", 
"sfc"), precision = 0, bbox = structure(c(xmin = 6.1444416896701, 
ymin = 43.3338302197931, xmax = 6.17837798466806, ymax = 43.3486637269545
), class = "bbox"))), row.names = c(NA, 4L), class = c("sf", 
"data.frame"), sf_column = "geometry", agr = structure(integer(0), class = "factor", .Label = c("constant", 
"aggregate", "identity"), .Names = character(0)))

ADM are the two orange squares, and SHAPES are the black polygons.

I expected this simple pipe to work, I remember seing somewhere that with dyplr, it's not mandatory to use intersect functions from the sf package.

ADM %>% 
  group_by(ID) %>% 
  sum(st_area(SHAPES))

Yet, it's not working.

Answer 1

I see some issues on the approach:

• SHAPES does not have ID column as per your example.
• Personally, I prefer to extract the area to a column and summarise as in a regular data frame.

See an example:

library(sf)

ADM <- structure(list(geometry = structure(list(structure(list(structure(c(6.14080991669186, 
                                                                    6.13947697270692, 6.16645297360113, 6.16780165721562, 6.14080991669186, 
                                                                    43.3518254218006, 43.3317355383674, 43.3309488103689, 43.351153846062, 
                                                                    43.3518254218006), .Dim = c(5L, 2L))), class = c("XY", "POLYGON", 
                                                                                                                     "sfg")), structure(list(structure(c(6.16890637059613, 6.16740543717258, 
                                                                                                                                                         6.19429519052069, 6.19526756384816, 6.16890637059613, 43.3510636393805, 
                                                                                                                                                         43.3309787504817, 43.3300719009479, 43.3504615102039, 43.3510636393805
                                                                                                                     ), .Dim = c(5L, 2L))), class = c("XY", "POLYGON", "sfg"))), n_empty = 0L, crs = structure(list(
                                                                                                                       epsg = 4326L, proj4string = "+proj=longlat +datum=WGS84 +no_defs"), class = "crs"), class = c("sfc_POLYGON", 
                                                                                                                                                                                                                     "sfc"), precision = 0, bbox = structure(c(xmin = 6.13947697270692, 
                                                                                                                                                                                                                                                               ymin = 43.3300719009479, xmax = 6.19526756384816, ymax = 43.3518254218006
                                                                                                                                                                                                                     ), class = "bbox")), ID = 1:2), row.names = 1:2, sf_column = "geometry", agr = structure(c(ID = NA_integer_), class = "factor", .Label = c("constant", 
                                                                                                                                                                                                                                                                                                                                                                "aggregate", "identity")), class = c("sf", "data.frame"))

SHAPE <- structure(list(geometry = structure(list(structure(list(structure(c(6.14610631720318, 
                                                                    6.1444416896701, 6.14698673458045, 6.14872736274246, 6.14610631720318, 
                                                                    43.3486637269545, 43.3472082709604, 43.3459195725539, 43.3473149345871, 
                                                                    43.3486637269545), .Dim = c(5L, 2L))), class = c("XY", "POLYGON", 
                                                                                                                     "sfg")), structure(list(structure(c(6.14685762308362, 6.14559823469476, 
                                                                                                                                                         6.14883926645652, 6.15082025971966, 6.14685762308362, 43.3421058187015, 
                                                                                                                                                         43.3408121210182, 43.3392139218913, 43.3406600795448, 43.3421058187015
                                                                                                                     ), .Dim = c(5L, 2L))), class = c("XY", "POLYGON", "sfg")), structure(list(
                                                                                                                       structure(c(6.15692525266947, 6.15575738447906, 6.15884632257924, 
                                                                                                                                   6.16050796769147, 6.15692525266947, 43.3436054388965, 43.3425405435208, 
                                                                                                                                   43.3410621584899, 43.3424596372303, 43.3436054388965), .Dim = c(5L, 
                                                                                                                                                                                                   2L))), class = c("XY", "POLYGON", "sfg")), structure(list(
                                                                                                                                                                                                     structure(c(6.17287816794718, 6.17142829167592, 6.17667469587113, 
                                                                                                                                                                                                                 6.17837798466806, 6.17287816794718, 43.3368337767714, 43.3349675471364, 
                                                                                                                                                                                                                 43.3338302197931, 43.3359781405233, 43.3368337767714), .Dim = c(5L, 
                                                                                                                                                                                                                                                                                 2L))), class = c("XY", "POLYGON", "sfg"))), n_empty = 0L, crs = structure(list(
                                                                                                                                                                                                                                                                                   epsg = 4326L, proj4string = "+proj=longlat +datum=WGS84 +no_defs"), class = "crs"), class = c("sfc_POLYGON", 
                                                                                                                                                                                                                                                                                                                                                                                 "sfc"), precision = 0, bbox = structure(c(xmin = 6.1444416896701, 
                                                                                                                                                                                                                                                                                                                                                                                                                           ymin = 43.3338302197931, xmax = 6.17837798466806, ymax = 43.3486637269545
                                                                                                                                                                                                                                                                                                                                                                                 ), class = "bbox"))), row.names = c(NA, 4L), class = c("sf", 
                                                                                                                                                                                                                                                                                                                                                                                                                                        "data.frame"), sf_column = "geometry", agr = structure(integer(0), class = "factor", .Label = c("constant", 
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                        "aggregate", "identity"), .Names = character(0)))

Let's have a look to the objects:

ADM
#> Simple feature collection with 2 features and 1 field
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 6.139477 ymin: 43.33007 xmax: 6.195268 ymax: 43.35183
#> Geodetic CRS:  WGS 84
#>                         geometry ID
#> 1 POLYGON ((6.14081 43.35183,...  1
#> 2 POLYGON ((6.168906 43.35106...  2

SHAPE
#> Simple feature collection with 4 features and 0 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 6.144442 ymin: 43.33383 xmax: 6.178378 ymax: 43.34866
#> Geodetic CRS:  WGS 84
#>                         geometry
#> 1 POLYGON ((6.146106 43.34866...
#> 2 POLYGON ((6.146858 43.34211...
#> 3 POLYGON ((6.156925 43.34361...
#> 4 POLYGON ((6.172878 43.33683...



plot(st_geometry(ADM))
plot(st_geometry(SHAPE), add=TRUE)

We can confirm that SHAPE does not have ID. We can get that using st_intersection:

# Join to get ID on SHAPE

SHAPE_INT <- st_intersection(SHAPE,ADM)
#> although coordinates are longitude/latitude, st_intersection assumes that they are planar
#> Warning: attribute variables are assumed to be spatially constant throughout all
#> geometries

SHAPE_INT
#> Simple feature collection with 4 features and 1 field
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 6.144442 ymin: 43.33383 xmax: 6.178378 ymax: 43.34866
#> Geodetic CRS:  WGS 84
#>   ID                       geometry
#> 1  1 POLYGON ((6.148727 43.34731...
#> 2  1 POLYGON ((6.15082 43.34066,...
#> 3  1 POLYGON ((6.160508 43.34246...
#> 4  2 POLYGON ((6.178378 43.33598...

plot(SHAPE_INT)

Fine! Now we have IDs on SHAPE_INT. We can proceed to compute the areas on a new column, remove the geometry with st_drop_geometry and then just summarise the areas by ID:

library(dplyr)

SHAPE_INT$area <- st_area(SHAPE_INT)

SHAPE_INT %>% 
  st_drop_geometry() %>%
  group_by(ID) %>% 
  summarise(area = sum(area))
#> # A tibble: 2 x 2
#>      ID     area
#>   <int>    [m^2]
#> 1     1 173781.5
#> 2     2 111317.5

^{Created on 2021-05-27 by the reprex package (v2.0.0)}