CODEWITHSUNDEEP
rplotgeometrygisr-sf
adl
adl

Reputation: 1441

Make many circles with st_buffer at multiple geographic locations

I'm trying to create circles from a dataframe containing longitudes and latitudes, taking into consideration that the radius of every circle needs to be 400 Nautical Miles precisely. There are many points (approx. 8000) corresponding to aiports all over the world. I have been using the sf package to create circles with st_buffer, but this function takes distance argument in degrees. I`ve seen resolved questions that do not solve my problem completely since they involve single circle at a specific place for which a grid can be used. Here is the code I use: 

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)

I need these circles as sf objects in a dataframe because I will use them as polygons to find intersections between them and spatial lines (also sf) in other dataframe.

Results below - three circles are different size, two of them are distorted, one is missing:

enter image description here

Upvotes: 1

Views: 1680

Answers (1)

sebdalgarno
sebdalgarno

Reputation: 3197

you could transform each point according to it's respective utm zone and then get a buffer for each one separately. First, a function to find the utm zone proj4string from lat/lon for each point (requires purrr package):

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
}

now, convert each to utm and get buffer (this should be done on the df with all crs = 4326, not after transforming to Irish 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()

enter image description here

Upvotes: 3

Related Questions