DEV Community: Jakub Nowosad

Inset maps with ggplot2

Jakub Nowosad — Sun, 08 Dec 2019 00:00:00 +0000

Inset maps enable multiple places to be shown in the same geographic data visualisation, as described in the Inset maps section (8.2.7) of our open source book Geocomputation with R. The topic of inset maps has gained attention and recently Enrico Spinielli asked inset maps could be created for data in unusual coordinate systems:

Speaking of insets, do you know of any ggplot2 examples with an inset for showing where the (bbox of the) map is in an orthographic/satellite proj? It is usually done to provide geographic context... NYT/WP have fantastic maps like that

— espinielli (@espinielli) November 4, 2019

R’s flexibility allows inset maps to be created in various ways, using different approaches and packages. However, the main idea stays the same: we need to create at least two maps: a larger one, called the main map, that shows the central story and a smaller one, called the inset map, that puts the main map in context.

This blog post shows how to create inset maps with ggplot2 for visualization. The approach also uses the sf package for spatial data reading and handling, cowplot to arrange inset maps, and rcartocolor for additional color palettes. To reproduce the results on your own computer, after installing them, these packages can be attached as follows:

library(sf)
library(ggplot2)
library(cowplot)
library(rcartocolor)

Basic inset map

Let’s start by creating a basic inset map.

Data preparation

The first step is to read and prepare the data we want to visualize. We use the us_states data from the spData package as the source of the inset map, and north_carolina from the sf package as the source of the main map.

library(spData)
data("us_states", package = "spData")
north_carolina = read_sf(system.file("shape/nc.shp", package = "sf"))

Both objects should have the same coordinate reference system (crs). Here, we use crs = 2163, which represents the US National Atlas Equal Area projection.

us_states_2163 = st_transform(us_states, crs = 2163)
north_carolina_2163 = st_transform(north_carolina, crs = 2163)

We also need to have the borders of the area we want to highlight (use in the main map). This can be done by extracting the bounding box of our north_carolina_2163 object.

north_carolina_2163_bb = st_as_sfc(st_bbox(north_carolina_2163))

Maps creation

The second step is to create both inset and main maps independently. The inset map should show the context (larger area) and highlight the area of interest.

ggm1 = ggplot() + 
  geom_sf(data = us_states_2163, fill = "white") + 
  geom_sf(data = north_carolina_2163_bb, fill = NA, color = "red", size = 1.2) +
  theme_void()

ggm1

The main map’s role is to tell the story. Here we show the number of births between 1974 and 1978 in the North Carolina counties (the BIR74 variable) using the Mint color palette from the rcartocolor palette. We also customize the legend position and size - this way, the legend is a part of the map, instead of being somewhere outside the map frame.

ggm2 = ggplot() + 
  geom_sf(data = north_carolina_2163, aes(fill = BIR74)) +
  scale_fill_carto_c(palette = "Mint") +
  theme_void() +
  theme(legend.position = c(0.4, 0.05),
        legend.direction = "horizontal",
        legend.key.width = unit(10, "mm"))

ggm2

Maps joining

The final step is to join two maps. This can be done using functions from the cowplot package. We create an empty ggplot layer using ggdraw(), fill it with out main map (draw_plot(ggm2)), and add an inset map by specifing its position and size:

gg_inset_map1 = ggdraw() +
  draw_plot(ggm2) +
  draw_plot(ggm1, x = 0.05, y = 0.65, width = 0.3, height = 0.3)

gg_inset_map1

The final map can be saved using the ggsave() function.

ggsave(filename = "01_gg_inset_map.png", 
       plot = gg_inset_map1,
       width = 8, 
       height = 4,
       dpi = 150)

Advanced inset map

Let’s expand the idea of the inset map in ggplot2 based on the previous example.

Data preparation

This map will use the US states borders (states()) as the source of the inset map and the Kentucky Senate legislative districts (state_legislative_districts()) as the main map.

library(tigris)
us_states = states(cb = FALSE, class = "sf")
ky_districts = state_legislative_districts("KY", house = "upper",
                                           cb = FALSE, class = "sf")

