terra nearest()设centroids=FALSE时to_id列全为NA的解决问询

问题背景 #

现有一批海上地理点数据和瑞典海域地图数据，需完成两个需求：

• 确定每个点所属的海域
• 对不在海域内的点，计算其到最近海域的距离并获取对应海域的HID

使用terra包执行nearest函数后，返回结果的to_id列全部为NA且类型为逻辑型，无法获取对应海域的HID值。

问题原因 #

当nearest函数设置centroids=FALSE时，若点直接落在海域多边形内部，函数会返回距离为0，但不会自动关联对应多边形的ID（即HID），导致to_id列出现NA。要实现需求，需先区分在海域内的点和不在海域内的点，分别处理后再合并结果。

解决方案 #

步骤1：判断点是否在海域内并获取所属HID #

使用extract函数直接获取点所在多边形的属性，同时标记出不在任何海域内的点。

步骤2：对不在海域内的点，计算到最近海域的距离与HID #

对标记为“不在海域内”的点，用nearest函数找到最近的海域多边形，再关联其HID值。

完整代码示例 #

library(terra)
#> terra 1.7.71

# 加载海域地图数据并转换坐标系
map <- vect("/vsizip//vsicurl/https://www.smhi.se/polopoly_fs/1.140307!/Havsomr_SVAR_2016_3b.zip")
map <- project(map, "EPSG:3006")

# 构建点数据并转换坐标系
pts <- data.frame(id = 11:25,
                  lon = c(15.29416656, 15.10333347, 15.27083302, 15.27250004, 15.1291666, 
                          15.31333351, 15.24166679, 15.30716705, 15.32333374, 15.38833332, 
                          15.41416645, 15.40666676, 15.27166653, 15.26083374, 15.1916666),
                  lat = c(56.1566658, 56.15000153, 56.15250015, 56.13999939, 56.16383362, 
                          56.15333176, 56.1558342, 56.17416763, 56.12366486, 56.14250183, 
                          56.15916824, 56.17666626, 56.17666626, 56.17300034, 56.1570015)) |> 
  vect(geom = c("lon", "lat"), crs = "EPSG:4326") |> 
  project("EPSG:3006")

# 步骤1：提取点所在海域的HID，标记是否在海域内
inside <- extract(map, pts, "HID")
pts$inside_HID <- inside$HID
pts$is_inside <- !is.na(pts$inside_HID)

# 步骤2：处理不在海域内的点
pts_out <- pts[!pts$is_inside, ]
if (nrow(pts_out) > 0) {
  nearest_res <- nearest(pts_out, map, centroids = FALSE)
  # 通过to_id索引关联对应海域的HID
  nearest_res$to_HID <- map$HID[nearest_res$to_id]
  # 合并结果到原始点数据
  pts$nearest_HID <- NA
  pts$distance_to_sea <- 0
  pts[!pts$is_inside, "nearest_HID"] <- nearest_res$to_HID
  pts[!pts$is_inside, "distance_to_sea"] <- nearest_res$distance
} else {
  pts$nearest_HID <- NA
  pts$distance_to_sea <- 0
}

# 整理最终的海域ID：内部点用所属HID，外部点用最近海域HID
pts$final_HID <- ifelse(pts$is_inside, pts$inside_HID, pts$nearest_HID)

# 查看结果
head(pts)

结果说明 #

• final_HID列：每个点对应的海域ID（在海域内则为所属海域HID，不在则为最近海域HID）
• distance_to_sea列：点到最近海域的距离（在海域内则为0）
• is_inside列：标记点是否在海域内部

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