terra::distance计算点与多边形距离返回度数而非千米的问题

这不是bug，核心是你忽略了地理坐标系与投影坐标系的处理差异：

问题原因 #

EPSG:4326是地理坐标系（经纬度），而非平面投影坐标系。terra::distance的unit参数仅对平面投影坐标系生效，当输入为地理坐标系时，函数默认按平面坐标（度数）计算距离，此时unit="km"不会触发球面/椭球面距离计算，返回结果自然是度数单位（你看到的1.118度，对应赤道上约123.7km，和perim计算结果匹配）。

而perim(as.lines(...))能得到正确公里数，是因为perim处理地理坐标系的线对象时，会自动采用地理距离算法计算球面周长，除以1000后得到千米数。

正确解决方案 #

有两种可靠方法计算地理坐标系下的点与多边形顶点的真实球面距离：

方法1：转换到投影坐标系计算 #

将数据转换为平面投影坐标系（如UTM分区投影，需匹配数据所在地理区域），再用distance计算：

library(terra)

pol <- vect("POLYGON ((1 1, 3 1, 3 2, 4 3, 3 4, 2 4, 1 1))",crs="EPSG:4326")
pt <- vect(cbind(3.5,5),crs="EPSG:4326")

# 转换至UTM 32N投影（EPSG:32632，匹配示例数据的纬度范围）
pol_utm <- project(pol, "EPSG:32632")
pt_utm <- project(pt, "EPSG:32632")

# 计算距离（单位米），转换为千米
distance(pt_utm, pol_utm)/1000

方法2：使用专门的地理距离函数terra::geodist #

直接用geodist计算地理坐标系下的球面/椭球面距离，无需转换投影：

library(terra)

pol <- vect("POLYGON ((1 1, 3 1, 3 2, 4 3, 3 4, 2 4, 1 1))",crs="EPSG:4326")
pt <- vect(cbind(3.5,5),crs="EPSG:4326")

# 提取多边形第5个顶点
pol_pt5 <- as.points(pol)[5,]

# 计算两点间地理距离，单位千米
geodist(pt, pol_pt5, unit="km")

两种方法的计算结果都会和perim方法的123.7km左右一致。

内容的提问来源于stack exchange，提问作者schosse-sitzer