Skip to content

Instantly share code, notes, and snippets.

@clauswilke
Created November 24, 2018 07:17
Show Gist options
  • Save clauswilke/c288b94ebc44a41437f8dce874508051 to your computer and use it in GitHub Desktop.
Save clauswilke/c288b94ebc44a41437f8dce874508051 to your computer and use it in GitHub Desktop.
library(tidyverse)
library(sf)
# Development version of colorspace is needed. Install via:
# install.packages("colorspace", repos = "http://R-Forge.R-project.org")
library(colorspace)
crs_goode <- "+proj=igh"
world_goode <- st_as_sf(rworldmap::getMap(resolution = "low")) %>%
st_transform(crs = crs_goode) %>%
mutate(
area = st_area(geometry), # area in m^2
popdens = POP_EST/area,
logdens = log(as.numeric(popdens)*1e6)
)
# projection outline in long-lat coordinates
lats <- c(
90:-90, # right side down
-90:0, 0:-90, # third cut bottom
-90:0, 0:-90, # second cut bottom
-90:0, 0:-90, # first cut bottom
-90:90, # left side up
90:0, 0:90, # cut top
90 # close
)
longs <- c(
rep(180, 181), # right side down
rep(c(80.01, 79.99), each = 91), # third cut bottom
rep(c(-19.99, -20.01), each = 91), # second cut bottom
rep(c(-99.99, -100.01), each = 91), # first cut bottom
rep(-180, 181), # left side up
rep(c(-40.01, -39.99), each = 91), # cut top
180 # close
)
goode_outline <-
list(cbind(longs, lats)) %>%
st_polygon() %>%
st_sfc(
crs = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
) %>%
st_transform(crs = crs_goode)
# bounding box in transformed coordinates
xlim <- c(-21945470, 21963330)
ylim <- c(-9538022, 9266738)
goode_bbox <-
list(
cbind(
c(xlim[1], xlim[2], xlim[2], xlim[1], xlim[1]),
c(ylim[1], ylim[1], ylim[2], ylim[2], ylim[1])
)
) %>%
st_polygon() %>%
st_sfc(crs = crs_goode)
# area outside the earth outline
goode_without <- st_difference(goode_bbox, goode_outline)
limit_range <- function(x, range = c(0, 1)) {
x[x < range[1]] <- range[1]
x[x > range[2]] <- range[2]
x
}
ggplot(world_goode, aes(fill = logdens)) +
geom_sf(color = "black", size = 0.5/.pt) +
geom_sf(data = goode_without, fill = "white", color = NA) +
geom_sf(data = goode_outline, fill = NA, color = "grey30", size = 0.5/.pt) +
scale_x_continuous(name = NULL) +
scale_y_continuous(name = NULL) +
scale_fill_viridis_c(option = "B", oob = limit_range, limits = c(-1, 10)) +
scale_fill_continuous_sequential(
palette = "Lajolla", rev = TRUE, p1 = 2, p2 = 1.3,
name = "population density\n(persons / square km)",
limits = c(-1, 8),
breaks = log(c(1, 10, 100, 1000, 10000)),
labels = parse(text = c("10^0", "10^1", "10^2", "10^3", "10^4")),
guide = guide_colorbar(
title.position = "top",
title.hjust = 0.5,
frame.colour = "black",
ticks.colour = "white",
barwidth = grid::unit(120, "pt"),
barheight = grid::unit(20, "pt")
)
) +
coord_sf(xlim = 0.95*xlim, ylim = 0.95*ylim, expand = FALSE, ndiscr = 1000) +
ggtitle("World-wide population density") +
theme_minimal() +
theme(
panel.background = element_rect(fill = "#56B4E950", color = "white", size = 1),
panel.grid.major = element_line(color = "gray30", size = 0.25),
legend.position = "bottom"
)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment