Static maps with mapbox background images (snapbox)

maps Austria sf snapbox mapbox

How to create ggplot maps with mapbox images as background

true
11-25-2020

CONTEXT

Partly triggered by my last post on replicating NYT’s electoral maps, I came across the nascent/experimental snapbox package (link) by Anthony North and Miles McBain. The package provides a convenient way to add mapbox tiles as basemaps to maps created with ggplot (did I mention ‘maps’?`).

The package is quite straightforward and not much to explain about. What caused me some hassles though was to figure out how to combine the basemap with my data, plotted with geom_segment or geom_sf. The critical issue here is to ensure that the maps projection of the coordiante system are the same. Hence, this post is somewhat a note to myself, graphs are not polished and I don’t dig into substantive matters. For the purpose of demonstration, I’ll use electoral data from the city of Vienna and want to highlight how the number of eligible voters per precinct (‘Wahlsprengel’) changed between 2015 and 2020.

SETUP

Code: Load packages
library(snapbox)
library(sf)
library(tidyverse)
#library(magick)
knitr::knit_hooks$set(crop = knitr::hook_pdfcrop) #allows to crop white frame from ggplots.
#for details see https://bookdown.org/yihui/rmarkdown-cookbook/crop-plot.html
plot_bg_color <- readr::read_file(file=here::here("theme.css")) %>% 
  str_extract(., regex("(?<=background-color:).*?(?=;)")) %>% str_trim()

ELECTORAL DATA

Get the election results from 2015 and 2020, retrieve the number of eligible voters per precinct (`wber``), and calculate the change.

Code: Get electoral data + change in number of elegible voters
# https://www.data.gv.at/katalog/dataset/stadt-wien_wahleninwiengemeinderatswahlen/resource/a526df98-c307-4dee-8233-43b02856d1d7
# https://www.wien.gv.at/politik/wahlen/hilfetext/rohdaten.html
res20 <- readr::read_csv2(file = "https://www.wien.gv.at/politik/wahlen/ogd/gr201_99999999_9999_spr.csv")
res15 <- readr::read_csv2(file = "https://www.wien.gv.at/politik/wahlen/ogd/gr151_99999999_9999_spr.csv")

res <- bind_rows(elections=list('2020'=res20, '2015'=res15), .id="election") %>% 
  janitor::clean_names() %>% 
  janitor::remove_empty("cols") %>% 
  mutate(election=as.numeric(election)) %>% 
  mutate(indicator_sprengel=bz*1000+spr, .after=spr) %>% 
  select(election:wber) 


# change of eligible voters ----------------------------------------------------------

df_wber_diff <- res %>% 
  filter(wber>0) %>% #remove 'technical' precincts for mail & absentee voting 
  select(-contains("nuts")) %>% 
  group_by(indicator_sprengel) %>% 
  arrange(desc(election), .by_group=T) %>% 
  mutate(wber_diff=wber-lead(wber)) %>% 
  filter(election==2020) %>% 
  mutate(wber_diff_rel=wber_diff/wber) %>% 
  ungroup() %>% 
  mutate(change_indicator=case_when(wber_diff_rel>0 ~ "increase",
                                    wber_diff_rel<0 ~ "decrease",
                                    TRUE ~ as.character("stable")))
# A tibble: 1,483 x 4
   indicator_sprengel wber_diff wber_diff_rel change_indicator
                <dbl>     <dbl>         <dbl> <chr>           
 1               1001       -43       -0.0640 decrease        
 2               1002       -34       -0.0521 decrease        
 3               1003       -36       -0.0574 decrease        
 4               1004        -8       -0.0247 decrease        
 5               1005       -76       -0.104  decrease        
 6               1006       -43       -0.0814 decrease        
 7               1007       -74       -0.138  decrease        
 8               1008       -17       -0.0297 decrease        
 9               1009       -48       -0.0650 decrease        
10               1010       -31       -0.0442 decrease        
# ... with 1,473 more rows

Let’s have a quick look at the distribution of the change rates:

