ggplot, ggsave & coord_map/quickmap: how to save large spatial objects and get the projection right?

Question

I have a largish polyline shapefile (Bavarian rivers, which can be accessed here) which I would like to plot and save via ggplot. This can easily be done via e.g. this code:

library(ggplot2)
library(rgdal)
library(sp)
library(rgeos)

riv <- readOGR(paste0(getwd(),"\\rivers_bavaria","rivers_bavaria"))
riv1 <- subset(riv,WDM=="1310"|WDM=="1320")
riv2 <- subset(riv,WDM=="1330")

p <- ggplot() +
  geom_line(data=riv1, aes(x=long, y=lat, group=group), color="dodgerblue", size=1) +
  geom_line(data=riv2, aes(x=long, y=lat, group=group), color="dodgerblue")

ggsave(paste0(getwd(),"\\riv.tiff",p,device="tiff",units="cm",dpi=300)

This is not exactly efficient, due to the large file size, but it works. However, without further specifying aspect ratio or projection, the dimensions of the output file are defined by the plot window - not desirable for maps. This can be remedied by using coord_quickmap().

p1 <- ggplot() +
  geom_line(data=riv1, aes(x=long, y=lat, group=group), color="dodgerblue", size=1) +
  geom_line(data=riv2, aes(x=long, y=lat, group=group), color="dodgerblue") +
  coord_quickmap()
p1

Unfortunately, the projection is completely off. I have tried coord_map() for a better result, but due to the large file size, it takes forever and is therefore not a realistic option. Simplifying the polyline via gLinemerge() produces a much smaller object, but cannot be handled by ggplot, as it is a SpatialLines object. Using fortify() or data.frame() to coerce it into a ggplot-friendly data frame format also produces Error: ggplot2 doesn't know how to deal with data of class SpatialLines.

I'm therefore desperately looking for a workflow that will allow me to plot and save this kind of spatial data in good quality with ggplot. Any suggestions will be much appreciated!

Answer 1

Here's a quick walkthrough with sf. I recommend the sf vignettes and docs to see more details of any of the functions. I'm first reading the shapefile in as an sf object using sf::st_read, then filtering, mutating, and selecting the same as you would in dplyr to get a smaller version of the shape.

library(tidyverse)
library(sf)

rivers_sf <- st_read("rivers_bavaria/rivers_bavaria.shp") %>%
  filter(WDM %in% c("1310", "1320", "1330")) %>%
  mutate(name2 = ifelse(WDM == "1330", "river 2", "river 1")) %>%
  select(name2, NAM, geometry)

The object is pretty big, and will be very slow to plot, so I simplified it by uniting the geometries by name, then using st_simplify. There's also rmapshaper::ms_simplify, which uses Mapshaper and which I prefer for better control over how much information you keep. Then to show a CRS transformation, I picked a projection from Spatial Reference for Germany.

riv_simple <- rivers_sf %>%
  group_by(name2, NAM) %>%
  summarise(geometry = st_union(geometry)) %>%
  ungroup() %>%
  st_simplify(preserveTopology = T, dTolerance = 1e6) %>%
  st_transform(31493)

The dev version of ggplot2 on GitHub has a function geom_sf for plotting different types of sf objects. To get this version, run devtools::install_github("tidyverse/ggplot2").

geom_sf has some quirks, and works a little differently from other geoms, but it's pretty versatile. I believe it's being included in the next CRAN release. geom_sf has corresponding stat_sf and coord_sf. By default, it plots graticule lines; to turn those off, add coord_sf(ndiscr = F).

ggplot(riv_simple) +
  geom_sf(aes(size = name2), color = "dodgerblue", show.legend = "line") +
  scale_size_manual(values = c("river 1" = 1, "river 2" = 0.5)) +
  theme_minimal() +
  coord_sf(ndiscr = F)

Hope that helps you get started!