如何用spatialEco包计算多POINT要素加权质心？投影问题求解

1. 将投影转换为距离单位 #

wt.centroid要求输入的几何数据必须是投影坐标系（距离单位，如米），而你的数据当前是地理坐标系（经纬度，角度单位）。可以用sf包的st_transform()函数完成转换：

步骤： #

• 选择适合数据区域的投影坐标系：你的数据坐标位于欧洲西北部（法国附近），推荐使用UTM 30N投影（EPSG:32630），该投影以米为单位，覆盖该区域。
• 转换投影后再调用wt.centroid。

代码示例： #

library(sf)
library(spatialEco)

# 加载数据（从dput输出）
cities <- structure(list(weight = c(12, 8, 3, 15), geometry = structure(list(
    structure(c(-1.75245754632, 49.5393529841), class = c("XY", 
    "POINT", "sfg")), structure(c(1.09509946531, 49.4104887444
    ), class = c("XY", "POINT", "sfg")), structure(c(1.11637971624, 
    49.3172603768), class = c("XY", "POINT", "sfg")), structure(c(0.827643384885, 
    49.8731818649), class = c("XY", "POINT", "sfg"))), class = c("sfc_POINT", 
"sfc"), precision = 0, bbox = structure(c(xmin = -1.75245754632, 
ymin = 49.3172603768, xmax = 1.11637971624, ymax = 49.8731818649
), class = "bbox"), crs = structure(list(input = "EPSG:4236", 
    wkt = "GEOGCRS[\"Hu Tzu Shan 1950\",\n    DATUM[\"Hu Tzu Shan 1950\",\n        ELLIPSOID[\"International 1924\",6378388,297,\n            LENGTHUNIT[\"metre\",1]]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    USAGE[\n        SCOPE[\"Geodesy.\"],\n        AREA[\"Taiwan, Republic of China - onshore - Taiwan Island, Penghu (Pescadores) Islands.\"],\n        BBOX[21.87,119.25,25.34,122.06]],\n    ID[\"EPSG\",4236]]"), class = "crs"), n_empty = 0L)), row.names = c(NA, 
-4L), sf_column = "geometry", agr = structure(c(weight = NA_integer_), levels = c("constant", 
"aggregate", "identity"), class = "factor"), class = c("sf", 
"tbl_df", "tbl", "data.frame"))

# 转换为UTM 30N投影（米为单位）
cities_proj <- st_transform(cities, crs = 32630)

# 计算加权质心
weighted_centroid <- wt.centroid(cities_proj$geometry, cities_proj$weight)

# 可选：将结果转回原地理坐标系查看经纬度
weighted_centroid_latlon <- st_transform(weighted_centroid, crs = 4236)

2. 无需投影转换的替代方法 #

如果不想转换投影，可以基于球面坐标系直接计算加权质心，避免角度单位带来的偏差：

方法1：使用geosphere包的weightedMean #

geosphere包的weightedMean函数支持直接对经纬度点进行加权平均，基于球面几何计算：

library(geosphere)

# 提取经纬度坐标矩阵
coords <- st_coordinates(cities$geometry)

# 计算球面加权质心
weighted_centroid_geo <- weightedMean(coords, w = cities$weight)

# 转换为sf点对象（可选）
weighted_centroid_geo_sf <- st_sfc(st_point(weighted_centroid_geo), crs = st_crs(cities))

方法2：手动实现球面加权质心 #

通过将经纬度转换为三维笛卡尔坐标，加权平均后再转回经纬度：

# 定义经纬度转笛卡尔坐标函数
latlon_to_cartesian <- function(lat, lon) {
  lat_rad <- deg2rad(lat)
  lon_rad <- deg2rad(lon)
  x <- cos(lat_rad) * cos(lon_rad)
  y <- cos(lat_rad) * sin(lon_rad)
  z <- sin(lat_rad)
  c(x, y, z)
}

# 计算每个点的笛卡尔坐标并加权
cartesian_coords <- t(apply(st_coordinates(cities$geometry), 1, function(p) {
  latlon_to_cartesian(p[2], p[1]) * cities$weight[which(st_coordinates(cities$geometry) == p, arr.ind = TRUE)[1,1]]
}))

# 加权平均
mean_cartesian <- colSums(cartesian_coords) / sum(cities$weight)

# 笛卡尔坐标转回经纬度
lon_rad <- atan2(mean_cartesian[2], mean_cartesian[1])
lat_rad <- atan2(mean_cartesian[3], sqrt(mean_cartesian[1]^2 + mean_cartesian[2]^2))
weighted_centroid_manual <- c(rad2deg(lon_rad), rad2deg(lat_rad))

# 转换为sf点对象
weighted_centroid_manual_sf <- st_sfc(st_point(weighted_centroid_manual), crs = st_crs(cities))

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