9

I would like to plot a circle 110 NM (nautical miles) around Dublin airport using sf package. (Later on I will intersect via st_intersect that with flight position reports from ADS-B.)

I have defined a new unit for NM as follows:

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

Then defined Dublin airport coordinates:

# 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)

Hence defined the radius of the circle around the airport (in meters):

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

Now when I try st_buffer things are not working:

> 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.

If I try to pass a numeric value (203.72, these are km) as distance at least I get only a warning:

> 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.

But plotting it shows quite a too big circle

library(mapview)
mapview(dub110)

enter image description here

What are the units for dist I should enter in st_buffer? I read the documentation but didn't really find out what to do...

Any hints/helps really appreciated!

espinielli
  • 666
  • 1
  • 6
  • 19
  • You're using WGS84 for airport coordinates - with unit of degrees - and so you're asking for circle of radius 110°, not km. Convert to Irish Grid http://spatialreference.org/ref/epsg/29902/ which uses metres. – Phil Oct 12 '17 at 08:42

3 Answers3

9

Thanks to Phil and Jul the complete solution to the initial question is as follows:

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)

enter image description here

espinielli
  • 666
  • 1
  • 6
  • 19
5

Here's a pure sf answer if you prefer, but @Jul 's is perfectly serviceable.

Set up as your example:

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)

Then transform your coordinate to Irish Grid:

dub = st_transform(dub, 29902)

Create your buffer in metres around this point:

dub_buffer = st_buffer(dub, 110000)

Plot the result:

plot(dub_buffer)
plot(dub, add = TRUE)
Phil
  • 4,344
  • 2
  • 23
  • 33
  • Perfect, Thanks. Just to be complete I did add/modify slightly ```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) ``` – espinielli Oct 12 '17 at 09:20
2

As mentioned in Phil's, you need to transform your coordinates to a metres/'distance' projection rather than a degree-based projection.

I'm not familiar with sf, but with sp...

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

...before you run the buffer command should suffice. You may then want to convert the buffered object back to epsg:4326 for plotting/additional processing. e.g. 

dub110 <- spTransform(dub110,CRS("+init=epsg:4326"))
Jul
  • 1,129
  • 6
  • 12