R语言如何实现两地理焦点指定距离下的椭圆区域绘制

椭圆的核心几何定义就是平面内到两个焦点的距离之和等于固定常数的所有点的集合，完全匹配你要的「两点间总行程固定的可达范围」需求，不需要依赖专门的椭圆绘制包，基于平面几何参数直接生成轮廓转sf对象即可，全程可以和你现有sf工作流无缝衔接。

核心计算逻辑 #

所有计算放在平面投影坐标系下完成（你当前用的27700英国国家网格就很合适，单位为米，距离计算无变形），参数计算规则：

• 记两个地理点为椭圆的两个焦点，先计算两点直线间距、椭圆中心点位置
• 自定义的总行程长度必须大于两焦点直线间距，否则不存在符合条件的连续范围
• 长半轴长度 = 总行程 / 2
• 焦点到椭圆中心的距离 = 两焦点间距 / 2
• 短半轴长度 = 根号下(长半轴² - 焦点到中心距离²)
• 按两个焦点的连线方位角计算椭圆的旋转角度
• 按角度步长生成椭圆轮廓点，闭合为多边形即可

可直接运行的实现代码 #

你原来的双100km缓冲区并集和目标椭圆范围不是同一个概念：缓冲区并集只包含「到任意一个点距离≤100km」的位置，而椭圆包含所有「到两个点的距离之和≤200km」的位置，在两点连线的垂直方向覆盖范围更大。

library(tidyverse)
library(sf)

# 输入点位坐标与自定义参数
citylocations <- tibble::tribble(
  ~city,     ~lon,      ~lat,
  "London",  -0.1276, 51.5072,
  "Birmingham", -1.8904, 52.4862,
)
total_travel_dist <- 200000 # 总行程约束，单位：米，此处设为200km

# 转换为平面投影坐标系（英国国家网格EPSG:27700）
citysflocations <- st_as_sf(citylocations, coords = c("lon","lat" ), crs = 4326)
cityBNGsflocations <- st_transform(citysflocations, crs = 27700)
coords_mat <- st_coordinates(cityBNGsflocations)

# 计算椭圆几何参数
f1 <- coords_mat[1, ]
f2 <- coords_mat[2, ]
center <- (f1 + f2)/2 # 椭圆几何中心
d_foci <- sqrt(sum((f1 - f2)^2)) # 两焦点直线距离

# 参数合法性校验
if(total_travel_dist <= d_foci){
  stop("设定的总行程必须大于两点间的直线距离，否则无有效范围")
}
a <- total_travel_dist / 2 # 长半轴
c <- d_foci / 2 # 焦点到中心的距离
b <- sqrt(a^2 - c^2) # 短半轴
rot_angle <- atan2(f1[2] - f2[2], f1[1] - f2[1]) # 长轴旋转角度（弧度）

# 生成椭圆轮廓点（默认1度步长，精度足够常规可视化）
theta_seq <- seq(0, 2*pi, by = pi/180)
ellipse_points <- cbind(
  x = center[1] + a*cos(theta_seq)*cos(rot_angle) - b*sin(theta_seq)*sin(rot_angle),
  y = center[2] + a*cos(theta_seq)*sin(rot_angle) + b*sin(theta_seq)*cos(rot_angle)
)

# 转换为标准sf多边形，转回WGS84经纬度坐标系
ellipse_poly <- st_polygon(list(rbind(ellipse_points, ellipse_points[1,]))) |> 
  st_sfc(crs = 27700) |> 
  st_transform(4326)

# 生成原有双100km缓冲区并集用于对比
dat_circles <- st_buffer(cityBNGsflocations, dist = 100000)
join_circles <- st_union(dat_circles) |> st_transform(4326)

# 可视化对比
plot(ellipse_poly, col = rgb(0.2,0.4,0.8,0.3), border = "darkblue", lwd = 2)
plot(join_circles, col = rgb(0.8,0.4,0.2,0.3), border = "darkred", lwd = 1, add = T)
plot(st_geometry(citysflocations), pch = 21, bg = "red", cex = 1.5, add = T)
text(st_coordinates(citysflocations)[,1], st_coordinates(citysflocations)[,2],
     labels = citylocations$city, pos = 4, cex = 0.9)
legend("topright", 
       legend = c("200km总行程可达范围(椭圆)", "双100km缓冲区并集"),
       fill = c(rgb(0.2,0.4,0.8,0.3), rgb(0.8,0.4,0.2,0.3)),
       border = c("darkblue", "darkred"))

使用说明 #

• 如果需要更高的轮廓精度，把生成点的角度步长改小即可，比如将by = pi/180调整为by = pi/360就是0.5度间隔取点
• 更换研究区域时，只要替换为当地适合的平面投影坐标系（保证距离单位为米、投影变形小），核心计算代码不需要修改
• 生成的ellipse_poly是标准sf多边形对象，支持直接做空间相交、裁剪、导出shp/geojson等所有常规空间操作

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