为何rgeos/sp与terra计算多边形面积的结果存在差异？

我正在将多个函数从sp包迁移至terra包，过程中发现多边形面积的计算结果存在不一致性，以下代码示例可展示该问题：

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

我认为poly.1与poly.2描述的是同一个多边形，但二者的面积计算结果并不相同（虽相近，但差异足以破坏现有代码的兼容性）。请问我是否操作有误？或者有没有办法让terra计算出与area.1相同的结果？

解答 #

• 操作本身没有错误，结果差异来源于两个工具依赖的底层面积计算算法不同：

  • rgeos::gArea基于GEOS地理空间库实现面积计算，采用针对多边形环的精确几何计算逻辑
  • terra::expanse在transform=F参数下，使用terra自身的坐标面积计算逻辑，与GEOS的实现存在细微差异

• 要让terra输出与gArea完全一致的结果，直接调用terra的geos_area()函数即可——该函数同样依赖GEOS库，和rgeos::gArea使用完全相同的底层计算逻辑：

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::geos_area(poly.2)

area.1 - area.2 # 结果趋近于0，可认为完全一致

补充说明：

• 若多边形使用地理坐标系（经纬度），geos_area()会自动处理球面面积计算，结果与rgeos::gArea(poly, spheroid=TRUE)匹配
• 若多边形使用平面坐标系，geos_area()的结果与默认的rgeos::gArea()完全一致

内容的提问来源于stack exchange，提问作者randr