如何用R语言及sf包基于两点构建符合方位角的矩形？

问题背景 #

需要基于给定的两个WGS84坐标系下的点，构建一个以这两点为其中一条边起止点的矩形，但现有代码生成的是轴对齐矩形，无法遵循两点的方位角方向。

用户提供的初始代码及点数据：

library(sf)
library(units)

# 点数据
point1 <- st_sfc(st_point(c(-73.51295, 45.529)), crs = 4326)
point2 <- st_sfc(st_point(c(-73.51347, 45.52816)), crs = 4326)

# 初始尝试的代码
center <- st_point(c((st_coordinates(point1)[1] + st_coordinates(point2)[1]) / 2,
                     (st_coordinates(point1)[2] + st_coordinates(point2)[2]) / 2))
half_width <- as.numeric(st_distance(point1, point2) / 2)
constant_half_height <- 100 / 111000  # 近似1度=111000米
half_width <- half_width / 111000
constant_half_height <- set_units(100, "meters")

corners <- rbind(
  st_coordinates(center) + c(-half_width, -constant_half_height),
  st_coordinates(center) + c(half_width, -constant_half_height),
  st_coordinates(center) + c(half_width, constant_half_height),
  st_coordinates(center) + c(-half_width, constant_half_height),
  st_coordinates(center) + c(-half_width, -constant_half_height)
)

rectangle <- st_polygon(list(st_linestring(corners)))

问题原因 #

初始代码直接在WGS84地理坐标系（经纬度）下进行坐标加减，生成的是轴对齐矩形，没有考虑两点连线的方位角，因此无法满足沿两点方向构建矩形的需求。另外用111000米/度的近似值会带来误差，尤其是高纬度地区。

解决方案 #

核心思路 #

1. 将地理坐标系转换为投影坐标系（如UTM），用米作为单位，避免经纬度的距离近似误差；
2. 计算两点连线的方位角，用于后续矩形的旋转变换；
3. 生成轴对齐的矩形顶点，再围绕中心点旋转对应角度，得到沿两点方向的矩形；
4. 转回WGS84地理坐标系。

完整代码 #

library(sf)
library(units)

# 定义原始点（WGS84，EPSG:4326）
point1 <- st_sfc(st_point(c(-73.51295, 45.529)), crs = 4326)
point2 <- st_sfc(st_point(c(-73.51347, 45.52816)), crs = 4326)

# 转换到UTM投影坐标系（自动匹配对应带号，单位为米）
utm_crs <- st_crs(st_transform(point1, crs = 32618)) # 该区域对应UTM 18N，EPSG:32618
point1_utm <- st_transform(point1, crs = utm_crs)
point2_utm <- st_transform(point2, crs = utm_crs)

# 计算中心点（UTM坐标系）
center_utm <- st_centroid(st_union(point1_utm, point2_utm))

# 计算两点距离的一半（矩形的半长）
half_length <- as.numeric(st_distance(point1_utm, point2_utm)) / 2

# 矩形的半宽（用户需求的100米的一半）
half_width <- 50 # 单位：米

# 计算两点连线的方位角（弧度制）
coord1 <- st_coordinates(point1_utm)
coord2 <- st_coordinates(point2_utm)
dx <- coord2[1] - coord1[1]
dy <- coord2[2] - coord1[2]
azimuth <- atan2(dy, dx)

# 生成轴对齐的矩形顶点（相对于中心点）
base_corners <- rbind(
  c(-half_length, -half_width),
  c(half_length, -half_width),
  c(half_length, half_width),
  c(-half_length, half_width),
  c(-half_length, -half_width) # 闭合多边形
)

# 旋转变换矩阵
rotation_matrix <- matrix(c(cos(azimuth), -sin(azimuth),
                            sin(azimuth), cos(azimuth)),
                          nrow = 2, ncol = 2)

# 应用旋转并平移到中心点
rotated_corners <- t(rotation_matrix %*% t(base_corners)) + st_coordinates(center_utm)

# 构建矩形多边形并转换回WGS84
rectangle_utm <- st_polygon(list(st_linestring(rotated_corners)))
rectangle_wgs84 <- st_sfc(rectangle_utm, crs = utm_crs) %>% st_transform(crs = 4326)

# 查看结果
print(rectangle_wgs84)

代码说明 #

• 投影转换：UTM坐标系以米为单位，距离计算更准确，避免经纬度近似误差；
• 方位角计算：用atan2(dy, dx)得到两点连线的弧度角度，用于旋转矩阵；
• 旋转变换：通过旋转矩阵将轴对齐的矩形顶点旋转到目标方位角，再平移到中心点位置；
• 坐标系转回：最后将矩形转换回原始的WGS84地理坐标系，方便后续使用。

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