Code: Distribution of changes
pl_dist <- df_wber_diff %>% 
  ggplot()+
  ggridges::geom_density_line(aes(x=wber_diff_rel))+
  labs(x="% change number of eligible voters 2020 vs 2015")+
  theme(axis.text.y=element_blank())

table(df_wber_diff$change_indicator)

decrease increase   stable 
     978      451       54 

decrease increase   stable 
     978      451       54 

So, while the purpose of this post is not to dig into substantive matters, I was a bit surprised to see that a large majority of precincts features a decreased number of eligible voters.

MAP

Get boundaries of precincts

Get the shapefile of the precincts’ boundaries.

Code: Get shape with boundaries of precincts
# https://www.data.gv.at/katalog/dataset/79c1030d-5cf6-4d58-ade6-02f66fb4dffb
precincts <- st_read(here::here("_blog_data","static_map_snapbox","shapefile","WAHLSPRGR2020OGDPolygon.shp"),
              quiet=T) %>% 
  janitor::clean_names() %>% 
  mutate(sprengel=as.numeric(sprengel))

ggplot()+
  geom_sf(data=precincts)

Now let’s calculate the center of each precinct. Note that I transform the coordination projection from WGS84 to 3857. The latter is ‘a Spherical Mercator projection coordinate system popularized by web services such as Google and later OpenStreetMap.’1 After the transformation I take the x and y values from the centroids’ geometry. The function st_centroid gives us the centre of each polygon.