The states() function, in addition to the 50 states, also returns the District of Columbia, Puerto Rico, American Samoa, the Commonwealth of the Northern Mariana Islands, Guam, and the US Virgin Islands. For our purpose, we are interested in the continental 48 states and the District of Columbia only; therefore, we remove the rest of the divisions using subset().

us_states = subset(us_states, 
                   !NAME %in% c(
                     "United States Virgin Islands",
                     "Commonwealth of the Northern Mariana Islands",
                     "Guam",
                     "American Samoa",
                     "Puerto Rico",
                     "Alaska",
                     "Hawaii"
                   ))

The same as in the example above, we transform both objects to have the same projection.

ky_districts_2163 = st_transform(ky_districts, crs = 2163)
us_states_2163 = st_transform(us_states, crs = 2163)

We also extract the bounding box of the main object here. However, instead of using it directly, we add a buffer of 10,000 meters around it. This output will be handy in both inset and main maps.

ky_districts_2163_bb = st_as_sfc(st_bbox(ky_districts_2163))
ky_districts_2163_bb = st_buffer(ky_districts_2163_bb, dist = 10000)

The ky_districts_2163 object does not have any interesting variables to visualize, so we create some random values here. However, we could also join the districts’ data with another dataset in this step.

ky_districts_2163$values = runif(nrow(ky_districts_2163))

Map creation

The inset map should be as clear and simple as possible.

ggm3 = ggplot() + 
  geom_sf(data = us_states_2163, fill = "white", size = 0.2) + 
  geom_sf(data = ky_districts_2163_bb, fill = NA, color = "blue", size = 1.2) +
  theme_void()

ggm3

On the other hand, the main map looks better when we provide some additional context to our data. One of the ways to achieve it is to add the borders of the neighboring states.

Importantly, we also need to limit the extent of our main map to the range of the frame in the inset map. This can be done with the coord_sf() function.

ggm4 = ggplot() + 
  geom_sf(data = us_states_2163, fill = "#F5F5DC") +
  geom_sf(data = ky_districts_2163, aes(fill = values)) +
  scale_fill_carto_c(palette = "Sunset") +
  theme_void() +
  theme(legend.position = c(0.5, 0.07),
        legend.direction = "horizontal",
        legend.key.width = unit(10, "mm"),
        plot.background = element_rect(fill = "#BFD5E3")) +
  coord_sf(xlim = st_bbox(ky_districts_2163_bb)[c(1, 3)],
           ylim = st_bbox(ky_districts_2163_bb)[c(2, 4)])

ggm4

Finally, we draw two maps together, trying to find the best location and size for the inset map.

gg_inset_map2 = ggdraw() +
  draw_plot(ggm4) +
  draw_plot(ggm3, x = 0.02, y = 0.65, width = 0.35, height = 0.35)

gg_inset_map2

The final map can be saved using the ggsave() function.

ggsave(filename = "02_gg_inset_map.png", 
       plot = gg_inset_map2,
       width = 7.05, 
       height = 4,
       dpi = 150)

Summary

The above examples can be adjusted to any spatial data and location. It is also possible to put more context on the map, including adding main cities’ names, neighboring states’ names, and annotations (using geom_text(), geom_label()). The main map can also be enhanced with the north arrow and scale bar using the ggsn package.

As always with R, there are many possible options to create inset maps. You can find two examples of inset maps created using the tmap package in the Geocomputation with R book. The second example is a classic map of the United States, which consists of the contiguous United States, Hawaii, and Alaska. However, Hawaii and Alaska are displayed at different geographic scales than the main map there. This problem can also be solved in R, which you can see in the Making maps of the USA with R: alternative layout blogpost and the Alternative layout for maps of the United States repository.

The presented approaches also apply to other areas. For example, you can find three ways on how to create an inset map of Spain in the Alternative layout for maps of Spain repository. Other examples of inset maps with ggplot2 can be found in the Inset Maps vignette by Ryan Peek and the blog post Drawing beautiful maps programmatically with R, sf and ggplot2 by Mel Moreno and Mathieu Basille.

The decision which option to use depends on the expected map type preferred R packages, etc. Try different approaches on your own data and decide what works best for you!

Map coloring: the color scale styles available in the tmap package

