Visualizing buffers that cross the dateline

Viewed 41

I'm trying to create and visualize buffers around point locations with the sf package in R. An initial attempt looked like this:

library(sf)
library(dplyr)
library(mapview)

sf_use_s2(TRUE)

coord <- c(178.4, -80.1)
point <- st_sfc(st_point(coord), crs = 4326)
buffer <-  st_buffer(point, 2000000, max_cells = 10000)
buffer %>%
  st_wrap_dateline(options = c("WRAPDATELINE=YES", "DATELINEOFFSET=180")) %>%
  mapview() + mapview(point)

enter image description here

I was able to fix this using st_shift_longitude() (sort of, latitude doesn't stretch to -90):

buffer %>%
  st_shift_longitude() %>%
  st_wrap_dateline(options = c("WRAPDATELINE=YES", "DATELINEOFFSET=180")) %>%
  mapview() + mapview(point)

enter image description here

However, this approach fails for other points:

coord <- c(78.4, -80.1)
point <- st_sfc(st_point(coord), crs = 4326)
buffer <-  st_buffer(point, 2000000, max_cells = 10000)
buffer %>%
  st_shift_longitude() %>%
  st_wrap_dateline(options = c("WRAPDATELINE=YES", "DATELINEOFFSET=180")) %>%
  mapview() + mapview(point)

enter image description here

Is there a surefire way to produce buffers like this?

1 Answers

If your coordinates and polygons are near the poles as above, mapview's default projections (4326 I think 3857 web mercator) probably won't work well. You can use other projections (with the native.crs = T argument), but you'll have to supply polygon data for the landmass as well. The default 'background' of Earth's landmasses won't automatically appear.

Below I've used the crs and antarctic polygon from a github issue thread found here: https://github.com/r-spatial/mapview/issues/298. You might be able to find some other tips in the thread as well.

library(mapview)
library(sf)
library(leaflet)
library(dplyr)

#Loading data 
# steal the crs & antarctic polygon (SFPoly2) for the points
load(url("https://github.com/elgabbas/Misc/blob/master/Data.RData?raw=true"))

# Your data
coord <- c(178.4, -80.1)
point <- st_sfc(st_point(coord), crs = 4326)
# Transform to use the polar crs
point <- st_transform(point, st_crs(SFPoint))
buffer <-  st_buffer(point, 2000000, max_cells = 10000)

# Use mapview with the 'native.crs = T' argument
#  There will be many warnings about old-style crs & not using long-lat data
mapview(SFPoly2, native.crs = T) + # Antarctic landmass
  mapview(point, fill = 'red', native.crs = T) + 
  mapview(buffer, native.crs = T)

enter image description here

Related