Code: Get centroids for precincts
st_crs(precincts)
Coordinate Reference System:
  User input: WGS 84 
  wkt:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ID["EPSG",6326]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8901]],
    CS[ellipsoidal,2],
        AXIS["geodetic longitude",east,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic latitude",north,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]]]
cent_precincts <- precincts %>% 
  st_transform(., crs=3857) %>% 
  st_centroid() %>% 
  mutate(point_x=map_dbl(geometry, pluck, 1),
         point_y=map_dbl(geometry, pluck, 2))

ggplot()+
  geom_sf(data=precincts)+
  geom_sf(data=cent_precincts,
          size=0.1)

Add change indicators (arrows)

Calculate the length of the arrow indicating the change: Here the arrows should have an angle of 20 degrees; positive changes point towards the right, negative changes towards the left. With the sin and cos function we calculate the horizontal and vertical distance from the centroids. In addition, to make the arrows sufficiently long to be visible, I multiply them with a scale factor. Note that I add coord_sf() to the plot which ‘ensures that all layers use a common CRS.’2

Code: Calculate length of arrows indicating change
scale_factor <- 10^3

df_change <- cent_precincts %>% 
  left_join(., df_wber_diff,
            by=c("sprengel"="indicator_sprengel")) %>% 
  mutate(x_change=wber_diff_rel*sin(20),
         y_change=abs(wber_diff_rel)*cos(20)) %>% 
  mutate(point_x_end=point_x+(x_change*scale_factor), .after=point_x) %>% 
  mutate(point_y_end=(y_change*scale_factor)+point_y, .after=point_y) 

pl_change <- ggplot()+
  geom_sf(data=precincts,
          size=0.1,
          fill="transparent")+
  geom_segment(data=df_change,
               aes(x=point_x,
                   xend=point_x_end,
                   y=point_y,
                   yend=point_y_end,
                   color=change_indicator),
               linejoin="mitre",
               arrow=arrow(length=unit(0.05, "cm"),
                           angle=15,
                           type="closed"),
               size=0.5)+
  scale_color_manual(values=c("increase"="darkgreen",
                              "decrease"="maroon",
                              "stable"="orange"))+
  coord_sf(crs=3857)+
  theme_void()

Add mapbox basemap

Now, let’s add a mapbox basemap by using the layer_mapbox function of snapbox. To so, we first have to define the scope or boundary of the basemap. I do this by extracting the minimum and maximum x and y values from our map of the precincts. And then create a bbox object. Note that this boundary box is has the 4326 projection. The basemap for this boundary box is then added to the plot by using the layer_mapbox function. Note that you need an API-key to access basemaps of mapbox.

Et voila.

Code: Add mapbox basemap
#get boundaries of plot

#map boundaries: get max/min x, y values
d <- precincts$geometry
d <- attributes(d)$bbox
my_xmin <- d[[1]]
my_ymin <- d[[2]]
my_xmax <- d[[3]]
my_ymax <- d[[4]]
my_xmin
[1] 16.18183
my_ymin
[1] 48.1179
my_xmax
[1] 16.57751
my_ymax
[1] 48.32267
#get basemap
background_image <- st_bbox(
  c(xmin = my_xmin, ymin = my_ymin, xmax = my_xmax, ymax = my_ymax),
  crs=4326) 

#combine
pl_combined <- ggplot()+
  layer_mapbox(background_image,
               map_style = mapbox_galley_pencil(),
               mapbox_api_access_token = Sys.getenv("MAPBOX_ACCESS_TOKEN"),
               scale_ratio = 1)+
  geom_sf(data=precincts,
          size=0.2,
          fill="transparent")+
  geom_segment(data=df_change,
               aes(x=point_x,
                   xend=point_x_end,
                   y=point_y,
                   yend=point_y_end,
                   color=change_indicator),
               linejoin="mitre",
               arrow=arrow(length=unit(0.05, "cm"),
                           angle=15,
                           type="closed"),
               size=0.5)+
  scale_color_manual(values=c("increase"="darkgreen",
                              "decrease"="maroon",
                              "stable"="orange"))+
  coord_sf(crs=3857)+
  theme_void()+
  theme(legend.position = "none")

You can quickly change the background image by specifying a different map_style value, e.g. below I use the decimal style. And also combine with a different aesthetic, e.g. using fill to color the precincts according to their respective change.

Code: Gradient fill of precincts
df_fill <- precincts %>% 
  left_join(., df_change %>% select(sprengel, wber_diff_rel) %>% as.data.frame(),
            by=c("sprengel"="sprengel"))


pl_combined2<- ggplot()+
  layer_mapbox(background_image,
               map_style = mapbox_gallery_north_star(),
               mapbox_api_access_token = Sys.getenv("MAPBOX_ACCESS_TOKEN"),
               scale_ratio = 1)+
  geom_sf(data=df_fill,
          aes(fill=wber_diff_rel),
          size=0.2)+

  scale_fill_gradient2(low="red", high="green", mid="white",
                       breaks=scales::pretty_breaks(n=4)(min(df_fill$wber_diff_rel, na.rm = T):
                                                   abs(min(df_fill$wber_diff_rel, na.rm=T))),
                                                 guide = "colourbar",
                       name="% change eligible voters 2015-2020",
                       limits=c(min(df_fill$wber_diff_rel, na.rm=T), abs(min(df_fill$wber_diff_rel, na.rm=T))))+

  coord_sf(crs=3857)+
  theme_void()+
  theme(plot.background = element_rect(fill=plot_bg_color, color=plot_bg_color),
        legend.title = element_text(size=7, vjust=1),
        legend.text = element_text(size=5),
        legend.position = "bottom",
        legend.justification = "left")+
  guides(fill = guide_colourbar(barwidth = 3, barheight = .2, title.position = "top"))


  1. See here.↩︎

  2. See here.↩︎

Reuse

Text and figures are licensed under Creative Commons Attribution CC BY-NC-SA 4.0. The figures that have been reused from other sources don't fall under this license and can be recognized by a note in their caption: "Figure from ...".

Citation

For attribution, please cite this work as

Schmidt (2020, Nov. 25). werk.statt.codes: Static maps with mapbox background images (snapbox). Retrieved from https://werk.statt.codes/post/2020-11-25-static-maps-with-mapbox-background-images-snapbox/

BibTeX citation

@misc{schmidt2020static,
  author = {Schmidt, Roland},
  title = {werk.statt.codes: Static maps with mapbox background images (snapbox)},
  url = {https://werk.statt.codes/post/2020-11-25-static-maps-with-mapbox-background-images-snapbox/},
  year = {2020}
}