在多个地理位置用 st_buffer 圈出许多圈子

Make many circles with st_buffer at multiple geographic locations

我正在尝试根据包含经度和纬度的数据框创建圆,考虑到每个圆的半径需要精确到 400 海里。有很多点(约8000个)对应世界各地的机场。我一直在使用 sf 包来创建带有 st_buffer 的圆,但此函数采用以度为单位的距离参数。我看到已解决的问题并不能完全解决我的问题,因为它们涉及在可以使用网格的特定位置的单个圆圈。这是我使用的代码:

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)

# Sample data: 
df <- data.frame(lon = c(45,47,1, -109), lat = c(7, 10, 59, 30))

# Creating simple features with sf:
df <- df %>% st_as_sf(coords = c("lon", "lat"), dim = "XY")

# Applying Coordinate reference system WGS84:
df <- df %>% st_set_crs(4326)

# Transform to Irish grid - I know this step should be different for 
# different parts of the world but I don`t know how to make universal       
# solution
df <- st_transform(df$geometry, 29902)

# define radius of interest which is 400 NM
rad <- set_units(400, NM) %>% set_units(km) %>% set_units(m)

# make circles
df_buffer <- st_buffer(df, rad)

# visualise using mapview
mapview(df_buffer)

我需要这些圆圈作为数据框中的 sf 对象,因为我会将它们用作多边形来查找它们与其他数据框中的空间线(也称为 sf)之间的交点。

下面的结果 - 三个圆圈大小不一,其中两个变形,一个丢失:

您可以根据各自的 utm 区域转换每个点,然后分别为每个点获取一个缓冲区。 首先,一个从 lat/lon 中为每个点查找 utm zone proj4string 的函数(需要 purrr 包):

library(purrr)
utm_prj4 <- function(x) {

  coords <- cbind(x, st_coordinates(x))

  long <- coords$X
  lat <- coords$Y

  zone <- if(lat >= 56 && lat < 64 && long >= 3 && long < 12){x <- 32} else if(
    lat >= 72 && lat < 84 && long >= 0 && long < 9) {x <- 31} else if(
      lat >= 72 && lat < 84 && long >= 9 && long < 21) {x <- 33} else if(
        lat >= 72 && lat < 84 && long >= 21 && long < 33) {x <- 35} else if(
          lat >= 72 && lat < 84 && long >= 33 && long < 42) {x <- 37} else{
            x <- (floor((long + 180)/6) %% 60) + 1
          }
  prj <- purrr::map2_chr(zone, lat, function(y, z){
    if (z >= 0){
      paste0("+proj=utm +zone=", y, " +datum=WGS84 +units=m +no_defs")
    } else{
      paste0("+proj=utm +zone=", y, " +south", " +datum=WGS84 +units=m +no_defs")
    }})
  prj
}

现在,将每个转换为 utm 并获取缓冲区(这应该在所有 crs = 4326 的 df 上完成,而不是在转换为爱尔兰 crs 之后):

# creates a list of data.frames, each with different crs

dfs <- map2(1:4, utm_prj4(df), function(x, y){
  st_transform(df[x,], y)
})
map(dfs, ~ st_buffer(., rad)) %>% mapview()