@Chris wow, this looks really fast, but the intersect with [ does not provide the same result as as terra::intersect. What does the is.related do? Given that these data are in lon/lat WGS 84, to calculate the area of all polygons resulting from intersect, should the data be projected using an equal area projection?