When a web map is no good

Sometimes a web map is not a viable option. We have a number of static map options available. Before looking at them, I will switch back to static mapping:

tmap_mode("plot")
## tmap mode set to plotting

Bounding box

We can limit the extent of a static map with a bbox option in the first layer. By default, tmap uses the extent of the first layer to set the bounding box. So

tm_shape(welly) + 
  tm_polygons() +
  tm_shape(abb) + 
  tm_dots()

makes a map that only extends as far as that first layer. If we want it to extend to the full extent of the abb data then we can do

tm_shape(welly, bbox = abb) +
  tm_polygons() +
  tm_shape(abb) + 
  tm_dots()

We’ll worry about the missing ‘base’ layer later!

We can also use a bounding box to zoom in on an area. For example, if we make a bounding box that is from only the Lower Hutt City TA:

lower_hutt <- welly %>%
  filter(ta_name == "Lower Hutt City")
tm_shape(welly, bbox = lower_hutt) +
  tm_polygons() +
  tm_shape(abb) + 
  tm_dots()

Again, there is missing basemap, but we’ll get to that in a bit.

Using a bounding box for a detailed local map

If we want to make a specific bounding box, then we need to figure out the relevant coordinates. For this purpose, it is useful to make a temporary map with a grid.

tm_shape(welly) +
  tm_polygons() +
  tm_grid()

Now say we wanted a map showing only the central city. From the grid, we can see that this has an extent from 1,745,000 to 1,750,000 EW and 5,425,000 to 5,430,000 NS. We can use this information to make a bounding box using the bb() function from tmaptools:

cbd <- bb(xlim = c(1745000, 1750000), ylim = c(5425000, 5430000))

and now we can do

tm_shape(welly, bbox = cbd) +
  tm_polygons() +
  tm_shape(abb) + 
  tm_dots()

There are other options in the bb function, that allow you to specify it by using cx and cy (x and y centre coordinates) and width and height. The units are always the map projection units, which here are metres, so for example another way to make the above map would be

cbd <- bb(cx = 1747500, cy = 5427500, width = 5000, height = 5000)
tm_shape(welly, bbox = cbd) +
  tm_polygons() +
  tm_shape(abb) + 
  tm_dots()

This approach is a bit easier to recalculate if you want to zoom in closer or zoom out a bit, by just changing the width and height of the bounding box. There is also an ext option in the bb function to increase or shrink the bounding box relative to the requested width and height.

But what I really want is a basemap

Often what we really want a web map for is the basemap.

We can get basemaps in static maps using the get_tiles() function in the package maptiles and including it as a tm_rgb layer in the map. Here’s how we use it. First we grab the map tiles:

library(maptiles) # you may have to install this
basemap <- get_tiles(abb)

Then we make a map

tm_shape(basemap) + 
  tm_rgb() + 
  tm_shape(abb) + 
  tm_dots()

If the basemap is a bit blurry, then you can zoom in a bit. You have to be careful about this, as if you increase the zoom too much it will take ages to download (and you won’t be able to read the detail anyway). By experiment zoom levels in the range 9 to 11 are about right for the Wellington region, and get_tiles() is smart enough to get it about right most of the time, without you specifying a zoom level.

The only other thing to watch out for is that get_tiles() delivers whole web map tiles, so the cropping might seem a bit strange. You can fix this one of two ways, either: using bbox in tm_shape() as we have seen already:

basemap <- get_tiles(abb)
tm_shape(basemap, bbox = abb) + 
  tm_rgb() + 
  tm_shape(abb) + 
  tm_dots()

or, by setting the crop = TRUE parameter in the get_tiles() function:

basemap <- get_tiles(abb, crop = TRUE)
tm_shape(basemap) + 
  tm_rgb() + 
  tm_shape(abb) + 
  tm_dots()

Adjusting the zoom

Both the above examples suffer from a problem that get_tiles() has downloaded at lower zoom level than we’d like (zoom level 9 as it turns out), resulting in blurry labels. If we specify a higher zoom level, we get a finer resolution, although the download will take a bit longer:

basemap <- get_tiles(abb, crop = TRUE, zoom = 10)
tm_shape(basemap) + 
  tm_rgb() + 
  tm_shape(abb) + 
  tm_dots()

Removing that weird white margin

For some reason, get_tiles() adds a white margin around the downloaded image. You can get rid of this by buffering the crop area. The buffer distance is in the units of the data set (here metres):

basemap <- get_tiles(abb %>% st_buffer(2500), crop = TRUE, zoom = 10)
tm_shape(basemap, bbox = abb) + 
  tm_rgb() + 
  tm_shape(abb) + 
  tm_dots()

Some other refinements

Desaturate the basemap

If the colours of the basemap are distracting, you can wash them out with the saturation option in tm_rgb

tm_shape(basemap, bbox = abb) + 
  tm_rgb(saturation = 0.35) + 
  tm_shape(abb) + 
  tm_dots()

Other basemaps

Other basemaps are available than OpenStreetMap. You can see the options in the get_tiles() function help. For example, this is a nice one (if a bit distracting):

basemap <- get_tiles(abb, provider = "Stamen.TerrainBackground")
tm_shape(basemap, bbox = abb) + 
  tm_rgb() + 
  tm_shape(abb) + 
  tm_dots()

Some providers have labels-only tiles available so you can do this kind of thing:

labels <- get_tiles(abb, provider = "Stamen.TerrainLabels")
tm_shape(basemap, bbox = abb) + 
  tm_rgb() + 
  tm_shape(abb) + 
  tm_dots() +
  tm_shape(labels) + 
  tm_rgb()

Partial transparency of overlaid polygons

Remember that if you are overlaying polygons on a basemap you will probably want to set alpha (the opacity) to less than 1. So here is a final example.

basemap <- get_tiles(welly, crop = TRUE)
tm_shape(basemap, bbox = welly) +
  tm_rgb(saturation = 0.3) +
  tm_shape(welly) + 
  tm_polygons(col = "pop", palette = "viridis", alpha = 0.5) +
  tm_layout(legend.outside = TRUE)

Finally… here’s a little extra. The cbd bounding box from up above, might be something we want to use in a map like this. A bb object cannot be used directly to specify the limits of a web map, so we need to make it into a proper sf dataset for that, which you can do like this:

cbd_sf <- cbd %>%
  st_as_sfc() %>%
  st_set_crs(2193) # this is the known projection in this case
basemap <- get_tiles(cbd_sf, zoom = 13) # I experimented to get best zoom
tm_shape(basemap, bbox = cbd) +
  tm_rgb(saturation = 0.3) +
  tm_shape(welly) + 
  tm_polygons(col = "pop", palette = "viridis", alpha = 0.5) +
  tm_shape(abb) + 
  tm_dots() +
  tm_layout(legend.outside = TRUE)