3

I am converting several functions from the sp package to terra and I have encountered an inconsistency in the calculated area of polygons. The code sample below demonstrates my issue.

x = seq(0, 2 * pi, length.out = 100)
circle = cbind(sin(x), cos(x))

# sp/rgeos version
poly.1 = sp::SpatialPolygons(list(sp::Polygons( list(sp::Polygon(circle) ), "circle")))
area.1 = rgeos::gArea(poly.1)

# terra version
poly.2 = terra::vect(circle, type = "polygons", crs = "local")
area.2 = terra::expanse(poly.2, transform = F)

area.1 - area.2 # != 0

I have understood that poly.1 and poly.2 should be describing the same polygon. The area calculations for each are not the same though (similar, but different enough to break compatibility with existing code).

Can anyone piont out whether I am doing something wrong? Or is there an approach that could get terra to deliver the same value as area.1?

randr
  • 255
  • 1
  • 7

1 Answers1

0

The difference between the two values is -2.479456e-09. That is minuscule, and probably due to floating point imprecision. I would not say that the results are similar; they are essentially the same.

To get the same results, you can round the area to 6 (or in this case 8) decimals.

round(area.1, 6) == round(area.2, 6)
#[1] TRUE
Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63
  • After some testing (and familiarising myself with the terra package), it appears that all my differences are indeed on the scale of floating point precision. – randr Apr 12 '23 at 14:03