如何确定polygon_sf2最近点相对polygon_sf1的8方位方向

问题描述 #

如何找出polygon_sf2的最近点相对于polygon_sf1的方向，并将其归类为以下八个方向之一：north（北）、north-east（东北）、east（东）、south-east（东南）、south（南）、south-west（西南）、west（西）、north-west（西北）？在本次案例中，该方向应为东北或东。

初始代码 #

library(sf)

# 创建多边形
coords <- matrix(c(
  -1.6522, 55.571,
  -1.6487, 55.568,
  -1.6522, 55.565,
  -1.6550, 55.568,
  -1.6522, 55.571
), ncol = 2, byrow = TRUE)

polygon <- st_sfc(st_polygon(list(coords)), crs = 4326)

polygon_sf1 <- st_sf(geometry = polygon)

coords <- matrix(c(
  -1.645, 55.575,
  -1.6449, 55.571,
  -1.6472, 55.569,
  -1.649, 55.571,
  -1.645, 55.575
), ncol = 2, byrow = TRUE)

polygon2 <- st_sfc(st_polygon(list(coords)), crs = 4326)

polygon_sf2 <- st_sf(geometry = polygon2)

# 计算两个多边形间的最近点距离
st_distance(polygon_sf1, polygon_sf2)

# 需确定polygon_sf2最近点相对polygon_sf1的8方位方向，本次案例应为东北或东

解决方案 #

可以通过获取最近点对→计算方位角→映射8方向三步实现需求：

1. 获取两个多边形的最近点对

用st_nearest_points提取最近点连线，再拆分出两个多边形的最近点坐标：

# 获取最近点连线
nearest_line <- st_nearest_points(polygon_sf1, polygon_sf2)

# 提取端点坐标：pt1是polygon_sf1的最近点，pt2是polygon_sf2的最近点
pt_coords <- st_coordinates(nearest_line)[, c("X", "Y")]
pt1 <- pt_coords[1, ]
pt2 <- pt_coords[2, ]

2. 计算从pt1指向pt2的方位角

基于经纬度差计算方位角（从正北顺时针为正，范围0-360度）：

# 计算经纬度差值
delta_lon <- pt2[1] - pt1[1]
delta_lat <- pt2[2] - pt1[2]

# 计算方位角并转换为0-360度范围
bearing <- atan2(delta_lon, delta_lat) * (180 / pi)
bearing <- ifelse(bearing < 0, bearing + 360, bearing)

3. 将方位角映射到8个方向

根据方位角区间匹配对应的方向：

# 定义各方向的方位角区间
direction_ranges <- list(
  north = c(337.5, 22.5),
  north_east = c(22.5, 67.5),
  east = c(67.5, 112.5),
  south_east = c(112.5, 157.5),
  south = c(157.5, 202.5),
  south_west = c(202.5, 247.5),
  west = c(247.5, 292.5),
  north_west = c(292.5, 337.5)
)

# 匹配方位角对应的方向
get_direction <- function(bearing) {
  for (dir in names(direction_ranges)) {
    range <- direction_ranges[[dir]]
    if (range[1] < range[2]) {
      if (bearing > range[1] && bearing <= range[2]) return(dir)
    } else {
      # 处理跨0度的north区间
      if (bearing > range[1] || bearing <= range[2]) return(dir)
    }
  }
}

# 获取最终方向
direction <- get_direction(bearing)
cat("polygon_sf2最近点相对于polygon_sf1的方向：", direction, "\n")

运行结果

整合代码运行后，输出为：

polygon_sf2最近点相对于polygon_sf1的方向： north_east

符合案例预期的东北方向。

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