Search code examples
rggplot2coordinatesspatialggspatial

How to use `ggplot2::geom_segment()` or `ggspatial::geom_spatial_segment()` for sf objects not centered on Greenwich?


I want to add paths into a map generated using geom_sf in ggplot2.

I am familiar with using the standard map centred on Greenwich. But when I change the map design(? crs), the usual methods fail.

The 'normal' map working fine

point_df=data.frame(
  location=c('Tokyo','Brasilia','France'),
  lat=c(35.65,-15.79,48.86),
  long=c(139.84,-47.88,2.35)
)
point_df2=data.frame(
  from_location=c('Tokyo','Brasilia'),
  from_lat=c(35.65,-15.79),
  from_long=c(139.84,-47.88),
  to_location=c('France'),
  to_lat=c(48.86),
  to_long=c(2.35)
)

library("rnaturalearth")
world <- ne_countries(scale = "medium", returnclass = "sf")
transpoint = st_as_sf(point_df,coords=c("long","lat"),crs=4326)
world_sf <- st_transform(world, crs=4326)

ggplot(world_sf)+geom_sf()+
  geom_point(data=transpoint,aes(geometry=geometry,color=location),
stat="sf_coordinates")+
  geom_text_repel(data=transpoint,aes(geometry=geometry,label=location),
stat="sf_coordinates")+
  geom_curve(data=point_df2,
             aes(x=from_long,y=from_lat,xend=to_long,yend=to_lat),
colour='blue')

Greenwich centered map

Not working: map with different CRS or not-centred in Greenwich

However, if I prefer a map with the Pacific Ocean at the center, or an elliptical map, they fail.

pacific centred map with sf::st_break_antimeridian()

For the pacific centred map, I used a code adapted from Closed boundaries on a map when using geom_sf to cross dateline Note the geom_point() works fine. The geom_segment() fails completely...


robinson <- "+proj=robin +lon_0=-90 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"

world_pac <- world_sf %>%
  st_break_antimeridian(lon_0 = -90) %>% 
  st_transform(crs = robinson)

transpoint = st_as_sf(point_df,coords=c("long","lat"),crs=4326)
dtran = st_transform(transpoint,robinson)
library(ggrepel)
ggplot(world_pac)+geom_sf()+
  geom_point(data=dtran,aes(geometry=geometry,color=location),stat="sf_coordinates")+
  geom_text_repel(data=dtran,aes(geometry=geometry,label=location),stat="sf_coordinates")+
  geom_segment(data=point_df2, aes(x=from_long,y=from_lat,xend=to_long,yend=to_lat),
                       colour='blue')


Pacific centered geom_segment

The ggspatial::geom_spatial_segment() show curves/segments, but they are weird.

ggplot(world_pac)+geom_sf()+
  geom_point(data=dtran,aes(geometry=geometry,color=location),stat="sf_coordinates")+
  geom_text_repel(data=dtran,aes(geometry=geometry,label=location),stat="sf_coordinates")+
  geom_spatial_segment(data=point_df2, aes(x=from_long,y=from_lat,xend=to_long,yend=to_lat),crs=4326,
                       wrap_dateline = FALSE,
                       colour='blue')

Pacific centered with geom_spatial_segment

Ellipsis map

I am also unable to make geom_segment() or geom_spatial_segment() work if I transform to crs=2163.


world_ellipse <- world_sf %>%
  st_transform(crs=2163)

transpoint = st_as_sf(point_df,coords=c("long","lat"),crs=4326)
dtran = st_transform(transpoint,crs=2163)

ggplot(world_ellipse)+geom_sf()+
  geom_point(data=dtran,aes(geometry=geometry,color=location),stat="sf_coordinates")+
  geom_text_repel(data=dtran,aes(geometry=geometry,label=location),stat="sf_coordinates")
geom_spatial_segment(data=point_df2,
                     aes(from_long, from_lat, xend = to_long, yend = to_lat),crs = 2163,
                     wrap_dateline = FALSE)


Ellipsis with geom_spatial_segment

What am I missing here? How to work with geom_segment() and/or geom_spatial_segment() with ggplot2::geom_sf(), when we use different CRS?


Solution

  • The easiest way to get around this is to generate an sf version of your points of interest, then create the from/to variables in your sf for the geom_curve() .

    I have dropped ggrepel from this example as it is like playing Whack-A-Mole when it comes to getting satisfactory results. I also adjusted the curvature of the lines as the default value bends the France/Tokyo line off the custom Robinson map. If these parameters don't suit, you can tweak them where necessary. The example plot contains an arrow option (commented out) in case you want to demonstrate direction.

    This works for both your custom CRS, and for EPSG:9311 (the updated version of EPSG:2163, which has been deprecated). The only modification you'll need to make to generate an EPSG:9311 version of your map is to replace st_transform(robinson) with st_transform(9311).

    Note also that there is an error in your ellipsis map code. You have omitted the "+" after the geom_text_repel() function. Perhaps a copy/paste error, but geom_spatial_segment() will still not plot point_df2. However, it will work with the sf_points object created in this solution.

    Here is a full working repex:

    library(rnaturalearth)
    library(sf)
    library(dplyr)
    library(ggplot2)
    
    # Your point data
    point_df <- data.frame(
      location = c("Tokyo", "Brasilia", "France"),
      lat = c(35.65,-15.79,48.86),
      long = c(139.84,-47.88,2.35))
    
    # Custom CRS
    robinson <- "+proj=robin +lon_0=-90 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
    
    # Create 90W-centred version of ne_countries with custiom CRS
    world_pac <- ne_countries(scale = "medium", returnclass = "sf") %>%
      st_break_antimeridian(lon_0 = -90) %>% 
      st_transform(robinson)
    
    # Create sf of your point data
    transpoint <- st_as_sf(point_df, coords = c("long","lat"), crs = 4326) %>%
      st_transform(robinson)
    
    # Generate from/to variables between France and the rest
    sf_points <- transpoint %>%
      st_set_geometry("to") %>%
      mutate(from = to[location == "France"]) %>%
      filter(!location == "France") %>%
      mutate(x1 = st_coordinates(from)[,1],
             y1 = st_coordinates(from)[,2],
             x2 = st_coordinates(to)[,1],
             y2 = st_coordinates(to)[,2])
    
    #Plot
    ggplot() +
      geom_sf(data = world_pac, colour = "grey70") +
      geom_sf(data = transpoint, colour = "#DF536B", size = 3) +
      geom_curve(data = sf_points, 
                 aes(x = x1, y = y1, 
                     xend = x2, yend = y2), 
                 colour = "#0072B2",
                 # arrow = arrow(length = unit(0.03, "npc"), type="closed"),
                 curvature = 0.4) +
      geom_sf_text(data = transpoint,
                   aes(label = location),
                   size = 4,
                   fun.geometry = st_centroid,
                   vjust = 2,
                   colour = "black") +
      theme(axis.title = element_blank())
    

    Custom CRS result: result

    EPSG:9311 result: result1