I have loved the graphics from "A Portrait of the City as a Squiggly Line" from the from the "We Alone on Earth" blog for the past 15 years (!), and have long aspired to recreate it. With some extensive help from Claude, I have done so (R Code at the bottom).
The original (based on 2000 census data)
My recreation (2023 ACS data):
I suspect that the creator of the original did some additional (visual) smoothing for the second chart.
R code follows:
#-------------------------------------------------------------------------------
# Population density experienced by residents, compared across eight 2020 Census Urban Areas:
# Los Angeles--Long Beach--Anaheim, CA San Francisco--Oakland, CA
# Chicago, IL--IN New York--Jersey City--Newark, NY--NJ
# Boston, MA--NH Seattle--Tacoma, WA
# Las Vegas--Henderson--Paradise, NV Atlanta, GA
#
# Universe: census tracts whose interior point falls inside the urban area
# boundary. Whole tracts are used (their full population and land area), so
# tracts straddling the edge are in or out as a unit.
#
# Density = tract population / tract LAND area (people per km²).
#
# Two charts:
# 1. Share of residents living at AT LEAST each density (log x-axis).
# Read: "x% of the urban area's people live in tracts at least this dense."
# 2. Number of people living in tracts within +/-15% of each density.
# Read: "this many people live at about this density."
#
# Both are computed directly from tracts -- no smoothing.
#
# Population: most recent ACS 5-year (tries newer vintages first, falls back).
# Land area: TIGER/Line cartographic boundary tracts (ALAND) for the same year.
# Urban areas: TIGER/Line 2020 Urban Areas (2020 criteria).
#
# Requires: tidycensus, tigris, sf, dplyr, ggplot2, scales, patchwork
# install.packages(c("tidycensus", "tigris", "sf", "dplyr", "ggplot2", "scales", "patchwork"))
# One-time API key setup (free key: https://api.census.gov/data/key_signup.html):
# tidycensus::census_api_key("YOUR_KEY", install = TRUE); readRenviron("~/.Renviron")
# ---------------------------------------------------------------------------
suppressPackageStartupMessages({
library(tidycensus)
library(tigris)
library(sf)
library(dplyr)
library(ggplot2)
library(scales)
})
options(tigris_use_cache = TRUE)
# --- settings --------------------------------------------------------------
# One row per urban area. `ua_name` must match the Census name exactly --
# note the double hyphen in multi-city names. If a name is off, the script
# stops and suggests close matches. `state` = the state(s) whose tracts to
# search; for a UA crossing a state line, add a second row with the same
# label and ua_name and the other state.
areas <- data.frame(
label = c("Los Angeles", "Bay Area",
"Chicago", "Chicago",
"New York", "New York",
"Boston", "Boston",
"Seattle", "Las Vegas", "Atlanta"),
ua_name = c("Los Angeles--Long Beach--Anaheim, CA", "San Francisco--Oakland, CA",
"Chicago, IL--IN", "Chicago, IL--IN",
"New York--Jersey City--Newark, NY--NJ", "New York--Jersey City--Newark, NY--NJ",
"Boston, MA--NH", "Boston, MA--NH",
"Seattle--Tacoma, WA", "Las Vegas--Henderson--Paradise, NV", "Atlanta, GA"),
state = c("CA", "CA",
"IL", "IN",
"NY", "NJ",
"MA", "NH",
"WA", "NV", "GA"),
stringsAsFactors = FALSE
)
ua_year <- 2023 # fallback TIGER vintage if the 2020-criteria download fails
# Line colors, keyed by label -- matplotlib's classic cycle (b, g, r, c, m,
# y, k), matching the reference figure, except Atlanta, which is pale blue
# here instead of repeating Los Angeles's blue.
line_cols <- c("Los Angeles" = "#0000FF",
"Bay Area" = "#008000",
"Chicago" = "#FF0000",
"New York" = "#00BFBF",
"Boston" = "#BF00BF",
"Seattle" = "#BFBF00",
"Las Vegas" = "#000000",
"Atlanta" = "#7FB2E5") # pale blue, to separate it from LA's blue
x_range <- c(10^2.5, 1e5) # x-axis limits for both charts, people/km².
# NULL = fit the data.
chart1_xmin <- 10^2.5 # chart 1 x-axis starts here (~316 people/km²)
chart2_xmin <- 10^2.5 # chart 2 x-axis starts here (~316 people/km²)
window_pct <- 15 # chart 2: count people within +/- this % of each density
chart2_log <- TRUE # chart 2 x-axis on a log scale (FALSE = linear)
grid_n <- 400 # chart 2: number of density points evaluated
panel_wh <- 2 # each chart's plot area is this many times wider than tall
years_to_try <- 2025:2022 # newest first; first one that answers wins
out_png <- "density_profiles.png" # both charts stacked, one figure
out_png1 <- "density_share_at_least.png" # chart 1 alone
out_png2 <- "density_people_within_window.png" # chart 2 alone
# --- download --------------------------------------------------------------
if (!nzchar(Sys.getenv("CENSUS_API_KEY")))
stop("No Census API key found. Run once:\n",
' tidycensus::census_api_key("YOUR_KEY", install = TRUE)\n',
' readRenviron("~/.Renviron")\n',
"Free key: https://api.census.gov/data/key_signup.html", call. = FALSE)
states <- unique(areas$state)
# Statewide tract pulls; the urban area filter happens spatially below.
pull_pop <- function(year) {
bind_rows(lapply(states, function(s)
get_acs(geography = "tract", variables = c(pop = "B01003_001"),
state = s, survey = "acs5", year = year, cache_table = TRUE)))
}
acs <- NULL
errs <- character(0)
for (yr in years_to_try) {
res <- tryCatch(pull_pop(yr), error = function(e) e)
if (inherits(res, "error")) {
errs[as.character(yr)] <- conditionMessage(res)
message(sprintf(" %d: %s", yr, conditionMessage(res)))
next
}
if (nrow(res) > 0) { acs <- res; acs_year <- yr; break }
}
if (is.null(acs)) stop("Every year failed. Errors:\n",
paste0(" ", names(errs), ": ", errs, collapse = "\n"),
call. = FALSE)
message(sprintf("Using ACS %d-%d 5-year estimates", acs_year - 4, acs_year))
# --- urban area boundaries -------------------------------------------------
# tigris only builds a valid download URL for the 2020-criteria urban areas
# as the full TIGER/Line file with year = 2020 + criteria = "2020" (no
# cartographic-boundary version). It's a one-time ~100 MB download, cached
# afterwards. The second attempt is a fallback for other tigris versions.
ua_attempts <- list(
function() urban_areas(year = 2020, criteria = "2020", progress_bar = FALSE),
function() urban_areas(year = ua_year, progress_bar = FALSE)
)
ua_all <- NULL
for (f in ua_attempts) {
ua_all <- tryCatch(f(), error = function(e) { message(" urban areas: ", conditionMessage(e)); NULL })
if (!is.null(ua_all)) break
}
if (is.null(ua_all)) stop("Could not download urban area boundaries (see messages above).",
call. = FALSE)
name_col <- grep("^NAME", names(ua_all), value = TRUE)[1] # NAME20 in 2020-based files
missing <- setdiff(areas$ua_name, ua_all[[name_col]])
if (length(missing) > 0) {
hints <- vapply(missing, function(m) {
key <- sub("[-, ].*$", "", m) # first word of the name
hits <- grep(key, ua_all[[name_col]], value = TRUE, ignore.case = TRUE)
paste0(" '", m, "' -- did you mean: ",
if (length(hits)) paste(head(hits, 5), collapse = " | ") else "(no match)")
}, character(1))
stop("Urban area name(s) not found:\n", paste(hints, collapse = "\n"), call. = FALSE)
}
# Equal-area projection so interior points and joins are planar-correct.
ua <- ua_all[ua_all[[name_col]] %in% areas$ua_name, name_col]
names(ua)[names(ua) == name_col] <- "ua_name"
ua <- st_transform(ua, 5070)
# --- tracts: land area + urban area assignment ----------------------------
# Same vintage as the ACS so tract GEOIDs match. ALAND excludes water.
tr <- do.call(rbind, lapply(states, function(s)
tracts(state = s, year = acs_year, cb = TRUE, progress_bar = FALSE)[, c("GEOID", "ALAND")])) |>
st_transform(5070)
# A tract belongs to an urban area if its interior point (always inside the
# polygon, unlike a centroid) falls within the boundary.
tract_ua <- suppressWarnings(st_point_on_surface(tr)) |>
st_join(ua, join = st_within, left = FALSE) |>
st_drop_geometry() |>
inner_join(distinct(areas, ua_name, area = label), by = "ua_name") |>
select(GEOID, ALAND, area)
df <- acs |>
inner_join(tract_ua, by = "GEOID") |>
transmute(area, GEOID, pop = estimate,
land_km2 = as.numeric(ALAND) / 1e6,
density = pop / land_km2)
n_in_ua <- nrow(df)
df <- df |> filter(!is.na(pop), pop > 0, land_km2 > 0) # log axis needs density > 0
if (n_in_ua > nrow(df))
message(sprintf(" %d urban-area tracts dropped (no population or no land area)",
n_in_ua - nrow(df)))
# --- summaries -------------------------------------------------------------
# Population-weighted median: half of residents live at lower density, half
# at higher. On chart 1 it's where each curve crosses 50%.
wmedian <- function(x, w) { o <- order(x); x[o][which(cumsum(w[o]) / sum(w) >= 0.5)[1]] }
stats <- df |>
group_by(area) |>
# med must come before pop: summarise() evaluates in order, and pop = sum(pop)
# would otherwise replace the per-tract weights with a single total.
summarise(tracts = n(),
med = wmedian(density, pop),
pop = sum(pop), .groups = "drop")
print(stats)
# Plot order and colors follow `areas` and `line_cols`.
lvl <- unique(areas$label)
pal <- line_cols[lvl]
if (anyNA(pal))
stop("No color in line_cols for: ", paste(lvl[is.na(pal)], collapse = ", "), call. = FALSE)
df <- df |> mutate(area = factor(area, levels = lvl))
src <- sprintf("Census tracts, ACS %d-%d 5-year estimates; land area from TIGER/Line",
acs_year - 4, acs_year)
# --- shared styling: matplotlib "classic" look -----------------------------
# White panel, full black box, ticks pointing inward, no gridlines, thin lines.
base_theme <- theme_bw(base_size = 12) +
theme(panel.grid = element_blank(),
panel.border = element_rect(colour = "black", fill = NA, linewidth = 0.6),
axis.ticks = element_line(colour = "black", linewidth = 0.4),
axis.ticks.length = unit(-4, "pt"), # inward ticks
axis.text = element_text(colour = "black", size = 11),
axis.text.x = element_text(margin = margin(t = 7)),
axis.text.y = element_text(margin = margin(r = 7)),
axis.title = element_text(size = 11),
plot.margin = margin(6, 14, 6, 6),
aspect.ratio = 1 / panel_wh) # panel height = width / panel_wh
x_lim <- if (is.null(x_range)) range(df$density) else x_range
log_x <- scale_x_log10(labels = label_log(), breaks = 10^(0:6), expand = c(0, 0))
# Minor log ticks, drawn inside the panel along the bottom like matplotlib's.
log_ticks <- annotation_logticks(sides = "b", outside = FALSE, colour = "black",
short = unit(2, "pt"), mid = unit(3, "pt"),
long = unit(4, "pt"), linewidth = 0.3)
line_w <- 0.95 # line thickness for both charts
# --- chart 1: share of residents living at AT LEAST each density ----------
# For each tract (sorted by density), the share of the urban area's population
# in tracts at that density or higher. Starts at 100% and steps down.
cum <- df |>
group_by(area) |>
arrange(density, .by_group = TRUE) |>
mutate(share_at_least = 100 * rev(cumsum(rev(pop))) / sum(pop)) |>
ungroup()
p1 <- ggplot(cum, aes(x = density, y = share_at_least, colour = area)) +
geom_step(direction = "vh", linewidth = line_w) +
scale_colour_manual(values = pal, name = NULL) +
log_x + log_ticks +
scale_y_continuous(breaks = seq(0, 100, 20), expand = c(0, 0)) +
# Zoom only: tracts below chart1_xmin still count toward the shares, so the
# curves enter at below 100% if some residents live at lower density.
coord_cartesian(xlim = c(chart1_xmin, x_lim[2]), ylim = c(0, 100)) +
labs(x = "People per sq. km.", y = "Percentage living at at least that density") +
base_theme +
# Boxed legend inside the top-right corner.
theme(legend.position = c(0.99, 0.98),
legend.justification = c(1, 1),
legend.background = element_rect(colour = "black", fill = "white", linewidth = 0.4),
legend.key = element_blank(),
legend.key.width = unit(2.4, "lines"),
legend.text = element_text(size = 13),
legend.margin = margin(4, 8, 4, 4))
# --- chart 2: people living within +/- window_pct of each density ---------
# Evaluate a grid of densities; at each one, sum the population of tracts
# whose density is within +/- window_pct of it. The window is relative, so
# it widens as density rises (at 15%: +/-150 at 1,000; +/-1,500 at 10,000).
rng <- x_lim
rng[1] <- chart2_xmin # no need to evaluate densities left of the axis
grid <- if (chart2_log) 10^seq(log10(rng[1]), log10(rng[2]), length.out = grid_n) else
seq(rng[1], rng[2], length.out = grid_n)
f <- window_pct / 100
win <- df |>
group_by(area) |>
group_modify(function(d, ...) {
inside <- outer(d$density, grid, function(x, g) abs(x - g) <= f * g)
data.frame(density = grid, people = colSums(inside * d$pop))
}) |>
ungroup()
p2 <- ggplot(win, aes(x = density, y = people, colour = area)) +
geom_line(linewidth = line_w) +
scale_colour_manual(values = pal, guide = "none") + # legend lives on chart 1
{if (chart2_log) log_x else scale_x_continuous(labels = label_comma(), expand = c(0, 0))} +
{if (chart2_log) log_ticks} +
scale_y_continuous(labels = label_number(big.mark = ""), breaks = breaks_extended(8),
expand = expansion(mult = c(0, 0.02))) +
coord_cartesian(xlim = rng) +
labs(x = "People per sq. km.",
y = sprintf("People living within %.1f%% of that density", window_pct)) +
base_theme
# --- output ----------------------------------------------------------------
fig_title <- "Density profiles of cities"
if (requireNamespace("patchwork", quietly = TRUE)) {
# Stack the two charts in one figure, like the reference.
fig <- patchwork::wrap_plots(p1, p2, ncol = 1) +
patchwork::plot_annotation(
title = fig_title,
caption = src,
theme = theme(plot.title = element_text(hjust = 0.5, size = 15),
plot.caption = element_text(colour = "grey40", size = 8)))
print(fig)
ggsave(out_png, fig, width = 10, height = 10.5, dpi = 300, bg = "white")
message("Saved ", out_png)
} else {
message("Install patchwork for the stacked figure: install.packages(\"patchwork\")")
print(p1); print(p2)
}
ggsave(out_png1, p1 + ggtitle(fig_title) + theme(plot.title = element_text(hjust = 0.5)),
width = 10, height = 5.8, dpi = 300, bg = "white")
ggsave(out_png2, p2 + ggtitle(fig_title) + theme(plot.title = element_text(hjust = 0.5)),
width = 10, height = 5.8, dpi = 300, bg = "white")
message("Saved ", out_png1, " and ", out_png2)
# ---------------------------------------------------------------------------
# Notes
# * Chart 1 is the complementary cumulative distribution: every curve starts
# at 100% on the left and falls to 0 at the urban area's densest tract. A curve
# further right means more residents live at higher density.
# * Chart 2 is raw counts, so bigger urban areas sit higher (New York's
# ~20 million people dwarf Las Vegas's ~2.2 million). For shape
# comparison, divide `people` by each urban area's total population.
# * Multi-state urban areas (Chicago, New York, Boston) have one row per
# state in `areas`; all their tracts land in the same curve.
# * Chart 2 is jagged by design: each tract enters and leaves the window as a
# block of people. Raise window_pct (e.g. 15) for a smoother line.
# * Edge tracts: a tract whose interior point is inside the urban area counts
# in full, including any rural land it holds outside the boundary, so the
# low-density tail is slightly overstated. To check the boundary match, map
# it: plot(st_geometry(ua)); plot(st_geometry(tr[tr$GEOID %in% df$GEOID, ]), add = TRUE)
# * People per mi² = people per km² x 2.59.
# ---------------------------------------------------------------------------
# Chicago, IL--IN New York--Jersey City--Newark, NY--NJ
# Boston, MA--NH Seattle--Tacoma, WA
# Las Vegas--Henderson--Paradise, NV Atlanta, GA
#
# Universe: census tracts whose interior point falls inside the urban area
# boundary. Whole tracts are used (their full population and land area), so
# tracts straddling the edge are in or out as a unit.
#
# Density = tract population / tract LAND area (people per km²).
#
# Two charts:
# 1. Share of residents living at AT LEAST each density (log x-axis).
# Read: "x% of the urban area's people live in tracts at least this dense."
# 2. Number of people living in tracts within +/-15% of each density.
# Read: "this many people live at about this density."
#
# Both are computed directly from tracts -- no smoothing.
#
# Population: most recent ACS 5-year (tries newer vintages first, falls back).
# Land area: TIGER/Line cartographic boundary tracts (ALAND) for the same year.
# Urban areas: TIGER/Line 2020 Urban Areas (2020 criteria).
#
# Requires: tidycensus, tigris, sf, dplyr, ggplot2, scales, patchwork
# install.packages(c("tidycensus", "tigris", "sf", "dplyr", "ggplot2", "scales", "patchwork"))
# One-time API key setup (free key: https://api.census.gov/data/key_signup.html):
# tidycensus::census_api_key("YOUR_KEY", install = TRUE); readRenviron("~/.Renviron")
# ---------------------------------------------------------------------------
suppressPackageStartupMessages({
library(tidycensus)
library(tigris)
library(sf)
library(dplyr)
library(ggplot2)
library(scales)
})
options(tigris_use_cache = TRUE)
# --- settings --------------------------------------------------------------
# One row per urban area. `ua_name` must match the Census name exactly --
# note the double hyphen in multi-city names. If a name is off, the script
# stops and suggests close matches. `state` = the state(s) whose tracts to
# search; for a UA crossing a state line, add a second row with the same
# label and ua_name and the other state.
areas <- data.frame(
label = c("Los Angeles", "Bay Area",
"Chicago", "Chicago",
"New York", "New York",
"Boston", "Boston",
"Seattle", "Las Vegas", "Atlanta"),
ua_name = c("Los Angeles--Long Beach--Anaheim, CA", "San Francisco--Oakland, CA",
"Chicago, IL--IN", "Chicago, IL--IN",
"New York--Jersey City--Newark, NY--NJ", "New York--Jersey City--Newark, NY--NJ",
"Boston, MA--NH", "Boston, MA--NH",
"Seattle--Tacoma, WA", "Las Vegas--Henderson--Paradise, NV", "Atlanta, GA"),
state = c("CA", "CA",
"IL", "IN",
"NY", "NJ",
"MA", "NH",
"WA", "NV", "GA"),
stringsAsFactors = FALSE
)
ua_year <- 2023 # fallback TIGER vintage if the 2020-criteria download fails
# Line colors, keyed by label -- matplotlib's classic cycle (b, g, r, c, m,
# y, k), matching the reference figure, except Atlanta, which is pale blue
# here instead of repeating Los Angeles's blue.
line_cols <- c("Los Angeles" = "#0000FF",
"Bay Area" = "#008000",
"Chicago" = "#FF0000",
"New York" = "#00BFBF",
"Boston" = "#BF00BF",
"Seattle" = "#BFBF00",
"Las Vegas" = "#000000",
"Atlanta" = "#7FB2E5") # pale blue, to separate it from LA's blue
x_range <- c(10^2.5, 1e5) # x-axis limits for both charts, people/km².
# NULL = fit the data.
chart1_xmin <- 10^2.5 # chart 1 x-axis starts here (~316 people/km²)
chart2_xmin <- 10^2.5 # chart 2 x-axis starts here (~316 people/km²)
window_pct <- 15 # chart 2: count people within +/- this % of each density
chart2_log <- TRUE # chart 2 x-axis on a log scale (FALSE = linear)
grid_n <- 400 # chart 2: number of density points evaluated
panel_wh <- 2 # each chart's plot area is this many times wider than tall
years_to_try <- 2025:2022 # newest first; first one that answers wins
out_png <- "density_profiles.png" # both charts stacked, one figure
out_png1 <- "density_share_at_least.png" # chart 1 alone
out_png2 <- "density_people_within_window.png" # chart 2 alone
# --- download --------------------------------------------------------------
if (!nzchar(Sys.getenv("CENSUS_API_KEY")))
stop("No Census API key found. Run once:\n",
' tidycensus::census_api_key("YOUR_KEY", install = TRUE)\n',
' readRenviron("~/.Renviron")\n',
"Free key: https://api.census.gov/data/key_signup.html", call. = FALSE)
states <- unique(areas$state)
# Statewide tract pulls; the urban area filter happens spatially below.
pull_pop <- function(year) {
bind_rows(lapply(states, function(s)
get_acs(geography = "tract", variables = c(pop = "B01003_001"),
state = s, survey = "acs5", year = year, cache_table = TRUE)))
}
acs <- NULL
errs <- character(0)
for (yr in years_to_try) {
res <- tryCatch(pull_pop(yr), error = function(e) e)
if (inherits(res, "error")) {
errs[as.character(yr)] <- conditionMessage(res)
message(sprintf(" %d: %s", yr, conditionMessage(res)))
next
}
if (nrow(res) > 0) { acs <- res; acs_year <- yr; break }
}
if (is.null(acs)) stop("Every year failed. Errors:\n",
paste0(" ", names(errs), ": ", errs, collapse = "\n"),
call. = FALSE)
message(sprintf("Using ACS %d-%d 5-year estimates", acs_year - 4, acs_year))
# --- urban area boundaries -------------------------------------------------
# tigris only builds a valid download URL for the 2020-criteria urban areas
# as the full TIGER/Line file with year = 2020 + criteria = "2020" (no
# cartographic-boundary version). It's a one-time ~100 MB download, cached
# afterwards. The second attempt is a fallback for other tigris versions.
ua_attempts <- list(
function() urban_areas(year = 2020, criteria = "2020", progress_bar = FALSE),
function() urban_areas(year = ua_year, progress_bar = FALSE)
)
ua_all <- NULL
for (f in ua_attempts) {
ua_all <- tryCatch(f(), error = function(e) { message(" urban areas: ", conditionMessage(e)); NULL })
if (!is.null(ua_all)) break
}
if (is.null(ua_all)) stop("Could not download urban area boundaries (see messages above).",
call. = FALSE)
name_col <- grep("^NAME", names(ua_all), value = TRUE)[1] # NAME20 in 2020-based files
missing <- setdiff(areas$ua_name, ua_all[[name_col]])
if (length(missing) > 0) {
hints <- vapply(missing, function(m) {
key <- sub("[-, ].*$", "", m) # first word of the name
hits <- grep(key, ua_all[[name_col]], value = TRUE, ignore.case = TRUE)
paste0(" '", m, "' -- did you mean: ",
if (length(hits)) paste(head(hits, 5), collapse = " | ") else "(no match)")
}, character(1))
stop("Urban area name(s) not found:\n", paste(hints, collapse = "\n"), call. = FALSE)
}
# Equal-area projection so interior points and joins are planar-correct.
ua <- ua_all[ua_all[[name_col]] %in% areas$ua_name, name_col]
names(ua)[names(ua) == name_col] <- "ua_name"
ua <- st_transform(ua, 5070)
# --- tracts: land area + urban area assignment ----------------------------
# Same vintage as the ACS so tract GEOIDs match. ALAND excludes water.
tr <- do.call(rbind, lapply(states, function(s)
tracts(state = s, year = acs_year, cb = TRUE, progress_bar = FALSE)[, c("GEOID", "ALAND")])) |>
st_transform(5070)
# A tract belongs to an urban area if its interior point (always inside the
# polygon, unlike a centroid) falls within the boundary.
tract_ua <- suppressWarnings(st_point_on_surface(tr)) |>
st_join(ua, join = st_within, left = FALSE) |>
st_drop_geometry() |>
inner_join(distinct(areas, ua_name, area = label), by = "ua_name") |>
select(GEOID, ALAND, area)
df <- acs |>
inner_join(tract_ua, by = "GEOID") |>
transmute(area, GEOID, pop = estimate,
land_km2 = as.numeric(ALAND) / 1e6,
density = pop / land_km2)
n_in_ua <- nrow(df)
df <- df |> filter(!is.na(pop), pop > 0, land_km2 > 0) # log axis needs density > 0
if (n_in_ua > nrow(df))
message(sprintf(" %d urban-area tracts dropped (no population or no land area)",
n_in_ua - nrow(df)))
# --- summaries -------------------------------------------------------------
# Population-weighted median: half of residents live at lower density, half
# at higher. On chart 1 it's where each curve crosses 50%.
wmedian <- function(x, w) { o <- order(x); x[o][which(cumsum(w[o]) / sum(w) >= 0.5)[1]] }
stats <- df |>
group_by(area) |>
# med must come before pop: summarise() evaluates in order, and pop = sum(pop)
# would otherwise replace the per-tract weights with a single total.
summarise(tracts = n(),
med = wmedian(density, pop),
pop = sum(pop), .groups = "drop")
print(stats)
# Plot order and colors follow `areas` and `line_cols`.
lvl <- unique(areas$label)
pal <- line_cols[lvl]
if (anyNA(pal))
stop("No color in line_cols for: ", paste(lvl[is.na(pal)], collapse = ", "), call. = FALSE)
df <- df |> mutate(area = factor(area, levels = lvl))
src <- sprintf("Census tracts, ACS %d-%d 5-year estimates; land area from TIGER/Line",
acs_year - 4, acs_year)
# --- shared styling: matplotlib "classic" look -----------------------------
# White panel, full black box, ticks pointing inward, no gridlines, thin lines.
base_theme <- theme_bw(base_size = 12) +
theme(panel.grid = element_blank(),
panel.border = element_rect(colour = "black", fill = NA, linewidth = 0.6),
axis.ticks = element_line(colour = "black", linewidth = 0.4),
axis.ticks.length = unit(-4, "pt"), # inward ticks
axis.text = element_text(colour = "black", size = 11),
axis.text.x = element_text(margin = margin(t = 7)),
axis.text.y = element_text(margin = margin(r = 7)),
axis.title = element_text(size = 11),
plot.margin = margin(6, 14, 6, 6),
aspect.ratio = 1 / panel_wh) # panel height = width / panel_wh
x_lim <- if (is.null(x_range)) range(df$density) else x_range
log_x <- scale_x_log10(labels = label_log(), breaks = 10^(0:6), expand = c(0, 0))
# Minor log ticks, drawn inside the panel along the bottom like matplotlib's.
log_ticks <- annotation_logticks(sides = "b", outside = FALSE, colour = "black",
short = unit(2, "pt"), mid = unit(3, "pt"),
long = unit(4, "pt"), linewidth = 0.3)
line_w <- 0.95 # line thickness for both charts
# --- chart 1: share of residents living at AT LEAST each density ----------
# For each tract (sorted by density), the share of the urban area's population
# in tracts at that density or higher. Starts at 100% and steps down.
cum <- df |>
group_by(area) |>
arrange(density, .by_group = TRUE) |>
mutate(share_at_least = 100 * rev(cumsum(rev(pop))) / sum(pop)) |>
ungroup()
p1 <- ggplot(cum, aes(x = density, y = share_at_least, colour = area)) +
geom_step(direction = "vh", linewidth = line_w) +
scale_colour_manual(values = pal, name = NULL) +
log_x + log_ticks +
scale_y_continuous(breaks = seq(0, 100, 20), expand = c(0, 0)) +
# Zoom only: tracts below chart1_xmin still count toward the shares, so the
# curves enter at below 100% if some residents live at lower density.
coord_cartesian(xlim = c(chart1_xmin, x_lim[2]), ylim = c(0, 100)) +
labs(x = "People per sq. km.", y = "Percentage living at at least that density") +
base_theme +
# Boxed legend inside the top-right corner.
theme(legend.position = c(0.99, 0.98),
legend.justification = c(1, 1),
legend.background = element_rect(colour = "black", fill = "white", linewidth = 0.4),
legend.key = element_blank(),
legend.key.width = unit(2.4, "lines"),
legend.text = element_text(size = 13),
legend.margin = margin(4, 8, 4, 4))
# --- chart 2: people living within +/- window_pct of each density ---------
# Evaluate a grid of densities; at each one, sum the population of tracts
# whose density is within +/- window_pct of it. The window is relative, so
# it widens as density rises (at 15%: +/-150 at 1,000; +/-1,500 at 10,000).
rng <- x_lim
rng[1] <- chart2_xmin # no need to evaluate densities left of the axis
grid <- if (chart2_log) 10^seq(log10(rng[1]), log10(rng[2]), length.out = grid_n) else
seq(rng[1], rng[2], length.out = grid_n)
f <- window_pct / 100
win <- df |>
group_by(area) |>
group_modify(function(d, ...) {
inside <- outer(d$density, grid, function(x, g) abs(x - g) <= f * g)
data.frame(density = grid, people = colSums(inside * d$pop))
}) |>
ungroup()
p2 <- ggplot(win, aes(x = density, y = people, colour = area)) +
geom_line(linewidth = line_w) +
scale_colour_manual(values = pal, guide = "none") + # legend lives on chart 1
{if (chart2_log) log_x else scale_x_continuous(labels = label_comma(), expand = c(0, 0))} +
{if (chart2_log) log_ticks} +
scale_y_continuous(labels = label_number(big.mark = ""), breaks = breaks_extended(8),
expand = expansion(mult = c(0, 0.02))) +
coord_cartesian(xlim = rng) +
labs(x = "People per sq. km.",
y = sprintf("People living within %.1f%% of that density", window_pct)) +
base_theme
# --- output ----------------------------------------------------------------
fig_title <- "Density profiles of cities"
if (requireNamespace("patchwork", quietly = TRUE)) {
# Stack the two charts in one figure, like the reference.
fig <- patchwork::wrap_plots(p1, p2, ncol = 1) +
patchwork::plot_annotation(
title = fig_title,
caption = src,
theme = theme(plot.title = element_text(hjust = 0.5, size = 15),
plot.caption = element_text(colour = "grey40", size = 8)))
print(fig)
ggsave(out_png, fig, width = 10, height = 10.5, dpi = 300, bg = "white")
message("Saved ", out_png)
} else {
message("Install patchwork for the stacked figure: install.packages(\"patchwork\")")
print(p1); print(p2)
}
ggsave(out_png1, p1 + ggtitle(fig_title) + theme(plot.title = element_text(hjust = 0.5)),
width = 10, height = 5.8, dpi = 300, bg = "white")
ggsave(out_png2, p2 + ggtitle(fig_title) + theme(plot.title = element_text(hjust = 0.5)),
width = 10, height = 5.8, dpi = 300, bg = "white")
message("Saved ", out_png1, " and ", out_png2)
# ---------------------------------------------------------------------------
# Notes
# * Chart 1 is the complementary cumulative distribution: every curve starts
# at 100% on the left and falls to 0 at the urban area's densest tract. A curve
# further right means more residents live at higher density.
# * Chart 2 is raw counts, so bigger urban areas sit higher (New York's
# ~20 million people dwarf Las Vegas's ~2.2 million). For shape
# comparison, divide `people` by each urban area's total population.
# * Multi-state urban areas (Chicago, New York, Boston) have one row per
# state in `areas`; all their tracts land in the same curve.
# * Chart 2 is jagged by design: each tract enters and leaves the window as a
# block of people. Raise window_pct (e.g. 15) for a smoother line.
# * Edge tracts: a tract whose interior point is inside the urban area counts
# in full, including any rural land it holds outside the boundary, so the
# low-density tail is slightly overstated. To check the boundary match, map
# it: plot(st_geometry(ua)); plot(st_geometry(tr[tr$GEOID %in% df$GEOID, ]), add = TRUE)
# * People per mi² = people per km² x 2.59.
# ---------------------------------------------------------------------------
No comments:
Post a Comment
And your thoughts on the matter?