计算两个数据框的平均纬度差并转换为千米

核心思路 #

你的数据覆盖北美大范围，采用EPSG:4326（WGS84经纬度）。要将纬度均值的差值转换为千米，有两种可靠路径：直接计算球面距离（无需转投影，精度更高），或选择适配北美的投影系统转换后计算直线距离。

方法一：直接计算球面距离（推荐） #

经纬度是球面坐标，直接用椭球体距离公式计算两个均值点的距离，避免投影变形带来的误差。

用terra包实现 #

library(terra)
# 创建两个均值点的空间矢量对象
point1 <- vect(cbind(-104.81, 51.95062), crs = "EPSG:4326")
point2 <- vect(cbind(-104.70, 52.03755), crs = "EPSG:4326")
# 计算距离（单位米，转换为千米）
distance_km <- distance(point1, point2) / 1000
distance_km

用geosphere包实现（哈维正弦公式） #

library(geosphere)
point1 <- c(-104.81, 51.95062)
point2 <- c(-104.70, 52.03755)
# 计算距离并转千米
distance_km <- distHaversine(point1, point2) / 1000
distance_km

方法二：使用北美专用投影转换 #

如果必须通过投影转换计算，推荐选择等面积/等角的区域投影，减少大范围数据的变形：

• EPSG:3347（Canada Atlas Lambert）：适配北美大范围的等面积投影
• EPSG:26910（UTM Zone 10N）：适配北美西部的等角投影

代码示例 #

library(terra)
# 创建点对象
point1 <- vect(cbind(-104.81, 51.95062), crs = "EPSG:4326")
point2 <- vect(cbind(-104.70, 52.03755), crs = "EPSG:4326")
# 转换到Canada Atlas Lambert投影
point1_proj <- project(point1, "EPSG:3347")
point2_proj <- project(point2, "EPSG:3347")
# 计算距离并转千米
distance_km <- distance(point1_proj, point2_proj) / 1000
distance_km

关于terra::expanse()的补充 #

terra::expanse()基于WGS84椭球体的梯形面积公式计算栅格面积，和点距离的逻辑不同，但核心都依赖椭球体参数。你可以用近似值快速估算：在52°纬度附近，1°纬度对应的地面距离约为111.19千米，因此0.0869°的差值约为 0.0869 * 111.19 ≈ 9.66千米，和球面距离计算结果基本一致。

内容的提问来源于stack exchange，提问作者Sophie Père