使用R的sf与st_distance()匹配最近经纬度点结果不一致问题

问题背景 #

拥有大型站点数据集，每行对应一个站点，需求是每年范围内，为每个站点匹配使用不同设备类型（gear）的最近站点，并将站点对信息合并或关联。现有实现结果中部分站点关联了非最近的站点（绘图可见：顶部站点本该关联右侧站点，却被关联到下方站点），同时需要更高效的处理方法。

错误原因分析 #

核心错误在于sf坐标顺序搞反：

• sf包要求坐标输入必须是先经度（longitude），后纬度（latitude）
• 你创建sf对象时用了coords=5:6，对应lat和lon列（即纬度在前，经度在后），这导致距离计算时坐标错位，最终匹配出错误的最近邻站点。

修正方案 + 高效实现 #

1. 修正坐标顺序的基础版本 #

先修正sf对象的坐标顺序，确保距离计算正确：

library(ggplot2)
library(sf)
library(dplyr)
library(purrr)

# 生成数据
set.seed(123)
latitude <- runif(100, 72, 81)
longitude <- runif(100, 20, 60)
gear <- factor(sample(1:2, 100, replace = TRUE))
year <- factor(sample(c(2020, 2021), 100, replace = TRUE))
orig.data <- data.frame(latitude, longitude, gear, year)

# 关键修正：sf坐标必须是【经度，纬度】，直接指定列名避免顺序错误
df <- st_as_sf(orig.data, coords = c("longitude", "latitude"), crs = 4326) # WGS84坐标系

# 定义分组内的最近邻匹配函数
match_nearest <- function(group_df) {
  # 拆分不同gear的站点
  gear1 <- group_df %>% filter(gear == 1)
  gear2 <- group_df %>% filter(gear == 2)
  
  if(nrow(gear1) == 0 || nrow(gear2) == 0) return(NULL)
  
  # 给gear1找最近的gear2站点
  nn_1to2 <- st_nearest_feature(gear1, gear2)
  gear1_matched <- gear1 %>%
    mutate(
      match_id = row.names(gear2)[nn_1to2],
      match_lon = st_coordinates(gear2)[nn_1to2, 1],
      match_lat = st_coordinates(gear2)[nn_1to2, 2],
      match_gear = gear2$gear[nn_1to2],
      distance = st_distance(gear1, gear2[nn_1to2, ], by_element = TRUE)
    )
  
  # 给gear2找最近的gear1站点
  nn_2to1 <- st_nearest_feature(gear2, gear1)
  gear2_matched <- gear2 %>%
    mutate(
      match_id = row.names(gear1)[nn_2to1],
      match_lon = st_coordinates(gear1)[nn_2to1, 1],
      match_lat = st_coordinates(gear1)[nn_2to1, 2],
      match_gear = gear1$gear[nn_2to1],
      distance = st_distance(gear2, gear1[nn_2to1, ], by_element = TRUE)
    )
  
  # 合并结果
  bind_rows(gear1_matched, gear2_matched) %>%
    st_drop_geometry() # 移除sf几何列，转为普通数据框
}

# 按年份分组执行匹配
result <- orig.data %>%
  group_split(year) %>%
  map_dfr(match_nearest)

# 筛选gear=1的匹配结果用于绘图
nnij2 <- result %>% filter(gear == 1)

# 绘图验证
ggplot(data = nnij2, aes(x = longitude, y = latitude, shape = gear)) +
  geom_point(size = 3) +
  geom_point(aes(x = match_lon, y = match_lat, shape = match_gear), color = "red", size = 3) +
  geom_segment(aes(x = longitude, y = latitude, xend = match_lon, yend = match_lat, colour = distance)) +
  facet_wrap(~year)

2. 超大型数据集的高效实现 #

如果数据量极大（10万+站点），生成全距离矩阵会导致内存溢出，推荐使用nngeo包的空间索引方法，速度更快：

library(nngeo)
library(data.table)

# 转为data.table提升处理速度
dt <- as.data.table(orig.data)
dt[, geometry := st_sfc(st_point(c(longitude, latitude))), by = 1:nrow(dt)]
dt <- st_as_sf(dt, crs = 4326)

# 按年份分组匹配
result_dt <- dt[, {
  current_year <- .SD
  gear1 <- current_year[gear == 1]
  gear2 <- current_year[gear == 2]
  
  if(nrow(gear1) == 0 || nrow(gear2) == 0) return(NULL)
  
  # 找gear1到gear2的最近邻
  nn1 <- st_nn(gear1, gear2, k = 1, returnDist = TRUE)
  gear1[, `:=`(
    match_id = gear2$.[nn1$nn],
    match_lon = gear2$longitude[nn1$nn],
    match_lat = gear2$latitude[nn1$nn],
    match_gear = gear2$gear[nn1$nn],
    distance = nn1$dist
  )]
  
  # 找gear2到gear1的最近邻
  nn2 <- st_nn(gear2, gear1, k = 1, returnDist = TRUE)
  gear2[, `:=`(
    match_id = gear1$.[nn2$nn],
    match_lon = gear1$longitude[nn2$nn],
    match_lat = gear1$latitude[nn2$nn],
    match_gear = gear1$gear[nn2$nn],
    distance = nn2$dist
  )]
  
  rbind(gear1, gear2)
}, by = year] %>% st_drop_geometry()

关键优化点 #

• 避免全距离矩阵：原方法生成n×n的距离矩阵，数据量稍大就会内存溢出；分组+空间索引的方法只计算必要的最近邻距离
• 明确坐标系：指定crs=4326（WGS84），如需精确球面距离，可转换为UTM等投影坐标系后再计算
• 分组处理：按年份分组后匹配，确保不会跨年份关联站点

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