如何用terra包计算点到多边形最近边的距离？

数据准备 #

• 谷歌地图经纬度坐标数据：

tpts <- data.frame(lat = c(40.91947468397188, 39.13618031492076,
                           44.54912813480803, 38.55433576050368,
                           36.833190022953396, 41.111132432433234), 
                   lon = c(-2.341558113881383, -6.780034564585769,
                           -7.219487678516898, 2.602289417843798,
                           -10.339604787427902, -4.234324))

• tidyterra包提供的多边形示例数据：

extfile <- system.file("extdata/cyl.gpkg", package = "tidyterra")
cyl <- terra::vect(extfile)

尝试的代码 #

目标是计算tpts中每个点到cyl中每个多边形最近边的距离，编写了以下代码：

# 将点转换为SpatVector格式，适配distance()的输入要求
pts.sv <- terra::vect(tpts, geom = c("lon","lat"), crs = crs("+proj=longlat +datum=WGS84"))
# 将点转换到和多边形相同的CRS（不确定是否必要？）
pts.sv <- terra::project(pts.sv, crs(cyl))

# 可视化验证图层匹配情况
ggplot() + geom_spatvector(data = cyl) +
  geom_spatvector(data = pts.sv)

# 计算每个点到每个多边形最近边的距离
terra::distance(pts.sv, cyl, unit = "m")

疑问解答 #

1. lonlat和haversine参数的作用及报错原因 #

• lonlat=TRUE：仅适用于输入数据为地理坐标系（经纬度）且未做投影转换的场景。你已经把点数据转到了cyl的投影坐标系（UTM平面坐标系），此时用这个参数会报错，因为它只针对地理坐标系下的球面距离计算。
• haversine：这是旧版terra（或raster包）的参数，当前版本的terra::distance()已经移除了该参数，所以添加后会报错。如果要在地理坐标系下算球面距离，直接用lonlat=TRUE即可（前提是数据没投影）。

2. 确认terra::distance()的计算逻辑 #

当x是点、y是多边形时，terra::distance(x, y)计算的是点到多边形最近边的距离，不是到质心或其他位置的距离：

• 点在多边形内部时，距离为0；
• 点在外部时，是点到多边形边界上最近点的直线距离（平面坐标系）或球面距离（地理坐标系+lonlat=TRUE）。
  你对比质心距离结果不同，正好验证了这一点——质心距离是点到多边形中心的距离，和到最近边的距离逻辑完全不同。

3. 投影转换的必要性 #

这一步是必须的：

• 如果直接在地理坐标系下算距离，需要开lonlat=TRUE来算球面距离；
• 如果要算平面直线距离，必须把点和多边形转到同一投影坐标系（比如cyl的UTM投影），否则经纬度的差值没法直接转换成准确的米制距离。你的代码里这一步是正确的。

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