Plotting around the +/- 180 dateline using sf and ggplot2

Viewed 732

I am using sf to plot a series of shapefiles and point data. This data spans the dateline at =/-180 longitude.

sf::st_shift_longitude() handles this situation and plots points as expected. However, ggplot is behaving strangely when assigning the longitude tick marks. Code below provides one example with two points on one side of dateline- note ggplot adds logical tick marks for longitude. In second example, points span dateline.

Edit: Added third case, with manually expanded x limits. When these get large enough, the graticule plots as expected. However, I only found 'large enough' by experimentation with xlim.

library(sf)
library(rnaturalearth)
library(ggplot2)

#get basic world map 
world = ne_coastline(scale = "medium", returnclass = "sf")

#example 1: without dateline
#minimal data- two points on one side of dateline
sites = data.frame(longitude = c(-173.9793, -177.7405), latitude = c(52.21415, 51.98994))

#convert to sf object
sites = st_as_sf(sites, coords = c("longitude", "latitude"), crs = 4326)

#plot with ggplot
ggplot()+
  geom_sf(data = world, fill = 'transparent')+
  geom_sf(data = sites)+
  #set the limits for the plot
  coord_sf(crs = 4326,
           xlim = c(min(st_coordinates(sites)[,1]) -1, max(st_coordinates(sites)[,1])+1),
           ylim = c(min(st_coordinates(sites)[,2]) -1, max(st_coordinates(sites)[,2])+1))+
  labs(title = 'data on one side of dateline only- looks good')+
  theme_bw()


#example 2: with dateline
#deal with dateline using st_shift_longitude 
world_2 = st_shift_longitude(world)

#minimal data- a point on each side of dateline
sites_2 = data.frame(longitude = c(-173.9793, 177.7405), latitude = c(52.21415, 51.98994))

#convert to sf object
sites_2 = st_as_sf(sites_2, coords = c("longitude", "latitude"), crs = 4326)
#and deal with dateline using st_shift_longitude 
sites_2 = st_shift_longitude(sites_2)

#plot with ggplot
ggplot()+
  geom_sf(data = world_2, fill = 'transparent')+
  geom_sf(data = sites_2)+
  #set the limits for the plot
  coord_sf(crs = 4326,
           xlim = c(min(st_coordinates(sites_2)[,1]) -1, max(st_coordinates(sites_2)[,1])+1),
           ylim = c(min(st_coordinates(sites_2)[,2]) -1, max(st_coordinates(sites_2)[,2])+1))+
  labs(title = 'data on both sides of dateline - grid wrong')+
  theme_bw()

#plot with manually expanded limits- graticule works
ggplot()+
  geom_sf(data = world_2, fill = 'transparent')+
  geom_sf(data = sites_2)+
  #set the limits for the plot
  coord_sf(crs = 4326,
           xlim = c(175, 195),
           ylim = c(min(st_coordinates(sites_2)[,2]) -1, max(st_coordinates(sites_2)[,2])+1))+
  labs(title = 'data on both sides of dateline - manually expand x lims')+
  theme_bw()



enter image description here enter image description here enter image description here

Do you have any ideas on working around the dateline to not get this behavior when plotting with ggplot?

Thanks!

1 Answers

I don't know enough of ggplot, but my guess is that internally the graticule is restricted to -180,180, and event though you shifted the sf object ggplot doesn't recognise that on creating the axis and the grid.

I created an example with base plot that seems to work as you expect. The grid is created with graphics::grid()



library(sf)
#> Warning: package 'sf' was built under R version 3.5.3
#> Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
library(rnaturalearth)
#> Warning: package 'rnaturalearth' was built under R version 3.5.3

#get basic world map
world = ne_coastline(scale = "medium", returnclass = "sf")

#example 1: without dateline
#minimal data- two points on one side of dateline
sites = data.frame(
  longitude = c(-173.9793,-177.7405),
  latitude = c(52.21415, 51.98994)
)

#convert to sf object
sites = st_as_sf(sites,
                 coords = c("longitude", "latitude"),
                 crs = 4326)

#deal with dateline using st_shift_longitude
world_2 = st_shift_longitude(world)

#minimal data- a point on each side of dateline
sites_2 = data.frame(
  longitude = c(-173.9793, 177.7405),
  latitude = c(52.21415, 51.98994)
)

#convert to sf object
sites_2 = st_as_sf(sites_2,
                   coords = c("longitude", "latitude"),
                   crs = 4326)
#and deal with dateline using st_shift_longitude
sites_2 = st_shift_longitude(sites_2)

#Plot
plot(
  st_geometry(world_2),
  axes = TRUE,
  xlim = c(min(st_coordinates(sites_2)[, 1]) - 1, max(st_coordinates(sites_2)[, 1]) +
             1),
  ylim = c(min(st_coordinates(sites_2)[, 2]) - 1, max(st_coordinates(sites_2)[, 2]) +
             1)
)
grid()
plot(st_geometry(sites_2), pch = 20, add = TRUE)

Created on 2020-03-28 by the reprex package (v0.3.0)

Related