如何让gCentroid在地球极点区域正确计算地表空间质心？

问题原因 #

gCentroid 函数本身仅支持平面坐标系下的质心计算，你直接传入WGS84（EPSG:4326）地理坐标时，它会把经度、纬度视为普通XY平面坐标做算术平均，完全不考虑地球球面曲率，因此高纬度区域的计算结果不符合球面质心预期。你示例中四个点的经度算术平均为0、纬度算术平均为80，所以输出了(0,80)的结果。

解决方法 #

方法1：投影转换后用gCentroid计算 #

将地理坐标投影到适配极区的方位投影下，完成质心计算后再转换回WGS84坐标系，代码示例如下：

library(rgeos)
library(rgdal)

# 构造原始点数据
pts <- SpatialPoints(
  coords = data.frame(
    longitude = c(-135,-45,45,135),
    latitude = c(80,80,80,80)
  ),
  proj4string = CRS('EPSG:4326')
)

# 定义以北极为中心的方位等距投影，适配极区计算场景
polar_proj <- CRS("+proj=aeqd +lat_0=90 +lon_0=0 +datum=WGS84 +units=m +no_defs")

# 原始点坐标投影转换
pts_polar <- spTransform(pts, polar_proj)

# 计算投影后平面坐标的质心
centroid_polar <- gCentroid(pts_polar)

# 质心坐标转换回WGS84经纬度
centroid_wgs84 <- spTransform(centroid_polar, CRS('EPSG:4326'))
centroid_wgs84

运行后输出的质心坐标接近北极点(0,90)，符合预期。

方法2：使用支持球面计算的工具包 #

直接用原生支持球面地理计算的包（如sf）完成计算，无需手动转换投影，代码示例如下：

library(sf)

# 构造sf格式点数据
pts_sf <- st_as_sf(
  data.frame(
    longitude = c(-135,-45,45,135),
    latitude = c(80,80,80,80)
  ),
  coords = c("longitude", "latitude"),
  crs = 4326
)

# 开启s2球面计算支持，直接计算球面质心
st_centroid(pts_sf, use_s2 = TRUE)

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