Jakub Nowosad — Thu, 17 Oct 2019 00:00:00 +0000

This vignette builds on the making maps chapter of the Geocomputation with R book. Its goal is to demonstrate all possible map styles available in the tmap package.

Prerequisites

The examples below assume the following packages are attached:

library(spData) # example datasets
library(tmap) # map creation
library(sf) # spatial data reprojection

The world object containing a world map data from Natural Earth and information about countries’ names, regions, and subregions they belong to, areas, life expectancies, and populations. This object is in geographical coordinates using the WGS84 datum, however, for mapping purposes, the Mollweide projection is a better alternative (learn more in the modifying map projections section). The st_tranform function from the sf package allows for quick reprojection to the selected coordinate reference system (e.g., "+proj=moll" represents the Mollweide projection).

world_moll = st_transform(world, crs = "+proj=moll")

One color

Let’s start with the basics. To create a simple world map, we need to specify the data object (world_moll) inside the tm_shape() function, and the way we want to visualize it. The tmap package offers several visualisation possibilities for polygons, including tm_borders(), tm_fill(), and tm_polygons(). The last one draws the filled polygons with borders, where the fill color can be specified with the col argument:

tm_shape(world_moll) +
  tm_polygons(col = "lightblue")

The output is a map of world countries, where each country is filled with a light blue color.

Coloring of adjacent polygons

The col argument is very flexible, and its action depends on the value provided. In the previous example, we provided a single color value resulting in a map with one color. To create a map, where adjacent polygons do not get the same color, we need to provide a keyword "MAP_COLORS".

tm_shape(world_moll) +
  tm_polygons(col = "MAP_COLORS")

The default color can be changed using the palette argument - run the tmaptools::palette_explorer() function to see possible palettes’ names.

tm_shape(world_moll) +
  tm_polygons(col = "MAP_COLORS",
              palette = "Pastel1")

Additionally, in this case, it is possible to use the minimize argument, which triggers the internal algorithm to search for a minimal number of colors for visualization.

tm_shape(world_moll) +
  tm_polygons(col = "MAP_COLORS",
              minimize = TRUE)

The new map uses five colors. On a side note, in theory, no more than four colors are required to color the polygons of the map so that no two adjacent polygons have the same color (learn more about the four color map theorem on Wikipedia).

Categorical maps

The third use of the col argument is by providing the variable (column) name. In this case, the map will represent the given variable. By default, tmap behaves differently depending on the input variable type. For example, it will create a categorical map when the provided variable contains characters or factors. The tm_polygons(col = "subregion", style = "cat") code will be run automatically in this case.

tm_shape(world_moll) +
  tm_polygons(col = "subregion")+
  tm_layout(legend.outside = TRUE)

Discrete maps

Discrete maps represents continuous numerical variables using discrete class intervals. There are several ways to convert continuous variables to discrete ones implemented in tmap.

Pretty

When the variable provided as the col argument is numeric, tmap will use the "pretty" style as a default. In other words, it runs tm_polygons(col = "lifeExp", style = "pretty") invisibly to the user. This style rounds breaks into whole numbers where possible and spaces them evenly.

tm_shape(world_moll) +
  tm_polygons(col = "lifeExp",
              legend.hist = TRUE) +
  tm_layout(legend.outside = TRUE)

A histogram is added using legend.hist = TRUE in this and several next examples to show how the selected map style relates to the distribution of values.

It is possible to indicate a preferred number of classes using the n argument. Importantly, not every n is possible depending on the range of the values in the data.

tm_shape(world_moll) +
  tm_polygons(col = "lifeExp",
              legend.hist = TRUE,
              n = 4) +
  tm_layout(legend.outside = TRUE)

Fixed

The "jenks" style allows for a manual selection of the breaks in conjunction with the breaks argument.

tm_shape(world_moll) +
  tm_polygons(col = "lifeExp", 
              style = "fixed",
              breaks = c(45, 60, 75, 90),
              legend.hist = TRUE) +
  tm_layout(legend.outside = TRUE)

Additionally, the default labels can be overwritten using the labels argument.

