Search code examples
rgoogle-mapsggplot2facet

R: Facet a map with ggplot2

I am trying to facet a map with ggplot2 because is faster than tmap package. Lamentably the map that I got is not that I have in mind.

library(ggplot2)
library(raster)
library(ggspatial)


chile  <- getData("GADM",country="Chile",level=1)

chile2= chile[c(2,4:5,7,8,12:16),]
chile2$grupo=1
chile3= chile[c(1,3,6,9:11),]
chile3$grupo=2
mapa=rbind(chile2, chile3)

ggplot() +
  layer_spatial(mapa) +
  lims(x = c( -77.1,-65), y = c(-57, -15))+
  facet_wrap(~grupo)

ggplot map

Lamentably the map above is not what I need. With tmap I got the map below and is the map that really I need:

tm map

Do you know how solve this using ggplot2?

Solution

  • This can be done, though I'm not sure if it's really worth the effort.

    There are two issues here:

    1. Usually, you can get different scales (latitudes / longitudes in this case) in each facet of facet_wrap by setting scales = "free", but layer_spatial only works with the coord_sf coordinate system, and coord_sf is hard-coded to only work with fixed scales (there's a discussion on the ggplot GH page that covers how this came about);

    2. Setting lims() explicitly defeats the purpose of free scales anyway.

    For the first issue, we can go against the package developers' intention to define an alternative version of the coordinate system that accepts free scales. (I'm not saying that it should be done, but it can be. Caveat emptor.)

    CoordSf2 <- ggproto("CoordSf2",
                    CoordSf,
                    is_free = function() TRUE)
trace(coord_sf, edit = TRUE)

    Running the trace(...) line will result in a pop-up window with the code for coord_sf. Modify the last section from ggproto(NULL, CoordSf, ...) to ggproto(NULL, CoordSf2, ...) will point coord_sf to our modified CoordSf2 instead of the original CoordSf. This effect will remain in place until the end of the current R session, or you can terminate it earlier by running untrace(coord_sf).

    For the second issue, I suppose the limits were set to show only polygons within the limits. We can perform this filtering step in the data frame, before passing it to ggplot:

    library(dplyr)

keep <- lapply(mapa@polygons, bbox) %>%                    # retrieve bounding box for each polygon
  lapply(function(x) c(x) <= c(-77.1, -57, -65, -15)) %>%  # compare each polygon's bbox against
  lapply(function(x) x == c(F, F, T, T)) %>%               # desired limits, if within, 
  sapply(all)                                              # results should be FFTT

keep # only the 10th polygon is outside the limits
 [1]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE
[11]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE

    Now plot with the subset of mapa polygons that fall within the desired limits, with facet scales set to "free" & coord_sf pointing to the modified version that allows it:

    ggplot() +
  layer_spatial(mapa[keep, ]) +
  facet_wrap(~grupo, scales = "free")

    result

    Personally, though? I'd probably just make separate plots & stitch them together, like the cowplot method demonstrated in this question.