Global population density

population_density_map.R

library(rgdal)
library(maps)
library(viridis)
library(extrafont)
font_import(pattern = "GIL", prompt = FALSE)  # Import Gill family
loadfonts(device="win")  # Load them all
fonts()  # See what fonts are available
library(tidyverse)
#devtools::install_github("hadley/ggplot2", force = TRUE)
1

# Read population data, in grid form
pop <- read_csv("SEDAC_POP_2000-01-01_rgb_3600x1800.SS.CSV")

# Make it tall
longPop <- gather(pop[, -1])
longPop$Lat <- rep(pop$`lat/lon`, ncol(pop)-1)

# Omit empty coordinates
longPop <- longPop %>%
  filter(value != 99999) %>%
  mutate(Lon = as.numeric(Lat),
         Lat = as.numeric(key))

# Map projection
places_robin_df <- project(cbind(longPop$Lat, longPop$Lon), proj="+init=ESRI:54030") 
places_robin_df <- project(cbind(longPop$Lat, longPop$Lon), proj="+proj=wintri") 
longPop$transLat <- places_robin_df[, 1]
longPop$transLon <- places_robin_df[, 2]

# Marginal totals
latMargin <- longPop %>%
  group_by(Lat) %>%
  dplyr::summarise(Pop = sum(value))
lonMargin <- longPop %>%
  group_by(Lon) %>%
  dplyr::summarise(Pop = sum(value))

# Population by Lat, Lon, and Lat/Lon
zp1 <- ggplot(longPop) +
  geom_tile(aes(x = Lat, y = Lon, fill = value),
            colour = "transparent") +
  theme_void() +
  geom_rug(data = lonMargin,
           aes(y = Lon, colour = Pop), sides = "l") +
  geom_rug(data = latMargin,
           aes(x = Lat, colour = Pop), sides = "b") +
  coord_equal() +
  annotate(geom = "text", x = 125, y = -78, label = "@dsparks",
           family = "Gill Sans MT", size = 7,
           colour = gray(1), alpha = 1/2) +
  scale_fill_gradientn(guide = "none",
                       "Pop/Sq Mi",
                       breaks = 10 ^ (1:6),
                       colours = inferno(100), trans = "log") +
  scale_colour_gradientn(guide = "none",
                         "Aggregated by\nLat/Lon",
                         breaks = 10 ^ (1:6),
                         colours = inferno(100), trans = "log") +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  labs(title = "World population density",
       subtitle = "Logged, with marginal totals aggregated along latitude and longitude",
       caption = "NASA SEDAC, 2000") +
  theme(plot.title=element_text(family="Gill Sans MT", size = 27),
        plot.subtitle=element_text(family="Gill Sans MT", size = 18),
        plot.caption=element_text(family="Gill Sans MT", size = 15))
ggsave(plot = zp1, "Logged population density map.png", h = 10.80, w = 19.20, type = "cairo-png")


# A nicer projection
tX <- 6176370 * 0.9  # Apologies to Antarctica for the #branding
tY <- -8957861 # Apologies also for the hard-coded magic numbers

zp1 <- ggplot(longPop) +
  geom_point(aes(x = transLat, y = transLon, colour = value),
             fill = "transparent", shape = ".") +
  theme_void() +
  coord_equal() +
  annotate(geom = "text", x = tX, y = tY, label = "@dsparks",
           family = "Gill Sans MT", size = 7,
           colour = gray(1), alpha = 1/2) +
  scale_fill_gradientn(guide = "none",
                       "Pop/Sq Mi",
                       breaks = 10 ^ (1:6),
                       colours = inferno(100), trans = "log") +
  scale_colour_gradientn(guide = "none",
                         "Aggregated by\nLat/Lon",
                         breaks = 10 ^ (1:6),
                         colours = inferno(100), trans = "log") +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  labs(title = "World population density",
       caption = "NASA SEDAC, 2000") +
  theme(plot.title=element_text(family="Gill Sans MT", size = 27),
        plot.subtitle=element_text(family="Gill Sans MT", size = 18),
        plot.caption=element_text(family="Gill Sans MT", size = 15))
ggsave(plot = zp1, "Logged population density map (winkel tripel).png", h = 10.80, w = 19.20, type = "cairo-png")

# Find densest regions by sorting...
totalPop <- sum(longPop$value)
longPop <- longPop %>%
  arrange(desc(value)) %>%
  mutate(cumePop = cumsum(value),
         hiDensityRegion = cumePop / totalPop)

zp1 <- ggplot(longPop) +
  geom_point(aes(x = transLat, y = transLon, colour = hiDensityRegion <= 0.5),
             fill = "transparent", shape = ".") +
  annotate(geom = "text", x = tX, y = tY, label = "@dsparks",
           family = "Gill Sans MT", size = 7,
           colour = gray(1), alpha = 1/2) +
  theme_void() +
  coord_equal() +
  scale_colour_manual(guide = "none", values = inferno(7)[c(2, 6)]) +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  labs(title = "Highest population density regions",
       subtitle = "Half of the world's population lives in the yellow area.",
       caption = "NASA SEDAC, 2000") +
  theme(plot.title=element_text(family="Gill Sans MT", size = 27),
        plot.subtitle=element_text(family="Gill Sans MT", size = 18),
        plot.caption=element_text(family="Gill Sans MT", size = 15))
ggsave(plot = zp1, "Logged population density map (hi-density 50).png", h = 10.80, w = 19.20, type = "cairo-png")