Inconsistency between raster and polygon area

[AMBarbosa]

Hi,

I expected these expanse() outputs to be the same, but they vary slightly:

r <- rast(system.file("ex/elev.tif", package="terra"))
plot(r)
r_cellsize <- cellSize(r, mask = TRUE)
plot(r_cellsize)
r_pol <- as.polygons(r*0)
plot(r_pol, add = TRUE)

expanse(r)  # 2563610102
sum(data.frame(r_cellsize))  # same
expanse(r_pol)  # 2563614794

[rhijmans]

You get the same number as for the raster data if you do

as.polygons(r, aggregate=FALSE) |> expanse() |> sum()
#[1] 2563610102

And that is not unexpected given how the raster algorithm works.

So the question is why a single polygon covering the same area does not return the same number.

as.polygons(r*0, aggregate=TRUE) |> expanse() |> sum()
#[1] 2563614794

This difference only shows up when computing the area from angular coordinates. There is no difference if you use planar coordinates (but these numbers are less precise, unless the original data source had a planar crs).

r <- rast(system.file("ex/elev.tif", package="terra"))
r <- project(r, "+proj=utm +zone=32")
r_cellsize <- cellSize(r, mask = TRUE)
r_pol <- as.polygons(r*0)
expanse(r, transform=FALSE)[1,2]
#[1] 2567718739
expanse(r_pol, transform=FALSE)
#[1] 2567718739

So there can be a difference in precision in the lon/lat area computation when you compare the area for a single polygon versus the area of same polygon subdivided into many small cells. I suspect that the division into many small cells introduces additional floating point imprecision and that the number for the undivided polygon is closer to the truth.

The difference is very small, so there does not seem to be a reason to be worried about it.