tm_shape(world_moll) +
  tm_polygons(col = "lifeExp", 
              style = "fixed",
              breaks = c(45, 60, 75, 90),
              labels = c("low", "medium", "high"),
              legend.hist = TRUE) +
  tm_layout(legend.outside = TRUE)

Breaks based on the standard deviation value

The "sd" style calculates a standard deviation of a given variable, and next use this value as the break width.

tm_shape(world_moll) +
  tm_polygons(col = "lifeExp", 
              style = "sd",
              legend.hist = TRUE) +
  tm_layout(legend.outside = TRUE)

Fisher algorithm

The "fisher" style creates groups with maximalized homogeneity.¹

tm_shape(world_moll) +
  tm_polygons(col = "lifeExp",
              style = "fisher",
              legend.hist = TRUE) +
  tm_layout(legend.outside = TRUE)

Jenks natural breaks

The "jenks" style identifies groups of similar values in the data and maximizes the differences between categories.²

tm_shape(world_moll) +
  tm_polygons(col = "lifeExp",
              style = "jenks",
              legend.hist = TRUE) +
  tm_layout(legend.outside = TRUE)

Hierarchical clustering

In the "hclust" style, breaks are created using hierarchical clustering.³

tm_shape(world_moll) +
  tm_polygons(col = "lifeExp",
              style = "hclust",
              legend.hist = TRUE) +
  tm_layout(legend.outside = TRUE)

Bagged clustering

The "bclust" style uses the bclust function to generate the breaks using bagged clustering.⁴

tm_shape(world_moll) +
  tm_polygons(col = "lifeExp",
              style = "bclust",
              legend.hist = TRUE) +
  tm_layout(legend.outside = TRUE)

## Committee Member: 1(1) 2(1) 3(1) 4(1) 5(1) 6(1) 7(1) 8(1) 9(1) 10(1)
## Computing Hierarchical Clustering

k-means clustering

The "kmeans" style uses the kmeans function to generate the breaks.⁵

tm_shape(world_moll) +
  tm_polygons(col = "lifeExp", 
              style = "kmeans",
              legend.hist = TRUE) +
  tm_layout(legend.outside = TRUE)

Quantile breaks

The "quantile" style creates breaks with an equal number of features (polygons).

tm_shape(world_moll) +
  tm_polygons(col = "lifeExp", 
              style = "quantile",
              legend.hist = TRUE) +
  tm_layout(legend.outside = TRUE)

Equal breaks

The "equal" style divides input values into bins of equal range and is appropriate for variables with a uniform distribution. It is not recommended for variables with a skewed distribution as the resulting map may end-up having little color diversity.

tm_shape(world_moll) +
  tm_polygons(col = "lifeExp", 
              style = "equal",
              legend.hist = TRUE) +
  tm_layout(legend.outside = TRUE)

Learn more about the implementation of discrete scales in the classInt package’s documentation.

Continuous maps

The tmap package also allows for creating continuous maps.

Continuous

The "cont" style presents a large number of colors over the continuous color field.

tm_shape(world_moll) +
  tm_polygons(col = "lifeExp",
              style = "cont") +
  tm_layout(legend.outside = TRUE)

Order

The "order" style also presents a large number of colors over the continuous color field. However, this style is suited to visualize skewed distributions; notice that the values on the legend do not change linearly.

tm_shape(world_moll) +
  tm_polygons(col = "lifeExp",
              style = "order") +
  tm_layout(legend.outside = TRUE)

Logarithmic scales

The default numeric style, pretty, is easy to understand, but it is not proper for maps of variables with skewed distributions.

tm_shape(world_moll) +
  tm_polygons(col = "pop") +
  tm_layout(legend.outside = TRUE)

Another possible style, order works better in this case; however, it is not easy to interpret.

tm_shape(world_moll) +
  tm_polygons(col = "pop", 
              style = "order") +
  tm_layout(legend.outside = TRUE)

A better alternative, in this case, is to use a common logarithm (the logarithm to base 10) scale. The tmap package gives two possibilities in this case - "log10_pretty" and "log10". The "log10_pretty" style is a common logarithmic version of the regular pretty style.

tm_shape(world_moll) +
  tm_polygons(col = "pop", 
              style = "log10_pretty") +
  tm_layout(legend.outside = TRUE)

