Search code examples
rterra

Summing polygons by value using terra in R

I have a terra Spatvector with many polygons and I would like to sum the polygons, so that the new set of polygons is the sum of all the overlapping values. This is the same question as this, but I would like to use the terra package in R.

Here is a reprex to create suitable data, and I tried using aggregate() but that doesn't quite give what I'm after. I could rasterize all the polygons and then sum the rasters, but I feel like there must be a simple vector math type solution!

library(terra)
#> terra 1.8.15

v <- vect(system.file("ex/lux.shp", package="terra"))

#create two circles overlaying vector v
circles <- vect(rbind(c(6.05, 49.6), c(6.1, 49.8)), crs = "epsg:4326") |>
  buffer(10e3)

circles$ID_1 <- 4:5

#bind everything into a single vector
v2 <- rbind(v, circles)

plot(v2, "ID_1", alpha = 0.8)



#I want to sum the polygons using "ID_1"
#aggregating doesn't sum overlapping areas: the bottom circle should have a value of 3+4 = 7, but it is just overlaying the aggregated polygon with value (ID_1) = 3
v2_agg <- aggregate(v2, by = "ID_1", fun = "sum") 
  
plot(v2_agg, "ID_1", alpha = 0.8)

Created on 2025-02-11 with reprex v2.1.1

Solution

  • With your example data

    library(terra)
v <- vect(system.file("ex/lux.shp", package="terra"))
circles <- vect(rbind(c(6.05, 49.6), c(6.1, 49.8)), crs = "epsg:4326") |> buffer(10e3)
circles$ID_1 <- 4:5
v2 <- rbind(v, circles)

    You can do

    u <- union(v2)
s <- colSums(v2$ID_1 * t(values(u)))
values(u) <- data.frame(sum=s)

    Below I show the rasterization approach. That is not the ideal solution unless you wanted to end up with a raster anyway. But there is no need to rasterize the polygons separately as you can provide argument fun="sum".

    r <- rast(v2, res=.005)
x <- rasterize(v2, r, "ID_1", fun="sum")
p <- as.polygons(x)

    Compare the two approaches by rasterizing the first one

    y <- rasterize(u, r, "sum", wopt=list(names="ID_1"))
all.equal(x, y)
#[1] TRUE