R terra overlap vectors and return land areas

Question

I am following the example for terra::intersect

library(terra)
f <- system.file("ex/lux.shp", package="terra")
v <- vect(f)
e <- ext(5.6, 6, 49.55, 49.7)
x <- intersect(v, e)
p <- vect(c("POLYGON ((5.8 49.8, 6 49.9, 6.15 49.8, 6 49.6, 5.8 49.8))",
            "POLYGON ((6.3 49.9, 6.2 49.7, 6.3 49.6, 6.5 49.8, 6.3 49.9))"), crs=crs(v))
values(p) <- data.frame(pid=1:2, area=expanse(p))
y <- intersect(v, p)

and what I ultimately want it so summarise the land area of all the polygons created by joining p and v. What I mean is, for each polygon in p, how many ha of polygon is in Diekirch, Grevenmacher, and Luxembourg.

I try this:

lapply(y, FUN=expanse, unit="ha")
expanse(v)

yy <- union(v,p)
lapply(yy, FUN=expanse, unit="ha")

but in both cases, expanse returns the same values as

expanse(v)

Answer 1

You can do

e <- expanse(y, unit="km")
aggregate(cbind(area=e), y[,c("pid", "NAME_1"), drop=TRUE], sum)

#  pid       NAME_1       area
#1   1     Diekirch 226.271445
#2   2     Diekirch   8.147181
#3   2 Grevenmacher 286.019968
#4   1   Luxembourg 194.340843
#5   2   Luxembourg  32.281530

Or like this (I would prefer the approach above because it avoids the need to aggregate polygons)

z <- aggregate(y, c("pid", "NAME_1"))
cbind(z[, c("pid", "NAME_1"), drop=TRUE], 
                 area=expanse(z, unit="km"))

#  pid       NAME_1       area
#1   1     Diekirch 226.271445
#2   1   Luxembourg 194.340843
#3   2     Diekirch   8.147181
#4   2 Grevenmacher 286.019968
#5   2   Luxembourg  32.281530