Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Summing polygons by value using terra in R

Tags:

r

terra

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

like image 561
Jason Avatar asked Nov 21 '25 21:11

Jason


1 Answers

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
like image 65
Robert Hijmans Avatar answered Nov 24 '25 12:11

Robert Hijmans



Donate For Us

If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!