terra::buffer()功能异常：生成指定面积圆形多边形结果不符

核心问题是**terra::expanse()在EPSG:2150这类投影坐标系下，实际计算的是椭球面积而非文档描述的平面面积**，而terra::buffer()在投影坐标系下创建的是平面圆，两者面积计算逻辑不匹配，导致结果偏差。

具体细节：

• 你用expanse(test[whichpoly])得到的是原多边形的真实地理椭球面积。
• 基于这个面积计算的radius2，对应的是平面圆的半径（该平面圆的平面面积等于原多边形的椭球面积）。
• EPSG:2150是Lambert等角投影，卢森堡所在区域的面积变形因子约为1.7，因此平面圆的真实椭球面积会是平面面积的1/1.7，而expanse(buff2)返回的正是这个椭球面积，最终导致两者比值约为1.7。

解决方案 #

根据你的需求选择以下两种方式：

方式一：创建与原多边形真实地理面积相等的圆 #

这种方式得到的圆在真实地理空间中的面积和原多边形一致：

library(terra)

whichpoly <- 8

# 导入数据并转换为WGS84地理坐标系
v <- vect(system.file("ex/lux.shp", package="terra"))
test_wgs84 <- project(v, "EPSG:4326")
cent_wgs84 <- centroids(test_wgs84[whichpoly])

# 获取原多边形的真实地理面积
target_area <- expanse(test_wgs84[whichpoly])

# 迭代法找到对应测地线圆的半径
find_radius <- function(target, cent, tol=1e-6) {
  r <- sqrt(target/pi) # 初始值
  current_area <- expanse(buffer(cent, r))
  while(abs(current_area - target) > tol) {
    r <- r * sqrt(target/current_area)
    current_area <- expanse(buffer(cent, r))
  }
  return(r)
}

radius_geo <- find_radius(target_area, cent_wgs84)
buff_geo <- buffer(cent_wgs84, radius_geo)

# 转换回EPSG:2150可视化
buff_proj <- project(buff_geo, "EPSG:2150")

# 验证面积
expanse(test_wgs84[whichpoly])
expanse(buff_geo)
expanse(test_wgs84[whichpoly])/expanse(buff_geo) # 结果接近1

方式二：创建与原多边形投影平面面积相等的圆 #

这种方式得到的圆在EPSG:2150投影平面中的面积和原多边形一致：

library(terra)

whichpoly <- 8

# 导入数据并投影到EPSG:2150
v <- vect(system.file("ex/lux.shp", package="terra"))
test <- project(v, "EPSG:2150")
cent <- centroids(test[whichpoly])

# 用鞋带公式计算投影平面面积
shoelace_area <- function(poly) {
  coords <- crds(poly)[[1]]
  n <- nrow(coords)
  area <- abs(sum(coords[1:n-1,1]*coords[2:n,2]) - sum(coords[1:n-1,2]*coords[2:n,1]))/2
  return(area)
}

target_plane_area <- shoelace_area(test[whichpoly])

# 计算平面圆的半径
radius_plane <- sqrt(target_plane_area/pi)
buff_plane <- buffer(cent, radius_plane)

# 验证平面面积
shoelace_area(buff_plane)
target_plane_area
shoelace_area(buff_plane)/target_plane_area # 结果接近1

额外说明 #

当前terra版本（包括你使用的1.7-3）在处理部分投影坐标系（如EPSG:2150）的expanse()计算时，存在与文档描述不符的情况，实际返回椭球面积而非平面面积。若需精确平面面积，建议使用鞋带公式自行计算。

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