On the other hand, the "log10" style is a version of a continuous scale.

tm_shape(world_moll) +
  tm_polygons(col = "pop", 
              style = "log10") +
  tm_layout(legend.outside = TRUE)

Conclusions

Selecting a color scale style is not an easy task. It depends on the type of input variable and its distribution, but also the intended audience. Therefore, it is worth to spend some time and think about your readers (e.g., would they be able to understand the logarithmic scale or should you use the manual breaks instead?) and your data (e.g., how many breaks should there be to show different subgroups?). Now you know different color scale styles implemented in tmap , so let’s try using them for your own projects!

https://www.tandfonline.com/doi/abs/10.1080/01621459.1958.10501479↩
https://en.wikipedia.org/wiki/Jenks_natural_breaks_optimization↩
See the ?hclust documentation for more details.↩
See the ?bclust documentation for more details.↩
See the ?kmeans documentation for more details.↩

Grids and graticules in the tmap package

Jakub Nowosad — Wed, 04 Sep 2019 00:00:00 +0000

This vignette builds on the making maps chapter of the Geocomputation with R book. Its goal is to demonstrate how to set and modify grids and graticules in the tmap package.

Prerequisites

The examples below assume the following packages are attached:

library(spData) # example datasets
library(tmap) # map creation (>=2.3)
library(sf) # spatial data classes

Grids and graticules

The tmap package offers two ways to draws coordinate lines - tm_grid() and tm_graticules(). The role of tm_grid() is to represent the input data’s coordinates. For example, the nz object uses the New Zealand Transverse Mercator projection, with meters as its units.

tm_shape(nz) + 
  tm_polygons() +
  tm_grid()

tm_graticules() shows longitude lines (meridians) and latitude lines (parallels), with degrees as units (note the degree sign in the example below).

tm_shape(nz) + 
  tm_polygons() +
  tm_graticules()

Layers order

Both, tm_grid() and tm_graticules() could be placed above or below the main spatial data. Its position on the map depends on its place in the code. When tm_grid() or tm_graticules() are placed after the code drawing geometry (e.g. tm_polygons()), the grids or graticules are ploted on the top of the map. On the other hand, when tm_grid() or tm_graticules() are placed before the code drawing geometry (e.g. tm_polygons()), the grids or graticules are plotted behind the spatial data.

tm_shape(nz) +
  tm_graticules() + 
  tm_polygons()

Customization

Grids and graticules can be easily customized in tmap using several arguments. The first one, labels.inside.frame moves the labels inside the map grid (it is set to FALSE as the default).

tm_shape(nz) +
  tm_grid(labels.inside.frame = TRUE) + 
  tm_polygons()

The number of horizontal (x) and vertical (y) lines can be set using the n.x and n.y arguments. Importantly, tmap rounds coordinate values to equally spaced “round” values, so the number of actual labels may be slightly different than set with n.x and n.y.

tm_shape(nz) +
  tm_grid(n.x = 4, n.y = 3) + 
  tm_polygons()

By default, tm_grid() and tm_graticules() shows ticks and lines. They can be disabled using ticks = FALSE and lines = FALSE.

tm_shape(nz) +
  tm_grid(ticks = FALSE) +
  tm_polygons()

Especially, lines = FALSE could be useful when presenting raster data.

tm_shape(nz) +
  tm_grid(lines = FALSE) +
  tm_polygons()

It is also possible to customize tm_grid() and tm_graticules() apperance, for example by chaning the lines colors (col), width (lwd) or labels size (labels.size).

tm_shape(nz) +
  tm_grid(col = "red", lwd = 3, labels.size = 0.4) +
  tm_polygons()

The above examples uses tm_grid(), but the same arguments apply to the tm_graticules().

Layout settings

By default, tmap adds small inner margins between the presented data and the map frame. It works well in many cases, for example, see the map of New Zealand above. However, it does not look perfect for world maps.

tm_shape(world) + 
  tm_graticules() + 
  tm_polygons()

The way to fix this is to use the tm_layout() function and set its inner.margins argument to 0.

tm_shape(world) + 
  tm_graticules() + 
  tm_polygons() +
  tm_layout(inner.margins = 0)