World population growth through time

Scaled cartogram of the world population changes between 1800 and 2100.

Starting

Let’s start with the packages. If you are new to R, you may want to take a read of this first, which points to various resources for setting-up R for geographic data. When you have a recent R version and the appropriate packages installed (e.g. by executing devtools::install_github("geocompr/geocompkg")) the packages can be attached as follows:

library(sf)             # spatial data classes
library(rnaturalearth)  # world map data
library(readxl)         # reading excel files
library(dplyr)          # data manipulation
library(stringr)        # data manipulation
library(tidyr)          # data manipulation
library(purrr)          # data manipulation
library(cartogram)      # cartograms creation
library(gganimate)      # animation creation

Getting data

To create the maps of the world population for each year we will need two datasets – one containing spatial data of the world’s countries and one non-spatial with information about the annual population in the world’s countries. The first one can be easily downloaded from the Natural Earth website, for example using the rnaturalearth package:

world_map = ne_countries(returnclass = "sf")

The second one is available from the Gapminder foundation. Gapminder provides a dataset with population data for all countries and world regions from 1800 to 2100. We can download and read the dataset using the code below:

if(!dir.exists("data")) dir.create("data")
download.file("http://gapm.io/dl_pop", destfile = "data/pop1800_2100.xlsx")
world_pop = read_xlsx("data/pop1800_2100.xlsx", sheet = 7)

Cleaning

As always when working with multiple datasets – some data cleaning is necessary. Our world_map dataset has many columns irrelevant for cartograms creation and we do not need spatial data of Antarctica. We can also transform our data into a more appropriate projection².

world_map = world_map %>% 
  filter(str_detect(type, "country|Country")) %>% 
  filter(sovereignt != "Antarctica") %>% 
  select(sovereignt) %>% 
  st_transform(world_map, crs = "+proj=moll")

We need to have a common identifier to combine our spatial and non-spatial datasets, for example, names of the countries. However, there are inconsistencies between some of the names. We can fix it manually:

world_pop = world_pop %>% 
  mutate(sovereignt = name) %>% 
  mutate(sovereignt = replace(sovereignt, sovereignt == "Tanzania", "United Republic of Tanzania")) %>% 
  mutate(sovereignt = replace(sovereignt, sovereignt == "United States", "United States of America")) %>% 
  mutate(sovereignt = replace(sovereignt, sovereignt == "Congo, Dem. Rep.", "Democratic Republic of the Congo")) %>% 
  mutate(sovereignt = replace(sovereignt, sovereignt == "Bahamas", "The Bahamas")) %>% 
  mutate(sovereignt = replace(sovereignt, sovereignt == "Serbia", "Republic of Serbia")) %>% 
  mutate(sovereignt = replace(sovereignt, sovereignt == "Macedonia, FYR", "Macedonia")) %>% 
  mutate(sovereignt = replace(sovereignt, sovereignt == "Slovak Republic", "Slovakia")) %>% 
  mutate(sovereignt = replace(sovereignt, sovereignt == "Congo, Rep.", "Republic of Congo")) %>% 
  mutate(sovereignt = replace(sovereignt, sovereignt == "Kyrgyz Republic", "Kyrgyzstan")) %>% 
  mutate(sovereignt = replace(sovereignt, sovereignt == "Lao", "Laos")) %>% 
  mutate(sovereignt = replace(sovereignt, sovereignt == "Cote d'Ivoire", "Ivory Coast")) %>% 
  mutate(sovereignt = replace(sovereignt, sovereignt == "Timor-Leste", "East Timor")) %>% 
  mutate(sovereignt = replace(sovereignt, sovereignt == "Guinea-Bissau", "Guinea Bissau"))

Preparing

Now we can join our two datasets, remove missing values and unimportant columns. We also should transform the data into a long format³, which is accepted by the plotting packages (such as tmap or ggplot2).

world_data = left_join(world_map, world_pop, by = "sovereignt") %>% 
  na.omit() %>% 
  select(-geo, -name, -indicator) %>% 
  gather(key = "year", value = "population", `1800.0`:`2100.0`,
         convert = TRUE)

Additionally, we can calculate total global populations in each year:

world_data = world_data %>% 
  group_by(year) %>% 
  mutate(total_pop = sum(as.numeric(population), na.rm = TRUE)) %>% 
  mutate(title = paste0("Year: ", year, "\nTotal population (billions): ", 
                        round(total_pop / 1e9, 2))) %>% 
  ungroup()

Subsetting

Now, our data contains information about the world population for each year between 1800 and 2100. It is possible to use it to create cartograms, however, to reduce calculation time and simplification of the results we will only use data for every 25 years:

world_data = world_data %>% 
  filter(year %in% seq(1800, 2100, by = 25))

Transforming

With newly created data, we are able to create our cartograms. We need to split our data into independent annual datasets, create cartograms based on the population variable, and combine all of the cartograms back into one object.

world_data_carto = world_data %>% 
  split(.$year) %>%
  map(cartogram_cont, "population", maxSizeError = 2, threshold = 0.1) %>%
  do.call(rbind, .)

Scaling

Predicted global population is the largest for the year 2100 and consists of more than 11 billion human beings. We can resize a map for each year using this value as our reference point, where the rest of the maps would be proportionally smaller than the one for 2100. To scale the areas of the other maps we need to divide a total population in each year by the maximum value and compute the square root of the result. The last step of calculations is to resize the geometries (coordinates) by multiplying it with the scales values⁴:

world_data_scaled = world_data_carto %>% 
  mutate(scales = sqrt(total_pop / max(total_pop))) %>% 
  group_by(year) %>% 
  mutate(geometry = geometry * scales)

For more on the topic, we recommend checking out the Affine transformations section of Geocomputation with R.

Visualizing

Now when we have all of the pieces needed, the animated map can be created using either tmap⁵ or gganimate. It consists of three steps:

1. Creating maps for each year using a combination of the ggplot(), geom_sf(), and transition_states() functions.
2. Rendering an animation using animate(). It allows for changing the size of a resulting image, and the speed of animation – you can control it with any two combinations of nframes, fps, and duration.
3. Saving an animation with anim_save().

worlds_anim = ggplot() +
  geom_sf(data = world_data_scaled, aes(fill = population / 1000000), lwd = 0.1) +
  coord_sf(datum = NA) +
  scale_fill_viridis_c(name = "Population (mln)") +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 22),
    plot.background = element_rect(fill = "#f5f5f4", color = NA), 
    panel.background = element_rect(fill = "#f5f5f4", color = NA), 
    legend.background = element_rect(fill = "#f5f5f4", color = NA),
    legend.position = c(0.1, 0.2),
    legend.title = element_text(size = 18)
  ) +
  labs(title = "{closest_state}") +
  transition_states(title, transition_length = 2, state_length = 6) + 
  ease_aes("cubic-in-out")

worlds_animate = animate(worlds_anim, width = 1200, height = 550,
                         duration = 30, fps = 4)
anim_save("worlds_animate.gif", animation = worlds_animate)

1. You can learn more about it at https://doi.org/10.1016/j.jag.2018.09.013 or by reading its preprint at https://eartharxiv.org/k3rmn/.↩

2. You can read more about projections in the Reprojecting geographic data chapter of Geocomputation with R.↩

3. This process is called gathering.↩

4. Mutating of geometry columns is possible with the recent improvements in the sf package.↩

5. Visit the animated maps chapter of Geocomputation with R or my previous blog post to learn how to make animations with tmap.↩