用 st_buffer 围绕一个地理点

Posted 2023-02-16

【中文标题】用 st_buffer 围绕一个地理点

【英文标题】：circle around a geographic point with st_buffer

【发布时间】：2018-03-24 02:41:26

【问题描述】：

我想使用sf 包在都柏林机场周围绘制一个 110 NM（海里）的圆圈。 （稍后我将通过st_intersect 与来自 ADS-B 的飞行位置报告相交。）

我已经为 NM 定义了一个新的单位如下：

library(units)
library(tidyverse)
library(sf)
NM <- make_unit("NM")
install_conversion_constant("NM", "km", 1.852)

然后定义都柏林机场坐标：

# DUB/EIDW location, see 
# https://skyvector.com/airport/EIDW/Dublin-Airport
# Coordinates:
#   N53°25.28' / W6°16.20' (Degrees Decimal Minutes (DDM) format)
#   (-6.27, 53.421333) (lon/lat Decimal Degrees (DD))
# Elevation: 242.0 feet (MSL)
dub_lon <- -6.27
dub_lat <- 53.421333
dub_elv <- set_units(242.0, ft)

dub <- st_point( x = c(dub_lon, dub_lat, dub_elv), dim = "XYZ")
dub <- dub %>% st_sfc(crs = 4326)

因此定义了机场周围圆的半径（以米为单位）：

r110 <- set_units(110, NM) %>% set_units(km)

现在当我尝试st_buffer 时，一切都无法正常工作：

> r110 <- set_units(110, NM) %>% set_units(km)
Error: cannot convert km into °
In addition: Warning message:
In st_buffer.sfc(dub, dist = r110) :
  st_buffer does not correctly buffer longitude/latitude data, dist needs to be in decimal degrees.

如果我尝试传递一个数值（203.72，这些是公里）作为距离，至少我只会收到一个警告：

> dub110 <- st_buffer(dub, dist = 203.72)
Warning message:
In st_buffer.sfc(dub, dist = 203.72) :
  st_buffer does not correctly buffer longitude/latitude data, dist needs to be in decimal degrees.

但是绘制它显示了一个太大的圆圈

library(mapview)
mapview(dub110)

我应该在st_buffer 中输入dist 的单位是什么？ 我阅读了文档，但并没有真正找到该怎么做......

非常感谢任何提示/帮助！

【问题讨论】：

您使用 WGS84 作为机场坐标 - 以度为单位 - 因此您要求半径为 110° 的圆，而不是公里。转换为使用米的爱尔兰网格spatialreference.org/ref/epsg/29902。

【参考方案1】：

感谢Phil 和Jul 最初问题的完整解决方案如下：

library(units)
library(tidyverse)
library(sf)
library(mapview)
library(units)

# define nautical miles (as per ICAO notation)
NM <- make_unit("NM")
install_conversion_constant("NM", "km", 1.852)

# DUB/EIDW location, see
# https://skyvector.com/airport/EIDW/Dublin-Airport
# Coordinates:
#   N53°25.28' / W6°16.20' (Degrees Decimal Minutes (DDM) format)
#   (-6.27, 53.421333) (lon/lat Decimal Degrees (DD))
# Elevation: 242.0 feet (MSL)
dub_lon <- -6.27
dub_lat <- 53.421333
dub_elv <- set_units(242.0, ft)
dub <- st_point(x = c(dub_lon, dub_lat, dub_elv), dim = "XYZ")
dub <- dub %>% st_sfc(crs = 4326)

# define radious of interest, i.e. 110 NM
r110 <- set_units(110, NM) %>% set_units(km) %>% set_units(m)

# change to Irish grid, which uses meters
dub <- st_transform(dub, 29902)
dub_buffer <-  st_buffer(dub, r110)

# eventually convert back to WSG84 if needed for other purposes
dub <- st_transform(dub, 4326)
dub_buffer <- st_transform(dub_buffer, 4326)
mapview(dub_buffer)

【参考方案2】：

如果您愿意，这里有一个纯粹的 sf 答案，但 @Jul 的答案非常有用。

设置为您的示例：

library(units)
library(tidyverse)
library(sf)
NM <- make_unit("NM")
install_conversion_constant("NM", "km", 1.852)

# DUB/EIDW location, see 
# https://skyvector.com/airport/EIDW/Dublin-Airport
# Coordinates:
#   N53°25.28' / W6°16.20' (Degrees Decimal Minutes (DDM) format)
#   (-6.27, 53.421333) (lon/lat Decimal Degrees (DD))
# Elevation: 242.0 feet (MSL)
dub_lon <- -6.27
dub_lat <- 53.421333
dub_elv <- set_units(242.0, ft)

dub <- st_point(x = c(dub_lon, dub_lat, dub_elv), dim = "XYZ")
dub <- dub %>% st_sfc(crs = 4326)

然后将你的坐标转换为Irish Grid：

dub = st_transform(dub, 29902)

围绕这一点创建以米为单位的缓冲区：

dub_buffer = st_buffer(dub, 110000)

绘制结果：

plot(dub_buffer)
plot(dub, add = TRUE)

【讨论】：

完美，谢谢。为了完整起见，我确实添加/修改了一点library(units) NM <- make_unit("NM") install_conversion_constant("NM", "km", 1.852) r110 <- set_units(110, NM) %>% set_units(km) %>% set_units(m) dub_buffer <- st_buffer(dub, r110) library(mapview) mapview(dub_buffer)

【参考方案3】：

如 Phil 所述，您需要将坐标转换为米/“距离”投影，而不是基于度数的投影。

我不熟悉sf，但是sp...

library(sp) 
dub_transformed <- spTransform(dub,CRS("+init=epsg:29902")) 

...在运行缓冲区命令之前就足够了。 然后，您可能希望将缓冲对象转换回 epsg:4326 以进行绘图/附加处理。例如

dub110 <- spTransform(dub110,CRS("+init=epsg:4326")) 

【讨论】：

非常感谢您的精彩回复！我授予菲尔作为第一个暗示解决方案的人......无论如何感谢一百万。