R语言UTM32N坐标系下点到多边形最近距离异常问题求助

问题背景 #

使用R与RStudio开发，全局环境中有德国土壤多边形数据buek250_shapefile_all，可通过leaflet正常绘制。编写lowestdistancebuek函数接收经纬度参数，期望返回与该点距离最近的多边形信息，但无论输入德国境内哪个经纬度值，函数始终返回同一个多边形，且计算距离在5000-7000km之间（远超德国范围）。

错误根源 #

1. 坐标顺序颠倒：sf包中st_point遵循GIS标准的x（经度）在前，y（纬度）在后顺序，但原代码中传入的是c(inputlat, inputlon)，把纬度放在了经度的位置，导致点的位置完全错误。
2. CRS不匹配：输入的经纬度属于WGS84坐标系（EPSG:4326），但原代码直接将其赋值为UTM 32N坐标系（EPSG:25832），相当于把WGS84的经纬度数值当作UTM坐标使用，这会让点的位置偏离实际区域数千公里，自然会计算出异常距离，且匹配到错误的多边形。

修正后的代码 #

lowestdistancebuek <- function(inputlat, inputlon) {
  # 1. 按正确顺序创建WGS84坐标系的点（先经度，后纬度）
  point <- st_point(c(inputlon, inputlat)) %>%
    st_sfc(crs = 4326) %>%
    # 2. 将点转换为与多边形数据一致的CRS（确保和目标数据坐标系匹配）
    st_transform(crs = st_crs(buek250_shapefile_all))
  
  # 检查点是否在多边形内部
  nearestshape <- st_intersects(point, buek250_shapefile_all, sparse = FALSE)
  nearestshape <- which(nearestshape, arr.ind = TRUE)[,2]

  if(length(nearestshape) == 1) {
    nearestshape <- buek250_shapefile_all[nearestshape, ]
  } else if (length(nearestshape) == 0){
    # 使用by_element=TRUE返回距离向量，更便捷获取最小值索引
    distances <- st_distance(point, buek250_shapefile_all,
                             by_element = TRUE,
                             tolerance = TRUE)
    nearestshape <- buek250_shapefile_all[which.min(distances), ]
  }
  return(nearestshape)
}

验证步骤 #

• 确认多边形数据的坐标系：运行st_crs(buek250_shapefile_all)，确保输出为EPSG:25832，与数据来源的UTM 32N坐标系一致。
• 测试时可打印转换后的点坐标：print(point)，确认点的UTM坐标在德国境内的合理范围（大致x: 320000-900000，y: 5000000-6200000）。

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