Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Global population density
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")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.