Andrew Heiss's blog

Generating universes within universes with a single seed

Andrew Heiss — Mon, 13 Apr 2026 04:00:00 GMT

“…as an ook cometh of a litel spyr…”

Geoffrey Chaucer (1343–1400),
Troilus and Criseyde Book II, line 1335.
In Chaucer’s Middle English, a “spyr” is a seedling.

Mighty oaks from little acorns grow.

Traditional English idiom, perhaps based on Chaucer.

If you’ve ever done anything involving randomness with statistical computing, you’ve likely come across the idea of a “seed.”

Let’s say you make a plot like this to show distribution of the counts of different species of penguins across different islands:

library(ggplot2)

ggplot(penguins, aes(x = species, y = island, color = island)) +
  geom_point() +
  guides(color = "none")

This plot is kind of helpful—it shows that Torgersen Island only has Adélie penguins, Dream Island has Adélies and Chinstraps, and Biscoe Island has Adélies and Gentoos. HOWEVER, there are actually 300+ points here that are all stacked on top of each other. We can jitter the points a little to see them:

ggplot(penguins, aes(x = species, y = island, color = island)) +
  geom_point(position = position_jitter()) +
  guides(color = "none")

If we make the plot again, the points will again be randomly shuffled, and it’ll look slightly different, which can be annoying if you want to get the same consistent plot:

ggplot(penguins, aes(x = species, y = island, color = island)) +
  geom_point(position = position_jitter()) +
  guides(color = "none")

To make the points randomly shuffle the same way every time, we can set a seed:

ggplot(penguins, aes(x = species, y = island, color = island)) +
  geom_point(position = position_jitter(seed = 1234)) +
  guides(color = "none")

Now every time we run the code for that plot, it’ll put those points in the same shuffled locations every time.

Seeds aren’t just for plotting. If you’ve ever done things with Bayesian modeling or other types of simulation, you’ve also likely seen seeds. Modern Bayesian methods estimate posterior parameters with fancy Markov chain Monte Carlo (MCMC) simulations, which randomly walk through possible parameter values to figure out integrals. It’s magical.¹ To make MCMC chains take the same random walk on different computers, you can set a seed, and you’ll get the same results every time. Like, run this {brms} regression model on your computer and you should get exactly the same values:

¹ See chapter 7 in Bayes Rules! (Johnson, Ott, and Dogucu 2022) for the most intuitive explanation of MCMC I’ve come across.

library(brms)
library(parameters)

model <- brm(
  bf(body_mass ~ flipper_len + species),
  data = penguins,
  family = gaussian(),
  backend = "cmdstanr",
  seed = 1234, refresh = 0
)

model_parameters(model, component = "location", verbose = FALSE)
#> Parameter        |   Median |               95% CI |     pd |  Rhat |  ESS
#> --------------------------------------------------------------------------
#> (Intercept)      | -4036.54 | [-5198.18, -2897.46] |   100% | 0.999 | 2318
#> flipper_len      |    40.73 | [   34.81,    46.92] |   100% | 0.999 | 2268
#> speciesChinstrap |  -207.92 | [ -326.51,   -90.19] |   100% | 1.000 | 2862
#> speciesGentoo    |   264.31 | [   77.93,   450.38] | 99.55% | 1.000 | 2229

Even non-statistical things use seeds. Every Minecraft world is procedurally generated based on random numbers, and when you create a world, you have an option to specify a seed so that you can create specific random worlds. Later on, you can find the seed for a Minecraft world with the /seed command and recreate that exact world again using that seed. There’s even a whole subreddit dedicated to sharing weird and interesting things that get generated when using specific seeds. Since these random worlds are reproducible, there are sites like Chunkbase that let you see the whole world map for any given seed.

For fun, I spun up Minecraft and created a new world. It happened to use the seed -6863317686781453459. I took a screenshot of my view right at the world spawn, and that exact spot is visible in Chunkbase’s map:

Here’s the map for this seed. That compass is where I’m standing; I’m looking across that river at the area that’s in the bottom left corner of that red rectangle.

I’ve been using and thinking about seeds for more than a decade, but I have no idea how they actually work. All I know is that setting a seed lets you get the same random numbers again. Somehow.

But I want to know why! Let’s figure it out.

Seeds and reproducible randomness

Put simply, a seed lets you generate random numbers in a reproducible way. Any time you do something random in R—use a function like sample() or rnorm() or shuffle points in a plot with position_jitter() or move labels around with geom_label_repel()—it uses some sort of starting point to make up those numbers. If you set a seed, your jittered points will be jittered the same way every time you make the plot; your random simulation will be the same every time you run it. That’s great for reproducibility and consistency. Seeds create reproducible randomness.

A seed is an arbitrary integer that R’s internal random number generator (RNG) algorithm needs to initialize its number-generating process. If you don’t set a seed yourself, R will make one for itself through a combination of (1) the current time on your computer and (2) the internal ID of the R process that your computer is running. Seeds can be any number you want. People choose things like 1, 1234, 12345, 42, 8675309, the date like 20260413, or whatever. (Though, we’ll see later in this post that limiting yourself to these super common seeds is actually bad!)

I’ll talk about R throughout this post, but the same principle applies everywhere. Every computer-based process for generating random numbers has to start with some initial number or seed. Often you have the ability to set that seed yourself—Python, Julia, Stata, and even Minecraft give you control over the starting value. Often you can’t set the seed. Excel has a RAND() function, but there’s no way to set a seed (though people try!). Javascript does random things, but cannot set a seed (though there are workarounds!). Switch games like Mario Kart and the Zelda franchise all use randomness, but there’s no way to set a seed (I’d guess there’s some secret developer mode to do it for testing though).

My (somewhat incorrect) mental model of how seeds work

When I teach about seeds, I’ll typically have everyone in the class run this to generate 5 random uniformly distributed numbers:

runif(5)
## [1] 0.99034 0.17013 0.08634 0.98750 0.09786

I then ask if anyone got the same five numbers. Nobody ever has.

I then have them run this:

set.seed(1234)
runif(5)
## [1] 0.1137 0.6223 0.6093 0.6234 0.8609

Students are generally surprised to see that everyone gets the same five numbers every time. They all get the same thing because we set the seed to 1234, so the random number generation process all starts at the same point and ends up with the same values.

I then have students run this:

set.seed(1234)
runif(5)
## [1] 0.1137 0.6223 0.6093 0.6234 0.8609
runif(5)
## [1] 0.640311 0.009496 0.232551 0.666084 0.514251

They’ll get the same 5 numbers as before, and then a different new set of 5 numbers. They all get the same new set of 5 numbers.

To explain why the second runif() produced different random numbers, I would talk about “rounds”, which generally made sense to me conceptually. In my mind, it felt like R was kind of like bumping the seed up each time a random function was used. It almost like the first random function used 1234, then the next time R did something random it used something analogous to 1234.1, then 1234.2, and so on. I know that this is 100% wrong, but that’s roughly how I’ve imagined it—a stack of seed-like “rounds” that R works through every time it does something random:

set.seed(1234)
runif(5)  # This uses "round 1" of the seed. The seed is 1234
## [1] 0.1137 0.6223 0.6093 0.6234 0.8609
runif(5)  # This uses "round 2" of the seed. Pretend it's like the seed is 1234.1 (it's not!)
## [1] 0.640311 0.009496 0.232551 0.666084 0.514251

That also fits with this idea of resetting the “rounds” when resetting the seed:

set.seed(1234)
runif(5)  # This uses "round 1" of the seed at 1234
## [1] 0.1137 0.6223 0.6093 0.6234 0.8609

set.seed(1234)
runif(5)  # This uses "round 1" of the seed again because it was reset to 1234
## [1] 0.1137 0.6223 0.6093 0.6234 0.8609

But then a couple weeks ago when teaching about seeds again, I noticed something new that broke my mental model of how seeds actually work.

To show this, let’s set the seed to 1234 and then create 5 random uniformly distributed numbers:

set.seed(1234)
runif(5)
## [1] 0.1137 0.6223 0.6093 0.6234 0.8609

Cool. That’s all normal and expected. Let’s reset the seed to 1234 and generate just 2 numbers:

set.seed(1234)
runif(2)
## [1] 0.1137 0.6223

Those two are the same two numbers that appeared at the beginning of runif(5). That also makes sense to me, since we’re using “round 1” of the seed.

Let’s make 3 more random numbers. In my mental version of how seeds work, this is a new random-related function, so it should start from a kinda-sorta different seed, or the “round 2” version of the seed:

runif(3)
## [1] 0.6093 0.6234 0.8609

It doesn’t though! Those three numbers are the same as the numbers 3, 4, and 5 in runif(5)! Here’s everything all at once:

set.seed(1234)
runif(5)
## [1] 0.1137 0.6223 0.6093 0.6234 0.8609

set.seed(1234)
runif(2)
## [1] 0.1137 0.6223
runif(3)
## [1] 0.6093 0.6234 0.8609

# Or all in the same vector!
set.seed(1234)
c(runif(2), runif(3))
## [1] 0.1137 0.6223 0.6093 0.6234 0.8609

Using a second random-related function didn’t bump up to a new kinda-sorta seed. R isn’t generating a new set of random numbers each time it does something random.

Instead, R is continuing to use a single list of random numbers. We showed 5 of them the first time, then reset the seed, then showed 2, and then showed the next 3. Every function that uses randomness basically moves a pointer or a cursor through the same list of random values.

That’s wild.

So where does this list come from? How does a single seed create one huge list of random values? Will it run out of random numbers? If I jitter a scatterplot with a million points (probably a bad idea), will that consume all the random values in the list and mess things up if I want to use sample() or rnorm() later?

Making “random” numbers with an equation

Around the same time I discovered that runif(5) and c(runif(2), runif(3)) did the same thing, I was reading/listening to the audiobook of The Art of Uncertainty (Spiegelhalter 2025). Spiegelhalter has a chapter on the uncertainty associated with randomness and he spends some time explaining how equations let you generate a deterministic set of pseudorandom (or random-looking) values. He walks through a basic algorithm for generating random numbers called the linear congruential generator, or LCG. I’ve always wondered how these magic RNG algorithms actually worked, and I had to take a break from doing the dishes to run to Wikipedia to look up exactly how to use this equation.

The basic LCG algorithm looks like this:

I’m going to change it a little bit from what Wikipedia has because I don’t like thinking about calculating the next value of (or ) and instead want to think about calculating the current value of (or ). The only thing that’s different here is that the subscripts are and instead of and :

What is that “mod” thing?

The only potentially unfamiliar thing in this formula is that “mod” operator. This represents a modulo operation. That might sound intimidating and weird, but in practice it’s something you learned back in elementary school when you learned long division.

For example, what’s 17 ÷ 5? With long division, we can calculate that it’s 3, with a remainder of 2:

That remainder is the result of the modulus operator. That’s it. “17 mod 5” is 2, since 17 divided by 5 is 3 with a remainder of 2.

What about “19 mod 3”? That’s the remainder from 19 ÷ 3, or 1:

We can calculate the modulus directly in R with the %% mod operator:

17 %% 5
## [1] 2
19 %% 3
## [1] 1

Python uses a single %:

17 % 5
## 2
19 % 3
## 1

So if in the LCG algorithm is 5, we need to calculate the remainder when dividing by 5. If we have to do something like 24 mod 5, that’s 24 ÷ 5, or 4 remainder 4.

24 mod 5 = 4:

24 %% 5
## [1] 4

Fun unrelated side note: Odd and even numbers

The mod operator is a super common way to determine if a number is even or odd. If a number is even, it should have no remainder when you divide it by 2. 24 mod 2 should be 0:

## [1] 0

If there’s a remainder when dividing by 2, it means that the number is not even. 25 mod 2 is 1 (since 25 ÷ 2 = 12 remainder 1)

## [1] 1

You can use this mod operator to make your own function for checking if a number is even/odd:

## [1] "Nope! The number is odd!"
## [1] "Yep! The number is even!"

This algorithm generates a sequence of pseudorandom numbers based on a few different inputs:

, or the modulus. This should be bigger than 0, or more formally, .
, or the multiplier. This should be bigger than 0 and generally be smaller than , or .
, or the increment. This should be bigger or equal to 0 and generally be smaller than , or .
, or the previous pseudorandom value in the sequence.

Plug four numbers into the equation and you’ll get a pseudorandom number. Plug that number back into the equation and you’ll get a new pseudorandom number. Easy peasy.

You can generally use whatever numbers you want for , , and , but notice how you need to use the previous value from the sequence () to generate a new value (). That’s fine if you’re in the middle of a sequence, but if you’re at the beginning and want to generate , there is no previous value at . You need to provide your own. This is the seed. The seed is the initial value for the algorithm that generates the rest of the sequence of pseudorandom numbers.

Let’s make some random numbers one step at a time to illustrate how this works. We’ll arbitrarily set some values for , , and , and we’ll use a seed (or ) of 1.

Parameter		Value
Modulus	$m$	8
Multiplier	$a$	5
Increment	$c$	3
Seed	$X_0$	1

To generate the first number in the sequence, we’ll plug in all those values:

Our first “random” number is 0. Neat!

To make the next one in the sequence, we’ll use that 0 for :

The second number in the sequence is 3. Let’s do one more step to get :

is 2.

So far, our sequence of pseudorandom numbers looks like this:

The LCG algorithm produces a pseudorandom uniform distribution of integers. We can scale these down to a range between 0 and 1 by dividing by :²

² This works because we’ll never actually have any random number higher than . When dividing by 8, we’ll never have a remainder of 8 or 9 or anything, and 7 will be the biggest possible remainder.

Instead of doing this slowly step by step, here’s a table showing the sequence:

Table 1: Twelve pseudorandom numbers generated with the LCG algorithm: with a seed () of 1

Parameter		Value
Modulus	$m$	8
Multiplier	$a$	5
Increment	$c$	3
Seed	$X_0$	1

			The “random” number!
	$X_{n - 1}$	$(a X_{n - 1} + c)$	… mod $m$	$X_n / m$
$X_{1}$	1	8	0	0.000
$X_{2}$	0	3	3	0.375
$X_{3}$	3	18	2	0.250
$X_{4}$	2	13	5	0.625
$X_{5}$	5	28	4	0.500
$X_{6}$	4	23	7	0.875
$X_{7}$	7	38	6	0.750
$X_{8}$	6	33	1	0.125
$X_{9}$	1	8	0	0.000
$X_{10}$	0	3	3	0.375
$X_{11}$	3	18	2	0.250
$X_{12}$	2	13	5	0.625

We successfully created a sequence of pseudorandom numbers!

Since this series of numbers is uniformly distributed, we can pick numbers from the list to generate sets of random numbers. Let’s pretend we have our own function named simple_runif() that mimics R’s runif() but draws from this table that we just made.

If we want 5 random numbers from a uniform distribution, we can take the first 5 values:

If we wanted another set of 5 random numbers, we’d get the next 5:

If we set the seed back to 1 with an imaginary function like simple_set.seed(1) and then generate 2 random numbers, and then three random numbers (like I did earlier with the real runif(2) and runif(3)), we’d go back to the beginning of the series of numbers and get these:

We can adjust any of the and parameters, the modulus , or the starting seed to get a different set of numbers:

Parameter		Value
Modulus	$m$	8
Multiplier	$a$	5
Increment	$c$	3
Seed	$X_0$	1234

			The “random” number!
	$X_{n - 1}$	$(a X_{n - 1} + c)$	… mod $m$	$X_n / m$
$X_{1}$	1234	6173	5	0.625
$X_{2}$	5	28	4	0.500
$X_{3}$	4	23	7	0.875
$X_{4}$	7	38	6	0.750
$X_{5}$	6	33	1	0.125
$X_{6}$	1	8	0	0.000
$X_{7}$	0	3	3	0.375
$X_{8}$	3	18	2	0.250
$X_{9}$	2	13	5	0.625
$X_{10}$	5	28	4	0.500
$X_{11}$	4	23	7	0.875
$X_{12}$	7	38	6	0.750

Parameter		Value
Modulus	$m$	23
Multiplier	$a$	11
Increment	$c$	4
Seed	$X_0$	10

			The “random” number!
	$X_{n - 1}$	$(a X_{n - 1} + c)$	… mod $m$	$X_n / m$
$X_{1}$	10	114	22	0.95652
$X_{2}$	22	246	16	0.69565
$X_{3}$	16	180	19	0.82609
$X_{4}$	19	213	6	0.26087
$X_{5}$	6	70	1	0.04348
$X_{6}$	1	15	15	0.65217
$X_{7}$	15	169	8	0.34783
$X_{8}$	8	92	0	0.00000
$X_{9}$	0	4	4	0.17391
$X_{10}$	4	48	2	0.08696
$X_{11}$	2	26	3	0.13043
$X_{12}$	3	37	14	0.60870

Parameter		Value
Modulus	$m$	57
Multiplier	$a$	54
Increment	$c$	17
Seed	$X_0$	42

			The “random” number!
	$X_{n - 1}$	$(a X_{n - 1} + c)$	… mod $m$	$X_n / m$
$X_{1}$	42	2285	5	0.08772
$X_{2}$	5	287	2	0.03509
$X_{3}$	2	125	11	0.19298
$X_{4}$	11	611	41	0.71930
$X_{5}$	41	2231	8	0.14035
$X_{6}$	8	449	50	0.87719
$X_{7}$	50	2717	38	0.66667
$X_{8}$	38	2069	17	0.29825
$X_{9}$	17	935	23	0.40351
$X_{10}$	23	1259	5	0.08772
$X_{11}$	5	287	2	0.03509
$X_{12}$	2	125	11	0.19298

Live interactive playground

To help with the intuition behind all these moving parts, I find that it’s really useful to play with a live, interactive version of the LCG algorithm. When I cover this in class, I use Excel to build a notebook with adjustable parameters. You can download your own Excel simulation:

seed-playground.xlsx

Here’s an interactive in-browser OJS-based version of the same thing. Adjust , , , , and to see how parameters like and influence the calculations, how determines the length of a period, and how the seed leads to different series of pseudorandom numbers.

viewof m = Inputs.range([2, 500], {value: 8, step: 1, label: htl.html`Modulus (m)`})
viewof a = Inputs.range([1, 20], {value: 5, step: 1, label: htl.html`Multiplier (a)`})
viewof c = Inputs.range([0, 20], {value: 3, step: 1, label: htl.html`Increment (c)`})
viewof seed = Inputs.range([0, 2000], {value: 1, step: 1, label: htl.html`Seed (X₀)`})
viewof n = Inputs.range([1, 50], {value: 12, step: 1, label: htl.html`

Numbers to generate`})

lcg_data = {
  let rows = [];
  let x = seed;
  for (let i = 1; i <= n; i++) {
    const ax_c = a * x + c;
    const mod = ax_c % m;
    rows.push({
      n: i,
      x_prev: x,
      ax_plus_c: ax_c,
      mod_m: mod,
      x_over_m: mod / m
    });
    x = mod;
  }
  return rows;
}

htl.html`                                                                                ${lcg_data.map((row, i) => {      const colors = ["#CC66774C", "#DDCC774C", "#88CCEE4C"];      const cur = colors[i % 3];      const prev = colors[(i + 2) % 3];      return htl.html`                                              `;    })}  





The "random" number!



X_n−1
(aX_n−1 + c)
… mod m
X_n / m








X_${row.n}
${i === 0 ? '#DDAA33' : prev}; font-weight: ${i === 0 ? 'bold' : 'normal'}">${row.x_prev}
${row.ax_plus_c}
${cur}">${row.mod_m}
${row.x_over_m.toFixed(3)}



`

			The "random" number!
X_${row.n}	${i === 0 ? '#DDAA33' : prev}; font-weight: ${i === 0 ? 'bold' : 'normal'}">${row.x_prev}	${row.ax_plus_c}	${cur}">${row.mod_m}	${row.x_over_m.toFixed(3)}

Cycles and fancier algorithms

Look back at Table 1 where we we used m = 8, a = 5, c = 3, and a seed of 1. Something peculiar happens at the 9th number:

	$X_{1}$	$X_{2}$	$X_{3}$	$X_{4}$	$X_{5}$	$X_{6}$	$X_{7}$	$X_{8}$	$X_{9}$	$X_{10}$	$X_{11}$	$X_{12}$
$X_n$	0	3	2	5	4	7	6	1	0	3	2	5
$X_{n} / m$	0.000	0.375	0.250	0.625	0.500	0.875	0.750	0.125	0.000	0.375	0.250	0.625

It starts over! and are the same; and are the same; and so on. We have a repeating cycle and our series of pseudorandom numbers gets recycled.

This happens because of a severe limitation in the LCG algorithm—cycles, or periods, are inevitable and are related to the size of the modulus parameter . Since we’re dividing by 8 and getting the remainder, we only have the numbers 0–7 to work with for each subsequent number, and eventually we’ll run out and get into a cycle. We can somewhat avoid this by boosting the . If we increase the modulus to some big prime number like 499 (try it in the interactive simulation up above), we’ll have remainders like 403, 291, and so on. That many unique possibilities will make the period last longer, but eventually it will recycle.

So to actually get around this, no statistical or computational software really uses LCG.³ Instead, they use much fancier techniques with more complex algorithms that prevent cycling. Fundamentally all these algorithms do the same thing—they start with a seed number, do some math to it, and spit out a random number that is then used to create the next number in the sequence. They just do it with really fancy math.

³ Though, surprisingly, Java used an LCG up until 2020(!). It only switched to something more complex with Java 17.

For instance, by default R uses an algorithm called the Mersenne Twister—which was only invented in 1997! (Matsumoto and Nishimura 1998)—that takes a dozen+ parameters and uses a special seed that is actually a vector with 624 elements in it (!). See Danielle Navarro’s fascinating post about this here, where she looks into the guts of R’s random number generation system, including some of the individual pieces of the vector-based seed.

The Mersenne Twister algorithm will eventually recycle and restart its random numbers for a given seed, but its period length is massive: 2¹⁹⁹³⁷ − 1. That number has 6,002 digits, like 10⁶⁰⁰², but instead of all 0s like 10⁶⁰⁰², it’s a prime number. You can actually see all the digits here. Claus Wilke calls this period length “an unimaginably large number” and his description of what that actually means is poetic:

To give you a sense of how large it is, for comparison, there are approximately 10⁸⁰ atoms in the universe. This is tiny compared to 10⁶⁰⁰⁰. The period of the Mersenne twister has space for entire universes for every single atom in the universe, and then some. In fact, you could create an entire universe for every atom, and then create another entire universe for every atom in every of the universes you have created, and keep nesting 75 times, and still you wouldn’t run out of room in the period of the Mersenne twister. If you used the Mersenne twister to create nested universes 75 times deep, all these universes inside universes inside other universes would be different from each other.

So unless you’re creating deeply nested universes, in practice you’ll never recycle the numbers that are generated from a single seed based on the Mersenne Twister algorithm.

Again, even though the Mersenne Twister (and whatever other fancy algorithms other languages and programs use) is complex, the intuition from the LCG remains the same: generate a big long list of pseudorandom numbers by starting with an initial value—or a seed—and then plugging each subsequent random number into an equation to make more.

Why does it matter if “random” numbers aren’t actually random?

I’ve long known that random computer generated numbers aren’t officially random, but with my wrong idea of R setting a kinda-sorta seed each time it uses a new function, it still felt fairly random-ish, since I figured that new sets of numbers would get generated.

But that’s not the case at all!

When you do something like set.seed(1234), you essentially create an unfathomably huge list of all possible pseudorandom numbers arranged in a predetermined order. Setting a seed basically locks you into a predetermined universe devoid of true randomness. Each random thing you do—shuffling points in a plot, running MCMC simulations, randomly assigning rows to treatment or control conditions, and so on—just moves a pointer down through the list of already-known numbers.

This isn’t just a theoretical stoner shower thought (“whoa man, there’s no free will”). The absence of true randomness has actual consequences for analysis!

You’re limiting yourself to narrow, known universes

In this really great blog post, Claus Wilke makes a Jenny Bryan-esque threat: “If your random seed is 42 I will come to your office and set your computer on fire🔥.” He argues that using common seeds like 42, 1234, 1, and so on locks analyses into predetermined universes. If everyone uses 42 or 1234, all “random” sampling, MCMC simulation, and so on will happen with the same list of Mersenne Twister numbers. Again, there are 2¹⁹⁹³⁷ − 1 of those before things maybe start to recycle, so we’re not worried about running out of numbers. What we’re worried about is overusing the first chunk of those numbers. As Claus says:

The Mersenne Twister has a state space large enough for universes within universes, but every data scientist in the entire world is using the same 10,000 “random” numbers that you get when starting with seed 42.

Setting seeds is good for reproducibility! You often want people to recreate your exact results. Setting seeds is great for plotting things too—if you’re using {ggrepel} to push labels around randomly, you want every iteration of your plot to place your labels in the same position any time. You don’t have to limit yourself to common seeds like 0, 1, 1234, or 42.

R does have a maximum seed size—it’s the biggest value allowed for a 32-bit signed integer, or:

.Machine$integer.max
## [1] 2147483647

That means you can set the seed to any number between −2,147,483,647 and 2,147,483,647. You have 4,294,967,294 possible seeds to work with—quit using 1234 and explore new universes!

You can seed hack and get any values you want

Since massive lists of pseudorandom numbers are inherently deterministic, you can brute-force search for a seed that produces any random outcome you want.

This isn’t always bad! I do this all the time with plotting. Like, let’s say I make a plot like this with geom_text_repel():

library(ggrepel)

ggplot(...) + 
  geom_point() + 
  geom_label_repel(..., seed = 1234)

I don’t like where the labels ended up. Maybe some overlapped, or maybe some are covering up important points. To fix it, I’ll change the seed to 12345:

... +
  geom_label_repel(..., seed = 12345)

…or 123:

... + 
  geom_label_repel(..., seed = 123)

…or anything else until it looks good.

That’s a manual version of seed hacking (though it’s super benign in this case, since it’s just moving labels around).

But you can do other things with seed hacking too. Like, let’s say I want to make it so that sample() flips an imaginary coin a bunch of times and produces 10 heads in a row at the beginning. I can try it with some arbitrary seed and hope I get 10 in a row:

withr::with_seed(1234, {
  sample(c("H", "T"), 10, replace = TRUE)
})
##  [1] "T" "T" "T" "T" "H" "T" "H" "H" "H" "T"

That didn’t work. I could adjust it to 12345 or 123 or 42 or whatever and eventually find one that has 10 heads.

Or—even better—I can brute force simulate it.

Let’s take the numbers 1–10,000 and use each one as a seed for flipping a simulated coin 10 times. If it happens to result in all heads, we’ll mark it:

library(tidyverse)

possible_seeds <- data.frame(seed = 1:10000) |>
  mutate(
    all_heads = map_lgl(seed, \(s) {
      withr::with_seed(s, {
        flips <- sample(c("H", "T"), 10, replace = TRUE)
        all(flips == "H")
      })
    })
  )

possible_seeds |> filter(all_heads)
##    seed all_heads
## 1   614      TRUE
## 2  1667      TRUE
## 3  3212      TRUE
## 4  4166      TRUE
## 5  4580      TRUE
## 6  5527      TRUE
## 7  5824      TRUE
## 8  7365      TRUE
## 9  7468      TRUE
## 10 8975      TRUE

Of those 10,000 seeds, 10 of them happened to create universes where the predetermined list of pseudorandom numbers results in 10 heads at the beginning.

Let’s try one and confirm. Seed 614 actually creates 13 heads in a row (!):

withr::with_seed(614, {
  sample(c("H", "T"), 15, replace = TRUE)
})
##  [1] "H" "H" "H" "H" "H" "H" "H" "H" "H" "H" "H" "H" "H" "T" "T"

Real world bad things can happen because of pseudorandom numbers

All the examples of list recycling and seed hacking above are pretty benign, but we could theoretically use this same process of working through a ton of seeds to get to any difference in means or coefficient or p-value that we want from simulated data.

For example, Naimi, Yu, and Bodnar (2024) use several different machine learning techniques like random forests and neural networks—all of which incorporate some element of randomness—to analyze the same data and estimate the same average treatment effect. The estimated causal effect varies substantially based on the seed. Some estimates are positive and significant, some are positive and not significant, some are null, and so on. They don’t really have a solution to this variability, other than warning that researchers who use methods that depend on pseudorandom numbers need to be aware of seed dependence.⁴ They (and others) also caution that researchers should avoid seed hacking to find favorable results, but this kind of potential research misconduct is a lot less well known than things like p-hacking and HARKing, so people often forget about it.

⁴ Bao and Bindschaedler (2025) run their machine learning-based analyses across 500 randomly selected seeds, which sounds like a neat way to reduce variability.

In 2017, Canada’s government ministry responsible for immigration—Immigration, Refugees and Citizenship Canada (IRCC)—instituted a new lottery system to randomly offer permanent resident status to the parents and grandparents of Canadian citizens. However, IRCC decided to use an Excel document to allocate these slots. At the time, Excel (and Windows in general) didn’t use fancier methods like the Mersenne Twister, so its period length for pseudorandom numbers was shorter, and in 2007 researchers were able to reverse engineer Windows’s pseudorandom numbers after seeing enough values. That means that (1) people could potentially figure out all the IRCC lottery draws in advance, and, even worse, (2) due to the algorithm’s poor coverage of possible values, some people in the lottery might not ever even appear in the list of pseudorandom values due to early recycling.

In an extreme case, in the 2010s, Eddie Tipton—the IT security director of the Multi-State Lottery Association in the United States—rigged several state lotteries by modifying the pseudorandom algorithm. This allowed him to predict winning lottery numbers on specific days in a year, and he ended up winning millions of dollars across at least five different lotteries. In 2017, he was convicted and sentenced to 25 years in prison (but was released on parole after five years in 2022).

Can computers even create true randomness?

Every single “random” function in statistical software, programming languages in general, and even things like Minecraft technically generate pseudorandom values. If you know the algorithm and the current state (e.g. the value of in LCG), you can predict all the numbers in the entire series.

Internet security protocols like SSL rely on true randomness for encryption. Pseudorandom values from things like the Mersenne Twister are insecure for encryption. If dedicated hackers ever figure out or guess the current value in the algorithm (or the initial seed), they’d be able to decrypt future (and past!) internet traffic. Actual real life SSL algorithms prevent this—don’t worry!

So is there any way to make computers create true non-deterministic random values? Yes!

You just need something outside of the computer—some sort of physical process that creates lots of unpredictable entropy. There are lots of possible ways to do this:

Moving a mouse around

If you create an encryption key on your computer with something like PGP, the terminal interface will actually tell you to move your mouse around and type stuff on your computer. You’ll see a message like this:

❯ gpg --gen-key
gpg (GnuPG) 2.4.9; Copyright (C) 2025 g10 Code GmbH
This is free software: you are free to change and redistribute it.
There is NO WARRANTY, to the extent permitted by law.

...

We need to generate a lot of random bytes. It is a good idea to perform
some other action (type on the keyboard, move the mouse, utilize the
disks) during the prime generation; this gives the random number
generator a better chance to gain enough entropy.

The software records the time between mouse movements, keystrokes, files getting accessed, network traffic, and other events on the computer.⁵ It then uses that entropy to seed a pseudorandom algorithm to create public and private encryption keys. That randomness is not based solely on deterministic pseudorandom seeds—it’s based on you seeding that process by moving and doing stuff.

⁵ Actually, in practice, most operating systems maintain an “entropy pool” like Linux’s /dev/random that continuously collects time-based randomness from hardware-related events. GPG uses that pool of entropy when it creates a new key, and if it’s low, it has to fill it back up. Wild stuff.

Lava lamps

Cloudflare is one of the largest content delivery networks (CDNs) on the internet and it handles a ton of encryption-related tasks like SSL encryption (i.e. https stuff). For security reasons, they supplement their entropy sources with unpredictable physical processes, making it so their encryption isn’t reverse engineerable.

To generate actually random numbers, Cloudflare famously has a wall of ≈100 lava lamps in the lobby of their San Francisco office:

Lava Lamps at Cloudflare, licensed under Creative Commons by FASTILY

They have a camera trained on the wall that takes regular pictures of the globs of fluid. Those pictures are then turned into pixels, which are assigned numeric values, which are then used as a series of truly random numbers that are used as seeds for encryption algorithms.

Atmospheric noise

Random.org uses a set of three radios that capture atmospheric noise. It then uses the entropy of that noise to generate truly random numbers. It’s essentially like the Cloudflare lava lamps, but using radio static from the atmosphere instead.

It’s free too, and I use them all the time. They have all sorts of random-related services, like list randomizers, coin flippers, dice rollers, card shufflers, integer generators, and so on. When I sort students into teams for group projects, I use Random.org.

How I use true randomness in my own work

Using your mouse, a set of lava lamps, or the entire sky are all really neat ways to generate real randomness, but they’re slow, rate-limited,⁶ and resource-intensive. Even Cloudflare doesn’t use their lava lamps for every single random number—they use them to generate completely random seeds, which they then use in fancy cryptographically secure pseudorandom number generators.

⁶ Random.org actually has a daily quota system to make sure that people don’t use up its entropy.

R doesn’t have a way to use atmospheric noise for every random number—that would be way too resource intensive. Like, running runif(10000) would involve 10,000 API calls to Random.org, and they would probably not appreciate that.

So instead, I basically do what Cloudflare does. I use Random.org to generate a truly random value, then I use that number as my seed for whatever analysis or simulation I’m doing. I’ve been doing this for years (see here):

Search results for “From random.org” in my GitHub repositories. The seeds in that screenshot are all fairly small numbers because back when I wrote those files, I was generating only 4-digit random seeds, for whatever reason 🤷‍♂️

I actually use Raycast (more on Raycast here!) to run a little R script that grabs an integer between 100000 and 999999⁷ from Random.org and puts it in my clipboard. I then paste that into set.seed().

⁷ Visit this for an example! Though now that I’ve learned that R can use one of 4.2 billion possible seeds, I might expand this range.

Is it over the top and extra? Heck yeah. Is it better than using 1234 or 42? Yep.

get_seed.R

#!/usr/bin/env Rscript

# Required parameters:
# @raycast.schemaVersion 1
# @raycast.title Get Seed
# @raycast.mode silent

# Optional parameters:
# @raycast.packageName athcast
# @raycast.icon 🌱

# Documentation:
# @raycast.author Andrew Heiss
# @raycast.authorURL https://www.andrewheiss.com
# @raycast.description Generate a random seed from random.org

library(httr2)

seed <- request(
  "https://www.random.org/integers/?num=1&min=100000&max=999999&col=1&base=10&format=plain&rnd=new"
) |>
  req_method("GET") |>
  req_perform() |>
  resp_body_string(encoding = "UTF-8")

seed <- gsub("[\r\n]", "", seed)

clipr::write_clip(seed, object_type = "character", allow_non_interactive = TRUE)
glue::glue("{seed} copied to clipboard")

“…as an ook cometh of a litel spyr…”

I started this post with some pretentious Chaucerian epigraphs about seeds and oaks and acorns: mighty oaks from little acorns grow. Seeds contain all the genetic information and potential to grow massive trees.

Similarly, seeds for pseudorandom number generators do the same thing. Any “random” numbers that you use on a computer come from an astronomically massive deterministic list of values that are generated with some fancy math, which is kicked off with a single seed value. Each time you do something random, the computer moves along this list of theoretically already-known algorithmically generated pseudorandom values. Once you set a seed, every coin flip, every MCMC iteration, every random assignment to treatment status, every shuffled label, and every jittered point is already decided before you run a single line of code.⁸

Kind of. There are all sorts of caveats. For instance, the order of operations matters. Random functions will consume different values in the list at different rates: jittering 1,000 points and shuffling 20 labels would use the random numbers – and –, while shuffling 20 labels and jittering 1,000 points would use numbers – and –.

Like, see this for an example. Both of these end up at the same random number , but they use the first 1,020 random values in different ways, so any plots made prior to that point in the list will look different.

withr::with_seed(1234, {
  # Jitter 1000 points with X1-X1000
  # (but only show the first 5 numbers here)
  print(runif(1000)[1:5])
  # Shuffle 20 labels with X1001-X1020
  # (but only show the first 5 numbers here)
  print(runif(20)[1:5])
  # This should be X1021
  print(runif(1))
})
## [1] 0.1137 0.6223 0.6093 0.6234 0.8609
## [1] 0.8376 0.4875 0.1103 0.3514 0.7611
## [1] 0.1454

withr::with_seed(1234, {
  # Shuffle 20 labels with X1-X20
  # (but only show the first 5 numbers here)
  print(runif(20)[1:5])
  # Jitter 1000 points with X21-X1020
  # (but only show the first 5 numbers here)
  print(runif(1000)[1:5])
  # This should be X1021
  print(runif(1))
})
## [1] 0.1137 0.6223 0.6093 0.6234 0.8609
## [1] 0.3166 0.3027 0.1590 0.0400 0.2188
## [1] 0.1454

Practically speaking, does any of this matter? Yes! Limiting yourself to the same seeds everyone else is using (1, 1234, 42, etc) means that you’re exploring the same universes everyone else is exploring. Use big, uncommon seeds so you’re not sharing a universe with every other researcher.

The stakes of this go beyond jittered points and MCMC chains. As we saw, deterministic randomness has real world consequences like multimillion dollar fraud and faulty immigration lotteries. If you’re doing real world random work—assigning actual study participants to experimental treatments, running any sort of lottery, etc.—make sure you use good randomization with natural entropy.

Inject some entropy into your own analyses by grabbing a seed from Random.org—you’ll still have a big list of predetermined pseudorandom numbers, but that list starts with a process that you have no control over.

References

Bao, Wenxuan, and Vincent Bindschaedler. 2025. “Towards Reliable and Generalizable Differentially Private Machine Learning.” arXiv. https://doi.org/10.48550/arXiv.2508.15141.

Johnson, Alicia A., Miles Q. Ott, and Mine Dogucu. 2022. Bayes Rules!: An Introduction to Applied Bayesian Modeling. 1st ed. Boca Raton: Chapman and Hall/CRC. https://doi.org/10.1201/9780429288340.

Matsumoto, Makoto, and Takuji Nishimura. 1998. “Mersenne Twister: A 623-Dimensionally Equidistributed Uniform Pseudo-Random Number Generator.” ACM Transactions on Modeling and Computer Simulation 8 (1): 3–30. https://doi.org/10.1145/272991.272995.

Naimi, Ashley I., Ya-Hui Yu, and Lisa M. Bodnar. 2024. “Pseudo-Random Number Generator Influences on Average Treatment Effect Estimates Obtained with Machine Learning.” Epidemiology 35 (6): 779–86. https://doi.org/10.1097/EDE.0000000000001785.

Spiegelhalter, David J. 2025. The art of uncertainty: How to navigate chance, ignorance, risk and luck. New York, NY: W.W. Norton & Company.

Citation

BibTeX citation:

@online{heiss2026,
  author = {Heiss, Andrew},
  title = {Generating Universes Within Universes with a Single Seed},
  date = {2026-04-13},
  url = {https://www.andrewheiss.com/blog/2026/04/13/seeds-predetermined-universes/},
  doi = {10.59350/1z4vr-nss70},
  langid = {en}
}

For attribution, please cite this work as:

Heiss, Andrew. 2026. “Generating Universes Within Universes with a Single Seed.” April 13, 2026. https://doi.org/10.59350/1z4vr-nss70.

How to automatically convert TikZ images to SVG with Quarto

Andrew Heiss — Wed, 25 Feb 2026 05:00:00 GMT

Back in 2021 I wrote a little guide to how to use tikz chunks in R Markdown documents. {knitr} supports tikz natively, and with some some shenanigans and special settings—and after recompiling LaTeX with one extra feature and setting a new environment variable on macOS—it’s possible to have tikz chunks automatically knit to PDF in PDFs and SVG in HTML and Word. While convoluted, it’s a neat process and I’ve used it for several published papers (like, see Figures A2 and A3 here—those are built with tikz!).

While that process works, I’ve also found it really annoying, and for years I’ve looked for better ways. The very first pre-release version of Quarto was also in 2021, but it was too young to do much with languages beyond R and Python. Since then, though, Quarto has added native support for two different diagram languages: Mermaid and Graphviz. They work well for lots of diagrams, but not for most of the diagrams I make, where (1) I want exact control over where items go, and (2) I want them to look nice and visually fit with the rest of the document.

For instance, let’s say I want to include a minimalist little standard directed acyclic graph (DAG) that shows that some treatment causes some outcome in the subsequent time period, but the effect is confounded by a third variable . Something like this:

A simple DAG

Quick, easy, ugly diagrams with Mermaid and Graphviz

We can create this DAG with Mermaid, but there are limited styling options and there’s not really any good way to control the layout, so we’re pretty much stuck with this:

```{mermaid}
graph LR
  X([X]) --> Y([Y])
  Z([Z]) --> X([X])
  Z([Z]) --> Y([Y])
```

graph LR
  X([X]) --> Y([Y])
  Z([Z]) --> X([X])
  Z([Z]) --> Y([Y])

We can also use Graphviz to make the same DAG, and it has a few extra options for layout—like here, we tell it that the Z node is the focal point of the diagram and that X and Y should be pushed down:

```{dot}
digraph {
  node [shape=circle]
    
  { rank=source; Z }
  { rank=sink; X; Y }
    
  X -> Y
  Z -> X
  Z -> Y
}
```

That’s all fine, but what if we want, like, control over these diagrams? Or math? Or typographic things like subscripts or superscripts? Good luck.

Nicer diagrams with TikZ

TikZ is very good at making diagrams. It handles math. It handles fancy formatting. It’s also absolutely miserable to work with—you have to draw diagrams with a kinda-sorta-dialect of LaTeX. (Making initial TikZ diagrams is one area where LLMs are incredibly useful!). TikZ also converts exclusively to PDF,¹ but you can convert those diagrams from PDF to other formats like SVG or PNG.

¹ Technically DVI, but whatever.

Here’s that same DAG as TikZ:

\begin{tikzpicture}[>={stealth}]
  \node (x) at (0,0) {$X_{t}$};
  \node (y) at (2,0) {$Y_{t + 1}$};
  \node (z) at (1,1) {$Z$};
  \path[->] (z) edge (x);
  \path[->] (z) edge (y);
  \path[->] (x) edge (y);
\end{tikzpicture}

However, Quarto won’t render this as an image because it requires feeding it through LaTeX. As I show here, it’s possible to use R Markdown with {knitr} in a special tikz chunk, which would appear as a nice diagram when rendering to PDF. With a bunch of extra setup, you can also get tikz chunks to appear as SVG when rendering to non-PDF formats.

My guide still works with Quarto and you can just do that if you want to use TikZ chunks! I’ve been doing it for years! But Quarto re-renders and re-converts each diagram every time you render the document, which is excessive and slow. Getting all the special scaffolding set up is annoying too. It only works with R-based documents. You have to reinstall LaTeX with special extra settings on macOS. You have to create a special chunk template and then hook that into the tikz engine. And you have to do that every time. Oof.

There’s an easier way though!

Converting TikZ to both PDF and SVG with Quarto and `diagram`

The pandoc team has an extension called diagram that is essentially a fancy Lua filter that will take code for several different diagram languages (including TikZ) and convert them into images.² In 2023 they reorganized it to work as a Quarto extension, so installing it and using it is nice and straightforward. It’s language agnostic—you don’t have to use it only in an R document. You don’t have to recompile LaTeX on macOS or set any systemwide environment variables. It Just Works™.

² There are two other Quarto extensions/filters that do something similar (imagify and quarto_tikz). I’m sure they work great too—albeit differently from diagrams. I super trust the pandoc people, so I’m happy to stick with diagrams.

Here’s what you have to do:

Download and install Inkscape and tell your terminal that it exists. The diagram extension uses Inkscape, an open source vector editor (like Adobe Illustrator, but free) to handle the PDF → SVG conversion. This is the only external piece of software you need to get this process to work, and you never really need to open it. You’ll install it and set it up once—after that, things will work normally.
TipPutting inkscape in your PATH
Inkscape can be used from the terminal using the inkscape command, but that program lives in different places depending on your operating system. The diagram extension assumes that you can access it by running inkscape in your terminal, but this likely won’t work immediately.
- For macOS: On macOS, inkscape lives in /Applications/Inkscape.app/Contents/MacOS/inkscape. You can symlink it to a more accessible location with this terminal command:
  
  ln -s /Applications/Inkscape.app/Contents/MacOS/inkscape /usr/local/bin/inkscape
- For Windows: On Windows, inkscape.exe lives in the bin folder of wherever you installed Inkscape, likely somewhere like C:\Program Files\Inkscape\bin. You’ll need to add that to your PATH. You can do that by going to System Properties > Environment Variables > and append C:\Program Files\Inkscape\bin to the PATH variable, or you can run this in PowerShell:
  
  [Environment]::SetEnvironmentVariable("Path", $env:Path + ";C:\Program Files\Inkscape\bin", "Machine")
- For both: Once you’ve symlinked the program or added it to your path, you should be able to open a new terminal window and run this:
  
  inkscape --version
Install the diagram extension in a Quarto project. Create a new Quarto project. From the terminal, run this to install the diagram extension:
```
quarto install extension pandoc-ext/diagram
```
Tell your document to use the extension. Add this to the YAML header of a Quarto document:
```
---
filters:
  - diagram
---
```

Add a .tikz block to your document. Add a tikz code block to your document like this:

```{.tikz}
\begin{tikzpicture}[>={stealth}]
  \node (x) at (0,0) {$X_{t}$};
  \node (y) at (2,0) {$Y_{t + 1}$};
  \node (z) at (1,1) {$Z$};
  \path[->] (z) edge (x);
  \path[->] (z) edge (y);
  \path[->] (x) edge (y);
\end{tikzpicture}
```

Render! Render the document to PDF and you should see a diagram. Render it to HTML or .docx or Typst and it should use inkscape to convert it to SVG. Magic.

That’s it! No weird LaTeX recompilation. No special R-only knitr templates and engine overrides. Quarto, pandoc, Lua, and Inkscape take care of everything automatically.

Extra settings and complete example

You can configure lots of other settings too—see the TikZ section of the documentation at sample.md here or the extension’s test file for tikz for more details.

For example, you can control the LaTeX engine that gets used, add LaTeX and TikZ packages and other preamble-y things, and enable caching so that figures don’t need to be recompiled every time. You can even use native Quarto cross referencing and captions, but you need to use its cross-reference div syntax.

Here’s a complete mini example. Because we enabled caching, it won’t rebuild and re-convert each time you render, and because we gave it a filename, there will be a standalone SVG file in example_files/mediabag/a-neat-dag.svg. It also defines some custom colors and fonts.

example.qmd

---
title: "Some neat document"

filters:
  - diagram

diagram:
  cache: true
  cache-dir: ./tikz-cache
  engine:
    tikz:
      execpath: lualatex
      header-includes:
        - '\usepackage{libertine}'
        - '\usepackage{libertinust1math}'
        - '\definecolor{xcolor}{HTML}{d7433b}'
        - '\definecolor{ycolor}{HTML}{008d98}'
        - '\definecolor{zcolor}{HTML}{ffcc3d}'
---

See @fig-dag-thing for a neat example!

::: {#fig-dag-thing}

```{.tikz}
%%| filename: a-neat-dag
%%| width: 50%

\begin{tikzpicture}[>={stealth}]
  \node [xcolor] (x) at (0,0) {$X_{t}$};
  \node [ycolor] (y) at (2,0) {$Y_{t + 1}$};
  \node [zcolor] (z) at (1,1) {$Z$};
  \path[->] (z) edge (x);
  \path[->] (z) edge (y);
  \path[->] (x) edge (y);
\end{tikzpicture}
```

Caption for the DAG

:::

That will create a diagram like this:

Figure 1: Caption for the DAG

This process works with any TikZ diagram (these came from Christian Testa’s tikz examples). There’s no way you can do stuff like this with Mermaid or Graphviz :)

(a) A cone with a height and radius (original)

(b) Definition of an intersection (original)

Figure 2: Some example TikZ diagrams

Here’s the Markdown to make Figure 2 (a) and Figure 2 (b) in Figure 2:

example.qmd

Here's the Markdown to make @fig-cone and @fig-set in @fig-tikz-examples:

:::: {#fig-tikz-examples layout-ncol=2 layout-valign="bottom"}

::: {#fig-cone}

```{.tikz}
%%| filename: cone
%%| width: 100%

\usetikzlibrary{calc, backgrounds}

\begin{tikzpicture}[scale=.75,
  background rectangle/.style={fill=white},
  show background rectangle]
  \draw[dashed] (0,0) arc (170:10:2cm and 0.4cm)coordinate[pos=0] (a);
  \draw (0,0) arc (-170:-10:2cm and 0.4cm) coordinate (b);
  \draw[densely dashed] 
      ([yshift=4cm]$(a)!0.5!(b)$) -- 
      node[right,font=\footnotesize] {1} coordinate[pos=0.95] (aa)($(a)!0.5!(b)$) -- 
      node[above,font=\footnotesize] {1} coordinate[pos=0.1] (bb) (b);
  \draw (aa) -| (bb);
  \draw (a) -- ([yshift=4cm]$(a)!0.5!(b)$) -- (b);
\end{tikzpicture}
```

A cone with a height and radius ([original](https://github.com/ctesta01/tikz-examples/blob/main/2023/cone/cone.tex))

:::

::: {#fig-set}

```{.tikz}
%%| filename: set-belonging
%%| width: 100%

\usetikzlibrary{backgrounds}

\begin{tikzpicture}[scale=1.5, background rectangle/.style={fill=white}, show background rectangle]
  % Draw the universe rectangle
  \draw[thick] (-2,-1.75) rectangle (2,1.75) node[below left] {\(\Omega\)};

  % Draw the fill for circle A and label it
  \draw[thick, fill=blue!20, draw opacity = .65] (-0.5,0) circle (1);
  \node at (-1, 0.2) {\(A\)};

  % Draw circle B and label it
  \draw[thick] (0.5,0) circle (1);
  \draw[thick, fill=red!20, opacity = .65] (0.5,0) circle (1);

  % Draw the outer circle for A so it appears "on top of" the intersection
  \node at (1,0.2) {\(B\)};

  %
  \draw[thick] (-0.5,0) circle (1);

  % Label the intersection (A \cap B)
  \node at (0,0) {\(A \cap B\)};
\end{tikzpicture}
```

Definition of an intersection ([original](https://github.com/ctesta01/tikz-examples/blob/main/2023/intersection/intersection.tex))

:::

Some example TikZ diagrams

::::

Citation

BibTeX citation:

@online{heiss2026,
  author = {Heiss, Andrew},
  title = {How to Automatically Convert {TikZ} Images to {SVG} with
    {Quarto}},
  date = {2026-02-25},
  url = {https://www.andrewheiss.com/blog/2026/02/25/tikz-quarto-svg-fun/},
  doi = {10.59350/5m018-qax35},
  langid = {en}
}

For attribution, please cite this work as:

Heiss, Andrew. 2026. “How to Automatically Convert TikZ Images to SVG with Quarto.” February 25, 2026. https://doi.org/10.59350/5m018-qax35.

How to make your data analysis life easier using Positron, Raycast, and Espanso

Andrew Heiss — Tue, 13 Jan 2026 05:00:00 GMT

On January 13, 2026, I was on Posit’s Data Science Lab to talk about my Positron settings and other neat workflow-y things I use to make my life easier.

This isn’t a standard blog post—it’s mostly just a list of the stuff we talked about, with lots of different links to other resources.

Positron stuff

Favorite Positron extensions:

My custom Positron settings (these all live in settings.json):

Click to expand

How to open a folder as a Positron project with macOS Quick Actions

Andrew Heiss — Tue, 22 Jul 2025 04:00:00 GMT

tl;dr

Download and unzip this:

“Open in Positron” Quick Action

Put Open in Positron.workflow in ~/Library/Services.

Now you can right click on folders in Finder, go to Quick Actions, and select “Open in Positron” to open folders as projects in Positron.

One thing I’m still adjusting to in Positron is how it treats projects and folders. RStudio makes it easy to follow a project-oriented workflow, and project-specific settings get set in a special Rproj file. When you open an Rproj file from Finder or Windows Explorer, a new RStudio instance gets opened and the working directory is automatically set to the project directory. It’s portable across systems (and within your own computer), and you never have to worry about running setwd(). It’s the ideal workflow.

Projects in Positron

Positron doesn’t have the concept of an Rproj file, but it most definitely still works with project-based workflows. You can open a project by dragging a folder onto the Positron icon in the Dock (on macOS; Windows people—you do… something else?), and the whole folder opens in Positron, with the working directory automatically set there. You can have as many folders/projects open as you want, and instead of getting a bunch of separate RStudio icons in the Dock, you get one Positron icon with a ton of windows:

Lots of Positron projects open all at once

Positron also has a helpful project switcher menu in the top right corner, just like RStudio:

Positron’s project menu

You can also use ctrl+r in Positron to see a list of recently opened projects/folders:

Recent Positron projects with ctrl+r

That’s all great and wonderful and fine.

Quickly opening folders in Positron

My main source of friction with Positron-based projects (and this applies to Visual Studio Code too) is quickly opening a folder as a project from Finder. Since there’s no concept of an Rproj, I can’t just double click on a file to open Positron and have it pointed at the right place. I have to drag it down to the Dock.

Well, kind of. Positron has a few recommended ways to speed up the folder opening process.

Workspaces

Positron/Visual Studio Code has the idea of workspaces. You can save a folder (or a collection of folders) as a workspace to store extra project-specific settings. The Positron documentation has more details about this too.

These are a little weird compared to Rproj files. An Rproj file is a plain text file that contains project-specific settings and directs RStudio to set the working directory to its location. With Positron/Visual Studio Code workspaces, the project-specific settings are stored in .vscode/settings.json, but there’s nothing to tell Positron to open a project. The closest analogue is to go to File > Save Workspace As… in Positron and save a .code-workspace file. This is technically called a “Multi-root workspace” and allows you to include multiple folders in the same workspace, but it also works for a single folder. I do this for my website here—I have an ath-quarto.Rproj file for opening this project in RStudio and an ath-quarto.code-workspace file that I can double click on to open the project in Positron:

A code-workspace file

The one downside to this is that Positron treats multi-root workspaces differently visually—it includes an extra “(Workspace)” suffix in the window title, and it displays as a file and not a folder in the recently opened projects area in macOS’s Mission Control:

Different icon styles in Mission Control

Using multi-root workspaces for single-directory projects works fine—it’s just a little annoying to have these weirdly superpowered workspaces that get treated specially. But it’s also convenient if you want to have a file that you can open like an Rproj. I actually have a little folder in my Dock that lets me open my website either in RStudio or in Positron, which is only possible because of the code-workspace file:

Pinned Dock folder with shortcuts to an Rproj file and a code-workspace file

Project Manager + Raycast

Another alternative that Positron recommends is to use the Project Manager extension in conjunction with Raycast (see here for more). This is super neat—once you set it up, you can add a keyboard shortcut in Raycast to open a list of all your saved projects. For me, I can mash ⌘⌥^⇧P and open any of the projects that I’ve configured in Project Manager:¹

¹ Though note here that the saved code-workspace projects like “ath-quarto” and “2025 playground” get different icons, and macOS has weirdly decided to show those as Skype files???—yet another minorly annoying downside to the multi-root code-workspace approach.

Positron projects in Raycast

The one downside to this is that I have to remember to manually add a project to the Project Manager extension, which is extra friction and something I always forget to do.

Open a folder in Positron macOS Quick Actions

On Windows, the Positron installer provides an option to add an “Open with Positron” action to the Windows Explorer directory context menu so that when you right click on a folder, it gives an option to open the folder in Positron. That’s neat and easy.

macOS doesn’t have the exact equivalent and restricts what can be added to the right click context menu. But you can actually add special automated workflows as something called either Services or Quick Actions to Finder’s right click context menu.

Services vs. Quick Actions

Services have been part of macOS for a long time. They can be made with Automator, or incorporated directly in applications. Quick Actions are newer (though they’ve been around since at least Mojave in 2018) and are incorporated directly in Finder. I think they can only be made through Automator(??), though maybe applications can provide them directly as well(????)

It’s fairly straightforward to create your own Quick Action workflow:

Open Automator. It’s already installed on your computer; everybody always forgets this exists.
Go to File > New and select “Quick Action”
Create a workflow with all these settings (a screenshot of the finished workflow is included below):
1. Set “Workflow receives current” to “folders” in Finder²
2. Drag the “Run Shell Script” action from the Actions Library to the workflow
3. Set “Pass input” to “as arguments”
4. Paste this script into the script area:
```
for f in "$@"
do
    open -a "Positron" "$f"
done
```
Save it as “Open in Positron” (or whatever else you want)

² You could also set this to “files or folders” to add the Quick Action to all files too, but I find this excessive.

Here’s what the final thing should look like:

Finished Quick Action workflow

Once you’ve saved the workflow, it will live at ~/Library/Services/Open in Positron.workflow.

Now when you right click on a folder in Finder, you can go to Quick Actions > Open in Positron and the folder will open in Positron as a project!

Context menu for opening a folder in Positron

If you use the Gallery View in Finder (I never do—Column View forever), you’ll also get a neat little “Open in Positron” button when you look at the details of a folder:

“Open in Positron” button in Gallery View

This is so helpful! It stinks that “Open in Positron” is buried under one extra level of menus (Quick Actions), but having the ability to open a project from Finder without needing to drag things to the Dock is really great.

Making this native?

Theoretically, I think this extra Service/Quick Action thing could be bundled into Positron itself. Lots of macOS apps add stuff to the Service menu. Like, I can open a new iTerm2 tab to this folder, or open the folder in DiffMerge, or encrypt/decrypt files with OpenPGP. Those options are all available because I have iTerm2, DiffMerge, and GPG Suite installed.

Lots of services in Finder

These services are defined in each app’s respective .app/Contents/Info.plist file (see Apple’s official documentation here, or this example). For instance, DiffMerge has this entry in its Info.plist:

Info.plist








<key>NSServiceskey>
<array>
  <dict>
    <key>NSMenuItemkey>
    <dict>
      <key>defaultkey>
      <string>Open in DiffMergestring>
    dict>
    <key>NSMessagekey>
    <string>openFilesViaServicestring>
    <key>NSUserDatakey>
    <string>SGDMstring>
    <key>NSPortNamekey>
    <string>DiffMergestring>
    <key>NSRequiredContextkey>
    <dict>
      <key>NSTextContentkey>
      <string>FilePathstring>
    dict>
    <key>NSSendTypeskey>
    <array>
      <string>NSFilenamesPboardTypestring>
    array>

    <key>NSReturnTypeskey>
    <array/>
  dict>
array>

BUT I don’t know anything about official macOS development and how to create plist files like this. Theoretically it’s possible though!

Citation

BibTeX citation:

@online{heiss2025,
  author = {Heiss, Andrew},
  title = {How to Open a Folder as a {Positron} Project with {macOS}
    {Quick} {Actions}},
  date = {2025-07-22},
  url = {https://www.andrewheiss.com/blog/2025/07/22/positron-open-with-finder/},
  doi = {10.59350/pt66c-57w73},
  langid = {en}
}

For attribution, please cite this work as:

Heiss, Andrew. 2025. “How to Open a Folder as a Positron Project with macOS Quick Actions.” July 22, 2025. https://doi.org/10.59350/pt66c-57w73.

How to use Positron’s Connections Pane with DuckDB

Andrew Heiss — Thu, 10 Jul 2025 04:00:00 GMT

One more Positron-related post! It’s a quick one, just to highlight one feature I think is really neat and helpful: the Connections Pane.

In a newer research project I’m working on, I have geocoded data for every foreign aid project run by most donor countries since 1989. This involves millions of rows and nearly a thousand columns across three CSV files, totalling 3.3 GB. R can load this data, but it takes a long time and slows down the session substantially (and often crashes R entirely).

Fortunately there’s a good solution for this huge data. DuckDB is a lightning fast database system that works really well with R. The DuckDB file that I’ve made for the foreign aid data compresses things substantially—it only takes up 600 MB on the disk instead of 3.3 GB. It also lets me query the huge data using either SQL commands or with {dplyr} functions, and I can load subsets of the huge data into R really quickly.

DuckDB, {DBI}, and the difficulty of discerning data in a database

Here’s a quick little example showing the general process with some toy data. I’ll create a little .duckdb file on the disk named example_database.duckdb—I could also make a temporary in-memory database by naming the file ":memory:":

library(dplyr)

db_file <- "example_database.duckdb"

# Connect to the database
con <- DBI::dbConnect(duckdb::duckdb(), db_file)

Next I’ll add some stuff to it:

# Add stuff to the database
copy_to(con, penguins, name = "penguins", overwrite = TRUE, temporary = FALSE)
copy_to(con, gapminder::gapminder, name = "gapminder", overwrite = TRUE, temporary = FALSE)
copy_to(con, ggplot2::mpg, name = "mpg", overwrite = TRUE, temporary = FALSE)
copy_to(con, ggplot2::diamonds, name = "diamonds", overwrite = TRUE, temporary = FALSE)

And then I’ll close the connection:

DBI::dbDisconnect(con)

That worked, but it’s really hard to see if it did anything. The database file does contain stuff—it’s 1.6 MB now:

❯ ls -lh example_database.duckdb
# Permissions Size User   Date Modified    Name
# .rw-r--r--@ 1.6M andrew 2025-07-10 12:13 example_database.duckdb

…and I can connect to the database with R and get data out of it:

con <- DBI::dbConnect(duckdb::duckdb(), db_file)

# Get stuff out of the database
adelie_query <- tbl(con, I("penguins")) |> 
  filter(species == "Adelie")

# Check it out! It's SQL!
show_query(adelie_query)
## 
## SELECT q01.*
## FROM (FROM penguins) q01
## WHERE (species = 'Adelie')

# Actually run the query
penguins_from_db <- adelie_query |> 
  collect()
penguins_from_db
## # A tibble: 152 × 8
##    species island    bill_len bill_dep flipper_len body_mass sex     year
##                                  
##  1 Adelie  Torgersen     39.1     18.7         181      3750 male    2007
##  2 Adelie  Torgersen     39.5     17.4         186      3800 female  2007
##  3 Adelie  Torgersen     40.3     18           195      3250 female  2007
##  4 Adelie  Torgersen     NA       NA            NA        NA     2007
##  5 Adelie  Torgersen     36.7     19.3         193      3450 female  2007
##  6 Adelie  Torgersen     39.3     20.6         190      3650 male    2007
##  7 Adelie  Torgersen     38.9     17.8         181      3625 female  2007
##  8 Adelie  Torgersen     39.2     19.6         195      4675 male    2007
##  9 Adelie  Torgersen     34.1     18.1         193      3475     2007
## 10 Adelie  Torgersen     42       20.2         190      4250     2007
## # ℹ 142 more rows

DBI::dbDisconnect(con)

…but there’s not an easy way to check what’s actually going on in there. I only know there’s a table named penguins in there because I just made it.

DuckDB provides a helpful browser-based UI that you can load from the terminal like this:

❯ duckdb -ui example_database.duckdb

This gives you a notebook where you can write SQL commands, which I find less helpful because I typically interact with the database with {dplyr}—I really only use DuckDB for fast data storage and fast data access. Most importantly, though, it gives you a sidebar that shows all the tables inside the database, and you can preview the stuff inside:

The DuckDB browser-based UI

This all works fine for remembering what’s in the database, BUT it’s always a bit of a hassle because you can’t run commands from R as long as the the UI is open. You can only connect to a database once, so you’ll get this error:

con <- DBI::dbConnect(duckdb::duckdb(), db_file)
#> Error in `.local()`:
#> ! rapi_startup: Failed to open database: {"exception_type":"IO","exception_message":"Could not 
#> set lock on file \"example_database.duckdb\": Conflicting lock is held in /opt/homebrew/Cellar/
#> duckdb/1.2.2/bin/duckdb (PID 48019) by user andrew. 
#> See also https://duckdb.org/docs/stable/connect/concurrency","errno":"35"}

Disconnect when done!

When you’re done with the browser UI, you need to type .quit from the DuckDB terminal to close the connection:

❯ duckdb -ui example_database.duckdb
┌──────────────────────────────────────┐
│                result                │
│               varchar                │
├──────────────────────────────────────┤
│ UI started at http://localhost:4213/ │
└──────────────────────────────────────┘
v1.2.2 7c039464e4
Enter ".help" for usage hints.
D .quit

So you have to go back to the terminal, run .quit, and then go back to R and try again. If you want to check what’s in the database later, you have to disconnect from R, go to the terminal, run duckdb -ui, explore around in the browser, run .quit again, and go back to R.

DuckDB, {connections}, and the magical Connections Pane

There’s a better way to do this though! RStudio and Positron both have a Connections Pane that lets you save and explore database connections, and it magically fixes the problem with concurrent connections to a database so you can explore a database from inside your R session without needing to remember to turn off the in-browser GUI.

To make this work, you need to use the {connections} package, which is mostly a wrapper around standard {DBI} connection functions.

If you connect to the database with {connections}…

# Regular DBI way---don't do this
# con <- DBI::dbConnect(duckdb::duckdb(), db_file)

# Fancy {connections} way---do this. It connects *and* adds the connection to the sidebar
library(connections)
con <- connection_open(duckdb::duckdb(), db_file)

…you’ll automatically see a new connection in Positron’s Connections Pane. Click on the little → arrow icon to open it:

Basic DuckDB connection

…and you’ll see all the tables inside!

Exploring DuckDB tables

You can explore the data with Positron’s Data Explorer if you click on the little table icon:

Looking at the DuckDB table with Postrion’s data explorer

And you can work with the database with R too, without running into concurrency issues:

# Get stuff out of the database
tbl(con, I("penguins")) |> 
  filter(species == "Adelie") |> 
  collect()
## # A tibble: 152 × 8
##    species island    bill_len bill_dep flipper_len body_mass sex     year
##                                  
##  1 Adelie  Torgersen     39.1     18.7         181      3750 male    2007
##  2 Adelie  Torgersen     39.5     17.4         186      3800 female  2007
##  3 Adelie  Torgersen     40.3     18           195      3250 female  2007
##  4 Adelie  Torgersen     NA       NA            NA        NA     2007
##  5 Adelie  Torgersen     36.7     19.3         193      3450 female  2007
##  6 Adelie  Torgersen     39.3     20.6         190      3650 male    2007
##  7 Adelie  Torgersen     38.9     17.8         181      3625 female  2007
##  8 Adelie  Torgersen     39.2     19.6         195      4675 male    2007
##  9 Adelie  Torgersen     34.1     18.1         193      3475     2007
## 10 Adelie  Torgersen     42       20.2         190      4250     2007
## # ℹ 142 more rows

When you’re all done, disconnect from the database:

# Regular DBI way---don't do this
# DBI::dbDisconnect(con)

# Fancy {connections} way---do this
connections::connection_close(con)

Bonus: Better support for DuckDB in the Connections Pane

That Connections Pane connection is neat, but it’s missing a little bit of extra metadata (like the location of the file), and it shows the columns using generic output, not with DuckDB-specific column types.

DuckDB has fancier support for the Connections Pane, but for whatever reason, as of July 2025, it’s not enabled by default. If you set an option before connecting, though, you’ll get a prettier DuckDB-specific connection with a nice icon, the file location, and more specific column types.

options("duckdb.enable_rstudio_connection_pane" = TRUE)
con <- connections::connection_open(duckdb::duckdb(), db_file)

Nicer DuckDB connection

Nicer list of DuckDB tables

The whole game

Here’s what the overall process looks like for connecting to a database, adding data to it, querying it, closing the connection, and making a plot:

library(tidyverse)

# Use nicer DuckDB Connections Pane features
options("duckdb.enable_rstudio_connection_pane" = TRUE)

# Connect to an in-memory database, just for illustration
con <- connections::connection_open(duckdb::duckdb(), ":memory:")

# Add stuff to it
copy_to(
  con,
  gapminder::gapminder,
  name = "gapminder",
  overwrite = TRUE,
  temporary = FALSE
)
## # Source:   table [?? x 6]
## # Database: DuckDB v1.3.2 [root@Darwin 24.5.0:R 4.5.0/:memory:]
##    country     continent  year lifeExp      pop gdpPercap
##                            
##  1 Afghanistan Asia       1952    28.8  8425333      779.
##  2 Afghanistan Asia       1957    30.3  9240934      821.
##  3 Afghanistan Asia       1962    32.0 10267083      853.
##  4 Afghanistan Asia       1967    34.0 11537966      836.
##  5 Afghanistan Asia       1972    36.1 13079460      740.
##  6 Afghanistan Asia       1977    38.4 14880372      786.
##  7 Afghanistan Asia       1982    39.9 12881816      978.
##  8 Afghanistan Asia       1987    40.8 13867957      852.
##  9 Afghanistan Asia       1992    41.7 16317921      649.
## 10 Afghanistan Asia       1997    41.8 22227415      635.
## # ℹ more rows

# Get stuff out of it
gapminder_2007 <- tbl(con, I("gapminder")) |>
  filter(year == 2007) |>
  collect()

# All done
connections::connection_close(con)

# Make a pretty plot, just for fun
ggplot(gapminder_2007, aes(x = gdpPercap, y = lifeExp)) +
  geom_point(aes(color = continent)) +
  scale_x_log10(labels = scales::label_dollar(accuracy = 1)) +
  scale_color_brewer(palette = "Set1") +
  labs(
    x = "GDP per capita",
    y = "Life expectancy",
    color = NULL,
    title = "This data came from a DuckDB database!"
  ) +
  theme_minimal(base_family = "Roboto Condensed")

It works, and it shows the in-memory database in the Connections Pane!

Exploring tables in an in-memory database

Citation

BibTeX citation:

@online{heiss2025,
  author = {Heiss, Andrew},
  title = {How to Use {Positron’s} {Connections} {Pane} with {DuckDB}},
  date = {2025-07-10},
  url = {https://www.andrewheiss.com/blog/2025/07/10/positron-database-connections/},
  doi = {10.59350/w37d8-vj489},
  langid = {en}
}

For attribution, please cite this work as:

Heiss, Andrew. 2025. “How to Use Positron’s Connections Pane with DuckDB.” July 10, 2025. https://doi.org/10.59350/w37d8-vj489.

Use Positron to run R inside a Docker image through SSH

Andrew Heiss — Sat, 05 Jul 2025 04:00:00 GMT

I’ve long been a proponent of making quantitative research reproducible. It’s the main reason I do all my scientific writing in Quarto—I can mix code and text in the same document so I don’t need to copy/paste numbers, tables, and figures from some statistical program into a word processor. Everything automatically ends up one compiled document based on the most current data.

There’s a growing trend where academic journals require that they be able to reproduce your results on their end before publication. For most articles in these journals, this involves having a PhD student run a bunch of Stata .do files, look through the logs and output, and ensure all the numbers and figures match what’s in the manuscript.

With a couple of my recent papers where the journal required this, I handed them a {targets} pipeline. The research assistant ran one line of R code (targets::tar_make()) to compile all the raw data, clean it, run several {brms} Bayesian models, and generate HTML, Word, and PDF versions of the manuscript through Quarto with all the latest values and tables and figures. Overkill, yes. But it’s so cool.

But getting the computing environment set up for this to work is difficult. Despite progress with {renv}, it’s still fairly tricky to get project-specific environments with both R and the project packages installed at specific versions. Unless you use something like rig, R by default uses the same installation systemwide, so once you upgrade from from something like R 4.4.3 to 4.5.0, you’ll use 4.5.0 for everything. {renv} lets you install packages at specific versions within folders, but even then it can often run into dependency hell if the lockfile gets out of sync or packages are removed from CRAN or GitHub (though Posit Package Manager aleviates a lot of this). It’s even trickier if you’re using other tools beyond R in your analysis—try installing Stan and CmdStanR in a cross platform way, with all the extra gcc compilers and Makefile settings!

In one recent project, I wrote up detailed instructions for installing the computing environment with {renv}, but it took a ton of back-and-forth emailing with the editorial office because their RA used an older version of Windows and they couldn’t get several parts of the environment installed.

To make life easier for everyone, I decided the best option for maximum cross platform reproducibility was to essentially give the journal’s RA a whole virtual computer with everything installed in it, all through the magic of Docker. Anyone can reproduce the results by following these instructions—in short, this uses Docker Compose to (1) build a Linux computer with R locked at version 4.4.1, thanks to the Rocker Project, (2) install all the version-locked R packages, which takes ≈15–20 minutes, (3) run a {targets} pipeline, which takes ≈20 minutes, and (4) spit out a final HTML/PDF/Word version of the manuscript and an accompanying analysis notebook.

Again, probably (definitely) overkill. But the RA at that journal was able to reproduce everything right away with Docker, and RAs at other journals have since done it flawlessly without a ton of back-and-forth troubleshooting.

But this post isn’t about that. Installing and configuring everything to get Docker Compose + rocker/tidyverse + {renv} + CmdStanR + Quarto working nicely is beyond the scope of this. It’s (unfortunately) left as an exercise to the reader. But I have heavily commented and documented examples if you want to try it! See the Dockerfile and docker-compose.yml here and the README for the project for a complete example.

This post is a lot simpler. It’s about how to actually use R within these containers and how the advent of Positron makes it way easier to write code in a regular IDE instead of inside a browser-based instance of RStudio.

How to use Docker-based R normally

Many (most?) of the Rocker R images come with RStudio Server built-in so that you can work with the container-based version of R directly from your browser.

Download this repository for a quick little example. It runs this Docker Compose file, which grabs the R 4.5.0 rocker/tidyverse, makes and RStudio Server available at http://localhost:8787, and mounts the main project folder into the container so that it’s mirrored and any file changes are persistent.

docker-compose-basic.yml

services:
  rstudio:
    platform: linux/amd64  # Need to include this when building on Silicon macOS
    image: rocker/tidyverse:4.5.0
    ports:
      - "8787:8787"  # Make localhost:8787 in the container accessible as localhost:8787 on the computer
    environment:
      - DISABLE_AUTH=true  # Allow password-less login to RStudio server instance
    volumes:
      - ./project:/home/rstudio/project  # Mount the local project folder into the container

Here’s how to get started:

Install Docker Desktop on your computer (instructions for macOS or Windows).
Install either (1) Positron and the Container Tools extension or (2) Visual Studio Code and its Container Tools extension. This is optional, since you can do all this Docker work with the terminal, but it makes working with Docker and Docker Compose way easier.
Open the example folder in Positron or Visual Studio Code. In the File Explorer Pane, right click on docker-compose-basic.yaml and choose “Compose Up”. Or alternatively, using a terminal, navigate to the example folder and run docker compose -f docker-compose-basic.yml up -d
Wait for the container to build. It goes pretty quickly, depending on your internet speed.
Go to http://localhost:8787 in your browser and you’ll have RStudio in your browser!
Open project/ and click on example-project.Rproj to open an RStudio Project for the example R code.
Render analysis.qmd to run the analysis.

The project/ folder in the example will automatically be mounted into the Docker container, so any edits you make while in the container will be saved on your local computer. If you render the analysis file, it’ll create an HTML file or Typst-based PDF, and those will appear on your local computer too because of the mounted folder

The same analysis output in the container and locally

This process is super neat. You can access the container-based R from your browser and run and edit everything interactively there, and all changes and output will automatically show up on your local computer.

The main downside to this, though, is that you have to use RStudio in the browser. It doesn’t know about any custom settings and configuration you have in your regular RStudio installation on your computer. If you like to work in a dark theme with a ligature-based font like Fira Code and you have a bunch of custom keyboard shortcuts set up, or custom panel arrangements, or anything else, none of those changes will appear in the browser version because it’s a completely different version of RStudio. You can change settings for RStudio in the browser version, but unless you change the Docker Compose file to make it so that RStudio’s settings files/folders are mounted to your computer, any changes you make are ephemeral and will disappear the next time you run the Docker image. Also, your browser’s keyboard shortcuts will conflict with RStudio’s keyboard shortcuts. Like in Arc, ⌘S is used to toggle the tab sidebar. If I use it to save a file in RStudio, it’ll complain. And—maybe most annoyingly—if you use the ⌘W keyboard shortcut to close a file tab inside RStudio, it’ll actually close the whole browser tab.

Doing actual development in the container/browser version of RStudio is (1) important, since you get to use the version-locked packages and R installation, but (2) annoying because it feels like using someone else’s computer without all your tweaks and adjustments. Because of this, in practice, for projects where I care about exact replicability with version-locked packages and dependencies, I generally do all my coding and writing and analysis on my local computer and then try running it in the container at the end.

It would be really neat if there was some way to point RStudio or Positron to the container-based version of R and its package library so that I can use my regular local IDE, but backed by the container environment, but there’s no way to do that—both RStudio and Positron can only run local R installations.

How to use Docker-based R from Positron

But with Positron, there’s now (kind of) a way! Positron has native support for remote SSH sessions, where a whole local Positron window actually runs on a remote server. I use this (and Visual Studio Code’s SSH support) all the time for accessing files on my web server and doing anything else SSH-related.

The trick here is to use Positron to SSH into the container rather than into some remote machine.

The Rocker R images don’t have SSH enabled—or even installed—so we have to do a little extra work to do that. The docker-compose.yml and Dockerfile files in the example repository show all the changes that need to happen.

In short, we use a custom Dockerfile to pull rocker/tidyverse at R 4.5.0, and then we install and configure an SSH server on it. We can either set it up to require a password to connect (which is annoying) or use password-less public key authentication (which is not annoying once you have a public key set up).

Before trying this, make sure you have a public SSH key set up and accessible at ~/.ssh on your computer. Here’s a little guide.

Here’s the updated docker-compose.yml and the new Dockerfile:

docker-compose.yml

services:
  tidyverse:
    platform: linux/amd64  # Need to include this when building on M1 macOS

    build:
      # Build the Dockerfile
      context: "."
      dockerfile: "Dockerfile"

    ports:
      - "8787:8787"
      - "2222:22"  # localhost:2222 on the computer maps to localhost:22 in the container

    environment:
      - DISABLE_AUTH=true  # Allow password-less login to RStudio server instance

    volumes:
      - ./project:/home/rstudio/project
      # Mount the local public key into the container to SSH into it passwordlessly
      - ~/.ssh/id_rsa.pub:/home/rstudio/.ssh/authorized_keys:ro

Dockerfile

FROM rocker/tidyverse:4.5.0

# Install SSH
RUN apt-get update && \
    apt-get install -y --no-install-recommends openssh-server sudo && \
    rm -rf /var/lib/apt/lists/*

# Set up SSH
RUN mkdir -p /var/run/sshd && \
    echo "PasswordAuthentication no" >> /etc/ssh/sshd_config && \
    echo "PermitRootLogin no" >> /etc/ssh/sshd_config && \
    echo "PubkeyAuthentication yes" >> /etc/ssh/sshd_config

# Password-less sudo in case things need to be installed from the terminal
# I have no idea what the sudo password is for this Docker container ¯\_(ツ)_/¯
RUN echo "rstudio ALL=(ALL) NOPASSWD:ALL" >> /etc/sudoers

# Make sure there's a .ssh folder with the right permissions. The 
# local id_rsa.pub gets mounted to here through Docker Compose
RUN mkdir -p /home/rstudio/.ssh && \
    chown rstudio:rstudio /home/rstudio/.ssh && \
    chmod 700 /home/rstudio/.ssh

EXPOSE 22 8787

# Start both SSH and RStudio
CMD service ssh start && /init

Right click on docker-compose.yml in the File Explorer Pane in Positron or Visual Studio Code and choose “Compose Up”, or run docker compose -f docker-compose.yml up -d in a terminal.

Once it’s done building, the container will still have an RStudio Server instance available at http://localhost:8787, and it will also have an SSH server exposed at port 2222. You can connect to it from the terminal on your computer with:

ssh rstudio@localhost -p 2222

Or you can set up a saved connection (which Positron will be able to read) by adding this entry to ~/.ssh/config:

~/.ssh/config

Host positron-docker
    HostName localhost
    User rstudio
    Port 2222

With that shortcut, you can connect like this from the terminal:

ssh positron-docker

Even better, once you add the connection to ~/.ssh/config, you can use Positron’s Remote Explorer to connect to the not-actually-remote remote SSH server running inside the container:

Connecting to remote connections with Positron’s Remote Explorer

This gives you a(n almost) fully armed and operational Positron session that uses the version-locked installation of R inside the container. All your settings, preferences, themes, fonts, and keyboard shortcuts work here. It’ll feel (almost) just like a local Positron instance. The only thing that doesn’t transfer over to the containerized version of R (or any remote SSH instance) is extensions—Positron/Visual Studio Code will not automatically install your local extensions into the remote computer.

Positron connected to a not-really-remote remote session inside a Docker container

This makes it way easier to develop directly in the container. You get all your regular settings and you don’t get behavior that conflicts with the browser.

Citation

BibTeX citation:

@online{heiss2025,
  author = {Heiss, Andrew},
  title = {Use {Positron} to Run {R} Inside a {Docker} Image Through
    {SSH}},
  date = {2025-07-05},
  url = {https://www.andrewheiss.com/blog/2025/07/05/positron-ssh-docker/},
  doi = {10.59350/fredm-56671},
  langid = {en}
}

For attribution, please cite this work as:

Heiss, Andrew. 2025. “Use Positron to Run R Inside a Docker Image Through SSH.” July 5, 2025. https://doi.org/10.59350/fredm-56671.

Open files in external programs with Positron or Visual Studio Code

Andrew Heiss — Thu, 03 Jul 2025 04:00:00 GMT

Update!

In March/April 2026, Typora created its own VS Code extension (here’s the Open VSX version) for opening Markdown files in Typora, so you don’t need to set up a custom task for that any more. I now have this set as my keyboard shortcut:

keybindings.json

{
    "key": "ctrl+shift+alt+cmd+t",
    "command": "typora.openInTypora"
}

Tasks are still a useful way to open any other type of file in any other external program though!

Positron’s Quarto support is really quite robust and it’s full of nice little features like previewing a citation using whatever bibliography style your document is using (Chicago, APA, etc.) when you hover over a citation key:

Preview of a citation

…or knowing about the different headings inside your document for tab-completed cross references:

Tab completion list of sections in a document

…or even knowing about those headings in other documents in a project:

Tab completion list of sections in a different document

That’s all neat, but for long-form writing, like blog posts, articles, and book chapters, there are tools that make writing with Markdown a lot nicer, especially on macOS. Within most text editing areas, macOS has really neat writing tools and spellchecking and autocorrect features. For example, if you press ⌘^D while your cursor is in a word, you’ll get a little dictionary+thesaurus popup:

Dictionary and thesaurus popup

You also get macOS’s built-in autocorrect and spellchecking and autocomplete:

You can also access macOS’s Writing Tools, which are like a built-in Grammarly. I actually really don’t like these except for the proofread option, which is like an LLM-backed grammar/style suggestion generator.

However, lots of these tools don’t work in Positron (or Visual Studio Code, or really any Electron app). There’s a longstanding bug in Chromium/Electron-based apps (like Positron, Visual Studio Code, Obsidian, and others) where some of these writing tools don’t work (this issue is my least favorite thing about Obsidian). Some non-native apps support spellchecking, and most let you look up selected words in the dictionary or thesaurus, but they don’t use macOS’s autocorrect, autocomplete, or Writing Tools.

Positron and Visual Studio Code don’t support macOS’s built-in dictionary unless you use the vscode-spellright extension which uses macOS’s native spelling API. However, even then, it doesn’t work like it does in regular apps—misspelled words appear with a red underline and display in the Problems Pane, but you can’t right click on misspelled words to change them automatically.

Fixing spelling problems in Positron with vscode-spellright

There are a ton of Markdown editors for macOS, and I’ve probably used them all since 2008(!). Most of them fully support all of macOS’s writing features, and they’re generally designed to create delightful writing experiences. I’m particularly fond of iA Writer, Ulysses, and Typora.

Typora in particular is nice because it’s the only one I’ve found that’s based on pandoc-flavored Markdown, so it can properly preview things like subscripts (CO~2~ becomes CO₂), superscripts (x^2^ becomes x²), strikethroughs (~~deleted~~ becomes ~~deleted~~), inline footnotes (Blah.^[A note.] makes a footnote), and other pandoc-specific syntax. Typora can also display LaTeX math, and it can open Quarto Markdown files without any problems—it just can’t run any code chunks, but that’s not surprising, since it’s not an IDE.

Quarto and Typora

My normal workflow for writing documents in Quarto Markdown (and R Markdown before that) is to open a .qmd file in multiple apps. I open files in Positron or RStudio when working on code-related stuff, and I open files in Typora when focusing more on writing and editing. (I rarely actually ever put them side-by-side, and don’t typically always have the same document open in each simultaneously—it depends on what I’m doing: mainly writing or mainly coding).

The same document open in both Positron and Typora

I’m happy with this system and I’ve used it for years. If Positron / Visual Studio Code or RStudio had full access to all of the macOS spellchecking and autocorrect and writing features, I could do everything there and skip Typora, but IDEs have a different purpose than writing apps, so I don’t think they’ll ever be as polished and nice as Typora or other Markdown editors.

One point of friction in this workflow is getting a .qmd file in Positron open in Typora. I typically need to right click on the .qmd file in the File Explorer Pane, select “Reveal in Finder” to open it in macOS Finder, then right click on the .qmd file in Finder, select “Open With” and finally choose Typora:

That’s a lot of clicking and it’s often enough of a hurdle where I won’t open the file in Typora, and I’ll end up with lots of unnecessary typos or or [sic] repeated words from doing lots of writing in Positron.

It would be great if Positron / Visual Studio Code had an option to open a file in an external program from the File Explorer Pane, but it doesn’t.

But I recently discovered a workaround!

Running external processes with tasks

Both Visual Studio Code and Positron have support for Tasks, which let you run external tools and programs that you define in a tasks.json file.

You can edit the systemwide tasks.json if you go to Positron > Settings > Tasks (on macOS; on Windows, it’ll be somewhere in the settings menu):

Opening tasks.json

The full documentation for the Tasks syntax is online, but it’s fairly straightforward. We can add a task to open a file with Typora like this:

tasks.json

{
    "version": "2.0.0",
    "tasks": [
        {
            "label": "Open in Typora",
            "type": "shell",
            "command": "open -a Typora '${file}'"
        }
    ]
}

To run the task, open the Command Palette, search for “Tasks: Run Task”, and then choose the “Open in Typora” task.

Magical!

Launch tasks more quickly

One downside to this, though, is that it’s a two-step process with the Command Palette—you need to go into “Tasks: Run Task” first, then find the exact task there. We can speed this up in a couple ways.

First, we can assign a keyboard shortcut to the task (see this for complete documentation). In keybindings.json, add an entry for workbench.action.tasks.runTask with the args argument set to whatever you used as label in the task definition. With this, I can mash ⌘⌥⇧^T and open the current file in Typora:

keybindings.json

{
    "key": "shift+alt+cmd+ctrl+t",
    "command": "workbench.action.tasks.runTask",
    "args": "Open in Typora"
}

Second, there’s a neat tasks extension that lets you add little shortcut buttons to the status bar in Positron / Visual Studio Code. After installing the extension, modify tasks.json to look like this:

tasks.json

{
    "version": "2.0.0",
    "statusbar.default.hide": true,  // Hide all tasks by default
    "tasks": [
        {
            "label": "Open in Typora",
            "type": "shell",
            "command": "open -a Typora '${file}'",
            "options": {
                "statusbar": {
                    "hide": false,  // Show this task
                    "label" : "$(pencil) Open in Typora",
                    "detail" : "Open current Markdown-ish file in Typora",
                    "filePattern" : ".*\\.\\S*md$"  // Show this only for Markdown-flavored files
                }
            }
        }
    ]
}

Limit button to specific file types

The filaPattern entry here lets you control when the status button bar appears. The regular expression incantation ".*\\.\\S*md$" will make it appear on any file with a Markdown-like extension, like .md, .qmd, .Rmd, and so on.

There should be a new button down in the status bar that you can click to run the task and open the file in Typora:

An “Open in Typora” button in Positron’s status bar

Do more with tasks

This is just an example of one simple task that opens a file in a different program. You can modify it to any other program, or make even more complex tasks like testing and building packages, or building and deplying websites. You can even chain tasks together in a sequence. See this or this or this for other examples.

Citation

BibTeX citation:

@online{heiss2025,
  author = {Heiss, Andrew},
  title = {Open Files in External Programs with {Positron} or {Visual}
    {Studio} {Code}},
  date = {2025-07-03},
  url = {https://www.andrewheiss.com/blog/2025/07/03/open-files-external-programs-positron/},
  doi = {10.59350/87hpe-4ah24},
  langid = {en}
}

For attribution, please cite this work as:

Heiss, Andrew. 2025. “Open Files in External Programs with Positron or Visual Studio Code.” July 3, 2025. https://doi.org/10.59350/87hpe-4ah24.

How to use a histogram as a legend in {ggplot2}

Andrew Heiss — Wed, 19 Feb 2025 05:00:00 GMT

On Bluesky the other day, I came across this neat post that suggested using a histogram as a plot legend to provide additional context for the data being shown:

Joey Cherdarchuk’s original post

Here’s a closer comparison of those two maps (click to zoom):

Joey Cherdarchuk’s maps side-by-side

This histogram legend is especially useful for choropleth maps where units like counties are sized differently, which can create an illusion of a different distribution. For instance, in that original post, larger dark blue areas stand out a lot visually—like in Alaska, New Mexico, Arizona, and Central California—and make it seem like unemployment is fairly high.

But looking at the histogram that’s not actually the case. Most counties have an unemployment rate around 3–6%. This illusion is happening because land isn’t unemployed—people are.

I thought this was a cool approach, so I figured I’d try to replicate it with R. In the original post, the map was created with D3, the bar chart legend was created with Excel, and the two were combined with Figma. That process is a little too manual for me, but with the magic of R, {ggplot2}, and {patchwork}, we can create the same map completely programmatically.

Let’s do it!

Clean and join data

First, let’s load some packages and tweak some theme settings:

library(tidyverse)
library(readxl)
library(sf)
library(tigris)
library(patchwork)

# Add some font settings to theme_void()
theme_fancy_map <- function() {
  theme_void(base_family = "IBM Plex Sans") +
    theme(
      plot.title = element_text(face = "bold", hjust = 0.13, size = rel(1.4)),
      plot.subtitle = element_text(hjust = 0.13, size = rel(1.1)),
      plot.caption = element_text(hjust = 0.13, size = rel(0.8), color = "grey50"),
    )
}

BLS unemployment data

Next, we can get 2016 unemployment data from the Bureau of Labor Statistics. BLS offers county-level data on annual average labor force participation here, both as plain text and Excel files. The plain text data is structured a little goofily (it’s not comma-separated; it’s a fixed width format where column headings span multiple lines), but the Excel version is in nice columns and is easier to work with. Though even then, we need to skip the first few rows, and the last few rows, and specify column names ourselves.

Download this first from the BLS:

Labor force data by county, 2016 annual averages (XLS)

For the sake of mapping, we’ll truncate the unemployment rate at 9% and mark any counties with higher than 9% unemployment with 9.1 and modify the legend to show “>9%”:

# Load BLS data and clean it up
bls_2016 <- read_excel(
  "laucnty16.xlsx",
  skip = 5,
  col_names = c(
    "laus_code", "STATEFP", "COUNTYFP", "county_name_state",
    "year", "nothing", "labor_force", "employed", "unemployed", "unemp"
  )
) |> 
  # The last few rows in the Excel file aren't actually data, but extra notes,
  # so drop those rows here since they don't have a state FIPS code
  drop_na(STATEFP) |> 
  mutate(
    # Truncate the unemployment rate at 9
    unemp_truncated = ifelse(unemp > 9, 9.1, unemp),
    # Find difference from Fed target of 4%
    unemp_diff = unemp_truncated - 4
  )

bls_2016
## # A tibble: 3,219 × 12
##    laus_code       STATEFP COUNTYFP county_name_state   year  nothing labor_force employed unemployed unemp unemp_truncated unemp_diff
##                                                                           
##  1 CN0100100000000 01      001      Autauga County, AL  2016  NA            25710    24395       1315   5.1             5.1        1.1
##  2 CN0100300000000 01      003      Baldwin County, AL  2016  NA            89778    84972       4806   5.4             5.4        1.4
##  3 CN0100500000000 01      005      Barbour County, AL  2016  NA             8334     7638        696   8.4             8.4        4.4
##  4 CN0100700000000 01      007      Bibb County, AL     2016  NA             8539     7986        553   6.5             6.5        2.5
##  5 CN0100900000000 01      009      Blount County, AL   2016  NA            24380    23061       1319   5.4             5.4        1.4
##  6 CN0101100000000 01      011      Bullock County, AL  2016  NA             4785     4457        328   6.9             6.9        2.9
##  7 CN0101300000000 01      013      Butler County, AL   2016  NA             9116     8484        632   6.9             6.9        2.9
##  8 CN0101500000000 01      015      Calhoun County, AL  2016  NA            45450    42470       2980   6.6             6.6        2.6
##  9 CN0101700000000 01      017      Chambers County, AL 2016  NA            14858    14044        814   5.5             5.5        1.5
## 10 CN0101900000000 01      019      Cherokee County, AL 2016  NA            11241    10671        570   5.1             5.1        1.1
## # ℹ 3,209 more rows

Census geographic data

Next we’ll get geographic data from the US Census with {tigris}

Backup data source

At the time of this writing, {tigris} is working. It wasn’t working a couple weeks ago as the wildly illegal Department of Government Efficiency rampaged through different federal agencies—including the US Census—and shut down the Census’s GIS APIs. But it seems to be working for now?

If it’s not working, IPUMS’s NHGIS project offers the same shapefiles.

The BLS data and the Census data each have columns with state and county FIPS codes which we can use to join the two datasets:

# Get county and state shapefiles from Tigris
us_counties <- counties(year = 2016, cb = TRUE) |> 
  filter(as.numeric(STATEFP) <= 56) |> 
  shift_geometry()  # Move AK and HI

us_states <- states(year = 2016, cb = TRUE) |> 
  filter(as.numeric(STATEFP) <= 56) |> 
  shift_geometry()  # Move AK and HI

# Join BLS data to the map
counties_with_unemp <- us_counties |>
  left_join(bls_2016, by = join_by(STATEFP, COUNTYFP))

# Check out the joined data
counties_with_unemp |> 
  select(STATEFP, COUNTYFP, county_name_state, unemp_truncated, geometry)
## Simple feature collection with 3142 features and 4 fields
## Geometry type: GEOMETRY
## Dimension:     XY
## Bounding box:  xmin: -3112000 ymin: -1698000 xmax: 2258000 ymax: 1566000
## Projected CRS: USA_Contiguous_Albers_Equal_Area_Conic
## First 10 features:
##    STATEFP COUNTYFP    county_name_state unemp_truncated                       geometry
## 1       19      107    Keokuk County, IA             4.3 MULTIPOLYGON (((297173 4548...
## 2       19      189 Winnebago County, IA             3.4 MULTIPOLYGON (((163347 6734...
## 3       20      093    Kearny County, KS             3.1 MULTIPOLYGON (((-482328 605...
## 4       20      123  Mitchell County, KS             3.3 MULTIPOLYGON (((-212918 197...
## 5       20      187   Stanton County, KS             2.8 MULTIPOLYGON (((-528445 214...
## 6       21      005  Anderson County, KY             4.0 MULTIPOLYGON (((940067 1094...
## 7       21      029   Bullitt County, KY             4.1 MULTIPOLYGON (((873753 1022...
## 8       21      049     Clark County, KY             4.7 MULTIPOLYGON (((1012432 106...
## 9       21      059   Daviess County, KY             4.4 MULTIPOLYGON (((749702 5517...
## 10      21      063   Elliott County, KY             9.1 MULTIPOLYGON (((1102886 138...

The map works!

ggplot() +
  geom_sf(data = us_states, fill = "#0074D9", color = "white", linewidth = 0.25) +
  # Albers projection
  coord_sf(crs = st_crs("ESRI:102003"))

Map adjustments

We need to make a couple little adjustments to the map first. In the original image on Bluesky, there’s extra space on the right side of the map to allow for the legend. We can change the plot window by adding 10% of the width of the map to the right.

Technically we don’t have to work with percents here; the data is currently using the Albers projection, which works in meters, so we could add something like 500,000 meters / 500 km to the left. But this is a more general solution and also works if the map data is in decimal degrees instead of meters.

Also, the far western Aleutian islands mess with the visual balance of the map (and they don’t appear because they’re so small), so we’ll also subtract 10% of the map from the left.

# Get x-axis limits of the bounding box for the state data
xlim_current <- st_bbox(us_states)$xlim

# Add 540ish km (or 10% of the US) to the bounds (thus shifting the window over)
xlim_expanded <- c(
  xlim_current[1] + (0.1 * diff(xlim_current)), 
  xlim_current[2] + (0.1 * diff(xlim_current))
)

ggplot() +
  geom_sf(data = us_states, fill = "#0074D9", color = "white", linewidth = 0.25) +
  coord_sf(crs = st_crs("ESRI:102003"), xlim = xlim_expanded)

Extract interior state borders

Because we’re using color = "white", linewidth = 0.25, every state gets a thin white border. This causes some issues though. All the states that share borders actually get a thicker border, since a state’s western border joins up with its neighbor’s eastern border. Also, all the coastlines and islands get borders, which diminishes the landmass—especially on a white background.

Like, look at Alaska’s Aleutian Islands, or Hawai’i’s smaller islands, or Michigan’s Les Cheneaux Islands and Isle Royale, or California’s Channel Islands, or the Florida Keys, or North Carolina’s Outer Banks—they all basically disappear.

To fix this, we can use st_intersection() to identify the intersections of all the state shapes (see this and this for more details)

Now all the islands and coastlines have much better definition and the borders between states are truly sized at 0.25:

interior_state_borders <- st_intersection(us_states) |>
  filter(n.overlaps > 1) |> 
  # Remove weird points that st_intersection() adds
  filter(!(st_geometry_type(geometry) %in% c("POINT", "MULTIPOINT")))

ggplot() +
  geom_sf(data = us_states, fill = "#0074D9", linewidth = 0) +
  geom_sf(data = interior_state_borders, linewidth = 0.25, color = "white") +
  coord_sf(crs = st_crs("ESRI:102003"), xlim = xlim_expanded)

Map with horizontal gradient step legend

Now that we have cleaned and adjusted geographic and unemployment data, we can make a fancy map! Instead of building this sequentially, I’ve included all the code all at once, with lots of comments at each step.

A few things to note:

scale_fill_stepsn() lets you use distinct bins of color instead of a continuous gradient
We position the legend inside the plot with theme(legend.position = "inside", legend.position.inside = c(0.86, 0.32)). Those 0.86, 0.32 coordinates took a lot of tinkering to get! The units for legend.position.inside are based on percentages of the plot, so the legend appears where x is 86% across and 32% up. The position changes every time the plot dimensions change. To make life easier as I played with different values, I used {ggview} to specify and lock in exact dimensions of the plot:
```
library(ggview)

p <- ggplot(...) +
  geom_sf(...)

p + canvas(7, 5)
```
I’m not using ggview::canvas() here in the post because I’m specifying figure dimensions with Quarto chunk options instead (fig-width: 7 and fig-height: 5).

Here’s the map!

ggplot() +
  # Add counties filled with unemployment levels
  geom_sf(
    data = counties_with_unemp, aes(fill = unemp_truncated), linewidth = 0
  ) +
  # Add interior state boundaries
  geom_sf(
    data = interior_state_borders, color = "white", linewidth = 0.25
  ) +
  # Show the unemployment legend as steps instead of a standard gradient
  scale_fill_stepsn(
    colours = scales::brewer_pal(palette = "YlGnBu")(9),
    breaks = 1:10,
    limits = c(1, 10),
    # Change the label for >9%
    labels = case_match(
      1:10,
      1 ~ "1%",
      10 ~ ">9%",
      .default = as.character(1:10)
    )
  ) +
  # Yay labels
  labs(
    title = "US unemployment rates",
    subtitle = "2016 annual averages by county",
    caption = "Source: US Bureau of Labor Statistics",
    fill = "Unemployment rate"
  ) +
  # Use Albers projection and new x-axis limits
  coord_sf(crs = st_crs("ESRI:102003"), xlim = xlim_expanded) +
  # Theme adjustments
  theme_fancy_map() +
  theme(
    legend.position = "inside",
    legend.position.inside = c(0.86, 0.32),
    legend.direction = "horizontal",
    legend.text = element_text(size = rel(0.55)),
    legend.title = element_text(hjust = 0.5, face = "bold", size = rel(0.7), margin = margin(t = 3)),
    legend.title.position = "bottom",
    legend.key.width = unit(1.55, "lines"),
    legend.key.height = unit(0.7, "lines")
  )

Map with histogram legend

We can replace the step gradient legend with a histogram that is filled using the same colors as the step legend.

The easiest method that gives us the most control over the legend histogram is to create a separate plot object for the histogram and place it inside the map with {patchwork}’s inset_element().

Here’s the histogram, again with comments at each step. Only one neat trick to note here:

geom_histogram automatically determines the bin width for the variable assigned to the x aesthetic. In order to fill each bar by bin-specific color, we need to access information about those newly created bins. We can do this with after_stat()—here we fill each bar using the already-calculated x bin categories with fill = after_stat(factor(x))

hist_legend <- ggplot(bls_2016, aes(x = unemp_truncated)) +
  # Fill each histogram bar using the x axis category that ggplot creates
  geom_histogram(
    aes(fill = after_stat(factor(x))), 
    binwidth = 1, boundary = 0, color = "white"
  ) +
  # Fill with the same palette as the map
  scale_fill_brewer(palette = "YlGnBu", guide = "none") +
  # Modify the x-axis labels to use >9%
  scale_x_continuous(
    breaks = 2:10, 
    labels = case_match(
      2:10,
      2 ~ "2%",
      10 ~ ">9%",
      .default = as.character(2:10)
    )
  ) +
  # Just one label to replicate the legend title
  labs(x = "Unemployment rate") +
  # Theme adjustments
  theme_fancy_map() +
  theme(
    axis.text.x = element_text(size = rel(0.55)),
    axis.title.x = element_text(size = rel(0.68), margin = margin(t = 3, b = 3), face = "bold")
  )
hist_legend

Next, we’ll place that hist_legend plot inside a map with inset_element(). Like legend.position.inside = c(0.86, 0.32) in the previous map, the left = 0.75, bottom = 0.26, right = 0.98, top = 0.5 values here are percentages of the plot area and they’re fully dependent on the overall dimensions of the plot. Getting these exact numbers took a lot of manual adjusting, and ggview::canvas() was once again indispensable for keeping the plot dimensions constant.

unemp_map <- ggplot() +
  # Add counties filled with unemployment levels
  geom_sf(
    data = counties_with_unemp, aes(fill = unemp_truncated), color = NA, linewidth = 0
  ) +
  # Add interior state boundaries
  geom_sf(
    data = interior_state_borders, color = "white", linewidth = 0.25, fill = NA
  ) +
  # Show the unemployment legend as steps instead of a standard gradient, but
  # don't actually show the legend
  scale_fill_stepsn(
    colours = scales::brewer_pal(palette = "YlGnBu")(9),
    breaks = 1:10, 
    guide = "none"
  ) +
  # Yay labels
  labs(
    title = "US unemployment rates",
    subtitle = "2016 annual averages by county",
    caption = "Source: US Bureau of Labor Statistics"
  ) +
  # Use Albers projection and new x-axis limits
  coord_sf(crs = st_crs("ESRI:102003"), xlim = xlim_expanded) +
  # Theme stuff
  theme_fancy_map()

# Add the histogram to the map
combined_map_hist <- unemp_map + 
  inset_element(hist_legend, left = 0.75, bottom = 0.26, right = 0.98, top = 0.45)
combined_map_hist

Map with automatic histogram legend with {legendry}

Finally, the new {legendry} package makes it so we can create a custom histogram-based legend without needing to use {patchwork} with a separate histogram plot!

It doesn’t provide as much control over the resulting histogram. The gizmo_histogram() function uses base R’s hist() behind the scenes, so we have to specify bin widths and other settings in hist.arg as base R arguments, like breaks = 10 instead of ggplot’s binwidth = 10.

Not all of hist()’s options seem to work here. For instance, I get a warning if I use border = "white" to add a white border around each bar (argument ‘border’ is not made use of), since that border option is disabled when using base R’s hist() with plot = FALSE:

hist(counties_with_unemp$unemp_truncated, breaks = 10, border = "white", plot = FALSE)
## Warning in hist.default(counties_with_unemp$unemp_truncated, breaks = 10, : argument 'border' is not made use of
## $breaks
##  [1]  1  2  3  4  5  6  7  8  9 10
## 
## $counts
## [1]  12 244 589 821 644 410 205  97 119
## 
## $density
## [1] 0.00382 0.07768 0.18752 0.26138 0.20503 0.13053 0.06527 0.03088 0.03789
## 
## $mids
## [1] 1.5 2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5
## 
## $xname
## [1] "counties_with_unemp$unemp_truncated"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"

Also, it’s currently filling each histogram bar with the full gradient, not the 9 distinct steps, and I can’t figure out how to define custom colors for each bar—and it might not even be possible since color settings aren’t picked up anyway because of {legendry}’s use of plot = FALSE 🤷‍♂️.

But despite these downsides, this automatic histogram legend with {legendry} is really neat!

library(legendry)

# Create a custom histogram guide
histogram_guide <- compose_sandwich(
  middle = gizmo_histogram(just = 0, hist.arg = list(breaks = 10)),
  text = "axis_base"
)

ggplot() +
  # Add counties filled with unemployment levels
  geom_sf(
    data = counties_with_unemp, aes(fill = unemp_truncated), color = NA, linewidth = 0
  ) +
  # Add interior state boundaries
  geom_sf(
    data = interior_state_borders, color = "white", linewidth = 0.25, fill = NA
  ) +
  # Show the unemployment legend with a custom histogram guide
  scale_fill_stepsn(
    colours = scales::brewer_pal(palette = "YlGnBu")(9),
    breaks = 1:10,
    limits = c(1, 10),
    guide = histogram_guide,
    # Change the label for >9%
    labels = case_match(
      1:10,
      1 ~ "1%",
      10 ~ ">9%",
      .default = as.character(1:10)
    )
  ) +
  # Yay labels
  labs(
    title = "US unemployment rates",
    subtitle = "2016 annual averages by county",
    caption = "Source: US Bureau of Labor Statistics",
    fill = "Unemployment rate"
  ) +
  # Use Albers projection and new x-axis limits
  coord_sf(crs = st_crs("ESRI:102003"), xlim = xlim_expanded) +
  # Theme stuff
  theme_fancy_map() +
  theme(
    legend.position = "inside",
    legend.position.inside = c(0.86, 0.32),
    legend.direction = "horizontal",
    legend.text = element_text(size = rel(0.55)),
    legend.title = element_text(hjust = 0.5, face = "bold", size = rel(0.7), margin = margin(t = 3)),
    legend.title.position = "bottom"
  )

Bonus! Use points instead of choropleths

We’re still using choropleth maps here, which still isn’t ideal for showing the idea that “land isn’t unemployed”. One solution is to plot points that are sized by population. This is pretty straightforward with {sf}—we need to convert the county polygons into single points, which we can do with st_point_on_surface(). Then, after a bunch of tinkering with legend options, we’ll have this gorgeous map:

# Convert the county polygons into single points
counties_with_unemp_points <- counties_with_unemp |> 
  st_point_on_surface()

unemp_map_points <- ggplot() +
  # Use a gray background
  geom_sf(data = us_states, fill = "gray90", linewidth = 0) +
  geom_sf(data = interior_state_borders, linewidth = 0.25, color = "white") +
  # Include semi-transparent points with shape 21 (so there's a border)
  geom_sf(
    data = counties_with_unemp_points, 
    aes(size = labor_force, fill = unemp_truncated), 
    pch = 21, color = "white", stroke = 0.25, alpha = 0.8
  ) +
  # Control the size of the points in the legend
  scale_size_continuous(
    range = c(1, 9), labels = scales::label_comma(), 
    breaks = c(10000, 100000, 1000000),
    # Make the points black and not have a border
    guide = guide_legend(override.aes = list(pch = 19, color = "black"))
  ) +
  # Show the unemployment legend as steps instead of a standard gradient, but
  # don't actually show the legend
  scale_fill_stepsn(
    colours = scales::brewer_pal(palette = "YlGnBu")(9),
    breaks = 1:10, 
    guide = "none"
  ) +
  # Labels
  labs(
    title = "US unemployment rates",
    subtitle = "2016 annual averages by county",
    caption = "Source: US Bureau of Labor Statistics",
    fill = "Unemployment rate",
    size = "Labor force"
  ) +
  # Albers
  coord_sf(crs = st_crs("ESRI:102003"), xlim = xlim_expanded) +
  # Theme stuff
  theme_fancy_map() +
  theme(
    legend.position = "inside",
    legend.position.inside = c(0.837, 0.13),
    legend.text = element_text(size = rel(0.55)),
    legend.title = element_text(hjust = 0.5, face = "bold", size = rel(0.7), margin = margin(t = 3)),
    legend.title.position = "bottom"
  )

# Add the histogram to the map
combined_map_hist_points <- unemp_map_points + 
  inset_element(hist_legend, left = 0.75, bottom = 0.26, right = 0.98, top = 0.45)
combined_map_hist_points

Bonus #2! Use a diverging color scheme + nested legend circles

But wait, there’s more! Based on discussions with really smart dataviz people on Bluesky in the wake of me posting about this blog post there, we can make two additional tweaks:

While the different sizes for the points are neat, I’m not a fan of how big the vertical spacing is between the 10,000; 100,000; and 1,000,000. Unfortunately there’s no way to change it. Technically we can use legend.key.spacing.y in theme() to adjust it, but that doesn’t work as expected here because each of those legend entries is sized to match the largest point—i.e., the point for 1,000,000 is the biggest, so the legend entries for all the other values match its height, even if they don’t need all that space.

To fix this, we can use guide_circles() from {legendry} to show the different point sizes as coencentric circles, which is more compact (and just looks neat).
Instead of showing a range of low → high values, we can color these counties based on a meaningful midpoint to help highlight which counties are doing great (low unemployment! good!) and which aren’t (high unemployment! bad!). That might not always necessarily be the best approach—showing the full range of actual values like in the original map is a way of just describing the range and doesn’t inherently imply good or bad. But in other plots where data might be more actionable, divergences from some central value would be much more helpful.

In the United States, the Federal Reserve has a unique dual mandate to use macroeconomic policies to target both inflation and unemployment (most other countries’ central banks only target inflation). The Fed typically aims for an inflation rate of 2% and an unemployment rate of 4ish%. So in this new map, we’ll center each county’s unemployment rate around 4% and show the percentage point deviations from that Fed target. Counties colored in darker red have higher unemployment rates than the target; counties colored in blue have lower rates than the target.

We can then imagine that we’re a policymaker interested in unemployment trends—we can look at the map and quickly identify areas that are doing poorly and doing well.

Up at the beginning of the document where we loaded and cleaned the bls_2016 dataset, I’ve added a new variable that centers the unemployment rate at 4:

mutate(unemp_diff = unemp_truncated - 4)

We can then use this to create a new histogram and new map colored with the “vik” palette from the {scico} package, which has lots of neat diverging palettes. We’ll also create a fancy circle-based legend with {legendry}. Here’s the fully annotated code and final map:

library(ggtext)  # For Markdown-based text in ggplot
library(scico)   # For perceptually uniform colors

# Make new histogram legend
hist_legend_diffs <- ggplot(bls_2016, aes(x = unemp_diff)) +
  # Fill each histogram bar using the x axis category that ggplot creates
  # Use boundary = 0.5 to shift the bin ranges from things like 1-2 to 1.5-2.5
  geom_histogram(
    aes(fill = after_stat((x))), 
    binwidth = 1, boundary = 0.5, color = "white"
  ) +
  # Fill with the same palette as the map
  # scale_fill_brewer(palette = "YlGnBu", guide = "none") +
  scale_fill_scico(palette = "vik", midpoint = 0, guide = "none") +
  # Modify the x-axis labels to show perentage point values and format them with
  # markdown to get original unemployment values on separate lines
  scale_x_continuous(
    breaks = -2:5, 
    labels = case_match(
      -2:5,
      -2 ~ "**−2 pp.**
(2%)",
      0 ~ "**0**
(4% ±
0.5 pp.)",
      5 ~ "**>+4 pp.**
(>9%)",
      .default = glue::glue(
        "**{x}**", 
        x = scales::label_comma(
          style_positive = "plus", style_negative = "minus"
        )(-2:5))
    )
  ) +
  # Just one label to replicate the legend title
  labs(x = "Difference from Fed target") +
  # Theme adjustments
  theme_fancy_map() +
  theme(
    axis.text.x = element_markdown(size = rel(0.5), vjust = 1, lineheight = 1.3),
    axis.title.x = element_text(size = rel(0.68), margin = margin(t = 3, b = 3), face = "bold")
  )

unemp_map_points_diffs <- ggplot() +
  # Use a lighter gray background
  geom_sf(data = us_states, fill = "gray95", linewidth = 0) +
  # Use slightly darker state borders
  geom_sf(data = interior_state_borders, linewidth = 0.25, color = "grey60") +
  # Include semi-transparent points with shape 21 (so there's a border)
  geom_sf(
    data = counties_with_unemp_points, 
    aes(size = labor_force, fill = unemp_diff), 
    shape = 21, color = "white", stroke = 0.25, alpha = 0.8
  ) +
  # Control the size of the points in the legend
  scale_size_continuous(
    range = c(1, 11), labels = scales::label_comma(), 
    breaks = c(100000, 1000000, 5000000),
    # Make the points black and not have a border
    guide = guide_circles(
      text_position = "right",
      override.aes = list(
        fill = "grey30", alpha = 0.8
      )
    )
  ) +
  # This is tricky! We want to use the diverging vik palette but have it 
  # centered at 0. With scale_fill_scico(), there's a midpoint argument, like we 
  # used for the histogram. For generating regular lists of colors with scico(), 
  # though, there's no midpoint argument. Instead, we need to make a few 
  # specific adjustments: 
  #
  # 1. Generate 11 possible colors, since there are 5 colors above the 0 
  #    midpoint in the histogram and we need 5 parallel negative colors below 0 
  #    (even though we're only using 2)
  # 2. Set the limits of the legend to the symmetrical -5 to 5 range so that 
  #    it's centered at 0
  # 3. Set the breaks to go asymmetrically from -2:5. But actually set them 
  #    from -2.5 to 4.5 since that matches the shifted histogram, which uses a 
  #    boundary of 0.5 instead of 0 (so the histogram bins cover ranges like 
  #    0.5 to 1.5 instead of 0 to 1)
  scale_fill_stepsn(
    colours = scico::scico(11, palette = "vik"),
    limits = c(-5, 5),
    breaks = seq(-2.5, 4.5, by = 1),
    guide = "none"
  ) +
  # Labels
  labs(
    title = "US unemployment (2016)",
    subtitle = "Differences from the Federal Reserve's 4% target",
    caption = "Source: US Bureau of Labor Statistics",
    fill = "Unemployment rate",
    size = "County labor force"
  ) +
  # Albers
  coord_sf(crs = st_crs("ESRI:102003"), xlim = xlim_expanded) +
  # Theme adjustments
  theme_fancy_map() +
  theme(
    # {legendry} complains if there's no legend.margin setting; using 
    # theme_void() removes that setting and breaks the plot, so we specify 
    # some 0 values here
    legend.margin = margin(0, 0, 0, 0, "pt"),
    legendry.legend.key.margin = margin(0, 5, 0, 0, "pt"),
    legend.ticks = element_line(colour = "black", linetype = "22"),
    legend.position = "inside",
    legend.position.inside = c(0.87, 0.17),
    legend.text = element_text(size = rel(0.55)),
    legend.title = element_text(hjust = 0.5, face = "bold", size = rel(0.7), margin = margin(t = 3)),
    plot.subtitle = element_text(hjust = 0.18),
    legend.title.position = "bottom"
  )

# Add the histogram to the map
combined_map_hist_points_diffs <- unemp_map_points_diffs + 
  inset_element(hist_legend_diffs, left = 0.75, bottom = 0.26, right = 0.98, top = 0.45)
combined_map_hist_points_diffs

Citation

BibTeX citation:

@online{heiss2025,
  author = {Heiss, Andrew},
  title = {How to Use a Histogram as a Legend in \{Ggplot2\}},
  date = {2025-02-19},
  url = {https://www.andrewheiss.com/blog/2025/02/19/ggplot-histogram-legend/},
  doi = {10.59350/gt0nr-wct91},
  langid = {en}
}

For attribution, please cite this work as:

Heiss, Andrew. 2025. “How to Use a Histogram as a Legend in {Ggplot2}.” February 19, 2025. https://doi.org/10.59350/gt0nr-wct91.

How to move Crimea from Russia to Ukraine in maps with R

Andrew Heiss — Thu, 13 Feb 2025 05:00:00 GMT

The Natural Earth Project

The Natural Earth Project provides high quality public domain geographic data with all sorts of incredible detail, at three resolutions: high (1:10m), medium (1:50m), and low (1:110m). I use their data all the time in my own work and research, and the {rnaturalearth} package makes it really easy to get their data into R for immediate mapping. I mean, look at this!

library(tidyverse)
library(sf)
library(rnaturalearth)

# Set some colors
ukr_blue <- "#0057b7"  # Blue from the Ukrainian flag
ukr_yellow <- "#ffdd00"  # Yellow from the Ukrainian flag
rus_red <- "#d62718"  # Red from the Russian flag

clr_ocean <- "#d9f0ff"
clr_land <- "#facba6"

# CARTOColors Prism (https://carto.com/carto-colors/)
carto_prism = c(
  "#5F4690", "#1D6996", "#38A6A5", "#0F8554", "#73AF48", "#EDAD08", 
  "#E17C05", "#CC503E", "#94346E", "#6F4070", "#994E95", "#666666"
)

ne_countries(scale = 110) |> 
  filter(admin != "Antarctica") |> 
  ggplot() + 
  geom_sf(aes(fill = continent), color = "white", linewidth = 0.1) +
  scale_fill_manual(values = carto_prism, guide = "none") +
  coord_sf(crs = "+proj=robin") +
  theme_void()

Natural Earth’s de facto policy

Maps are intensely political things. There are dozens of disputes over maps and land and territories (e.g., Palestine, Western Sahara, Northern Cyprus, Taiwan, Kashmir, etc.), and many UN member states don’t recognize other UN member states (Israel isn’t recognized by many Arab states; Pakistan doesn’t recognize Armenia).

The Natural Earth Project’s official policy for disputed territories is to reflect on-the-ground de facto control over each piece of land¹:

¹ OpenStreetMap does this too.

Natural Earth Vector draws boundaries of sovereign states according to defacto status. We show who actually controls the situation on the ground. For instance, we show China and Taiwan as two separate states. But we show Palestine as part of Israel.

Though they claim that this de facto policy “is rigorous and self consistent”, it gets them in trouble a lot. For instance, there are nearly two dozen issues on GitHub about Crimea, which is illegally occupied by Russia but de jure part of Ukraine. There are huge debates over the ethics of the de facto policy.

Treating the Natural Earth de facto policy as a de facto policy

I’m not weighing in on that policy here! I don’t super like it—it makes it really hard to map Palestine, for instance—but it is what it is. In this post I’m treating the de facto policy as the de facto situation of the data.

Natural Earth de jure points of view

Natural Earth’s solution for disputed territories is to offer different options to reflect country-specific de jure points of view. They offer pre-built high resolution shapefiles for 31 different points of views, so it’s possible to download data that reflect de jure boundaries for a bunch of different countries. Their other shapefiles all have columns like fclass_us, fclass_ua, and so on for doing… something?… with the point of view. I can’t figure out how these columns work beyond localization stuff (i.e. changing place names based on the point of view). The documentation doesn’t say much about how to actually use these different points of view, and pre-built medium and low resolution maps don’t exist yet.

For example, the US doesn’t de jure-ily recognize the Russian occupation of Crimea, so if we download the pre-built high resolution (10m) version of the world from the US point of view, we can see Crimea as part of Ukraine (we have to download this manually—rnaturalearth::ne_countries() doesn’t support POV files):

world_10_us <- read_sf("ne_10m_admin_0_countries_usa/ne_10m_admin_0_countries_usa.shp")

world_10_us |> 
  filter(ADMIN == "Ukraine") |> 
  ggplot() +
  geom_sf(fill = ukr_blue) +
  theme_void()

Unfortunately, the pre-built point-of-view datasets only exist for the 10m high resolution data. If we want to show medium or low resolution maps, we’re stuck with the de facto version of the map, which means Crimea will be shown as part of Russia. Here’s the low resolution version of Ukraine, with Crimea in Russia:

world <- ne_countries(scale = 110, type = "map_units")

ukraine <- world |> filter(admin == "Ukraine")
russia <- world |> filter(admin == "Russia")

ukraine_bbox <- ukraine |> 
  st_buffer(dist = 100000) |>  # Add 100,000 meter buffer around the country 
  st_bbox()

ggplot() +
  geom_sf(data = world, fill = clr_land) +
  geom_sf(data = russia, fill = rus_red) + 
  geom_sf(data = ukraine, fill = ukr_blue, color = ukr_yellow, linewidth = 2) + 
  coord_sf(
    xlim = c(ukraine_bbox["xmin"], ukraine_bbox["xmax"]), 
    ylim = c(ukraine_bbox["ymin"], ukraine_bbox["ymax"])
  ) +
  theme_void() +
  theme(panel.background = element_rect(fill = clr_ocean))

Relocating Crimea manually with R and {sf}

Natural Earth’s recommendation is to “mashup our countries and disputed areas themes to match their particular political outlook”, so we’ll do that here. Though we won’t use any of the point-of-view themes or features because I have no idea how to get those to work.

Instead we’ll manipulate the geometry data directly and move Crimea from the Russia shape to the Ukraine shape by extracting the Crimea POLYGON from Russia and merging it with Ukraine.

The actual geometric shapes for all the countries in world are MULTIPOLYGONs, or collections of POLYGON geometric objects. For instance, Russia is defined as a single MULTIPOLYGON:

russia |> st_geometry()
## Geometry set for 1 feature 
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -180 ymin: 41.15 xmax: 180 ymax: 81.25
## Geodetic CRS:  WGS 84
## MULTIPOLYGON (((178.7 71.1, 180 71.52, 180 70.8...

We can split MULTIPOLYGONs into their component POLYGONs with st_cast(). Russia consists of 14 different shapes:

russia_polygons <- russia |> 
  st_geometry() |> 
  st_cast("POLYGON")

russia_polygons
## Geometry set for 14 features 
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -180 ymin: 41.15 xmax: 180 ymax: 81.25
## Geodetic CRS:  WGS 84
## First 5 geometries:
## POLYGON ((178.7 71.1, 180 71.52, 180 70.83, 178...
## POLYGON ((49.1 46.4, 48.65 45.81, 47.68 45.64, ...
## POLYGON ((93.78 81.02, 95.94 81.25, 97.88 80.75...
## POLYGON ((102.8 79.28, 105.4 78.71, 105.1 78.31...
## POLYGON ((138.8 76.14, 141.5 76.09, 145.1 75.56...

The second one is the main Russia landmass:

plot(russia_polygons[2])

The last one is the Crimean peninsula:

plot(russia_polygons[14])

Identifying the Crimea POLYGON from a POINT

The only way I figured out what of these POLYGONs were was to plot them individually until I saw a recognizable shape. And if I use a different map (like the 50m or 10m resolution maps), there’s no guarantee that Russia will have 14 POLYGONs or that the 14th one will be Crimea. We need a more reliable way to find the Crimea shape.

One way to do this is to create a POINT object based somewhere in Crimea and do some geometric set math to identify which Russian POLYGON contains it. The point 45°N 34°E happens to be in the middle of Crimea:

crimea_point <- st_sfc(st_point(c(34, 45)), crs = st_crs(world))

ggplot() +
  geom_sf(data = world, fill = clr_land) +
  geom_sf(data = russia, fill = rus_red) + 
  geom_sf(data = ukraine, fill = ukr_blue, color = ukr_yellow, linewidth = 2) + 
  geom_sf(data = crimea_point) +
  coord_sf(
    xlim = c(ukraine_bbox["xmin"], ukraine_bbox["xmax"]), 
    ylim = c(ukraine_bbox["ymin"], ukraine_bbox["ymax"])
  ) +
  theme_void() +
  theme(panel.background = element_rect(fill = clr_ocean))

We can use it with st_intersects() to identify the Russia POLYGON that contains it:

# Extract the Russia MULTIPOLYGON and convert it to polygons
russia_polygons <- world |> 
  filter(admin == "Russia") |> 
  st_geometry() |> 
  st_cast("POLYGON")

# Extract the Russia polygon that has Crimea in it
crimea_polygon <- russia_polygons |>
  keep(\(x) st_intersects(x, crimea_point, sparse = FALSE))

# This is the same as russia_polygons[14]
plot(crimea_polygon)

Extracting the Crimea POLYGON from Russia

We can then remove that polygon from Russia and recombine everything back into a MULTIPOLYGON. It works!

# Remove Crimea from Russia
new_russia <- russia_polygons |>
  discard(\(x) any(st_equals(x, crimea_polygon, sparse = FALSE))) |> 
  st_combine() |> 
  st_cast("MULTIPOLYGON")

ggplot() +
  geom_sf(data = world, fill = clr_land) +
  geom_sf(data = new_russia, fill = rus_red) + 
  geom_sf(data = ukraine, fill = ukr_blue, color = ukr_yellow, linewidth = 2) + 
  coord_sf(
    xlim = c(ukraine_bbox["xmin"], ukraine_bbox["xmax"]), 
    ylim = c(ukraine_bbox["ymin"], ukraine_bbox["ymax"])
  ) +
  theme_void() +
  theme(panel.background = element_rect(fill = clr_ocean))

Adding the Crimea POLYGON to Ukraine

Next we need to merge crimea_polygon with Ukraine. We’ll convert Ukraine to its component POLYGONs, combine those with Crimea, and recombine everything back to a MULTIPOLYGON. It also works!

# Extract the Ukraine MULTIPOLYGON and convert it to polygons
ukraine_polygons <- world |> 
  filter(admin == "Ukraine") |> 
  st_geometry() |> 
  st_cast("POLYGON")

# Add Crimea to Ukraine
new_ukraine <- st_union(c(ukraine_polygons, crimea_polygon)) |>
  st_cast("MULTIPOLYGON")

ggplot() +
  geom_sf(data = world, fill = clr_land) +
  geom_sf(data = new_russia, fill = rus_red) + 
  geom_sf(data = new_ukraine, fill = ukr_blue, color = ukr_yellow, linewidth = 2) + 
  coord_sf(
    xlim = c(ukraine_bbox["xmin"], ukraine_bbox["xmax"]), 
    ylim = c(ukraine_bbox["ymin"], ukraine_bbox["ymax"])
  ) +
  theme_void() +
  theme(panel.background = element_rect(fill = clr_ocean))

Updating Russia and Ukraine in the full data

The last step is to modify the full world dataset and replace the existing geometry values for the two countries with the updated boundaries:

world_un <- world |>
  mutate(geometry = case_when(
    admin == "Ukraine" ~ new_ukraine,
    admin == "Russia" ~ new_russia,
    .default = geometry
  ))

Now that world_un has the corrected boundaries in it, it works like normal. Here’s a map of Eastern Europe, colored by mapcolor9 (a column that comes with Natural Earth data that lets you use 9 distinct colors to fill all countries without having bordering countries share colors). Crimea is in Ukraine now:

eastern_eu_bbox <- ukraine |> 
  st_buffer(dist = 700000) |>  # Add 700,000 meter buffer around the country 
  st_bbox()

ggplot() +
  geom_sf(data = world_un, aes(fill = factor(mapcolor9)), linewidth = 0.25, color = "white") +
  scale_fill_manual(values = carto_prism, guide = "none") +
  coord_sf(
    xlim = c(eastern_eu_bbox["xmin"], eastern_eu_bbox["xmax"]), 
    ylim = c(eastern_eu_bbox["ymin"], eastern_eu_bbox["ymax"])
  ) +
  theme_void() +
  theme(panel.background = element_rect(fill = clr_ocean))

The whole game

Everything above was fairly didactic, with illustrations at each intermediate step. Here’s the whole process all in one place:

world_110 <- ne_countries(scale = 110, type = "map_units")

crimea_point_110 <- st_sfc(st_point(c(34, 45)), crs = st_crs(world_110))

# Extract the Russia MULTIPOLYGON and convert it to polygons
russia_polygons_110 <- world_110 |> 
  filter(admin == "Russia") |> 
  st_geometry() |> 
  st_cast("POLYGON")

# Extract the Russia polygon that has Crimea in it
crimea_polygon_110 <- russia_polygons_110 |>
  keep(\(x) st_intersects(x, crimea_point_110, sparse = FALSE))

# Remove Crimea from Russia
new_russia_110 <- russia_polygons_110 |>
  discard(\(x) any(st_equals(x, crimea_polygon_110, sparse = FALSE))) |> 
  st_combine() |> 
  st_cast("MULTIPOLYGON")

# Extract the Ukraine MULTIPOLYGON and convert it to polygons
ukraine_polygons_110 <- world_110 |> 
  filter(admin == "Ukraine") |> 
  st_geometry() |> 
  st_cast("POLYGON")

# Add Crimea to Ukraine
new_ukraine_110 <- st_union(c(ukraine_polygons_110, crimea_polygon_110)) |>
  st_cast("MULTIPOLYGON")

world_un_110 <- world_110 |>
  mutate(geometry = case_when(
    admin == "Ukraine" ~ new_ukraine_110,
    admin == "Russia" ~ new_russia_110,
    .default = geometry
  ))

Moving Crimea with medium resolution (50m) data

This same approach works for other map resolutions too, like 50m:

world_50 <- ne_countries(scale = 50, type = "map_units")

crimea_point_50 <- st_sfc(st_point(c(34, 45)), crs = st_crs(world_50))

# Extract the Russia MULTIPOLYGON and convert it to polygons
russia_polygons_50 <- world_50 |> 
  filter(admin == "Russia") |> 
  st_geometry() |> 
  st_cast("POLYGON")

# Extract the Russia polygon that has Crimea in it
crimea_polygon_50 <- russia_polygons_50 |>
  keep(\(x) st_intersects(x, crimea_point_50, sparse = FALSE))

# Remove Crimea from Russia
new_russia_50 <- russia_polygons_50 |>
  discard(\(x) any(st_equals(x, crimea_polygon_50, sparse = FALSE))) |> 
  st_combine() |> 
  st_cast("MULTIPOLYGON")

# Extract the Ukraine MULTIPOLYGON and convert it to polygons
ukraine_polygons_50 <- world_50 |> 
  filter(admin == "Ukraine") |> 
  st_geometry() |> 
  st_cast("POLYGON")

# Add Crimea to Ukraine
new_ukraine_50 <- st_union(c(ukraine_polygons_50, crimea_polygon_50)) |>
  st_cast("MULTIPOLYGON")

world_un_50 <- world_50 |>
  mutate(geometry = case_when(
    admin == "Ukraine" ~ new_ukraine_50,
    admin == "Russia" ~ new_russia_50,
    .default = geometry
  ))

Here’s a higher quality map of Eastern Europe with Crimea in Ukraine:

ggplot() +
  geom_sf(data = world_un_50, aes(fill = factor(mapcolor9)), linewidth = 0.25, color = "white") +
  scale_fill_manual(values = carto_prism, guide = "none") +
  coord_sf(
    xlim = c(eastern_eu_bbox["xmin"], eastern_eu_bbox["xmax"]), 
    ylim = c(eastern_eu_bbox["ymin"], eastern_eu_bbox["ymax"])
  ) +
  theme_void() +
  theme(panel.background = element_rect(fill = clr_ocean))

Using the adjusted Natural Earth data as GeoJSON in Observable JS

This updated shapefile works with Observable Plot too (see here for more about how to make nice maps with Observable), but requires one strange tweak because of weird behavior with the GeoJSON file format.

Broken GeoJSON

Let’s export the adjusted geographic data to GeoJSON:

# Save as geojson for Observable Plot
st_write(
  obj = world_un, 
  dsn = "ne_110m_admin_0_countries_un_BROKEN.geojson", 
  driver = "GeoJSON",
  quiet = TRUE,
  delete_dsn = TRUE  # Overwrite the existing .geojson if there is one
)

And then load it with Observable JS:

world_broken = FileAttachment("ne_110m_admin_0_countries_un_BROKEN.geojson").json()

clr_ocean = "#d9f0ff"
clr_land = "#facba6"
ukr_blue = "#0057b7"
ukr_yellow = "#ffdd00"
rus_red = "#d62718"

And then plot it:

Plot.plot({
  projection: "equal-earth",
  marks: [
    Plot.sphere({ fill: clr_ocean }),
    Plot.geo(world_broken, {
      stroke: "black",
      strokeWidth: 0.5,
      fill: clr_land
    }) 
  ]
})

lol what even. The new Ukraine shape seems to have broken boundaries that distort everything else in the map. Weirdly, Ukraine is filled with the ocean color while the rest of the globe—both the ocean and whatever countries didn’t have their borders erased—is the color of land.

Let’s zoom in on just Ukraine:

ukraine = world_broken.features.find(d => d.properties.name === "Ukraine")

Plot.plot({
  projection: { 
    type: "equal-earth", 
    domain: ukraine, 
    inset: 50 
  }, 
  width: 800, 
  marks: [
    Plot.sphere({ fill: clr_ocean }),
    Plot.geo(world_broken, {
      stroke: "black",
      strokeWidth: 0.5,
      fill: clr_land
    }),
    Plot.geo(ukraine, { fill: ukr_blue })
  ]
})

¯\_(ツ)_/¯. Now Ukraine is the color of the ocean and the whole rest of the world is the dark blue of the Ukrainian flag. And it didn’t zoom in at all.

GeoJSON and ↻ ↺ winding order ↻ ↺

This is a symptom of an issue with GeoJSON winding order. GeoJSON cares about the direction that country borders (and all LINESTRING elements) are drawn in. Exterior borders should be drawn counterclockwise; interior borders should be drawn clockwise. If a geographic shape doesn’t follow this winding order, bad things happen. Specifically:

a shape that represents a tiny speck of land becomes inflated to represent the whole globe minus that tiny speck of land, the map fills with a uniform color, the local projection explodes. (via @fil)

That’s exactly what’s happening here. Somehow the winding order is getting reversed when we combine Ukraine with Crimea. {sf} itself doesn’t care about winding order, so everything works fine within R; GeoJSON is picky about winding order, so things break.

Fixing it is tricky though! {sf} uses a bunch of different libraries behind the scenes to do its geographic calculations, including GEOS and S2, and they all have different approaches to polygon creation. Apparently GEOS goes clockwise by default while others go counterclockwise, or something. It should theoretically be possible to fix by adding st_sfc(check_ring_dir = TRUE) after making the new Ukraine shape:

# It would be cool if this worked but it doesn't :(
new_ukraine <- st_union(c(ukraine_polygons, crimea_polygon)) |>
  st_sfc(check_ring_dir = TRUE) |> 
  st_cast("MULTIPOLYGON")

But that doesn’t change anything (nor does it work for this person at GitHub).

Clean GeoJSON with correct winding order

BUT there’s another solution. We can force {sf} to not use the S2 library (which it uses by default, I guess?), since S2 seems go in the wrong direction. If we turn off S2 with sf_use_s2(FALSE), make the new Ukraine shape, and then turn S2 back on with sf_use_s2(TRUE), things work!

# Add Crimea to Ukraine
sf_use_s2(FALSE)
new_ukraine_110 <- st_union(c(ukraine_polygons_110, crimea_polygon_110)) |>
  st_cast("MULTIPOLYGON")
sf_use_s2(TRUE)

Here’s the full process with the 110m map:

world_110 <- ne_countries(scale = 110, type = "map_units")

crimea_point_110 <- st_sfc(st_point(c(34, 45)), crs = st_crs(world_110))

# Extract the Russia MULTIPOLYGON and convert it to polygons
russia_polygons_110 <- world_110 |> 
  filter(admin == "Russia") |> 
  st_geometry() |> 
  st_cast("POLYGON")

# Extract the Russia polygon that has Crimea in it
crimea_polygon_110 <- russia_polygons_110 |>
  keep(\(x) st_intersects(x, crimea_point_110, sparse = FALSE))

# Extract the Ukraine MULTIPOLYGON and convert it to polygons
ukraine_polygons_110 <- world_110 |> 
  filter(admin == "Ukraine") |> 
  st_geometry() |> 
  st_cast("POLYGON")

# Add Crimea to Ukraine
sf_use_s2(FALSE)
## Spherical geometry (s2) switched off
new_ukraine_110 <- st_union(c(ukraine_polygons_110, crimea_polygon_110)) |>
  st_cast("MULTIPOLYGON")
## although coordinates are longitude/latitude, st_union assumes that they are planar
sf_use_s2(TRUE)
## Spherical geometry (s2) switched on

# Remove Crimea from Russia
new_russia_110 <- russia_polygons_110 |>
  discard(\(x) any(st_equals(x, crimea_polygon_110, sparse = FALSE))) |> 
  st_combine() |> 
  st_cast("MULTIPOLYGON")

# Add the modified Russia and Ukraine to the main data
world_un_110_fixed <- world_110 |>
  mutate(geometry = case_when(
    admin == "Ukraine" ~ new_ukraine_110,
    admin == "Russia" ~ new_russia_110,
    .default = geometry
  ))

# Save as GeoJSON
st_write(
  obj = world_un_110_fixed, 
  dsn = "ne_110m_admin_0_countries_un.geojson", 
  driver = "GeoJSON",
  quiet = TRUE,
  delete_dsn = TRUE
)

Here’s the new world map with the the correct Ukraine:

world_fixed = FileAttachment("ne_110m_admin_0_countries_un.geojson").json()

Plot.plot({
  projection: "equal-earth",
  marks: [
    Plot.sphere({ fill: clr_ocean }),
    Plot.geo(world_fixed, {
      stroke: "black",
      strokeWidth: 0.5,
      fill: clr_land
    }) 
  ]
})

We can zoom in on Ukraine too:

ukraine_good = world_fixed.features.find(d => d.properties.name === "Ukraine")
russia = world_fixed.features.find(d => d.properties.name === "Russia")

Plot.plot({
  projection: { 
    type: "equal-earth", 
    domain: ukraine_good, 
    inset: 50 
  }, 
  width: 800, 
  marks: [
    Plot.sphere({ fill: clr_ocean }),
    Plot.geo(world_fixed, {
      stroke: "black",
      strokeWidth: 0.5,
      fill: clr_land
    }),
    Plot.geo(russia, { fill: rus_red }),
    Plot.geo(ukraine_good, { 
      fill: ukr_blue, 
      stroke: ukr_yellow, 
      strokeWidth: 3
    })
  ]
})

Alternative data sources

Natural Earth isn’t the only source of geographic data online, and other sources use de jure borders instead of de facto borders, like these:

GISCO

The European Commission’s Eurostat hosts the Geographic Information System of the Commission (GISCO), which provides GIS data for the EU. They offer global shapefiles that follow EU-based de jure borders. The {giscoR} package provides a nice frontend for getting that data into R at 5 different resolutions (1:60m, 1:20m, 1:10m, 1:3m, and super detailed 1:1m!). It does not come with additional metadata for each country, though (i.e. there are no regional divisions, population values, map colors, names in other languages, and so on), so it requires some extra cleaning work if you want those details. For example, can add region information with {countrycode}:

library(giscoR)
library(countrycode)

world_gisco <- gisco_get_countries(
  year = "2024",
  epsg = "4326",
  resolution = "60"
) |> 
  # Add World Bank regions
  mutate(region = countrycode(ISO3_CODE, origin = "iso3c", destination = "region"))

world_gisco |> 
  filter(NAME_ENGL != "Antarctica") |> 
  ggplot() + 
  geom_sf(aes(fill = region), color = "white", linewidth = 0.1) +
  scale_fill_manual(values = carto_prism, guide = "none") +
  coord_sf(crs = "+proj=robin") +
  theme_void()

Since the EU doesn’t de jure-ily recognize the Russian occupation of Crimea, Crimea is in Ukraine:

ukraine_gisco <- world_gisco |> filter(NAME_ENGL == "Ukraine")
russia_gisco <- world_gisco |> filter(NAME_ENGL == "Russian Federation")

ukraine_gisco_bbox <- ukraine_gisco |> 
  st_buffer(dist = 100000) |>  # Add 100,000 meter buffer around the country 
  st_bbox()

ggplot() +
  geom_sf(data = world_gisco, fill = clr_land) +
  geom_sf(data = russia_gisco, fill = rus_red) + 
  geom_sf(data = ukraine_gisco, fill = ukr_blue, color = ukr_yellow, linewidth = 2) + 
  coord_sf(
    xlim = c(ukraine_gisco_bbox["xmin"], ukraine_gisco_bbox["xmax"]), 
    ylim = c(ukraine_gisco_bbox["ymin"], ukraine_gisco_bbox["ymax"])
  ) +
  theme_void() +
  theme(panel.background = element_rect(fill = clr_ocean))

It’s possible to use GISCO data with Observable too. We could load the data into R and clean it up there (like adding regions and other details) and the save it as GeoJSON, like we did with the Natural Earth data.

Or we can grab the original raw GeoJSON from Eurostat directly. For example, here’s the raw 60M 2024 world map using the WGS84 (4326) projection that we grabbed with gisco_get_countries() earlier.

world_gisco = await FileAttachment("https://gisco-services.ec.europa.eu/distribution/v2/countries/geojson/CNTR_RG_60M_2024_4326.geojson").json()

ukraine_gisco = world_gisco.features.find(d => d.properties.NAME_ENGL === "Ukraine")
russia_gisco = world_gisco.features.find(d => d.properties.NAME_ENGL === "Russian Federation")

Plot.plot({
  projection: { 
    type: "equal-earth", 
    domain: ukraine_gisco, 
    inset: 50 
  }, 
  width: 800, 
  marks: [
    Plot.sphere({ fill: clr_ocean }),
    Plot.geo(world_gisco, {
      stroke: "black",
      strokeWidth: 0.5,
      fill: clr_land
    }),
    Plot.geo(russia_gisco, { fill: rus_red }),
    Plot.geo(ukraine_gisco, { 
      fill: ukr_blue, 
      stroke: ukr_yellow, 
      strokeWidth: 3
    })
  ]
})

Visionscarto

There are JSON-based map files created by @fil at Observable as part of the Visionscarto project.

They’re based on Natural Earth, but with some specific adjustments like adding Crimea to Ukraine, making Gaza a little bit bigger so that it doesn’t get dropped at lower resolutions like 110m, and using UN boundaries for Western Sahara.

Like GISCO, though, these don’t have the additional columns that Natural Earth comes with (country names in a bunch of languages, region and continent designations, map coloring schemes, population and GDP estimates, etc.), and those would need to be added manually in R or Observable or whatever.

import {world110m} from "@visionscarto/geo"

countries110m = topojson.feature(world110m, world110m.objects.countries)
ukraine_visionscarto = countries110m.features.find(d => d.properties.name === "Ukraine")
russia_visionscarto = countries110m.features.find(d => d.properties.name === "Russia")

Plot.plot({
  projection: { 
    type: "equal-earth", 
    domain: ukraine_visionscarto, 
    inset: 50 
  }, 
  width: 800, 
  marks: [
    Plot.sphere({ fill: clr_ocean }),
    Plot.geo(countries110m, {
      stroke: "black",
      strokeWidth: 0.5,
      fill: clr_land
    }),
    Plot.geo(russia_visionscarto, { fill: rus_red }),
    Plot.geo(ukraine_visionscarto, { 
      fill: ukr_blue, 
      stroke: ukr_yellow, 
      strokeWidth: 3
    })
  ]
})

Less automatic sources

There are other sources too, but they require manual downloading:

Citation

BibTeX citation:

@online{heiss2025,
  author = {Heiss, Andrew},
  title = {How to Move {Crimea} from {Russia} to {Ukraine} in Maps with
    {R}},
  date = {2025-02-13},
  url = {https://www.andrewheiss.com/blog/2025/02/13/natural-earth-crimea/},
  doi = {10.59350/28kp0-nbq92},
  langid = {en}
}

For attribution, please cite this work as:

Heiss, Andrew. 2025. “How to Move Crimea from Russia to Ukraine in Maps with R.” February 13, 2025. https://doi.org/10.59350/28kp0-nbq92.

Using USAID data to make fancy world maps with Observable Plot

Andrew Heiss — Mon, 10 Feb 2025 05:00:00 GMT

As part of Elon Musk’s weird Department of Government Efficiency’s unconstitutional rampage through the federal government, USAID’s ForeignAssistance.gov was taken offline on January 31, 2025. It reappeared on February 3, but it’s not clear how long it will be available, especially as USAID is gutted (despite court orders and injunctions to stop).

I study civil society, human rights, and foreign aid and rely on USAID aid data for several of my research projects, so as a backup, I used Datasette to create a mirror website/API of the entire ForeignAssistance.gov dataset at https://foreignassistance-data.andrewheiss.com/. Everything as of December 19, 2024 is available there, both as a queryable SQL database and as downloadable CSV files.

I also made a little frontend website with links to each individual dataset. As I built that website, I decided to try recreating the ForeignAssistance.gov dashboard, which had neat interactive maps and tables.

ForeignAssistance.gov dashboard

Since Quarto has native support for Observable JS for interactive work, and since I’ve meant to really dig into Observable and figure out how to make more interactive graphs, I figured I’d play around with the rescued USAID data.

So in this post, I show what I learned about working with geographic data and making pretty maps with Observable Plot,

Caveat!

I’m really bad at Javascript! The code here is probably wildly inefficient and feels R-flavored.

But it works, and that’s all that matters :)

Working with map data

Get map data

Observable Plot uses the d3-geo module behind the scenes to parse and work with map data, and D3 typically works with data formatted as GeoJSON. There are tons of high quality geographic data sources online, like the ~~US Census~~ (they’ve been removing those in the past few weeks), IPUMS NHGIS, IPUMS IHGIS, and the Natural Earth project, and cities and states typically offer GIS data for public sector-related data. These data sources tend to be stored as shapefiles, which are a fairly complex (but standard) format for geographic data that involve multiple files.

Observable Plot/D3 might be able to work with shapefiles directly, but it’s nowhere in the documentation. They seem to expect GeoJSON instead. We could hunt around online for GeoJSON data, but—even better—we can use the {sf} package in R to convert any shapefile-based data into GeoJSON by setting driver = "GeoJSON" in sf::st_write(). Here we’ll load two datasets from Natural Earth—(1) small scale low resolution 1:110m data for mapping the whole world and (2) medium scale 1:50m data for mapping specific regions and countries—and convert them to GeoJSON files.

library(sf)
library(rnaturalearth)

# Get low resolution Natural Earth data as map units instead of countries because of France
world <- ne_countries(scale = 110, type = "map_units")

# Save as geojson for Observable Plot
st_write(
  obj = world, 
  dsn = "ne_110m_admin_0_countries.geojson", 
  driver = "GeoJSON"
)

# Save medium resolution geojson
st_write(
  obj = ne_countries(scale = 50, type = "countries"), 
  dsn = "ne_50m_admin_0_countries.geojson", 
  driver = "GeoJSON"
)

Maybe skip intermediate saving?

We could probably use Quarto’s special R-to-OJS function ojs_define() and make these R objects directly accessible to OJS without needing to save intermediate files:

ojs_define(world = ne_countries(scale = 110, type = "map_units"))

…but geographic data is complex and I don’t know how things like Observable Plot’s Plot.geo() handle data that’s not read as GeoJSON. So to keep things simple, I ended up just saving these as GeoJSON. 🤷‍♂️

Maps and projections with Observable Plot

We can load these into our document with OJS with FileAttachment():

world = FileAttachment("ne_110m_admin_0_countries.geojson").json()
world_medium = FileAttachment("ne_50m_admin_0_countries.geojson").json()

Check out the structure of world. It’s a FeatureCollection with a slot named crs with the projection information and a slot named features with entries for each country. Each country Feature has a slot named properties with columns like name, iso_a3, formal_en, pop_est, and other details.

world

To plot it, we can use the Geo mark:

Plot.plot({
  marks: [
    Plot.geo(world)
  ]
})

To make things look nicer throughout this post, we’ll define some nicer colors for countries and land and ocean from CARTOColors:

carto_prism = [
  "#5F4690", "#1D6996", "#38A6A5", "#0F8554", "#73AF48", "#EDAD08", 
  "#E17C05", "#CC503E", "#94346E", "#6F4070", "#994E95", "#666666"
]

// From R:
// clr_ocean <- colorspace::lighten("#88CCEE", 0.7)
clr_ocean = "#D9F0FF"

// From CARTOColors Peach 2
clr_land = "#facba6"

We’ll make the land be orange-ish, add some thin black borders around the countries, and include a blue background color with Plot.frame():

Plot.plot({
  marks: [
    Plot.frame({ fill: clr_ocean }),  
    Plot.geo(world, { 
      stroke: "black", 
      strokeWidth: 0.5, 
      fill: clr_land 
    }) 
  ]
})

Built-in projections

Taking a round globe and smashing it on a two-dimensional surface always requires geometric shenanigans to get things flat. We can control how things get flattened by specifying the projection for the map. Here we’ll use the Equal Earth projection (invented in 2018 to show countries and continents at their true relative sizes to each other). Since projections contain relative height and width details, we need to specify a width for the plot now. I arbitrarily chose 1000 pixels here, which is the maximum width—it should autoshrink in smaller browser windows, and the height should be calculated automatically. Finally, instead of adding the background color with Plot.frame(), we can use Plot.sphere() to get a nicer background that uses the specified projection:

Plot.plot({
  projection: "equal-earth", 
  width: 1000, 
  marks: [
    Plot.sphere({ fill: clr_ocean }), 
    Plot.geo(world, {
      stroke: "black",
      strokeWidth: 0.5,
      fill: clr_land
    })
  ]
})

The Observable Plot library includes a bunch of common built-in projections:

viewof projection = Inputs.select(
  ["equirectangular", "equal-earth", "mercator", "transverse-mercator", "azimuthal-equal-area", "gnomonic"],
  {value: "azimuthal-equal-area", label: "Projection"}
)

Plot.plot({
  projection: projection,
  // width,
  marks: [
    Plot.sphere({ fill: clr_ocean }),
    Plot.geo(world, {
      stroke: "black",
      strokeWidth: 0.5,
      fill: clr_land
    })
  ]
})

Other projections

Observable Plot can support any other D3 projection too. There are a whole bunch of projections in the main d3-geo module, and there’s a separate d3-geo-projection module for dozens of others. My favorite global projection is Robinson (the foundation for Equal Earth), which lives in d3-geo-projection. To use it, we can import the module with require() and then access it with d3_geo_projection.geoRobinson():

d3_geo_projection = require("d3-geo-projection") 

Plot.plot({
  projection: d3_geo_projection.geoRobinson(), 
  width: 1000,
  height: 500, 
  marks: [
    Plot.sphere({ fill: clr_ocean }),
    Plot.geo(world, {
      stroke: "black",
      strokeWidth: 0.5,
      fill: clr_land
    })
  ]
})

Filtering map data and adjusting projections

Removing elements

Now that we have a nice projection, we can tweak the map a little. Antarctica is taking up a big proportion of the southern hemisphere, so we’ll filter it out. The world object that has all the map data keeps each country object inside a features slot:

world.features

We can filter it using Javascript’s .filter() function. To make sure that the resulting array keeps the geographic-ness of the data and is a FeatureCollection, we need to create a similarly structured object, with type and features slots:

world_sans_penguins = ({ 
  type: "FeatureCollection", 
  features: world.features.filter(d => d.properties.iso_a3 !== "ATA") 
}) 

Plot.plot({
  projection: d3_geo_projection.geoRobinson(),
  width: 1000,
  height: 500,
  marks: [
    Plot.sphere({ fill: clr_ocean }),
    Plot.geo(world_sans_penguins, { 
      stroke: "black",
      strokeWidth: 0.5,
      fill: clr_land
    })
  ]
})

That works and Antarctica is gone, as expected, but in reality the map didn’t actually change that much. Even if we stop using the sphere background and just fill the plot frame, we can see that the area where Antarctica was is still there, it’s just missing the land itself:

Plot.plot({
  projection: d3_geo_projection.geoRobinson(),
  width: 1000,
  height: 500,
  marks: [
    Plot.frame({ fill: clr_ocean, stroke: "black", strokeWidth: 1 }), 
    Plot.geo(world_sans_penguins, {
      stroke: "black",
      strokeWidth: 0.5,
      fill: clr_land
    })
  ]
})

Quick and dirty cheating method: change the width or height

One quick and dirty solution is to mess with the dimensions and shrink the height. After some trial and error, 430 pixels looks good:

Plot.plot({
  projection: d3_geo_projection.geoRobinson(),
  width: 1000,
  height: 430, 
  marks: [
    Plot.frame({ fill: clr_ocean, stroke: "black", strokeWidth: 1 }),
    Plot.geo(world_sans_penguins, {
      stroke: "black",
      strokeWidth: 0.5,
      fill: clr_land
    })
  ]
})

While this works in this case, it’s not a universal solution. The only reason this works is because Antarctica happens to be at the bottom of the map. When you adjust the height of the plot area, the map itself is anchored to the top. Like, if we set the height to 215, we’ll get just the northern hemisphere:

Plot.plot({
  projection: d3_geo_projection.geoRobinson(),
  width: 1000,
  height: 215, 
  marks: [
    Plot.frame({ fill: clr_ocean, stroke: "black", strokeWidth: 1 }),
    Plot.geo(world_sans_penguins, {
      stroke: "black",
      strokeWidth: 0.5,
      fill: clr_land
    })
  ]
})

As far as I can tell, there’s no way to anchor the map in any other position. If we filter the map data to only look at one continent, there’s no easy way to focus on just that continent by adjusting only the width or height options. Here’s Africa all by itself in a big empty plot area:

just_africa = ({ 
    type: "FeatureCollection", 
    features: world.features.filter(d => d.properties.continent == "Africa") 
}) 

Plot.plot({
  projection: d3_geo_projection.geoRobinson(),
  width: 1000,
  height: 430, 
  marks: [
    Plot.frame({ fill: clr_ocean, stroke: "black", strokeWidth: 1 }),
    Plot.geo(just_africa, {
      stroke: "black",
      strokeWidth: 0.5,
      fill: clr_land
    })
  ]
})

If we adjust the width or the height, the plot area will be resized with the map anchored in the top left corner so we’re left with just the northwestern part of Africa (and big empty areas where North America, South America, and Europe would be):

Plot.plot({
  projection: d3_geo_projection.geoRobinson(),
  width: 550, 
  height: 215, 
  marks: [
    Plot.frame({ fill: clr_ocean, stroke: "black", strokeWidth: 1 }),
    Plot.geo(just_africa, {
      stroke: "black",
      strokeWidth: 0.5,
      fill: clr_land
    })
  ]
})

That’s not great, but there are better ways!

Built-in projections and domain settings

The official Observable Plot method for fitting the plot window to a specific area of the map is to define a “domain” for one of the built-in projections to zoom in on specific areas. The documentation shows how to use special functions in d3-geo to create a circle around a point, but you can also pass a GeoJSON object and Plot will use its boundaries for the domain. The built-in projection options also let us control the outside margin of the domain with inset.

Here’s the world map without Antarctica with the Equal Earth projection, with the projection resized to fit within the bounds of world_sans_penguins, with 10 pixels of padding around the landmass. Antarctica is gone now and the rest of the map is vertically centered within the plot area:

Plot.plot({
  projection: { 
    type: "equal-earth", 
    domain: world_sans_penguins, 
    inset: 10 
  }, 
  width: 1000,
  marks: [
    Plot.frame({ fill: clr_ocean, stroke: "black", strokeWidth: 1 }),
    Plot.geo(world_sans_penguins, {
      stroke: "black",
      strokeWidth: 0.5,
      fill: clr_land
    })
  ]
})

We can see what’s happening behind the scenes if we add Plot.sphere() back in. The rounded globe area is still there, but it’s shifted down and out of the frame. We’re essentially panning around and zooming in on the Equal Earth projection:

Plot.plot({
  projection: {
    type: "equal-earth",
    domain: world_sans_penguins, 
    inset: 10
  },
  width: 1000,
  marks: [
    Plot.frame({ stroke: "black", strokeWidth: 1 }), 
    Plot.sphere({ fill: clr_ocean }), 
    Plot.geo(world_sans_penguins, {
      stroke: "black",
      strokeWidth: 0.5,
      fill: clr_land
    })
  ]
})

Passing a GeoJSON object as the domain is really neat because it makes it straightforward to zoom in on specific areas. For instance, here’s the complete medium resolution world map zoomed in around the just_africa object, which keeps non-African countries in the Middle East and southern Europe:

Plot.plot({
  projection: {
    type: "equal-earth",
    domain: just_africa, 
    inset: 10
  },
  width: 600, 
  height: 600, 
  marks: [
    Plot.frame({ fill: clr_ocean, stroke: "black", strokeWidth: 1 }),
    Plot.geo(world_medium, { 
      stroke: "black",
      strokeWidth: 0.5,
      fill: clr_land
    })
  ]
})

We could also extract Africa from the medium resolution world map and plot only that continent, omitting the Middle East and Europe:

just_africa_medium = ({ 
    type: "FeatureCollection", 
    features: world_medium.features.filter(d => d.properties.continent == "Africa") 
}) 

Plot.plot({
  projection: {
    type: "equal-earth",
    domain: just_africa,
    inset: 10
  },
  width: 600,
  height: 600,
  marks: [
    // Use a white background since we don't want to make it look like  
    // the Sinai peninsula has a coastline  
    Plot.frame({ fill: "white", stroke: "black", strokeWidth: 1 }), 
    Plot.geo(just_africa_medium, {  
      stroke: "black",
      strokeWidth: 0.5,
      fill: clr_land
    })
  ]
})

Other projections and `.fitExtent()`

Unfortunately, it’s a little bit trickier to set the domain and inset for projections that aren’t built in to Plot. We can’t do this:

Plot.plot({
  projection: {
    type: d3_geo_projection.geoRobinson(),
    domain: world_sans_penguins,
    inset: 10
  },
  ...
})

Instead, we need to adjust the size of the projection window itself and build in the inset with d3-geo’s .fitExtent(). This function takes four arguments in an array like [[x1, y1], [x2, y2]], defining the top left and bottom right corners (in pixels) of a window that is centered in the middle of a given GeoJSON object. Here, for instance, we create a copy of the Robinson projection that has a window around just_africa with a top left corner at (30, 30) and a bottom right corner at (570, 570):

inset_africa = 30
africa_map_width = 600
africa_map_height = 600

africa_robinson = d3_geo_projection.geoRobinson()
  .fitExtent(
    [[inset_africa, inset_africa],  // Top left
     [africa_map_width - inset_africa, africa_map_height - inset_africa]],  // Bottom right 
    just_africa
  )

Plot.plot({
  projection: africa_robinson,
  width: africa_map_width,
  height: africa_map_height,
  marks: [
    Plot.frame({ fill: clr_ocean, stroke: "black", strokeWidth: 1 }),
    Plot.geo(world_sans_penguins, {
      stroke: "black",
      strokeWidth: 0.5,
      fill: clr_land
    })
  ]
})

We can use the same approach with individual countries. For extra fun, we’ll fill these countries with distinct colors using Natural Earth’s mapcolor7 column, which assigns countries one of 7 different colors that don’t border other countries (so neighboring countries will never be the same color). We’ll also add some labels in the middle of each country.

egypt = world_medium.features.find(d => d.properties.name === "Egypt")

inset_egypt = 75
egypt_map_width = 600
egypt_map_height = 600

robinson_egypt = d3_geo_projection.geoRobinson()
  .fitExtent(
    [[inset_egypt, inset_egypt],  // Top left
     [egypt_map_width - inset_egypt, egypt_map_height - inset_egypt]],  // Bottom right
    egypt
  )

Plot.plot({
  projection: robinson_egypt,
  width: egypt_map_width,
  height: egypt_map_height,
  marks: [
    Plot.frame({ fill: clr_ocean, stroke: "black", strokeWidth: 1 }),
    Plot.geo(world_medium, Plot.centroid({
      fill: d => d.properties.mapcolor7,
      stroke: "black", 
      strokeWidth: 0.5
    })),
    Plot.geo(egypt, { stroke: "yellow", strokeWidth: 3 }),
    Plot.tip(world_medium.features, Plot.centroid({
      title: d => d.properties.name, 
      anchor: "top",
      fontSize: 13,
      fontWeight: "bold",
      textPadding: 3
    }))
  ],
  color: {
    range: carto_prism
  }
})

The approach works for the whole world_sans_penguins object as well. This addresses our original problem—here’s a world map with the Robinson projection without Antarctica that fills the plot area correctly:

inset_world = 10
world_map_width = 1000
world_map_height = 450

world_sans_penguins_robinson = d3_geo_projection.geoRobinson()
  .fitExtent(
    [[inset_world, inset_world],
     [world_map_width - inset_world, world_map_height - inset_world]],
    world_sans_penguins
  )

Plot.plot({
  projection: world_sans_penguins_robinson,
  width: world_map_width,
  height: world_map_height,
  marks: [
    Plot.frame({ fill: clr_ocean, stroke: "black", strokeWidth: 1 }),
    Plot.geo(world_sans_penguins, {
      stroke: "black",
      strokeWidth: 0.5,
      fill: clr_land
    })
  ]
})

Arbitrary areas and `.fitExtent()`

For bonus fun, this approach also works for any arbitrary rectangles. For example, we can use OpenStreetMap’s neat Export tool to pick the top, bottom, left, and right edges of a box that focuses on Western Europe.

Rectangle around western Europe with OpenStreetMap’s Export tool

We can then use those coordinates to create a MultiPoint geometric feature/object, which essentially acts like a rectangular fake country/region that can be used as the domain or extent of the map:

inset_europe = 10
europe_map_width = 800
europe_map_height = 800

europe_box = ({
  type: "Feature",
  geometry: {
    type: "MultiPoint",
    coordinates: [
      [-13, 35],  // [left/west, bottom/south] (or bottom left corner)
      [21, 60]    // [right/east, top/north]   (or top right corner)
    ]
  }
})

europe_robinson = d3_geo_projection.geoRobinson()
  .fitExtent(
    [[inset_europe, inset_europe],
     [europe_map_width - inset_europe, europe_map_height - inset_europe]], 
    europe_box
  )

Plot.plot({
  projection: europe_robinson,
  width: europe_map_width,
  height: europe_map_height,
  marks: [
    Plot.frame({ fill: clr_ocean, stroke: "black", strokeWidth: 1 }),
    Plot.geo(world_medium, Plot.centroid({
      fill: d => d.properties.mapcolor9,
      stroke: "white", 
      strokeWidth: 0.25
    })),
    Plot.tip(world.features, Plot.centroid({
      title: d => d.properties.name, 
      anchor: "bottom",
      fontSize: 13,
      fontWeight: "bold",
      textPadding: 3
    }))
  ],
  color: {
    range: carto_prism
  }
})

Working with USAID data

Get USAID data

To make it easier to access and filter and manipulate things, I put the rescued data on a Datasette instance, which is nice front-end for an SQLite database. This makes it possible to run SQL queries directly in the browser and generate custom datasets without needing to load the full massive CSV files into R or Python or Stata or whatever.

For example, one of the rescued USAID datasets is named us_foreign_aid_country and it contains 22,000+ rows, with data on aid obligations, appropriations, and disbursements starting in 1999.

If we want to get a total of all constant USD aid obligations by country in 2023, omitting regional and world totals, we could do something like this with R and {dplyr}:

library(tidyverse)

# Download the raw CSV and put it somewhere
us_foreign_aid_country <- read_csv("us_foreign_aid_country.csv")

us_foreign_aid_country |>
  filter(
    `Fiscal Year` == 2023, 
    `Transaction Type Name` == "Obligations",
    !str_detect(`Country Name`, "Region"),
    `Country Name` != "World"
  ) |>
  group_by(`Country Code`, `Country Name`, `Region Name`) |>
  summarize(total_constant_amount = sum(constant_amount)) |>
  arrange(desc(total_constant_amount))
#>  A tibble: 176 × 4
#>  Groups:   Country Code, Country Name [176]
#>   `Country Code` `Country Name`   `Region Name`                total_constant_amount
#>                                                                 
#> 1 UKR            Ukraine          Europe and Eurasia                     17193710403
#> 2 ISR            Israel           Middle East and North Africa            3302860882
#> 3 JOR            Jordan           Middle East and North Africa            1686862605
#> 4 EGY            Egypt            Middle East and North Africa            1503609426
#> 5 ETH            Ethiopia         Sub-Saharan Africa                      1457374911
#> 6 SOM            Somalia          Sub-Saharan Africa                      1181033990
#> 7 NGA            Nigeria          Sub-Saharan Africa                      1019947490
#> 8 COD            Congo (Kinshasa) Sub-Saharan Africa                       990456757
#> 9 AFG            Afghanistan      South and Central Asia                   886536741
#> 0 KEN            Kenya            Sub-Saharan Africa                       846303488
#>  ℹ 166 more rows
#>  ℹ Use `print(n = ...)` to see more rows

Or we could get that data extract directly from the database without needing to load the huge original CSV file. We can run an SQL query like this at the Datasette website:

SELECT "Country Code", "Country Name", "Region Name", SUM("constant_amount") AS total_constant_amount
  FROM "./us_foreign_aid_country"
  WHERE 
    "Fiscal Year" = '2023' 
    AND "Transaction Type Name" = 'Obligations' 
    AND "Country Name" NOT LIKE '%Region%' 
    AND "Country Name" != "World"
  GROUP BY "Country Code", "Country Name", "Region Name"
  ORDER BY total_constant_amount DESC;

SQL query and results

Since we’re working with interactive Observable Javascript, we can load that data directly into the browser instead of downloading intermediate CSV files. There’s a neat Datasette database client for Observable that lets us run SQL queries (there are lots of other clients too, if you want to connect to things like DuckDB, SQLite, MySQL, Snowflake, and so on).

import { DatasetteClient } from "@ambassadors/datasette-client"

aid_db = new DatasetteClient(
  "https://foreignassistance-data.andrewheiss.com/2025-02-03_foreign-assistance"
)

recipient_countries = await aid_db.sql`
  SELECT "Country Code", "Country Name", "Region Name", SUM("constant_amount") AS total_constant_amount
  FROM "./us_foreign_aid_country"
  WHERE 
    "Fiscal Year" = '2023' 
    AND "Transaction Type Name" = 'Obligations' 
    AND "Country Name" NOT LIKE '%Region%' 
    AND "Country Name" != "World"
  GROUP BY "Country Code", "Country Name", "Region Name"
  ORDER BY total_constant_amount DESC;
`

Through the magic of this Datasette client, we now have a pre-summarized dataset to work with!

// I don't want to keep hitting the Datasette server with requests, so I'm 
// cheating and loading a CSV extract instead. It comes from this query: 
// https://foreignassistance-data.andrewheiss.com/2025-02-03_foreign-assistance.csv?sql=SELECT+%22Country+Code%22%2C+%22Country+Name%22%2C+%22Region+Name%22%2C+SUM%28%22constant_amount%22%29+AS+total_constant_amount%0D%0A++FROM+%22.%2Fus_foreign_aid_country%22%0D%0A++WHERE+%22Fiscal+Year%22+%3D+%272023%27+%0D%0A++++AND+%22Transaction+Type+Name%22+%3D+%27Obligations%27+%0D%0A++++AND+%22Country+Name%22+NOT+LIKE+%27%25Region%25%27+%0D%0A++++AND+%22Country+Name%22+%21%3D+%22World%22%0D%0A++GROUP+BY+%22Country+Code%22%2C+%22Country+Name%22%2C+%22Region+Name%22%0D%0A++ORDER+BY+total_constant_amount+DESC%3B&_size=max
recipient_countries = await FileAttachment("recipient_countries.csv").csv({ typed: true })

recipient_countries

Connect USAID data to the map data

Following Observable Plot’s choropleth tutorial, to show these totals on a map, we need to create a Map object,¹ which is like a Python dictionary or an R data frame with two columns, where we have (1) a name that shares a name with something in the geographic data, like an ISO3 country code, and (2) a value with the thing we want to plot.

¹ This term is admittedly confusing because it has nothing to do with geographic maps and is instead related to functional programming.

country_totals = new Map(recipient_countries.map(d => [d["Country Code"], d.total_constant_amount]))
country_totals

This lets us get specific totals with the .get() method. Here’s Ukraine, for example:

country_totals.get("UKR")

We can feed the ISO3 code of each country-level geographic shape into this country_totals object to extract the total amount of aid for each country. We’ll use the Antarctica-free Robinson projection we made earlier, and we’ll remove the ocean fill since we’ll ultimately make this interactive and hoverable:

Plot.plot({
  projection: world_sans_penguins_robinson,
  width: world_map_width,
  height: world_map_height,
  marks: [
    Plot.frame({ stroke: "black", strokeWidth: 1} ),
    Plot.geo(world_sans_penguins, Plot.centroid({
      fill: d => country_totals.get(d.properties.iso_a3)
    }))
  ]
})

Improving the map

We have a choropleth! But this is hardly publication worthy. We need to fix a bunch of issues with it.

First, countries that don’t receive aid don’t appear in the map. Let’s add borders to all the countries:

Plot.plot({
  projection: world_sans_penguins_robinson,
  width: world_map_width,
  height: world_map_height,
  marks: [
    Plot.frame({ stroke: "black", strokeWidth: 1} ),
    Plot.geo(world_sans_penguins, Plot.centroid({
      fill: d => country_totals.get(d.properties.iso_a3)
    })),
    Plot.geo(world_sans_penguins, {  
      stroke: "black",  
      strokeWidth: 0.5  
    })  
  ]
})

The coloring here is gross because of some huge outliers (Ukraine) that make most countries black/dark blue. There’s also no legend to show what these values are. We can address all of this by adjusting the legend options. We’ll log total aid, include the legend, add a nice label, and use a single-hue coloring scheme with gray for countries without aid:

Plot.plot({
  projection: world_sans_penguins_robinson,
  width: world_map_width,
  height: world_map_height,
  marks: [
    Plot.frame({ stroke: "black", strokeWidth: 1} ),
    Plot.geo(world_sans_penguins, Plot.centroid({
      fill: d => country_totals.get(d.properties.iso_a3)
    })),
    Plot.geo(world_sans_penguins, {
      stroke: "black",
      strokeWidth: 0.5
    })
  ],
  color: {  
    scheme: "blues",  
    unknown: "#f2f2f2",  
    type: "log",   
    legend: true,  
    label: "Total obligations",  
  }  
})

Next, let’s make this interactive by turning on hovering tooltips:

Plot.plot({
  projection: world_sans_penguins_robinson,
  width: world_map_width,
  height: world_map_height,
  marks: [
    Plot.frame({ stroke: "black", strokeWidth: 1} ),
    Plot.geo(world_sans_penguins, Plot.centroid({
      fill: d => country_totals.get(d.properties.iso_a3),
      tip: true 
    })),
    Plot.geo(world_sans_penguins, {
      stroke: "black",
      strokeWidth: 0.5
    })
  ],
  color: {
    scheme: "blues",
    unknown: "#f2f2f2",
    type: "log", 
    legend: true,
    label: "Total obligations",
  }
})

That’s so cool. Hover over Mexico and you’ll see “Total obligations 232,214,023”.

We can make this tooltip more informative by including the country name and formatting the amount to show dollars. Instead of using tip: true, we can add the country name as a channel (Observable Plot’s version of a ggplot aesthetic), and format the tip so that the country name comes first and the total amount is formatted with d3.format():

Plot.plot({
  projection: world_sans_penguins_robinson,
  width: world_map_width,
  height: world_map_height,
  marks: [
    Plot.frame({ stroke: "black", strokeWidth: 1} ),
    Plot.geo(world_sans_penguins, Plot.centroid({
      fill: d => country_totals.get(d.properties.iso_a3),
      channels: { 
        Country: d => d.properties.name, 
      }, 
      tip: { 
        format: { 
          Country: true, 
          fill: d3.format("$,d") 
        } 
      } 
    })),
    Plot.geo(world_sans_penguins, {
      stroke: "black",
      strokeWidth: 0.5
    })
  ],
  color: {
    scheme: "blues",
    unknown: "#f2f2f2",
    type: "log", 
    legend: true,
    label: "Total obligations",
  }
})

Now hover over Mexico and you’ll see the country name and the amount of aid in dollars.

Fixing labelling issues

We have two final super minor issues to address.

First hover over a country that didn’t receive aid, like the United States or Australia. The total reported aid displays as “$NaN”. That’s gross. It’d be nicer if it said something else, like “$0” or “No aid” or something more informative.

To fix this, we can make a little function that formats the given value as a dollar amount if it’s an actual value, and formats it as something else if it’s missing or not a number (like log(0)):

function format_aid_total(value) {
  return value ? d3.format("$,d")(value) : "No aid";
}

That works nicely:

format_aid_total(394023)

format_aid_total(NaN)

The other problem is in the legend, which uses a logarithmic scale and includes breaks for 10k, 1M, 100M, and 10G, representing $10,000, $1 million, $100 million, and $10 billion in aid.

The issue is the $10 billion, which is abbreviated with “G”.

This is happening because d3.format() uses SI (Système international d’unités, or International System of Units) values for its numeric formats, which means that it uses SI metric prefixes. Those legend breaks, therefore, actually technically mean this:

10k: 10 kilodollars
1M: 1 megadollar
100M: 100 megadollars
10G: 10 gigadollars

lol, I should start talking about big dollar amounts with these values (“the 2022 US federal budget deficit was 1.4 teradollars”)

The first letters of many of these SI prefixes happen to line up with US-style large numbers:

In the US we already commonly use “k” for thousand
The initial “m” in “mega” aligns with “million”
The initial “t” in tera aligns with “trillion”

But “giga” doesn’t align with “billion”, hence the strange “G” here for dollar amounts.

People have requested that d3-format include an option for switching the abbreviation from G to B, but the developers haven’t added it (and probably won’t). Instead, a common recommended fix is to replace all “G”s with “B”s:

number_in_billions = 13840918291  // A big number I randomly typed

// Billions of dollars instead of SI-style gigadollars
d3.format("$.4s")(number_in_billions).replace("G", "B")

We can add format_aid_total() and the .replace("G", "B") tweak and fix the labels in our interactive map:

Plot.plot({
  projection: world_sans_penguins_robinson,
  width: world_map_width,
  height: world_map_height,
  marks: [
    Plot.frame({ stroke: "black", strokeWidth: 1} ),
    Plot.geo(world_sans_penguins, Plot.centroid({
      fill: d => country_totals.get(d.properties.iso_a3),
      channels: {
        Country: d => d.properties.name,
      },
      tip: {
        format: {
          Country: true,
          fill: d => format_aid_total(d) 
        }
      }
    })),
    Plot.geo(world_sans_penguins, {
      stroke: "black",
      strokeWidth: 0.5
    })
  ],
  color: {
    scheme: "blues",
    unknown: "#f2f2f2",
    type: "log", 
    legend: true,
    label: "Total obligations",
    tickFormat: d => d3.format("$0.2s")(d).replace("G", "B") 
  }
})

Some final tweaks

We’re so close! Just a couple final incredibly minor changes:

We’ll boost the font size of the tooltip a little and increase the font size of the legend
We’ll switch from the built-in ColorBrewer blues palette to show how to use custom gradients, like CARTOColors’s PurpOr sequential palette

Plot.plot({
  projection: world_sans_penguins_robinson,
  width: world_map_width,
  height: world_map_height,
  marks: [
    Plot.frame({ stroke: "black", strokeWidth: 1} ),
    Plot.geo(world_sans_penguins, Plot.centroid({
      fill: d => country_totals.get(d.properties.iso_a3),
      channels: {
        Country: d => d.properties.name,
      },
      tip: {
        fontSize: 12, 
        format: {
          Country: true,
          fill: d => format_aid_total(d)
        }
      }
    })),
    Plot.geo(world_sans_penguins, {
      stroke: "black",
      strokeWidth: 0.15
    })
  ],
  color: {
    // scheme: "blues", 
    range: ["#f9ddda", "#f2b9c4", "#e597b9", "#ce78b3", "#ad5fad", "#834ba0", "#573b88"], 
    unknown: "#f2f2f2",
    type: "log", 
    legend: true,
    label: "Total obligations",
    tickFormat: d => d3.format("$0.2s")(d).replace("G", "B"),
    style: { 
      "font-size": "14px" 
    } 
  }
})

The full game: Complete final code

That final interactive map looks great! We could be even fancier with it by adding dropdowns for dynamically grabbing data for different years or different types of amounts (appropriations, allocations, etc.), or even filter by specific regions or countries. But we won’t.

The different colors and data sources we’ve used are scattered throughout this post. To simplify things, here’s the complete code all in one location. (This chunk doesn’t actually run, since Observable gets mad if you create a new variable with the same name as one that already exists.)

d3_geo = require("d3-geo")
d3_geo_projection = require("d3-geo-projection")

// ----------------------------------------------------------------------
// Map stuff
// ----------------------------------------------------------------------
world = FileAttachment("ne_110m_admin_0_countries.geojson").json()

// Antarctica's ISO3 code is ATA
world_sans_penguins = ({
  type: "FeatureCollection",
  features: world.features.filter(d => d.properties.iso_a3 !== "ATA")
})

inset_world = 10
world_map_width = 1000
world_map_height = 450

world_sans_penguins_robinson = d3_geo_projection.geoRobinson()
  .fitExtent(
    [[inset_world, inset_world],
     [world_map_width - inset_world, world_map_height - inset_world]],
    world_sans_penguins
  )

// ----------------------------------------------------------------------
// Data stuff
// ----------------------------------------------------------------------
import { DatasetteClient } from "@ambassadors/datasette-client"

aid_db = new DatasetteClient(
  "https://foreignassistance-data.andrewheiss.com/2025-02-03_foreign-assistance"
)

recipient_countries = await aid_db.sql`
  SELECT "Country Code", "Country Name", "Region Name", SUM("constant_amount") AS total_constant_amount
  FROM "./us_foreign_aid_country"
  WHERE 
    "Fiscal Year" = '2023' 
    AND "Transaction Type Name" = 'Obligations' 
    AND "Country Name" NOT LIKE '%Region%' 
    AND "Country Name" != "World"
  GROUP BY "Country Code", "Country Name", "Region Name"
  ORDER BY total_constant_amount DESC;
`

country_totals = new Map(recipient_countries.map(d => [d["Country Code"], d.total_constant_amount]))

function format_aid_total(value) {
  return value ? d3.format("$,d")(value) : "No aid";
}

// ----------------------------------------------------------------------
// Plot stuff
// ----------------------------------------------------------------------
Plot.plot({
  projection: world_sans_penguins_robinson,
  width: world_map_width,
  height: world_map_height,
  marks: [
    Plot.frame({ stroke: "black", strokeWidth: 1} ),
    Plot.geo(world_sans_penguins, Plot.centroid({
      fill: d => country_totals.get(d.properties.iso_a3),
      channels: {
        Country: d => d.properties.name,
      },
      tip: {
        fontSize: 12,
        format: {
          Country: true,
          fill: d => format_aid_total(d)
        }
      }
    })),
    Plot.geo(world_sans_penguins, {
      stroke: "black",
      strokeWidth: 0.15
    })
  ],
  color: {
    range: ["#f9ddda", "#f2b9c4", "#e597b9", "#ce78b3", "#ad5fad", "#834ba0", "#573b88"],
    unknown: "#f2f2f2",
    type: "log", 
    legend: true,
    label: "Total obligations",
    tickFormat: d => d3.format("$0.2s")(d).replace("G", "B"),
    style: {
      "font-size": "14px" 
    }
  }
})

Citation

BibTeX citation:

@online{heiss2025,
  author = {Heiss, Andrew},
  title = {Using {USAID} Data to Make Fancy World Maps with {Observable}
    {Plot}},
  date = {2025-02-10},
  url = {https://www.andrewheiss.com/blog/2025/02/10/usaid-ojs-maps/},
  doi = {10.59350/c0aep-hp989},
  langid = {en}
}

For attribution, please cite this work as:

Heiss, Andrew. 2025. “Using USAID Data to Make Fancy World Maps with Observable Plot.” February 10, 2025. https://doi.org/10.59350/c0aep-hp989.

Guide to comparing sample and population proportions with CPS data, both classically and Bayesianly

Andrew Heiss — Mon, 27 Jan 2025 05:00:00 GMT

Last week I was making some final revisions to a paper where we used a neat conjoint experiment to test the effect of a bunch of different treatments on nonprofit donor preference.

One of the peer reviewers asked us to compare the characteristics of our experimental sample with the general population so that we could speak a little to the experiment’s generalizability. This is a super common thing to do with survey research, and one of the main reasons survey researchers include demographic questions in their surveys.

Thanks to the wonders of the R community—and thanks to publicly accessible data—I was able to grab nationally representative demographic data, clean it up and summarize it, run some statistical tests, and make a table to meet the reviewer’s request, all in like 45 minutes.

It was a magically quick and easy process, so I figured I’d make a guide about it so that the rest of the world (but mostly future me) can see how to do it.

Nationally representative demographic data

Finding nationally representative demographic data (in the US, at least) is pretty easy, and there are two common sources for it:

The US Census’s American Community Survey (ACS) is a rolling monthly survey of ≈3.5 million (!!!) US households that’s compiled into an annual dataset.
The US Census’s Current Population Survey (CPS) is a monthly survey of ≈100,000 US individuals. A more comprehensive annual version—the Annual Social and Economic Supplement (ASEC)—is published every March.

The two surveys serve different purposes, and the Census has an FAQ fact sheet explaining the difference between the ACS and CPS. Notably, the ACS only surveys households, and it uses a shorter 8-question survey, while the CPS tries to reach the entire civilian noninstitutionalized population and uses a longer, more detailed survey.

Researchers use both surveys—I’ve used both in my own work. According to the Census, due to its detailed questionnaire and staff experience and regular frequency, the CPS ASEC is a “high quality source of information used to produce the official annual estimate of poverty, and estimates of a number of other socioeconomic and demographic characteristics”. The ACS also has demographic details and people use those instead too.

I’m not entirely sure which one is best—to me they’re both great 🤷‍♂️. Smarter people than me know and care about the difference.

Accessing US Census data

Getting data from the Census is a surprisingly complex process! There are websites and R packages that make it easier though.

ACS

For the ACS, the {tidycensus} R package provides an interface to the Census’s API, and its documentation is great and thorough. Working with the results is tricky though, and involves a lot of pivoting and reshaping and combining variables. I have a whole notebook showing how I access the ACS and create a bunch of variables, with little notes reminding myself how I constructed everything:

Explanation of how I calculated the proportion of households in a block group with a high school education

CPS (and others!)

{tidycensus} doesn’t provide Census API access to CPS data. Instead, IPUMS—a project housed at the University of Minnesota and supported by a consortium of other institutions and companies—provides easy access to all sorts of census and survey data, both through its website and through an API. It’s wild how much data they have. In addition to the CPS, they have the ACS, census microtata for 100+ countries, historical GIS shapefiles, time use surveys, and a ton of other things. It’s an incredible project.

Getting started

Who this guide is for

Here’s what I assume you know:

You’re familiar with R and the tidyverse (particularly {dplyr} and {ggplot2}).
You’re familiar with {brms} for running Bayesian regression models and {tidybayes} and {ggdist} for manipulating and plotting posterior draws.

In this guide, we’ll use IPUMS to get CPS data (monthly and ASEC) with R. Because their data explorer website takes a little while to get used to, I’ll show a bunch of step-by-step screenshots of how to navigate it. It’s possible to access the IPUMS API with the {ipumsr} package, and I also show how to do that in this guide. But in order to use the API, you still need to know how to use the website—you have to find variable names and figure out which samples include which variables. So the screenshots below are still important even if you’re using the API.

I’ll then show how to answer the question of whether a survey proportion is equivalent to a population proportion in a couple diffferent ways:

Frequentist/classical proportion tests with null hypothesis significance testing, and
Bayesian proportion tests and inference based on regions of practical equivalence, or ROPEs

Before getting started, let’s load all the packages we need and create some helpful functions and variables.

The experimental survey data here comes from Chaudhry, Dotson, and Heiss (2024). Since we haven’t published it yet (though it’s close—it’s under review post R&R now!), the original data isn’t quite public yet. So I used the {synthpop} R package to create a synthetic version of part of our data that has the same relationships and distributions as the real results, but is all fake. You can see the R code I used for that process here.

If you want to follow along, you can download this synthetic data here:

synthetic_data.rds: RDS version with all the variable attributes (i.e. factor levels and ordering) included
synthetic_data.csv: Plain-text version (but without variable attributes)

To reflect the fact that this is all public, government-created data, I’m using the Public Sans font, an open source font developed as part of the General Services Administration’s USWDS (US Web Design System) for making accessible federal government websites. I’m also using the USWDS’s basic color palette, developed by 18F.

Let’s get started!

library(tidyverse)   # {ggplot2}, {dplyr}, and friends
library(tinytable)   # Nice tables
library(brms)        # Best way to run Stan models
library(tidybayes)   # Manipulate Stan objects and draws
library(broom)       # Convert model objects to data frames
library(glue)        # Easier string construction
library(scales)      # Nicer labels
library(ggdist)      # Plot posterior distributions
library(ggforce)     # Extra ggplot things like facet_col()
library(patchwork)   # Combine ggplot plots

# Load the synthetic survey results
results <- readRDS("synthetic_data.rds")

# Use the cmdstanr backend for brms because it's faster and more modern than the
# default rstan backend. You need to install the cmdstanr package first
# (https://mc-stan.org/cmdstanr/) and then run cmdstanr::install_cmdstan() to
# install cmdstan on your computer.
options(
  mc.cores = 4,
  brms.backend = "cmdstanr"
)

# Set some global Stan options
CHAINS <- 4
ITER <- 2000
WARMUP <- 1000
BAYES_SEED <- 1234

# Nice ggplot theme
theme_public <- function() {
  theme_minimal(base_family = "Public Sans") +
    theme(
      panel.grid.minor = element_blank(),
      plot.title = element_text(family = "Public Sans", face = "bold", size = rel(1.25)),
      plot.subtitle = element_text(family = "Public Sans Light", face = "plain"),
      plot.caption = element_text(family = "Public Sans Light", face = "plain"),
      axis.title = element_text(family = "Public Sans Semibold", size = rel(0.8)),
      axis.title.x = element_text(hjust = 0),
      axis.title.y = element_text(hjust = 1),
      strip.text = element_text(
        family = "Public Sans Semibold", face = "plain",
        size = rel(0.8), hjust = 0
      ),
      strip.background = element_rect(fill = "grey90", color = NA),
      legend.title = element_text(family = "Public Sans Semibold", size = rel(0.8)),
      legend.text = element_text(size = rel(0.8)),
      legend.position = "bottom",
      legend.justification = "left",
      legend.title.position = "top",
      legend.margin = margin(l = 0, t = 0)
    )
}

theme_set(theme_public())
update_geom_defaults("text", list(family = "Public Sans"))
update_geom_defaults("label", list(family = "Public Sans"))

# USWDS basic palette
# https://designsystem.digital.gov/utilities/color/#basic-palette-2
clrs <- c(
  "#e52207", # .bg-red
  "#e66f0e", # .bg-orange
  "#ffbe2e", # .bg-gold
  "#fee685", # .bg-yellow
  "#538200", # .bg-green
  "#04c585", # .bg-mint
  "#009ec1", # .bg-cyan
  "#0076d6", # .bg-blue
  "#676cc8", # .bg-indigo
  "#8168b3", # .bg-violet
  "#d72d79" # .bg-magenta
)

# Some functions for creating percentage point labels
label_pp <- label_number(
  accuracy = 1, scale = 100, suffix = " pp.", style_negative = "minus"
)

label_pp_01 <- label_number(
  accuracy = 0.1, scale = 100, suffix = " pp.", style_negative = "minus"
)

Getting CPS data from the IPUMS website

Go to the IPUMS CPS website and create an account if you don’t already have one.

Once you’re logged in, go the “Select Data” page, where IPUMS tells you to do two things:

Select samples, or specific verisons of different surveys
Select variables, or specific columns in different surveys

Initial data extract page

Visually it looks like you should select samples first, but I actually find it easier to poke around for different variables first, since not all variables are recorded in every sample.

Finding variables

So first let’s look at a few variables to get a feel for the IPUMS data extract website. Click on the little “SEARCH 🔍” button in the “SELECT VARIABLES” section. We could do a bunch of fancy advanced search options, but for now, just search for “age”

Searching for age-related variables

There are 200+(!) age-related variables in the CPS data:

Search results related to age

This big list shows some useful information already. The first variable in the search results is AGE, and it’s the one we care about. It gets recorded in every CPS survey: each monthly one and in the annual ASEC one. Not all the age variables do this—notice that WHYSS1 only appears in the annual ASEC.

If you click on AGE in the “Variable” column, you can see detailed information about the variable, like how it’s coded, a description, and its availability. For example, prior to 1976, age was only available in the annual ASEC; starting in January 1976, it became a monthly thing.

AGE availability

Since we know we want this variable, we can add it to our “Data Cart”. IPUMS ues a shopping metaphor for building a data extract—we can add different variables and samples to a cart and then check it out (for free) once we’ve found everything we’re looking for.

Click on “Add to cart” to add it, then go back to the search page to look for more variables. You can also add variables to the cart without going to the variable details page—there’s a plus sign in the search results page next to each result that will add the variable for you.

Add AGE to cart

Now that we have age, we need to hunt around for other variables we care about, like sex, marital status, voting history, and so on. To speed things up, you can search for their official variable names and add each one to the cart:

Sex = SEX
Marital status = MARST
Education = EDUC
Donating = VLDONATE
Volunteering = VLSTATUS
Volunteering Supplement weight = VLSUPPWT
Voting = VOTED
Voter Supplement weight = VOSUPPWT

Pay attention to the details!

Looking at the details for these variables is helpful since they’re all categorical variables, unlike age. For instance, marital status has 9 different levels:

Nine different marital status (MARST) codes

Also, it’s important to check the variable details to check for availability. While basic demographic variables like age, sex, marital status, etc. are available in both the monthly surveys and in the annual ASEC, more specialized variables are not.

Variables related to philanthropy and volunteering are only available in September (since they’re part of a special CPS Volunteer Supplement), and only in some years:

Volunteer status (VLSTATUS) only available in September, and only in some years

Variables related to voting are only available in November in even-numbered years (since they’re part the CPS’s Voting and Registration Supplement)

Voting status (VOTED) only available in even-numbered years in November

Selecting samples

Great! If you check the cart, you’ll see all the variables we added at the bottom, along with a bunch of other pre-selected columns:

All variables added to the cart, but no samples are selected

We can’t download any data yet, though. We’ve selected the variables—now we need to select the samples. Go back to the “Select data” page and click on “Select samples” (or click on the “Add more samples” button at the top of the data cart page).

By default, IPUMS will have a bunch of different samples pre-checked. In my case, it grabbed all annual ASEC surveys from 2010–2024, and all monthly surveys from 2021–2024.

Pre-selected ASEC samples

Pre-selected monthly samples

Including all these samples would be useful if we were doing some sort of analysis of CPS trends over time, comparing changes in age or education or volunteering or whatever. But that doesn’t matter here—all we want to know is what age (and everything else) looked like at the time the survey was administered. That means we really just need one year.

However, we can’t just choose one sample. Things like demographics are availble in all annual and monthly samples, but volunteering is only available in September in specific years, and voting is only available in November in specific years.

This survey was administered in mid-2019, so we’ll choose samples that are as close to that as possible. Though demographics are available both monthly and annually, I like to use the annual versions because ASEC data is typically used to stand in for annual information—like if you were building a state-year panel dataset, you’d use ASEC data for each year. The ASEC occurs in March and actually overlaps with the monthly March data (IPUMS has a note about that), so the data technically is for March 2019, but whatever. We don’t have to be super precise here.

We’ll use the September 2019 sample for volunteering, even though that’s after the survey was administered. The next earliest volunteering data is the September 2017 sample, which is like 2 years before the survey. Things don’t line up precisely, but again, that’s fine.

Finally, we’ll use the November 2018 sample for voting. That’s before the survey, but it’s the closest we can get—the next alternative is November 2020, which is a year after the survey. Once again, nothing lines up exactly, but it’s fine.

In summary, here are the variables we want and the samples we’ll get them from:

Age (AGE): 2019 ASEC
Sex (SEX): 2019 ASEC
Marital status (MARST): 2019 ASEC
Education (EDUC): 2019 ASEC
Donating (VLDONATE): 2019-09 Monthly
Volunteering (VLSTATUS): 2019-09 Monthly
Volunteering Supplement weight (VLSUPPWT): 2019-09 Monthly
Voting (VOTED): 2018-11 Monthly
Voter Supplement weight (VOSUPPWT): 2018-11 Monthly

Select those three samples (2019 ASEC, September 2019, and November 2018) and click on “Submit sample selections” to add them to the cart.

March 2019 ASEC sample

September 2019 and November 2018 monthly samples

The cart should now have 9 variables and 3 samples. Conveniently, it has a little summary table showing which samples have which variables, where we can confirm that age, sex, marital status, and education are in all three, volunteering and donating are only in September 2019, and voting is only in November 2018.

All variables and samples selected and ready to go

Downloading the data

Now that we have all the variables and samples we care about in the data cart, we can create a data extract and download this stuff.

Click on the “Create data extract” button at the top of the data cart page, which will take you to the official Extract Request page. There are a bunch of extra options here, and you can optionally add a description to the extract, but we’ll ignore all those. Click on the “Submit Abstract” button and wait for the IPUMS server to compile it all.

Extract submission page

Once it’s ready, it’ll appear at your “My Data” page, which will have a list of all your past extracts.

To download the data, we actually need to download two things:

The data itself. Click on the big green “Download .DAT” button.

Button for download a .dat version of the data

.dat vs .dat.gz

Depending on your browser, the downloaded file will either end in .dat or .dat.gz. If it ends in .gz, it’ll be compressed and zipped (the compressed version of this extract is 7.6 MB); if it ends in .dat, it’ll be uncompressed and huge (≈55 MB in this case). Chrome and Firefox will keep the compressed .gz version; Safari will automatically unzip it and throw away the .gz version, which is annoying.

Try to keep the compressed version. You don’t even need to extract/unzip it—the {ipumsr} data loading functions will handle unzipping for you automatically behind the scenes.

The machine-readable XML codebook, or the DDI file. R uses to clean and relabel the raw data when you load it. If you click on the DDI link in the Codebook column, your browser will likely open a plain text XML file, which isn’t really what you want. Instead, right click on the DDI link and choose “Save file as…” or “Download file as…” or whatever your browser calls it. This will let you save the XML file to your computer.

Context menu for the DDI codebook link

There are some other helpful links there too:

If you click on the R link, it’ll give you a barebones R script for loading the data. It’ll look something like this:

# NOTE: To load data, you must download both the extract's data and the DDI
# and also set the working directory to the folder with these files (or change the path below).

if (!require("ipumsr")) stop("Reading IPUMS data into R requires the ipumsr package. It can be installed using the following command: install.packages('ipumsr')")

ddi <- read_ipums_ddi("cps_00001.xml")
data <- read_ipums_micro(ddi)

If you click on the basic codebook, you’ll get a short plain text version of the codebook, which I find really helpful for remembering which variables show up where and how each variable is coded.

Move the newly downloaded .dat and the .ddi files to the same folder somewhere on your computer (preferably in an RStudio Project or a Positron project/folder or wherever your R working directory is). I put mine in a folder named raw_data; you can put it wherever.

We’re finally ready to load this CPS data into R!

More reproducible alternative: using the IPUMS API

Alternatively, it’s possible to use the {ipumsr} package to access the IPUMS API directly and not need to manually download the data extract from the IPUMS website.

The {ipumsr} API functions essentially let you programmatically create and download a data extract cart. Unless you know the IPUMS CPS data really well, you’ll still likely need to hunt around the website for specific variables and their availabilities, so the whole previous section is still relevant.

The {ipumsr} vignette for working with the API is nice and complete—see that for full details.

Here’s an abbreviated example of how to get the same data extract we collected manually from the website:

Go to your IPUMS dashboard and create an API key. This needs to be stored as an environment variable named IPUMS_API_KEY. You can manually add this to your .Renviron file, or you can run this to make {ipumsr} do it for you:
```
ipumsr::set_ipums_api_key("BLAH", save = TRUE)
```

Build a data extract with define_extract_micro(). This is the equivalent of adding stuff to your data cart on the IPUMS website. You’ll need to know two things:

The variable names, which you can find by searching the IPUMS website

The sample IDs, which you can find by running get_sample_info():

library(tidyverse)
library(ipumsr)

all_cps_samples <- get_sample_info(collection = "cps")

# Find the names for 2019 samples
all_cps_samples |> 
  filter(str_detect(description, "ASEC 2019"))
#> # A tibble: 13 × 2
#>    name        description              
#>                               
#>  1 cps2019_01s IPUMS-CPS, January 2019  
#>  2 cps2019_02s IPUMS-CPS, February 2019 
#>  3 cps2019_03b IPUMS-CPS, March 2019    
#>  4 cps2019_04b IPUMS-CPS, April 2019    
#>  5 cps2019_05s IPUMS-CPS, May 2019      
#>  6 cps2019_06s IPUMS-CPS, June 2019     
#>  7 cps2019_03s IPUMS-CPS, ASEC 2019  
#>  ...

The code for creating the extract will look like this:

cps_extract_definition <- define_extract_micro(
  collection = "cps",
  description = "API extract for blog post",
  samples = c(
    "cps2019_03s",  # ASEC, March 2019
    "cps2019_09s",  # CPS, September 2019
    "cps2018_11s"   # CPS, November 2018
  ),
  variables = c(
    "AGE", "SEX", "MARST", "EDUC",
    "VLDONATE", "VLSTATUS", "VLSUPPWT",
    "VOTED", "VOSUPPWT"
  )
)

This is identical to what we had in the cart on the the website: 9 variables and 3 samples:

cps_extract_definition
#> Unsubmitted IPUMS CPS extract 
#> Description: API extract for blog post
#> 
#> Samples: (3 total) cps2019_03s, cps2019_09s, cps2018_11s
#> Variables: (9 total) AGE, SEX, MARST, EDUC, VLDONATE, VLSTATUS, VLSU...

Submit the request to the server to generate the extract. This is equivalent to checking out your data cart on the website.
```
cps_extract <- submit_extract(cps_extract_definition)
#> Successfully submitted IPUMS CPS extract number ZZZZ
```

Download the extract. The extract won’t be downloadable immediately—you need to wait for e-mail confirmation. Once it’s ready, you can download it with download_extract(), which will download both the .dat.gz data and the .xml codebook to your computer:

cps_downloaded <- download_extract(cps_extract, download_dir = "raw_data")
#>  |==================================================| 100%
#>  |==================================================| 100%
#> DDI codebook file saved to ~/blah/raw_data/cps_ZZZZ.xml
#> Data file saved to ~/blah/raw_data/cps_ZZZZ.dat.gz

Once it’s on your computer, you can load it with the standard {ipumsr} process shown in the next section.
```
cps_data <- read_ipums_micro(cps_downloaded)
```

The IPUMS API and literate programming

If you’re using a literate programming document Quarto or R Markdown, don’t include this API extraction process in your document. It will rerun every time you render your document and create a new IPUMS extract each time, which is excessive. It’s best to run this process in a separate R script or function (perhaps orchestrated with something like {targets}), and then load the DDI .xml and .dat data in the document.

Loading CPS data

Getting this data into R is easy thanks to the {ipumsr} package. We feed the XML DDI codebook into read_ipums_ddi() and then feed that into read_ipums_micro()

library(ipumsr)

ddi <- read_ipums_ddi("raw_data/cps_00001.xml")
cps_data <- read_ipums_micro(ddi, data_file = "raw_data/cps_00001.dat.gz", verbose = FALSE)

glimpse(cps_data)
## Rows: 421,402
## Columns: 21
## $ YEAR      2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018,…
## $ SERIAL    1, 1, 3, 4, 4, 4, 4, 5, 5, 6, 6, 7, 7, 7, 8, 8, 9, 9, 10, 10, 11, 12, 13, 13, 14, 14, 14, 15, 15, 16, 17, 19, 20, 20, 21, 21, 22, 23, 23, 23, 23, 24, 25, 26, 26, 26, 26, 26, 28, 28, 28, 28, 29, 30, 30, 30, 32, 36, 36, 37, 37, 37, 37, 39, 39, 39, 39, 39, 40, 40, 41, 41, 41, 41, 42,…
## $ MONTH     11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 1…
## $ HWTFINL   1704, 1704, 1957, 1688, 1688, 1688, 1688, 2090, 2090, 1832, 1832, 1779, 1779, 1779, 1853, 1853, 2077, 2077, 1427, 1427, 1611, 2044, 1738, 1738, 1690, 1690, 1690, 3135, 3135, 2679, 2253, 1639, 1615, 1615, 1515, 1515, 2254, 1459, 1459, 1459, 1459, 1960, 1942, 1701, 1701, 1701, 1701,…
## $ CPSID     2.017e+13, 2.017e+13, 2.018e+13, 2.017e+13, 2.017e+13, 2.017e+13, 2.017e+13, 2.018e+13, 2.018e+13, 2.017e+13, 2.017e+13, 2.018e+13, 2.018e+13, 2.018e+13, 2.017e+13, 2.017e+13, 2.018e+13, 2.018e+13, 2.018e+13, 2.018e+13, 2.017e+13, 2.017e+13, 2.017e+13, 2.017e+13, 2.017e+13, 2.017e…
## $ ASECFLAG  NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ ASECWTH   NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ PERNUM    1, 2, 1, 1, 2, 3, 4, 1, 2, 1, 2, 1, 2, 3, 1, 2, 1, 2, 1, 2, 1, 1, 1, 2, 1, 2, 3, 1, 2, 1, 1, 1, 1, 2, 1, 2, 1, 1, 2, 3, 4, 1, 1, 1, 2, 3, 4, 5, 1, 2, 3, 4, 1, 1, 2, 3, 1, 1, 2, 1, 2, 3, 4, 1, 2, 3, 4, 5, 1, 2, 1, 2, 3, 4, 1, 1, 1, 2, 1, 1, 2, 3, 1, 2, 1, 2, 3, 4, 1, 2, 1, 2, 3, 4,…
## $ WTFINL    1704, 1845, 1957, 1688, 2780, 2780, 2679, 2090, 2090, 1832, 2679, 1754, 1779, 2452, 1853, 1870, 2151, 2077, 2003, 1427, 1611, 2044, 1738, 1738, 1690, 2102, 2537, 3135, 4172, 2679, 2253, 1639, 1615, 2609, 1900, 1515, 2254, 1459, 1546, 1912, 1752, 1960, 1942, 1527, 1701, 2231, 2104,…
## $ CPSIDP    2.017e+13, 2.017e+13, 2.018e+13, 2.017e+13, 2.017e+13, 2.017e+13, 2.017e+13, 2.018e+13, 2.018e+13, 2.017e+13, 2.017e+13, 2.018e+13, 2.018e+13, 2.018e+13, 2.017e+13, 2.017e+13, 2.018e+13, 2.018e+13, 2.018e+13, 2.018e+13, 2.017e+13, 2.017e+13, 2.017e+13, 2.017e+13, 2.017e+13, 2.017e…
## $ CPSIDV    2.017e+14, 2.017e+14, 2.018e+14, 2.017e+14, 2.017e+14, 2.017e+14, 2.017e+14, 2.018e+14, 2.018e+14, 2.017e+14, 2.017e+14, 2.018e+14, 2.018e+14, 2.018e+14, 2.017e+14, 2.017e+14, 2.018e+14, 2.018e+14, 2.018e+14, 2.018e+14, 2.017e+14, 2.017e+14, 2.017e+14, 2.017e+14, 2.017e+14, 2.017e…
## $ ASECWT    NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ AGE       26, 26, 48, 53, 16, 16, 20, 22, 23, 57, 23, 61, 62, 39, 74, 49, 54, 52, 69, 76, 41, 56, 64, 62, 53, 13, 21, 28, 28, 40, 51, 78, 64, 40, 59, 36, 27, 35, 36, 5, 8, 76, 80, 36, 33, 3, 6, 11, 62, 80, 61, 61, 57, 24, 22, 26, 61, 24, 0, 37, 3, 5, 32, 33, 8, 9, 10, 15, 55, 61, 29, 29…
## $ SEX       2, 1, 2, 2, 1, 1, 2, 1, 1, 2, 2, 1, 2, 2, 2, 2, 1, 2, 1, 2, 1, 2, 1, 1, 1, 2, 2, 2, 1, 1, 2, 2, 2, 2, 1, 2, 2, 2, 1, 1, 1, 2, 1, 1, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 2, 1, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 1, 1, 1, 2, 2, 1, 2, 1, 1, 1, 1, 2, 2, 1, 2, 2, 2, 1, 1, 2, 2, 1, 2, 1, 2, 2, 2…
## $ MARST     6, 6, 4, 4, 6, 6, 6, 6, 6, 4, 6, 1, 1, 6, 5, 6, 1, 1, 1, 1, 6, 4, 6, 6, 5, 9, 6, 1, 1, 6, 3, 5, 6, 6, 1, 1, 6, 1, 1, 9, 9, 5, 4, 1, 1, 9, 9, 9, 6, 5, 4, 3, 4, 6, 6, 6, 3, 6, 9, 3, 9, 9, 6, 6, 9, 9, 9, 6, 1, 1, 1, 1, 9, 9, 6, 4, 1, 1, 4, 1, 1, 9, 1, 1, 1, 1, 9, 9, 1, 1, 4, 1, 5…
## $ EDUC      111, 123, 73, 81, 50, 50, 81, 81, 81, 111, 81, 81, 81, 92, 81, 81, 123, 111, 81, 60, 111, 73, 73, 73, 81, 1, 81, 81, 81, 91, 111, 60, 73, 73, 125, 81, 92, 124, 111, 1, 1, 60, 73, 123, 125, 1, 1, 1, 73, 30, 73, 73, 124, 81, 92, 73, 73, 73, 1, 20, 1, 1, 20, 73, 1, 1, 1, 30, 111,…
## $ VLSTATUS  NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ VLDONATE  NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ VLSUPPWT  NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ VOTED     98, 98, 2, 2, 99, 99, 2, 99, 99, 98, 98, 2, 2, 2, 1, 97, 2, 2, 1, 1, 2, 98, 2, 1, 2, 99, 1, 1, 2, 98, 98, 2, 98, 98, 98, 98, 2, 2, 2, 99, 99, 1, 2, 1, 1, 99, 99, 99, 2, 1, 1, 2, 2, 1, 1, 1, 2, 1, 99, 98, 99, 99, 98, 2, 99, 99, 99, 99, 1, 1, 2, 2, 99, 99, 2, 98, 2, 2, 96, 2, 2,…
## $ VOSUPPWT  1704, 1845, 1957, 1688, 2780, 2780, 2679, 2090, 2090, 1832, 2679, 1754, 1779, 2452, 1853, 1870, 2151, 2077, 2003, 1427, 1611, 2044, 1738, 1738, 1690, 2102, 2537, 3135, 4172, 2679, 2253, 1639, 1615, 2609, 1900, 1515, 2254, 1459, 1546, 1912, 1752, 1960, 1942, 1527, 1701, 2231, 2104,…

holy moly we have nearly half a million rows. That’s because we have three samples (2019 ASEC, September 2019, and November 2018) and they’re all stacked on top of each other in this data. We need to filter this huge data to extract the three samples. We’ll also remove rows with missing data.

cps_demographics <- cps_data |>
  # Only look at the 2019 ASEC data
  filter(YEAR == 2019, MONTH == 03, ASECFLAG == 1) |>
  # Remove rows that are missing or are "not in universe"
  mutate(
    SEX = ifelse(SEX == 9, NA, SEX),
    MARST = ifelse(MARST == 9, NA, MARST),
    EDUC = ifelse(EDUC < 1 | EDUC == 999, NA, EDUC)
  )

cps_volunteer <- cps_data |> 
  filter(YEAR == 2019, MONTH == 09) |> 
  # Remove rows that are missing or are "not in universe"
  mutate(across(c(VLSTATUS, VLDONATE), \(x) ifelse(x == 99, NA, x)))

cps_voting <- cps_data |> 
  filter(YEAR == 2018, MONTH == 11) |>
  # Remove rows that are missing or are "not in universe"
  mutate(VOTED = ifelse(VOTED == 99, NA, VOTED))

These counts are more reasonable (but still huge!)

nrow(cps_demographics)
## [1] 180101
nrow(cps_volunteer)
## [1] 118557
nrow(cps_voting)
## [1] 122744

Summarizing CPS data

Ultimately, our goal is to find the population-level average of a bunch of characteristics and see if our sample plausibly matches population averages.

Things get a little tricky and loosey-goosey here. The different levels measured by the CPS don’t always match what’s in the survey. For example, the CPS measures sex and provides only 2 levels (1 = male; 2 = female); the experiment called this construct gender and included male, female,¹ transgender, prefer not to say, and other.

¹ We should have called these “man” and “woman” since gender ≠ sex.

To make the survey question reasonably match what the CPS is capturing, I find that it’s easiest to collapse both the survey data and the CPS data to simpler constructs. Before we collapse things, though, we need to look at one statistical issue: weighting.

Weighting

In an effort to make the CPS nationally representative, every row is weighted—each individual does not represent the same number of persons in the population. The Census oversamples some subpopulations and shifts weights up and down to give individuals more or less statistical influence in the sample so that the survey results better approximate the characteristics of the general population. Any analysis we do with CPS data needs to take those weights into account.

The weights for ASEC variables are included in the ASECWT column; the weights for volunteering and voting variables are in the VLSUPPWT and VOSUPPWT columns

If we’re calculating basic averages, we can use weighted.mean() instead of mean(). Note the difference in average when we don’t weight!

cps_demographics |> 
  summarize(
    avg_age_weighted = weighted.mean(AGE, w = ASECWT),  # BAD
    avg_age_unweighted = mean(AGE)  # GOOD
  )
## # A tibble: 1 × 2
##   avg_age_weighted avg_age_unweighted
##                            
## 1             38.8               37.3

If we’re doing stuff with models, we can use the weights argument:

# BAD: Non-weighted intercept-only model
lm(AGE ~ 1, data = cps_demographics) |> 
  tidy(conf.int = TRUE)
## # A tibble: 1 × 7
##   term        estimate std.error statistic p.value conf.low conf.high
##                                   
## 1 (Intercept)     37.3    0.0540      690.       0     37.2      37.4

# GOOD: Weighted intercept-only model
lm(AGE ~ 1, data = cps_demographics, weights = ASECWT) |> 
  tidy(conf.int = TRUE)
## # A tibble: 1 × 7
##   term        estimate std.error statistic p.value conf.low conf.high
##                                   
## 1 (Intercept)     38.8    0.0542      715.       0     38.7      38.9

Base R only really has weighted.mean(). If we want other things, like a weighted variance, or weighted rank, or weighted table/crosstabs, we can use a bunch of different functions in the {Hmisc} package:

# Some Hmisc::wtd.*() things:
cps_demographics |> 
  summarize(
    avg_age = Hmisc::wtd.mean(AGE, weights = ASECWT),
    var_age = Hmisc::wtd.var(AGE, weights = ASECWT),
    sd_age = sqrt(var_age)
  )
## # A tibble: 1 × 3
##   avg_age var_age sd_age
##          
## 1    38.8    529.   23.0

Calculating population-level proportions

We’ll collapse these population-level CPS values into binary versions of each question so that we can look at things like the proportion of women, the proportion of people who volunteer, and so on. We’ll also collapse age into a binary above/below the median age—this isn’t necessary, and we could totally work with numeric age instead of proportions, but in our anonymized survey data, our age column is an indicator representing being above/below 36 (the median age at the time of the survey).

We’ll do some basic summarizing with weighted.mean() and calculate all these national proportions, along with the weighted standard deviations (which will be important for the Bayesian analysis later in this post).

national_demographics <- cps_demographics |> 
  summarize(
    # AGE is already numeric
    age = weighted.mean(AGE >= 36, ASECWT), 
    age_sd = sqrt(Hmisc::wtd.var(AGE >= 36, weights = ASECWT)),

    # 1 = Female
    female = weighted.mean(SEX == 2, ASECWT),
    female_sd = sqrt(Hmisc::wtd.var(SEX == 2, weights = ASECWT)),

    # 1 = Married, spouse present
    # 2 = Married, spouse absent
    married = weighted.mean(MARST %in% 1:2, na.rm = TRUE),
    married_sd = sqrt(Hmisc::wtd.var(MARST %in% 1:2, weights = ASECWT)),

    # 111 = Bachelor's degree
    college = weighted.mean(EDUC >= 111, ASECWT, na.rm = TRUE),
    college_sd = sqrt(Hmisc::wtd.var(EDUC >= 111, weights = ASECWT))
)

national_volunteer <- cps_volunteer |> 
  summarize(
    # 1 = Volunteer
    volunteering = weighted.mean(VLSTATUS == 1, VLSUPPWT, na.rm = TRUE),
    volunteering_sd = sqrt(Hmisc::wtd.var(VLSTATUS == 1, weights = VLSUPPWT)),

    # 2 = Yes, made a donation to charity in the past 12 months
    donating = weighted.mean(VLDONATE == 2, VLSUPPWT, na.rm = TRUE),
    donating_sd = sqrt(Hmisc::wtd.var(VLDONATE == 2, weights = VLSUPPWT))
  )

national_voting <- cps_voting |>
  summarize(
    # 2 = Voted in the most recent November election
    voting = weighted.mean(VOTED == 2, VOSUPPWT, na.rm = TRUE),
    voting_sd = sqrt(Hmisc::wtd.var(VOTED == 2, weights = VOSUPPWT))
  )

I like to store these in a little one-row data frame so that it’s easy to access invidiual values:

national_values <- bind_cols(
  national_demographics, national_volunteer, national_voting
)
national_values
## # A tibble: 1 × 14
##     age age_sd female female_sd married married_sd college college_sd volunteering volunteering_sd donating donating_sd voting voting_sd
##                                                                   
## 1 0.530  0.499  0.510     0.500   0.410      0.492   0.258      0.437        0.300           0.458    0.474       0.499  0.534     0.499

# Proportion of women
national_values$female
## [1] 0.5097

# Proportion that voted
national_values$voting
## [1] 0.5344

Summarizing sample proportions

We’re almost done! All that’s left is testing whether the demographic characteristics of the survey experiment respondents reasonably match their corresponding population proportions.

First, though, we need to make binary versions of the survey responses. To make life easier, we’ll use the same names as the CPS data:

results_to_test <- results |> 
  mutate(
    age = age == "More than median",
    female = gender == "Female",
    married = marital_status == "Married",
    college = education %in% c(
      "4 year degree", 
      "Graduate or professional degree", 
      "Doctorate"
    ),
    volunteering = volunteer_frequency != "Haven't volunteered in past 12 months",
    donating = donate_frequency == "More than once a month, less than once a year",
    voting = voted == "Yes"
  ) |> 
  select(female, age, married, college, volunteering, donating, voting)

glimpse(results_to_test)
## Rows: 1,300
## Columns: 7
## $ female        FALSE, FALSE, FALSE, FALSE, FALSE, TRUE, FALSE, FALSE, TRUE, FALSE, FALSE, FALSE, TRUE, TRUE, FALSE, FALSE, TRUE, FALSE, FALSE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, TRUE, FALSE, FALSE, FALSE, TRUE, TRUE, FALSE, TRUE, FALSE, FALSE, TRUE, FALSE, TRUE, …
## $ age           TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, TRUE, TRUE, TRUE, FALSE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, TRUE, FALSE, FALSE, FALSE, FALSE, TRUE, FALSE, TRUE, FALSE, TRUE, FALSE, FALSE, FALSE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, TRUE, TRUE, TRUE, TRUE, F…
## $ married       TRUE, FALSE, FALSE, FALSE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, TRUE, FALSE, FALSE, FALSE, TRUE, TRUE, FALSE, TRUE, FALSE, FALSE, TRUE, FALSE, TRUE, FALSE, TRUE, FALSE, FALSE, FALSE, TRUE, FALSE, FALSE, FALSE, TRUE, TRUE, FALSE, FALSE, FALSE, TRUE, TRUE, TRUE, FALSE, TRU…
## $ college       TRUE, FALSE, TRUE, FALSE, TRUE, FALSE, FALSE, FALSE, FALSE, TRUE, TRUE, TRUE, FALSE, FALSE, TRUE, TRUE, FALSE, TRUE, FALSE, FALSE, FALSE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, TRUE, TRUE, FALSE, TRUE, TRUE, FALSE, FALSE, TRUE, FALSE, TRUE, TRUE, …
## $ volunteering  TRUE, FALSE, TRUE, TRUE, FALSE, FALSE, TRUE, FALSE, TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, TRUE, TRUE, FALSE, TRUE, TRUE, FALSE, TRUE, FALSE, TRUE, FALSE, FALSE, TRUE, TRUE, TRUE, FALSE, TRUE, FALSE, TRUE, FALSE, TRUE, FALSE, TRUE, FALSE, TRUE, TRUE, FALSE, TRUE, TRUE, FALSE, TR…
## $ donating      FALSE, FALSE, FALSE, TRUE, FALSE, TRUE, TRUE, FALSE, FALSE, FALSE, TRUE, FALSE, FALSE, FALSE, TRUE, FALSE, TRUE, TRUE, TRUE, FALSE, TRUE, FALSE, FALSE, FALSE, TRUE, FALSE, TRUE, TRUE, TRUE, FALSE, TRUE, FALSE, TRUE, FALSE, TRUE, FALSE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, …
## $ voting        TRUE, TRUE, TRUE, FALSE, TRUE, TRUE, TRUE, FALSE, FALSE, TRUE, TRUE, TRUE, TRUE, FALSE, TRUE, TRUE, TRUE, TRUE, FALSE, TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, TRUE, TRUE, TRUE, TRUE, FALSE, TRUE, TRUE…

Testing sample vs. population proportions frequentist-ly

One-sample proportion test for age

As a quick and easy check, we can run a one-sample proportion test to see if the proportion of a variable is significantly different from a null value. We can do this with prop.test(), which works a bunch of different ways—with matrices, with vectors, and with single values (see this blog post for some other examples of prop.test()).

Let’s look at age first. 50.77% of people in the sample are older than 36; 53% of people in the population are older than 36:

# Proportion of sample older than 36
mean(results_to_test$age)
## [1] 0.5077

# CPS proportion older than 36
national_values$age
## [1] 0.5303

Is that an issue? Is the sample significantly younger than the rest of the country?

For this one-sample test, we need to feed prop.test() three things: (1) the number of “successes”, or the count of rows where the respondent is older than 36 (or where age is TRUE), (2) the number of rows in the sample, and (3) the null value, or the population-level CPS proportion:

prop_test_freq_age <- prop.test(
  x = sum(results_to_test$age),  # Number of "successes" (rows where age == TRUE)
  n = nrow(results_to_test),     # Sample size
  p = national_values$age        # Population-level proportion from the CPS
)

tidy(prop_test_freq_age)
## # A tibble: 1 × 8
##   estimate statistic p.value parameter conf.low conf.high method                                               alternative
##                                                                                   
## 1    0.508      2.59   0.108         1    0.480     0.535 1-sample proportions test with continuity correction two.sided

For fun, we can plot this too:

Code for making this plot ↓

Apple Music Wrapped with R

Andrew Heiss — Wed, 04 Dec 2024 05:00:00 GMT

’Tis the season for Spotify Wrapped stats and I love it, both for seeing what everyone listens to and because it’s such a cool way of presenting data. A few years ago on Twitter, Caitlin Hudon noted that

Spotify Wrapped is a great example of how you can build a fantastic data product without maching learning or AI. (@beeonaposy)

At its core, Spotify Wrapped is really just some grouped and summarized data—a PivotTable with some album cover art slapped on. And it’s fun and neat and everyone loves it!

I’ve always been jealous of everyone’s annual Spotify Wrapped reports, but since I don’t use Spotify, I’ve never gotten to see my own details.

Because I’m an Elder Millennial and started listening to music in the days of Napster, I prefer to control my music files rather than stream it Spotify-style, so I get all my stuff from either the Amazon Music store or Bandcamp since they both provide DRM-free MP3s. I listen to everything in the used-to-be-iTunes Music app (not to be confused with Apple’s music streaming service, Apple Music), and I use iTunes Match to access my library across all my devices.¹

¹ I also have it all backed up to a Plex server on a Synology NAS in my house, and my kids listen to music on it through the Plexamp app, but I don’t because I still prefer using the iTunes/Apple Music desktop app 🤷‍♂️.

iTunes/Music keeps track of some song metadata, like a count of the number of times a song has been played:

All that metadata is stored in a big ol’ gross XML file. In days of iTunes, you could find it at ~/Music/iTunes/iTunes Library.xml; with Apple Music, it’s hidden in ~/Music/Music/Music Library/Library.musicdb. The easiest way to access it is to export a copy of it from Music with File > Library > Export Library…. It has a bunch of neat details about each file in your library:

<key>34813key>
<dict>
    <key>Track IDkey><integer>34813integer>
    <key>Namekey><string>In Another Lifestring>
    <key>Artistkey><string>The Killersstring>
    <key>Album Artistkey><string>The Killersstring>
    <key>Albumkey><string>Pressure Machinestring>
    <key>Genrekey><string>Alternative Rockstring>
    <key>Kindkey><string>MPEG audio filestring>
    <key>Sizekey><integer>7632215integer>
    <key>Total Timekey><integer>225724integer>
    <key>Disc Numberkey><integer>1integer>
    <key>Disc Countkey><integer>1integer>
    <key>Track Numberkey><integer>8integer>
    <key>Track Countkey><integer>11integer>
    <key>Yearkey><integer>2021integer>
    <key>Date Modifiedkey><date>2021-08-13T17:38:22Zdate>
    <key>Date Addedkey><date>2021-08-13T13:38:36Zdate>
    <key>Bit Ratekey><integer>268integer>
    <key>Sample Ratekey><integer>44100integer>
    <key>Commentskey><string>Amazon.com Song ID: REDACTEDstring>
    <key>Play Countkey><integer>97integer>
    <key>Play Datekey><integer>3815892895integer>
    <key>Play Date UTCkey><date>2024-12-01T15:14:55Zdate>
    <key>Ratingkey><integer>100integer>
    <key>Album Ratingkey><integer>100integer>
    <key>Album Rating Computedkey><true/>
    <key>Normalizationkey><integer>6230integer>
    <key>Artwork Countkey><integer>1integer>
    <key>Sort Album Artistkey><string>Killersstring>
    <key>Sort Artistkey><string>Killersstring>
    <key>Persistent IDkey><string>211319FB11435185string>
    <key>Track Typekey><string>Filestring>
    <key>Locationkey><string>file:///Users/andrew/Music/iTunes/iTunes%20Music/Music/The%20Killers/Pressure%20Machine/08%20In%20Another%20Life.mp3string>
    <key>File Folder Countkey><integer>5integer>
    <key>Library Folder Countkey><integer>1integer>
dict>

It keeps track of play count…

    <key>Play Countkey><integer>97integer>
    <key>Play Datekey><integer>3815892895integer>
    <key>Play Date UTCkey><date>2024-12-01T15:14:55Zdate>

…but unfortunately for Spotify Wrapped purposes, it overwrites the count and date information when you listen to a track—it doesn’t keep track of individual play counts. Here’s what the XML for “In Another Life” looked like before I listened to the track while writing this post:

<key>Play Countkey><integer>96integer>
<key>Play Datekey><integer>3808562446integer>
<key>Play Date UTCkey><date>2024-09-07T18:00:46Zdate>

That September 7th listen was erased from history once I hit play in December :(

That means it’s impossible to figure out how many times you listen to a track during a given time period—the play count only shows the most recent listen. With one XML export, you can’t find Spotify Wrapped-like details about listening habits in a single year.

However, if you have two XML exports, you can!

Calculating 2024 play counts with R

I played the long game this year and exported a copy of my iTunes/Music library on the morning of January 1 and stored the XML file in a folder on my computer. I then exported a copy of the library as it stands today. With these two library files, I can subtract the play count from January 1 from the play count today and find how many times I listened to each track. It still doesn’t give me date information—there’s no way to see time trends like what I was listening to in March or whatever²—but it gives me good data to work with.

² If I were super on top of things and cared that much, I could set up a script to automatically export a copy of the library every day and then reverse engineer daily listening data, but that seems like an excessive amount of work.

In the spirit of Caitlin’s tweet, I’m going to keep the analysis of this data as simple and straightforward as possible—just filtering, grouping, and summarizing.

The only bit of fancy R work comes at the beginning with parsing and cleaning the Apple Music XML files. The track information is deeply nested inside a bunch of XML layers and untangling all that requires some data wrangling. Fortunatley Simon Couch already did it in his 2022 analysis of his music, and he even made an accompanying package {wrapped} for doing it yourself. His package is designed to extract the play counts of all the music added in a given year, while I want the counts for all years, so I modified his wrap_library() function slightly to ignore the year argument and just parse everything. The modified function, now read_itunes_library() is below, for the morbidly curious:

R code for read_itunes_library()

Guide to generating and rendering computational markdown content programmatically with Quarto

Andrew Heiss — Mon, 04 Nov 2024 05:00:00 GMT

This year, I’ve helped build the Idaho Secretary of State’s office’s election results website for both the primary and general elections. Working with election data is a complex process, with each precinct reporting results to their parent counties, which all use different systems and software and candidate identifiers. Those county results then go into a central state-level database that state officials have access to for analysis and reporting.

In 2024, Idaho used a Quarto website to present the results for each statewide, congressional, and legislative contest (the URL is results.voteidaho.gov, though the Quarto website part probably won’t live there forever):

I may write something up someday about the process of building the website, depending on NDAs and security and copyright arrangements and whatnot. There are some neat technical details involved in the whole process, like complex static {targets} branching, remote {targets} storage, and replicating the structure of the real, live results database with a local DuckDB database for testing things without connecting to the live database. The short, sanitized version is that it uses two {targets} pipelines:

An ETL pipeline connects to the central results database, retrieves the latest results for each contest, standardizes the idenfitiers and cleans up the results, generates summary tables and maps for each race, and pushes the cleaned results to a shared datastore.
A website pipeline grabs the latest cleaned results from the shared datastore and builds a Quarto website.

The two pipelines run independently of each other every few minutes, and thanks to the magic of {targets}, if there are no updates (i.e. if there’s a lull in the reporting on election night), nothing needs to rebuild. It’s really neat.

There’s one Quarto/R Markdown trick that I used extensively when building the site: it’s possible to use R to automatically generate Quarto markdown before the entire document runs, allowing you to create parameterized templates for repeated elements.

Each race uses a tabset to show panels for (1) a table and (2) an interactive map of the results, and the reporting status for all the counties involved in the race is included in a callout block.

The markdown for each of these race results sections looks something like this:

## Parent district

::: {.callout-note title="X/Y counties fully reported" collapse="true"}
**Complete**: list of counties

**In progress**: list of counties
:::

### Race name

::: {.panel-tabset}
#### Table

```{r}
# R code for creating the table
```

#### Map

```{r}
# R code for creating the map
```
:::

For pages where there’s only one race, like the presidential election and the state’s constitutional amendment election, it’s trivial enough to just copy/paste that general template and replace the corresponding R code. But for the state-level legislative page, there are dozens of races. Repeating and modifying all that markdown 100+ times would be miserable. So instead, we programmatically generate the markdown for each race before the site is rendered so that Quarto thinks it’s working with hand-typed markdown.

Generating big chunks of markdown like this is a really cool approach with all sorts of applications (generate sections of a website; generate panel tabsets; generate presentation slides; etc.), but it’s a little unwieldly at first. So in this post, I’ll (1) show why this is trickier than just using regular R chunks with results="asis", (2) present a detailed step-by-step explanation of how to pre-render generated computational chunks, (3) provide a shorter, simpler, less-annotated example, and (4) give a more complex, less-annotated example.

Why not just use `results="asis"`?

It’s easy to use R/Python chunks to generate HTML or LaTeX or Typst or markdown and have that output appear in the rendered document—this is essentially what table-making packages like {tinytable}, {gt}, and {kableExtra} all do. To illustrate this, let’s load {gapminder} data:

library(tidyverse)
library(glue)
library(gapminder)

gapminder_2007 <- gapminder |> filter(year == 2007)

We can create a markdown list of all the continents in the dataset. Here I do it with paste0(), but I could also use a nicer wrapper like {pander}, which includes all sorts of functions for generating markdown.

This correctly makes a list, but it doesn’t get rendered like a list—it’s displayed as code chunk output:

continents <- gapminder_2007 |> 
  distinct(continent) |> 
  pull(continent)

cat(paste0("- ", continents, collapse = "\n"))

- Asia
- Europe
- Africa
- Americas
- Oceania

We can tell Quarto to treat the output of that chunk as raw markdown instead by setting the results="asis" chunk option:

```{r}
#| results: asis
cat(paste0("- ", continents, collapse = "\n"))
```

Asia
Europe
Africa
Americas
Oceania

That’s great and normal and I use this approach all the time for generating non-computational markdown.

Where this doesn’t work is when you have R chunks that need to be computed. To show this, let’s make a list showing π rounded to different digits using inline code chunks. We can manually type it like this:

- `r round(pi, 1)`
- `r round(pi, 2)`
- `r round(pi, 3)`
- `r round(pi, 4)`
- `r round(pi, 5)`

…which renders like this:

3.1
3.14
3.142
3.1416
3.1416

But that’s a lot of typing. So let’s generate it automatically. I’ll do this in a data frame, just because I like working that way, but you could also use standalone vectors (or even—gasp—a loop!):

pi_stuff <- tibble(digits = 1:5) |> 
  mutate(list_element = paste0("- `r round(pi, ", digits, ")`"))

# Everything is in the list_element column:
pi_stuff
## # A tibble: 5 × 2
##   digits list_element      
##                  
## 1      1 - `r round(pi, 1)`
## 2      2 - `r round(pi, 2)`
## 3      3 - `r round(pi, 3)`
## 4      4 - `r round(pi, 4)`
## 5      5 - `r round(pi, 5)`

We can put that column in a results="asis" chunk…

```{r}
#| results: asis
cat(paste0(pi_stuff$list_element, collapse = "\n"))
```

r round(pi, 1)
r round(pi, 2)
r round(pi, 3)
r round(pi, 4)
r round(pi, 5)

…and it renders correctly as markdown, but it doesn’t run the inline chunks :(

This is because of an issue with ordering: Quarto renders the chunk with cat(paste0(...)) and then moves on to the next chunk in the document. It won’t render the R chunks that the pi_stuff$list_element object contains because they’re all nested inside the parent chunk, and Quarto’s rendering process has moved on by the time the newly generated R chunks appear.

The trick is to pre-render the chunks before they officially show up in the document.¹ We can feed the collapsed pi_stuff$list_element object to knitr::knit() in an inline chunk, which makes Quarto render all the R chunks inside the chunk first, then place the output in the document to be rendered in the correct order like normal chunks:

¹ I haven’t found this formally documented anywhere—I stumbled across this approach in this gist from 2015.

Here's some regular markdown text. Let's show a list of 
differently-rounded values of $\pi$ for fun:

`r knitr::knit(text = paste0(pi_stuff$list_element, collapse = "\n"))`

Isn't that neat?

That markdown will render to this:

Here’s some regular markdown text. Let’s show a list of differently-rounded values of for fun:

3.1
3.14
3.142
3.1416
3.1416

Isn’t that neat?

Building a panel tabset with an inline chunk

Technically there was no need to use knitr::knit() in an inline chunk for that previous example. It would be easier to generate the text output within the pi_stuff data frame instead of in a bunch of inline chunks, and then show the results like normal with results="asis"

```{r}
#| results: asis
pi_stuff_easier <- tibble(digits = 1:5) |> 
  mutate(list_element = paste0("- ", round(pi, digits)))

cat(paste0(pi_stuff_easier$list_element, collapse = "\n"))
```

3.1
3.14
3.142
3.1416
3.14159

However, using knitr::knit(text = BLAH) in an inline chunk like this is a powerful trick that lets you do all sorts of more complex document generation automation. Let’s make a more complicated example with real data instead of a bunch of π rounding.

For this example, let’s make a panel tabset with a plot for each continent in the gapminder_2007 dataset we made earlier.

First, we’ll make a list of plots. This can be done any number of ways—I like using group_by() |> nest() and {purrr} functions like map(), but any way will work as long as you have a list of ggplot objects in the end.

continents_plots <- gapminder_2007 |> 
  group_by(continent) |> 
  nest() |> 
  ungroup() |> 
  # We could use map2(), but I like using pmap() just in case I need to
  # expand it beyond 2 things
  mutate(plot = pmap(
    lst(data, continent), 
    \(data, continent) {
      plot_title <- paste0("Health and wealth in ", continent)
      ggplot(data, aes(x = gdpPercap, y = lifeExp)) +
        geom_point() +
        scale_x_log10(labels = scales::label_dollar(accuracy = 1)) +
        labs(title = plot_title)
    }))

continents_plots
## # A tibble: 5 × 3
##   continent data              plot  
##                    
## 1 Asia         
## 2 Europe       
## 3 Africa       
## 4 Americas     
## 5 Oceania

The plots are all in the plot column in continents_plots, which is a list of ggplot objects. Here’s one of them:

continents_plots$plot[[3]]

To make a tabset with a panel for each continent, we need to write markdown like this:

::: {.panel-tabset}
### Continent 1

```{r}
#| label: panel-continent-a
#| echo: false
continents_plots$plot[[1]]
```

### Continent 2

```{r}
#| label: panel-continent-b
#| echo: false
continents_plots$plot[[2]]
```

### (…and so on…)

:::

We could just copy/paste those continent sections over and over, but that’s tedious and not very dynamic. Instead, we can create a little markdown template for each panel and generate all these chunks. To do that, we’ll use {glue}, which is a lot nicer for building strings than using paste0(), since it uses Python-style string interpolation. glue::glue() replaces any text inside {}s with the corresponding variable value:

# Some values
animals <- "cats"
number <- 12
location <- "house"

# Ugly paste() way
paste0("There are ", number, " blue ", animals, " in the ", location)
## [1] "There are 12 blue cats in the house"

# Nice glue() way
glue("There are {number} blue {animals} in the {location}")
## There are 12 blue cats in the house

If you need to use literal curly braces in the text, you can either double them or change the delimiters:

pkg_name <- "ggplot2"

# The curly braces disappear
glue("The {pkg_name} package is delightful")
## The ggplot2 package is delightful

# Double them to keep them
glue("The {{{pkg_name}}} package is delightful")
## The {ggplot2} package is delightful

# Or change the delimiter
glue("The {<>} package is delightful", .open = "<<", .close = ">>")
## The {ggplot2} package is delightful

Being able to change the delimiter is useful since we’ll need to generate chunks that start with ```{r}.

We can use glue() in a function that takes a continent name and a row number and generates a markdown tabset panel:

build_panel <- function(panel_title, plot_index) {
  chunk_label <- glue("panel-continent-{title}", title = janitor::make_clean_names(panel_title))

  output <- glue("
  ### <>

  ```{r}
  #| label: <>
  #| echo: false
  continents_plots$plot[[<>]]
  ```", .open = "<<", .close = ">>")

  output
}

First let’s make sure it works by itself:

build_panel("Africa", 3)

### Africa

```{r}
#| label: panel-continent-africa
#| echo: false
continents_plots$plot[[3]]
```

Yep!

Now we can iterate through the data frame of all the continents (continents_plots) and make a column that contains the markdown panel text.

continents_plots_with_text <- continents_plots |> 
  mutate(row = row_number()) |> 
  mutate(markdown = pmap_chr(
    lst(continent, row), 
    \(continent, row) build_panel(panel_title = continent, plot_index = row)
  ))

continents_plots_with_text

# A tibble: 5 × 5
  continent data              plot     row markdown                                                                                                   
                                                                                                                           
1 Asia              1 "### Asia\n\n```{r}\n#| label: panel-continent-asia\n#| echo: false\ncontinents_plots$plot[[1]]\n```"      
2 Europe            2 "### Europe\n\n```{r}\n#| label: panel-continent-europe\n#| echo: false\ncontinents_plots$plot[[2]]\n```"  
3 Africa            3 "### Africa\n\n```{r}\n#| label: panel-continent-africa\n#| echo: false\ncontinents_plots$plot[[3]]\n```"  
4 Americas          4 "### Americas\n\n```{r}\n#| label: panel-continent-americas\n#| echo: false\ncontinents_plots$plot[[4]]\n`…
5 Oceania            5 "### Oceania\n\n```{r}\n#| label: panel-continent-oceania\n#| echo: false\ncontinents_plots$plot[[5]]\n```"

Check out that new markdown column—it has a third level heading with the continent name, followed by an R chunk that will display the corresponding plot object from continents_plots. Here’s what one panel looks like:

cat(continents_plots_with_text$markdown[[1]])

### Asia

```{r}
#| label: panel-continent-asia
#| echo: false
continents_plots$plot[[1]]
```

Finally, we need to concatenate that column into one big string and include it as an inline chunk inside Quarto’s syntax for tabsets:

Health and wealth are related in each continent.

::: {.panel-tabset}

`r knitr::knit(text = paste0(continents_plots_with_text$markdown, collapse = "\n\n"))`

:::

Automatic tabset panels!

Here’s what it looks like when rendered:

Health and wealth are related in each continent.

Automatic tabset panels!

Perfect! We successfully generated a bunch of R chunks, pre-rendered them with Quarto, and then rendered the rest of the document.

Condensed example showing the evolution of a ggplot plot

Now that we’ve walked through the general process in detail, we’ll look at a less didactic example. Suppose you want to show the step-by-step process of creating a ggplot plot in a tabset panel or in a Revealjs Quarto slideshow. You could manually copy/paste a bunch of markdown over and over (ew), or you could generate chunks and make Quarto make the panels or slides for you (yay).

First we’ll make a list of plots. There’s actually a neat new package—{ggreveal}—that can create a list of intermediate plots automatically, but I’ll just do it manually here (though this process will work the same with {ggreveal}).

p1 <- ggplot(gapminder_2007, aes(x = gdpPercap, y = lifeExp, color = continent)) +
  geom_point()
p1_text <- glue("
  ggplot(gapminder_2007, aes(x = gdpPercap, y = lifeExp, color = continent)) +
    geom_point()")

p2 <- p1 + 
  scale_x_log10(labels = scales::dollar_format(accuracy = 1))
p2_text <- glue("
  {p1_text} +
    scale_x_log10(labels = scales::dollar_format(accuracy = 1))
")

p3 <- p2 + 
  scale_color_viridis_d(option = "plasma", end = 0.9)
p3_text <- glue('
  {p2_text} +
    scale_color_viridis_d(option = "plasma", end = 0.9)
')

p4 <- p3 + 
  theme_minimal()
p4_text <- glue('
  {p3_text} +
    theme_minimal()
')

p5 <- p4 + 
  labs(x = "GDP per capita", y = "Life expectancy", color = "Continent")
p5_text <- glue('
  {p4_text} +
    labs(x = "GDP per capita", y = "Life expectancy", color = "Continent")
')

plot_list <- list(p1, p2, p3, p4, p5)
plot_text <- list(p1_text, p2_text, p3_text, p4_text, p5_text)

plot_list <- tribble(
  ~plot, ~code_text, ~description,
  p1, p1_text, "Start with the initial plot…",
  p2, p2_text, "…use a logarithmic x-axis…",
  p3, p3_text, "…change the color palette…",
  p4, p4_text, "…change the theme…",
  p5, p5_text, "…and change the default labels"
)

We want the overall tabset to look something like this:

::: {.panel-tabset}

### Step 1

Short description

```{r}
plot_list[[1]]
```

```r
# Code here
```

### Step 2

Short description

```{r}
plot_list[[2]]
```

```r
# Code here
```

### …and so on

:::

…so next we’ll generate markdown for each panel:

panels <- map_chr(
  seq_len(nrow(plot_list)), 
  \(i) {
    glue("
    ### Step <>

    `r plot_list$description[[<>]]`

    ```{r}
    #| label: plot-panel-<>
    #| echo: false
    plot_list$plot[[<>]]
    ```

    ```r
    `r plot_list$code_text[[<>]]`
    ```", .open = "<<", .close = ">>")
  }
)

Finally we’ll wrap Quarto’s special markdown syntax for tabsets around these panels and then include the combined text as an inline chunk in the document:

Let's slowly build up the plot: ::: {.panel-tabset} `r knitr::knit(text = paste0(panels, collapse = "\n\n"))` :::

Let’s slowly build up the plot:

Step 1

Step 2

Step 3

Step 4

Step 5

Start with the initial plot…

ggplot(gapminder_2007, aes(x = gdpPercap, y = lifeExp, color = continent)) + geom_point()

…use a logarithmic x-axis…

ggplot(gapminder_2007, aes(x = gdpPercap, y = lifeExp, color = continent)) + geom_point() + scale_x_log10(labels = scales::dollar_format(accuracy = 1))

…change the color palette…

ggplot(gapminder_2007, aes(x = gdpPercap, y = lifeExp, color = continent)) + geom_point() + scale_x_log10(labels = scales::dollar_format(accuracy = 1)) + scale_color_viridis_d(option = "plasma", end = 0.9)

…change the theme…

ggplot(gapminder_2007, aes(x = gdpPercap, y = lifeExp, color = continent)) + geom_point() + scale_x_log10(labels = scales::dollar_format(accuracy = 1)) + scale_color_viridis_d(option = "plasma", end = 0.9) + theme_minimal()

…and change the default labels

ggplot(gapminder_2007, aes(x = gdpPercap, y = lifeExp, color = continent)) + geom_point() + scale_x_log10(labels = scales::dollar_format(accuracy = 1)) + scale_color_viridis_d(option = "plasma", end = 0.9) + theme_minimal() + labs(x = "GDP per capita", y = "Life expectancy", color = "Continent")

Condensed example of continent-level mini reports

Finally, let’s look at one more example that’s similar to what I used for the Idaho election results website, making a sort of miniature report for each continent. This time, instead of creating a tabset panel for each continent, we’ll make a whole markdown section for each continent, with a tabset panel included in each.

Each continent section will look something like this:

### Continent name ::: {.callout-note icon="false" title="X countries" collapse="true"} Comma-separated list of countries ::: ::: {.panel-tabset} #### Details ```{r} #| label: table-summary-continent #| echo: false # A table showing average GDP per capita and average life expectancy ``` #### Plot ```{r} #| label: plot-summary-continent #| echo: false # A plot showing the relationship between GDP per capita and life expectancy ``` :::

As before, we’ll translate this template into a glue() string, feed some data into it, and generate a bunch of R chunks.

TipChild documents

This template is getting gnarly with so many moving parts and so many common {glue} delimiters like {}s and []s and <>s. An alternative approach is to put the template in a separate child Quarto document and then pre-render it with knitr::knit_child(), which behaves just like the inline knitr::knit(file = BLAH) approach we’ve been using. Quarto has an official example of how to do it here: example and code.

First, we’ll make a data frame with all the different pieces we want to include, using the same group_by(continent) |> nest() approach as before:

continent_report_items <- gapminder_2007 |> group_by(continent) |> nest() |> ungroup() |> mutate(country_list = map_chr(data, \(x) knitr::combine_words(x$country))) |> mutate(n_countries = map_int(data, \(x) nrow(x))) |> mutate(summary_details = map(data, \(x) { x |> summarize( `Average GDP per capita` = mean(gdpPercap), `Average life expectancy` = mean(lifeExp), `Average population` = mean(pop) ) |> pivot_longer(everything(), names_to = "Statistic", values_to = "Value") |> mutate(Value = scales::comma_format()(Value)) })) |> mutate(plot = map(data, \(x) { x |> ggplot(aes(x = gdpPercap, y = lifeExp)) + geom_point() + scale_x_log10(labels = scales::label_dollar(accuracy = 1)) + labs(title = glue("Health and wealth in {continent}")) }))

Next we’ll make a function that generates the markdown output:

build_continent_report <- function(i) { name_for_labels <- janitor::make_clean_names(continent_report_items$continent[[i]]) # Quarto and RStudio and Positron all really struggle with syntax highlighting # and parsing when there are multiple ```s inside a string, so we can make # life easier by splitting the output into a few parts here, ensuring that # there's a maximum of one set of triple backticks output_first_part <- glue(' ### `r continent_report_items$continent[[<>]]` ::: {.callout-note icon="false" title="`r continent_report_items$n_countries[[<>]]` countries" collapse="true"} `r continent_report_items$country_list[[<>]]` :::', .open = "<<", .close = ">>") output_panel_details <- glue(' #### Details ```{r} #| label: table-summary-<> #| echo: false continent_report_items$summary_details[[<>]] |> knitr::kable() ```', .open = "<<", .close = ">>") output_panel_plot <- glue(' #### Plot ```{r} #| label: plot-summary-<> #| echo: false continent_report_items$plot[[<>]] ```', .open = "<<", .close = ">>") # Combine all the pieces output <- glue(' {output_first_part} ::: {{.panel-tabset}} {output_panel_details} {output_panel_plot} ::: ') output }

Finally, we’ll loop through each row in continent_report_items and generate the markdown report using the template…

continent_reports <- map_chr( seq_len(nrow(continent_report_items)), \(i) build_continent_report(i) )

…and include all the generated markdown in an inline chunk:

## Continent reports Check out all these automatically generated continent reports! `r knitr::knit(text = paste0(continent_reports, collapse = "\n\n"))`

That single inline chunk automatically generates dozens of inline and block chunks of R code before the full document goes through Quarto, which means all this output gets included in the final rendered version:

TipGenerated output

Continent reports

Check out all these automatically generated continent reports!

Asia

Note33 countries

Afghanistan, Bahrain, Bangladesh, Cambodia, China, Hong Kong, China, India, Indonesia, Iran, Iraq, Israel, Japan, Jordan, Korea, Dem. Rep., Korea, Rep., Kuwait, Lebanon, Malaysia, Mongolia, Myanmar, Nepal, Oman, Pakistan, Philippines, Saudi Arabia, Singapore, Sri Lanka, Syria, Taiwan, Thailand, Vietnam, West Bank and Gaza, and Yemen, Rep.

Details

Plot

Statistic Value

Average GDP per capita 12,473

Average life expectancy 71

Average population 115,513,752

Europe

Note30 countries

Albania, Austria, Belgium, Bosnia and Herzegovina, Bulgaria, Croatia, Czech Republic, Denmark, Finland, France, Germany, Greece, Hungary, Iceland, Ireland, Italy, Montenegro, Netherlands, Norway, Poland, Portugal, Romania, Serbia, Slovak Republic, Slovenia, Spain, Sweden, Switzerland, Turkey, and United Kingdom

Details

Plot

Statistic Value

Average GDP per capita 25,054

Average life expectancy 78

Average population 19,536,618

Africa

Note52 countries

Algeria, Angola, Benin, Botswana, Burkina Faso, Burundi, Cameroon, Central African Republic, Chad, Comoros, Congo, Dem. Rep., Congo, Rep., Cote d’Ivoire, Djibouti, Egypt, Equatorial Guinea, Eritrea, Ethiopia, Gabon, Gambia, Ghana, Guinea, Guinea-Bissau, Kenya, Lesotho, Liberia, Libya, Madagascar, Malawi, Mali, Mauritania, Mauritius, Morocco, Mozambique, Namibia, Niger, Nigeria, Reunion, Rwanda, Sao Tome and Principe, Senegal, Sierra Leone, Somalia, South Africa, Sudan, Swaziland, Tanzania, Togo, Tunisia, Uganda, Zambia, and Zimbabwe

Details

Plot

Statistic Value

Average GDP per capita 3,089

Average life expectancy 55

Average population 17,875,763

Americas

Note25 countries

Argentina, Bolivia, Brazil, Canada, Chile, Colombia, Costa Rica, Cuba, Dominican Republic, Ecuador, El Salvador, Guatemala, Haiti, Honduras, Jamaica, Mexico, Nicaragua, Panama, Paraguay, Peru, Puerto Rico, Trinidad and Tobago, United States, Uruguay, and Venezuela

Details

Plot

Statistic Value

Average GDP per capita 11,003

Average life expectancy 74

Average population 35,954,847

Oceania

Note2 countries

Australia and New Zealand

Details

Plot

Statistic Value

Average GDP per capita 29,810

Average life expectancy 81

Average population 12,274,974

↑ that’s how we were able to generate race-specific output for 100+ individual contests for the 2024 Idaho elections. It looks messy at first, but it’s a billion times easier to work with than copying/pasting markdown text and manually modifying the code in the chunks.

Statistic	Value
Average GDP per capita	12,473
Average life expectancy	71
Average population	115,513,752

Statistic	Value
Average GDP per capita	25,054
Average life expectancy	78
Average population	19,536,618

Statistic	Value
Average GDP per capita	3,089
Average life expectancy	55
Average population	17,875,763

Statistic	Value
Average GDP per capita	11,003
Average life expectancy	74
Average population	35,954,847

Statistic	Value
Average GDP per capita	29,810
Average life expectancy	81
Average population	12,274,974

Citation

BibTeX citation:

@online{heiss2024,
  author = {Heiss, Andrew},
  title = {Guide to Generating and Rendering Computational Markdown
    Content Programmatically with {Quarto}},
  date = {2024-11-04},
  url = {https://www.andrewheiss.com/blog/2024/11/04/render-generated-r-chunks-quarto/},
  doi = {10.59350/pa44j-cc302},
  langid = {en}
}

For attribution, please cite this work as:

Heiss, Andrew. 2024. “Guide to Generating and Rendering Computational Markdown Content Programmatically with Quarto.” November 4, 2024. https://doi.org/10.59350/pa44j-cc302.

Fun with Positron

Andrew Heiss — Mon, 08 Jul 2024 04:00:00 GMT

At the end of June 2024, Posit released a beta version of its next-generation IDE for data science: Positron. This follows Posit’s general vision for language-agnostic data analysis software: RStudio PBC renamed itself to Posit PBC in 2022 to help move away from a pure R focus, and Quarto is pan-lingual successor to R Markdown. Having the name of the main programming language in the title of things is out—providing more general tools is in.

Positron is essentially a specialized version of Microsoft’s Visual Studio Code, and is a fork of the underlying Code - OSS that powers VS Code. I’m super excited about this—in my own work, I use RStudio for most things R-related and VS Code for everything else (Stan, Python, HTML, CSS, Lua, LaTeX, Typst, etc.). VS Code is phenomenal and I love using it. It’s the best way to edit files on a remote server. It’s the best way to interact with Docker containers and Docker Compose. GitHub Copilot Chat is fantastic.

But for me, it’s never quite been a replacement for RStudio. Every couple months, I play around with trying to use VS Code for R work full time, but the constellation of VS Code R extensions (like the R extension and Radian for the terminal) and general R support has never been what I want, and I always end up going back to RStudio. Which is fine! I adore RStudio too and have been using it since it first came out in beta in February 2011 (13 years!).

Positron brings pretty much all the little R-related things that I love from RStudio and have missed in Visual Studio Code. The regular collection of VS Code’s R extensions and add-ons is no longer necessary, since Posit has created a custom R kernel—Ark—for any text editor or IDE with Jupyter support. It’s still a beta product and a little rough around the edges, but I’ve found that it really is the perfect blend of the best parts of RStudio and VS Code.

Below, following the example of Marc Dotson and Christopher Kenny, I want to highlight some of the neat new things Positron can do and share some of the settings, extensions, and other customizations I’ve been using for the past couple weeks.

Some cool new things

Positron brings RStudio’s best features to other languages like Python, like line-by-line code execution:

Line-by-line execution in Python
Line-by-line execution in R

…and the variables panel (equivalent to the Environment panel in RStudio):

Variables panel in Python
Variables panel in R

Variables panel in Python

Variables panel in R

You can also switch between different R and Python installations and versions. If you have rig installed, you can switch between different R versions, and Positron scans your computer at startup to find all the different Python virtual environments you have. I think eventually that little menu might also have better support for {renv} too. Here I can switch between R 4.4.0 and 4.3.3, as well as a bunch of random Python virtual environments and installations (lol python installations are the worst):

Version switcher

Plus, since it’s just a fancy version of VS Code, Positron supports pretty much everything VS Code can do, including making complex layouts. Use the little menu in the top right corner to set up your workspace however you want:

Customize layout

For instance, here’s an example of a fully armed and operational layout I’ve been using for one project:

Full ultrawide workspace

My settings

TipHow to change settings in Positron / VS Code

Positron (and VS Code in general) stores all its settings in a JSON file named settings.json that’s stored somewhere on your computer. On macOS it’s in ~/Library/Application Support/Positron/User/settings.json (see here for other operating systems). But you don’t need to ever remember that!

When you open Positron’s settings (with ⌘, on a Mac; using… something… on Windows), Positron provides a nice frontend for searching, managing, and changing settings so you don’t need to edit raw JSON if you don’t want to.

If you click on the document button in the top right corner, you can open the actual settings.json file in the editor and make changes there. This is the easiest way to share settings with other people (like in this blog post) or with yourself (you can commit settings.json to a git repository, for instance).

Positron settings page
Positron
Configuring and customizing Positron involved basically copying most of my settings from VS Code’s settings.json into Positron’s settings.json. Here’s everything I have set up, with comments explaining stuff. A few things to note in particular:

I set rstudio.keymap.enable to true to enable most of RStudio’s R-related keyboard shortcuts (like ⌘⌥I for a new chunk, ⌥- to insert <-, etc.).

I’m using GitHub’s Monaspace font because it looks neat and it has excellent font ligatures. I’ve enabled a bunch of different stylistic sets for the ligatures.

I’m a big fan of the Monokai color theme and use it in RStudio and VS Code. It’s easy enough to set in Positron too, but for mysterious unknown reasons, it uses colors differently and is overly aggressive in what gets colorized. Compare this ggplot code across three different Monokais (VS Code, Positron, and RStudio). The Positron version is incredibly green and pink, while VS Code and RStudio use color more sparingly.

Monokai highlighting in VS Code, Positron, and RStudio + Positron Dark

So for now, I’m using the Positron Dark theme instead, which does the best job of highlighting the things that RStudio did. It’s nice enough.

Here’s my settings.json file. Adapt from it however you want. All these settings are also accessible in the GUI too.

GUI for enabling RStudio keymapping

There are some extension-specific options at the bottom that I’ll explain below too.

settings.json

{ // Positron-specific settings // ------------------------------------------------------------------------- "rstudio.keymap.enable": true, "python.defaultInterpreterPath": "/opt/homebrew/bin/python", // Editor settings // ------------------------------------------------------------------------- // Fonts // Use GitHub's Monaspace (https://github.com/githubnext/monaspace) and enable ligatures "editor.fontFamily": "'Monaspace Argon Var'", "editor.fontSize": 12.5, "editor.fontLigatures": "'ss01', 'ss02', 'ss03', 'ss04', 'ss05', 'ss06', 'ss07', 'ss08', 'calt', 'dlig', 'liga'", // Theme // Monakai would be nice, but it has issues in Positron // "workbench.colorTheme": "Monokai", "workbench.colorTheme": "Default Dark Modern", // Use nicer icons "workbench.productIconTheme": "fluent-icons", "workbench.iconTheme": "material-icon-theme", // Highlight modified/unsaved tabs "workbench.editor.highlightModifiedTabs": true, // Add some rulers "editor.rulers": [ 80, 100 ], // Indent with two spaces, but only for R "[r]": { "editor.tabSize": 2 }, // Nicer handling of end-of-document newlines, via // https://rfdonnelly.github.io/posts/sane-vscode-whitespace-settings/ "files.insertFinalNewline": true, "editor.renderFinalNewline": "dimmed", "editor.renderWhitespace": "trailing", "files.trimFinalNewlines": true, "files.trimTrailingWhitespace": true, // Various editor settings "editor.formatOnPaste": true, "editor.detectIndentation": false, "editor.showFoldingControls": "always", "window.newWindowDimensions": "inherit", "editor.scrollBeyondLastLine": false, "window.title": "${activeEditorFull}${separator}${rootName}", "editor.tabSize": 4, "editor.wordWrap": "on", "editor.multiCursorModifier": "ctrlCmd", "editor.snippetSuggestions": "top", // Hide things from the global search menu and watcher "files.exclude": { "**/.Rhistory": true, "**/.Rproj": true, "**/.Rproj.user": true, "**/renv/library": true, "**/renv/local": true, "**/renv/staging": true }, "files.watcherExclude": { "**/.Rproj/*": true, "**/renv/library": true, "**/renv/local": true, "**/renv/staging": true }, // Sign git commits "git.enableCommitSigning": true, // Extension-specific settings // ------------------------------------------------------------------------- // Markdown linting settings (idk if this stuff even works with Quarto though) "markdownlint.config": { "default": true, "MD012": { "maximum": 2 }, "MD025": false, "MD041": false }, // Wrap at 80 columns with the "Rewrap" extension "rewrap.wrappingColumn": 80, // Hacky "Open Remote - SSH" settings "remote.SSH.serverDownloadUrlTemplate": "https://github.com/gitpod-io/openvscode-server/releases/download/openvscode-server-v${version}/openvscode-server-v${version}-${os}-${arch}.tar.gz", "remote.SSH.experimental.serverBinaryName": "openvscode-server", // Don't phone home for the "YAML" extension "redhat.telemetry.enabled": false, }

My keyboard shortcuts

TipHow to change keyboard shortcuts in Positron / VS Code

Changing keyboard shortcuts is just like changing settings. All the settings are stored in a JSON file (keybindings.json) located in a special folder on your computer, but you don’t have to work with raw JSON if you don’t want to.

The easiest way to get to the keyboard shortcut settings page is to open the Command Palette (⌘⇧P on macOS; ctrl + shift + p on Windows) and search for “Open Keyboard Shortcuts”:

This will give you a nice page for changing different settings. There are hundreds of possible shortcuts, but there’s a nice filtering system you can use to narrow things down.

If you click on the little document icon at the top, it will open the actual JSON file, just like with settings.json:

Accessing keyboard shortcuts from the Command Palette
Keyboard shortcut editor
Keyboard shortcuts as JSON
Enabling Positron’s RStudio Keymap option with rstudio.keymap.enable takes care of like 90% of my keyboard customization needs. Years ago when I first switched to VS Code, I changed several of RStudio’s keyboard shortcuts to match VS Code’s like ⌘/ for toggling commented code instead of RStudio’s default ⌘⇧C. Positron uses ⌘/ by default for comment toggling too, but when you enable the RStudio Keymap option, that gets overridden with ⌘⇧C, so I disable that.

RStudio also uses ⌘D for deleting a line, while VS Code uses it for adding text to a selection (i.e. if I select the word “the” in this document and then press ⌘D a bunch of times, it’ll add all those “the”s to the selection). The RStudio Keymap option adds ⌘D to delete the current line, so I disable that shortcut too to bring things back in line with standard VS Code.

Finally, I use iTerm2 for macOS for my systemwide terminal, and I have it configured with a global hotkey ^` so I can access the terminal from everywhere. This conflicts with VS Code’s and Positron’s terminal toggling shortcut, which is the same, so I change it to be ^⇧`.

Here’s my keybindings.json file. Like with settings.json, these are also accessible in the GUI.

Custom keyboard shortcuts

keybindings.json

[ { "key": "ctrl+alt+`", "command": "workbench.action.terminal.new", "when": "terminalProcessSupported || terminalWebExtensionContributedProfile" }, { "key": "ctrl+shift+`", "command": "-workbench.action.terminal.new", "when": "terminalProcessSupported || terminalWebExtensionContributedProfile" }, { "key": "ctrl+shift+`", "command": "workbench.action.terminal.toggleTerminal", "when": "terminal.active" }, { "key": "ctrl+`", "command": "-workbench.action.terminal.toggleTerminal", "when": "terminal.active" }, { "key": "shift+cmd+c", "command": "-editor.action.commentLine", "when": "config.rstudio.keymap.enable && editorTextFocus" }, { "key": "alt+cmd+q", "command": "rewrap.rewrapComment", "when": "editorTextFocus" }, { "key": "alt+q", "command": "-rewrap.rewrapComment", "when": "editorTextFocus" }, { "key": "cmd+d", "command": "-editor.action.deleteLines", "when": "config.rstudio.keymap.enable && editorTextFocus" } ]

My extensions

TipHow to install extensions in Positron / VS Code

Installing extensions in Positron / VS Code is super straightforward (see here). Click on the Extensions icon in the main Activity Bar, search for an extension, and click on “Install”. You can also disable or uninstall existing extensions from here.

Extension page for Stan
One of the best things about Positron is that it has access to most of VS Code’s extensions. Positron is not allowed to access Microsoft’s Visual Studio Extension Marketplace, but it can access (and is a major sponsor of) the alternative Open VSX Registry. With the exception of Microsoft’s extensions like GitHub Copilot, Dev Containers, and Remote - SSH, Open VSX had pretty much all the extensions that I already regularly use in VS Code.

The only minor VS Code extension I normally use that I couldn’t install in Positron was Stata Enhanced (Not that I even ever use Stata—I don’t have it installed on my computer and don’t have a license, but it’s nice to be able to open .do files and see syntax highlighting). Stata Enhanced isn’t listed at Open VSX, but I’ve opened an issue requesting that it gets listed.

Here’s what I use:

Managing other environments

Docker: Manage Docker containers and volumes; right click on docker-compose.yml files to spin them up and shut them down; syntax highlighting for Dockerfiles and Docker Compose

Open Remote - SSH: Connect to remote servers with SSH. This is bundled with Positron and there’s no need to install anything.

Text editing

Rewrap: Automatically add line breaks in long comments or text (I have it set to wrap at 80 characters using ⌘⌥Q)

Better Comments: Add special syntax highlighting for some types of comments like TODO, ?, !, and so on

Shebang Snippets: Provides snippets for adding shebang directives (e.g. type #!python to get #!/usr/bin/env python)

Viewers and syntaxes

Excel Viewer: View .xlsx files

vscode-pdf: View PDFs

Rainbow CSV: Does neat syntax highlighting for CSV files (highlighting each column with specific colors)

Stan: Syntax highlighting for Stan

YAML: Syntax highlighting for YAML

Lua: Syntax highlighting for Lua

markdownlint: Linting and style suggestions for Markdown

Theme stuff

Material Icon Theme: Customize the icons associated with specific file types in the file explorer

Fluent Icons: Customize the icons in the general Positron app (primarily the icons in the Activity Bar, like Explorer, Search, Source Control, etc.)

Remote connections with SSH

One of the best features of VS Code is its ability to connect to remote servers through SSH, but because that’s enabled with a special closed source Microsoft extension, it doesn’t work in Positron.

The Open Remote - SSH extension replicates Microsoft’s remote SSH extension, and it’s available at Open VSX. ~~However, it doesn’t work with Positron immediately—you’ll get an error when connecting.~~ This now works and there’s no need to install anything! (See this for a previous partial workaround.)

ImportantAugust 2024 update!

As of August 2024, Positron now bundles an SSH extension that Just Works™. If you have R or Python installed on a remote server, you can connect to it and run code remotely and it’s all great and wonderful now.

Things I still wish Positron could do

Positron is still in beta and is undergoing rapid development, and that’s totally fine. Even though it’s not a finished product yet, it works really really well.

There are still some things I wish it could do though. Some of these will eventually be addressed; some can’t because of Microsoft.

Packages panel: I love RStudio’s Packages panel. It’s so helpful for seeing which packages are currently installed, which versions are installed, updating existing packages, and installing new ones.

RStudio’s Packages panel

Nothing like this exists for Positron right now, but there’s discussion about how to build something like it (that would also work for Python).

Plot dimensions: (This will hopefully be addressed someday). In RStudio, when working with Quarto and R Markdown documents, inline images use the dimensions that you set in the chunk options, which makes it really easy to tinker with plot dimensions (i.e. changing from fig-width: 2.5 to fig-width: 2.75 to make sure labels fit in the plot area). Current, plots in Positron show up in the plots panel and use whatever dimensions that panel is set to use, either by manually resizing it or by using a dropdown menu with specific sizes:

Positron’s plot panel

It would be cool if Positron’s plot panel could pick up the dimensions specified in a Quarto document and auto-resize to match. For now, I’ve just been using the “Custom Size” option. If I want to preview an image that’s 5 inches wide and 3.75 inches tall, I convert the ratio of width/height to pixels. It’s not exact—there are issues with different DPIs and retina screens—but it at least shows the correct proportion.

Custom sizes in Positron’s plot panel

Remote editing and execution with Open Remote - SSH: It would be incredible if (1) it were a lot easier to install Open Remote - SSH, and (2) it were possible to run code on remote servers. I think they’re working on supporting this.

Similar to this, but less important to me because I don’t use Docker containers this way, VS Code can work with Docker containers with the Dev Containers extension, similar to SSH, using Docker environments to run R/Python locally. They might be working on supporting this some day.

GitHub Copilot Chat: Being able to chat with GitHub Copilot in VS Code is fantastic and it’s like only LLM thing I use. But it only works through a closed source extension by Microsoft and probably won’t ever work outside of VS Code proper.

Citation
BibTeX citation:
@online{heiss2024, author = {Heiss, Andrew}, title = {Fun with {Positron}}, date = {2024-07-08}, url = {https://www.andrewheiss.com/blog/2024/07/08/fun-with-positron/}, doi = {10.59350/zs7da-17c67}, langid = {en} }
For attribution, please cite this work as:
Heiss, Andrew. 2024. “Fun with Positron.” July 8, 2024. https://doi.org/10.59350/zs7da-17c67.

Calculating the proportion of US state borders that are coastlines

Andrew Heiss — Wed, 08 May 2024 04:00:00 GMT
A few days ago, my wife, a bunch of my kids, and I were huddled around a big wall map of the United States, joking about the relative unimportance of Rhode Island, the smallest state in the US. It’s one of the states I never ever think about:

Tweet by me from May 6, 2020

…and it’s just so small.

Amid the joking, my wife came to Rhode Island’s defense by declaring that even though it’s so small, it has one of the highest proportions of coastline to land borders. We all gave it a metaphorical gold star for being so maritime-y and moved on with our days.

But as I thought about it later, I got curious about how much of Rhode Island’s border really is coastline and how that proportion compares to other states. New England in general has lots of inlets and islands; North Carolina has the complex Outer Banks; Louisiana has the Mississippi Delta; Michigan is split into two parts and surrounded by the Great Lakes; Florida is Florida. Lots of other states have lots of coastline.

Using R, the extremely powerful {sf} package for working with geospatial data, and some high quality public domain geographic data, we can find the actual answers about coastline proportions. Spoilers: Rhode Island does really well, as expected.

But doing this is a lot more complicated than you might think, both for technical reasons and for philosophical reasons.

Let’s explore this data and make some pretty maps!

library(tidyverse) library(tigris) library(sf) library(rnaturalearth) library(patchwork) library(gt) clrs <- rcartocolor::carto_pal(12, "Prism") clr_ocean <- colorspace::lighten("#88CCEE", 0.7) # Custom ggplot theme to make pretty plots # Get the font at https://fonts.google.com/specimen/Overpass theme_map <- function() { theme_void(base_family = "Overpass Light") + theme( plot.title = element_text(family = "Overpass", face = "bold", hjust = 0.5), plot.subtitle = element_text(family = "Overpass", face = "plain", hjust = 0.5) ) } update_geom_defaults("text", list(family = "Overpass"))

Use US Census data on states and coastlines

At first glance, calculating the proportion of coastline borders in states feels fairly straightforward. We take state boundaries, find where they intersect with coastline boundaries, extract those overlapping sections, and voilá—we’re done.

The US Census Bureau even has shapefiles ready to use, like US state borders and the national coastline for 2023, and the {tigris} package makes it really easy to load that data directly into R.

First, let’s grab US state data at medium resolution (1:5 million) and calculate the length of each state’s border. To make calculations easier, we’ll change the projection to Albers, which measures distances in meters instead of the default NAD 83 decimal degree system.

census_states <- states( cb = TRUE, resolution = "5m", year = 2023, progress_bar = FALSE, keep_zipped_shapefile = TRUE ) |> st_transform(crs = st_crs("ESRI:102003")) |> # Albers mutate( border_length = st_perimeter(geometry), border_length_miles = units::set_units(border_length, "miles") ) census_states |> select(NAME, border_length, border_length_miles) ## Simple feature collection with 56 features and 3 fields ## Geometry type: MULTIPOLYGON ## Dimension: XY ## Bounding box: xmin: -10430000 ymin: -1685000 xmax: 3408000 ymax: 5141000 ## Projected CRS: USA_Contiguous_Albers_Equal_Area_Conic ## First 10 features: ## NAME border_length border_length_miles geometry ## 1 New Mexico 2389194 [m] 1484.6 [miles] MULTIPOLYGON (((-1231344 -5... ## 2 Puerto Rico 683648 [m] 424.8 [miles] MULTIPOLYGON (((3306526 -15... ## 3 Texas 6218137 [m] 3863.8 [miles] MULTIPOLYGON (((-1e+06 -570... ## 4 Kentucky 2097216 [m] 1303.1 [miles] MULTIPOLYGON (((584560 -886... ## 5 Ohio 1614692 [m] 1003.3 [miles] MULTIPOLYGON (((1094061 536... ## 6 Georgia 1950980 [m] 1212.3 [miles] MULTIPOLYGON (((939223 -230... ## 7 Arkansas 2117255 [m] 1315.6 [miles] MULTIPOLYGON (((122656 -111... ## 8 Oregon 2301152 [m] 1429.9 [miles] MULTIPOLYGON (((-2285910 94... ## 9 Pennsylvania 1579842 [m] 981.7 [miles] MULTIPOLYGON (((1287712 486... ## 10 Missouri 2357068 [m] 1464.6 [miles] MULTIPOLYGON (((19009 34499...

Next we’ll grab coastline data and also convert it to the meter-based Albers projection:

census_coastline <- coastline( year = 2023, progress_bar = FALSE, keep_zipped_shapefile = TRUE ) |> st_transform(crs = st_crs("ESRI:102003")) census_coastline ## Simple feature collection with 4236 features and 2 fields ## Geometry type: LINESTRING ## Dimension: XY ## Bounding box: xmin: -10430000 ymin: -1685000 xmax: 3408000 ymax: 5140000 ## Projected CRS: USA_Contiguous_Albers_Equal_Area_Conic ## First 10 features: ## NAME MTFCC geometry ## 1 Atlántico L4150 LINESTRING (3232962 -158263... ## 2 Atlántico L4150 LINESTRING (3282570 -157424... ## 3 Atlántico L4150 LINESTRING (3278311 -157345... ## 4 Atlántico L4150 LINESTRING (3283967 -157426... ## 5 Atlántico L4150 LINESTRING (3282914 -157442... ## 6 Atlántico L4150 LINESTRING (3282240 -157431... ## 7 Atlántico L4150 LINESTRING (3280919 -157406... ## 8 Atlántico L4150 LINESTRING (3276748 -157350... ## 9 Atlántico L4150 LINESTRING (3276791 -157363... ## 10 Atlántico L4150 LINESTRING (3286366 -157464...

This coastline data isn’t state-based—instead, each row represents a segment of the US coastline. Here’s what it looks like when plotted:

ggplot() + geom_sf(data = census_coastline, linewidth = 0.1) + labs(title = "US Census coastline data") + coord_sf(crs = st_crs("ESRI:102003")) + theme_map()

To find where the two maps intersect, we can use st_intersection(), and then we can calculate the length of each combined segment with st_length():

census_combined <- census_states |> st_intersection(census_coastline) |> mutate(coastline_length = st_length(geometry)) census_combined |> select(NAME, border_length, coastline_length, geometry) ## Simple feature collection with 1276 features and 3 fields ## Geometry type: GEOMETRY ## Dimension: XY ## Bounding box: xmin: -10430000 ymin: -1685000 xmax: 3408000 ymax: 5140000 ## Projected CRS: USA_Contiguous_Albers_Equal_Area_Conic ## First 10 features: ## NAME border_length coastline_length geometry ## 2 Puerto Rico 683648 [m] 119.450 [m] LINESTRING (3282570 -157424... ## 2.1 Puerto Rico 683648 [m] 12.122 [m] LINESTRING (3282914 -157442... ## 2.2 Puerto Rico 683648 [m] 722.088 [m] LINESTRING (3282240 -157431... ## 2.3 Puerto Rico 683648 [m] 5.332 [m] LINESTRING (3280919 -157406... ## 2.4 Puerto Rico 683648 [m] 2014.901 [m] LINESTRING (3279028 -157391... ## 2.5 Puerto Rico 683648 [m] 2446.737 [m] LINESTRING (3081263 -163791... ## 2.6 Puerto Rico 683648 [m] 12209.219 [m] MULTILINESTRING ((3306976 -... ## 2.7 Puerto Rico 683648 [m] 2804.085 [m] MULTILINESTRING ((3178034 -... ## 2.8 Puerto Rico 683648 [m] 7978.399 [m] MULTILINESTRING ((3230856 -... ## 2.9 Puerto Rico 683648 [m] 6004.139 [m] MULTILINESTRING ((3172336 -...

Some of these segments are thousands of miles; some are only a few miles. We can do some grouping and summarizing to collapse these into single values for each state. border_length is a state-level variable, not a border-segment-level variable, so it’s repeated in each of the rows of the combined dataset, so we only need to keep one of the values—here I keep the max, but min would work (since they’re the same).

coastline_length_by_state <- census_combined |> st_drop_geometry() |> # Stop worrying about geographic stuff group_by(NAME) |> summarize( total_coastline_length = sum(coastline_length), total_perimeter = max(border_length), prop_coastline = as.numeric(total_coastline_length / total_perimeter) ) coastline_length_by_state ## # A tibble: 35 × 4 ## NAME total_coastline_length total_perimeter prop_coastline ## [m] [m] ## 1 Alabama 192355. 1907670. 0.101 ## 2 Alaska 19223761. 29671399. 0.648 ## 3 American Samoa 175827. 172059. 1.02 ## 4 California 1129738. 4191420. 0.270 ## 5 Commonwealth of the Northern Mariana Islands 261865. 327378. 0.800 ## 6 Connecticut 77127. 574613. 0.134 ## 7 Delaware 110370. 433573. 0.255 ## 8 Florida 1842297. 3795469. 0.485 ## 9 Georgia 114280. 1950980. 0.0586 ## 10 Guam 120734. 133668. 0.903 ## # ℹ 25 more rows

At first glance, this looks fine. We have each state’s total border length, coastline length, and a proportion. Neat.

But some of these proportions look wrong, like Hawaiʻi:

coastline_length_by_state |> filter(NAME %in% c("Hawaii", "North Carolina")) ## # A tibble: 2 × 4 ## NAME total_coastline_length total_perimeter prop_coastline ## [m] [m] ## 1 Hawaii 712716. 1413312. 0.504 ## 2 North Carolina 248590. 3498277. 0.0711

According to this, only 50% of Hawaiʻi’s borders are on the coast. That’s obviously wrong—that state shares no borders with other states and it’s in the middle of the Pacific Ocean.

If we plot the two datasets, we can see what’s going on. First, the resolutions of the state data and coastline data don’t quite match. Check out Oʻahu here—the purple coastline crosses some bays and misses some of the landmass:

# Extract the state boundaries census_hi <- census_states |> filter(NAME == "Hawaii") # Extract the bounding box around the state so we can zoom in on the coastline map bbox_hi <- census_hi |> st_transform(st_crs("EPSG:4269")) |> st_bbox() ggplot() + geom_sf(data = census_hi, linewidth = 0.1, fill = clrs[5]) + geom_sf(data = census_coastline, aes(color = "US Census coastline data")) + annotate(geom = "text", x = I(0.425), y = I(0.89), label = "Oʻahu") + annotate( geom = "rect", xmin = I(0.35), xmax = I(0.5), ymin = I(0.63), ymax = I(0.85), color = clrs[12], fill = NA, linetype = "21" ) + scale_color_manual(values = c(clrs[11])) + labs(title = "Hawaiʻi", color = NULL) + coord_sf( xlim = bbox_hi[c(1, 3)], ylim = bbox_hi[c(2, 4)], crs = st_crs("EPSG:4269") ) + theme_map() + theme(legend.position = "top")

This mismatch between datasets is even more obvious if we look at a state like North Carolina. The actual landmass in the Outer Banks along the east coast is complex, with all sorts of inlets and a big long barrier island, but the Census simplifies it down substantially (with good reason(!) as we’ll see later):

census_nc <- census_states |> filter(NAME == "North Carolina") bbox_nc <- census_nc |> st_transform(st_crs("EPSG:4269")) |> st_bbox() ggplot() + geom_sf(data = census_nc, linewidth = 0.1, fill = clrs[6]) + geom_sf(data = census_coastline, aes(color = "US Census coastline data")) + annotate( geom = "rect", xmin = I(0.57), xmax = I(0.97), ymin = I(0.005), ymax = I(0.995), color = clrs[12], fill = NA, linetype = "21" ) + scale_color_manual(values = c(clrs[11])) + labs(title = "North Carolina", color = NULL) + coord_sf( xlim = bbox_nc[c(1, 3)], ylim = bbox_nc[c(2, 4)], crs = st_crs("EPSG:4269") ) + theme_map() + theme(legend.position = "top")

This mismatch in borders and coastlines matter a lot for calculations. We’re using st_intersection() to find where the two maps overlap. If the two maps are misaligned, we can’t identify the correct overlaps, which means we can’t identify the full coastal border. Remember how we calculated that only 50% of Hawaiʻi’s borders are coastlines? If we plot the coastal borders of Hawaiʻi, we can see why:

ggplot() + geom_sf(data = census_combined, linewidth = 0.2) + labs(title = "Incomplete coastline overlaps in Hawaiʻi") + coord_sf( xlim = bbox_hi[c(1, 3)], ylim = bbox_hi[c(2, 4)], crs = st_crs("EPSG:4269") ) + theme_map()

There are a couple issues here. The dots scattered along the borders are a sign that there’s some misalignment between the two maps. The borders for the states sometimes cross the coastline borders in a single point instead of following the coast exactly. But the even bigger issue is that these coastlines look sketched out—there are so many major gaps that roughly 50% of the borders are missing. That’s again because of misalignment with the two maps. The borders don’t overlap exactly, so st_intersection() can’t pick them up.

So what do we do? We could try finding different coastline data that matches the same resolution as the Census data. The National Oceanic and Atmospheric Administration (NOAA) and the US Geological Survey (USGS) each have their own shoreline datasets, and we could download those and load them into R and hope that they align with one of the Census’s state maps. But they don’t.

So we give up.

Use the lack of land as the coastline

Just kidding. There’s a better solution.

We already kind of have coastline data embedded in the state map data. For states that border an ocean or lake, anywhere the blue of the water touches the land is technically a coastline. Take Hawaiʻi, for instance, where the whole state border is the coastline:

ggplot() + geom_sf(data = census_hi, linewidth = 0.1, fill = clrs[5]) + labs(title = "The islands of Hawaiʻi") + coord_sf(crs = st_crs("EPSG:4269")) + # NAD83 theme_map() + theme( panel.background = element_rect(fill = clr_ocean), plot.title = element_text(margin = margin(6.5, 0, 6.5, 0)) )

When looking at states with interior borders (like the other 49 states), though, the state map data has no way to distinguish which of those borders touch the ocean or touch other states. The east coast of North Carolina touches the ocean, but the northern, southern, and western borders do not. We need to somehow figure out which borders don’t touch other states.

ggplot() + geom_sf(data = census_nc, linewidth = 0.1, fill = clrs[6]) + annotate(geom = "text", x = I(0.6), y = I(0.89), label = "↑ Virginia ↑") + annotate(geom = "text", x = I(0.58), y = I(0.3), label = "↓ South Carolina ↓", angle = 313) + annotate(geom = "text", x = I(0.14), y = I(0.5), label = "↓ Georgia") + annotate(geom = "text", x = I(0.27), y = I(0.64), label = "← Tennessee") + labs(title = "The “island” of North Carolina") + coord_sf(crs = st_crs("EPSG:4269")) + # NAD83 theme_map() + theme( panel.background = element_rect(fill = clr_ocean), plot.title = element_text(margin = margin(6.5, 0, 6.5, 0)) )

We do have data about Virginia, South Carolina, Tennessee, and Georgia, though, so there is a way to know that those borders aren’t ocean borders.

One way to determine which borders are interior and which ones are coastal is to create a big unified shape of the United States and then use the shape of North Carolina as a kind of cookie cutter. We can assume that any exposed edges are coasts.

First, we’ll make our big unified country shape:

census_us_giant <- census_states |> st_union() |> st_transform(crs = st_crs("ESRI:102003")) ggplot() + geom_sf(data = census_us_giant, linewidth = 0, fill = clrs[1]) + labs(title = "One big US-shaped shape") + theme_map() + theme( panel.background = element_rect(fill = clr_ocean), plot.title = element_text(margin = margin(6.5, 0, 6.5, 0)) )

Before trying this with North Carolina, we’ll test the cookie cutter selection with Hawaiʻi, since we know that 100% of its borders are coastlines. If we use the Hawaiʻi shape to take a chunk out of the overall US shape, we get…

hi_ocean_border_census <- census_hi |> st_difference(census_us_giant) ggplot() + geom_sf(data = hi_ocean_border_census) + labs(title = "lol nothing") + theme_map() + theme( panel.border = element_rect(color = "black", fill = clr_ocean), plot.title = element_text(margin = margin(6.5, 0, 6.5, 0)) )

…nothing.

The state cookie cutter was too perfect and selected right up to the edge of the overall country, leaving nothing.

To fix this, we can expand the state shape just a tiiiiiiny bit. We can use st_buffer() to add a tiny amount of distance all around the shape—since we’re using the Albers projection, we’re working in meters, so let’s add just 1 millimeter around the border before finding the difference:

hi_ocean_border_census <- census_hi |> st_buffer(dist = 0.001) |> st_difference(census_us_giant) ggplot() + geom_sf(data = hi_ocean_border_census, linewidth = 0.2) + labs(title = "Hawaiʻi’s coastal borders extracted from US shape") + coord_sf( xlim = bbox_hi[c(1, 3)], ylim = bbox_hi[c(2, 4)], crs = st_crs("EPSG:4269") ) + theme_map()

Perfect.

Well, almost perfect. We’ve identified the coastal borders, but adding the buffer actually distorts things when we calculate the length of the border.

Let’s find the coast-to-border proportion for Hawaiʻi, which should be 100%. First we’ll find the perimeter of the state:

hi_border_perimeter <- census_hi |> st_perimeter() hi_border_perimeter ## 1413312 [m] units::set_units(hi_border_perimeter, "miles") ## 878.2 [miles]

Hawaiʻi’s border is 1,413,312 meters, or 878 miles. Sounds reasonable.

Next we’ll find the perimeter of the coastline:

hi_ocean_border_perimeter <- hi_ocean_border_census |> st_perimeter() hi_ocean_border_perimeter ## 2826623 [m] units::set_units(hi_ocean_border_perimeter, "miles") ## 1756 [miles]

Hrm. Hawaiʻi’s coastline is 2,826,623 meters, or 1,756 miles. That’s actually exactly twice the correct border length:

as.numeric(hi_ocean_border_perimeter / hi_border_perimeter) ## [1] 2

This happened because of the 1 mm buffer that we added around the state shape. Adding the buffer transformed the line into a polygon—a super tiny 1 mm-narrow polygon, but a polygon nonetheless. That means the border technically has a top and sides and a bottom. When calculating the length or perimeter of the border, we’re double counting because we’re getting both the top and the bottom of the hyper-thin polygon.

To better illustrate what’s going on, let’s add a three kilometer buffer around the borders.

hi_ocean_border_census_huge_buffer <- census_hi |> st_buffer(dist = 3000) |> st_difference(census_us_giant) ggplot() + geom_sf( data = hi_ocean_border_census_huge_buffer, linewidth = 0.2, fill = colorspace::lighten(clrs[5], 0.5) ) + labs(title = "Hawaiʻi’s borders with a 3 km buffer") + coord_sf( xlim = bbox_hi[c(1, 3)], ylim = bbox_hi[c(2, 4)], crs = st_crs("EPSG:4269") ) + theme_map()

Calculating the perimeter of these borders will add the inside ring and the outside ring, effectively doubling the distance.

It’s less obvious that this doubling is happening when we add just 1 mm, but it is, and it’s leading to incorrect calculations. Fortunately it’s easy to adjust—we can halve the ocean border distance. Coastal borders comprise 100% of Hawaiʻi’s state borders, as expected:

as.numeric((hi_ocean_border_perimeter / 2) / hi_border_perimeter) |> scales::label_percent()() ## [1] "100%"

Let’s do the same thing with North Carolina:

nc_ocean_border_census <- census_nc |> st_buffer(dist = 0.001) |> st_difference(census_us_giant) ggplot() + geom_sf(data = nc_ocean_border_census, linewidth = 0.2) + labs(title = "North Carolina’s coastal borders extracted from US shape") + coord_sf(crs = st_crs("EPSG:4269")) + theme_map()

Perfect! Those are all the state borders that don’t touch other states. That’s the coastline.

We can visualize this better by plotting the larger country shape, the North Carolina shape, and the coastal border:

ggplot() + geom_sf(data = census_us_giant, fill = clrs[1], alpha = 0.4) + geom_sf(data = census_nc, linewidth = 0.1, fill = clrs[6]) + geom_sf( data = nc_ocean_border_census, linewidth = 0.4, aes(color = "Coastal border extracted from US shape"), key_glyph = draw_key_path ) + annotate(geom = "text", x = I(0.6), y = I(0.89), label = "↑ Virginia ↑") + annotate(geom = "text", x = I(0.58), y = I(0.3), label = "↓ South Carolina ↓", angle = 313) + annotate(geom = "text", x = I(0.14), y = I(0.5), label = "↓ Georgia") + annotate(geom = "text", x = I(0.27), y = I(0.64), label = "← Tennessee") + scale_color_manual(values = c(clrs[8])) + labs(title = "North Carolina and its coastal border", color = NULL) + coord_sf( xlim = bbox_nc[c(1, 3)], ylim = bbox_nc[c(2, 4)], crs = st_crs("EPSG:4269") ) + theme_map() + theme( panel.background = element_rect(fill = clr_ocean), plot.title = element_text(margin = margin(6.5, 0, 6.5, 0)), legend.key.height = unit(0.5, "lines"), legend.key = element_rect(fill = NA, color = NA), legend.position = "top" )

And we can find the proportion of the state’s borders that are coastline:

nc_border_perimeter <- census_nc |> st_perimeter() nc_border_perimeter ## 3498277 [m] units::set_units(nc_border_perimeter, "miles") ## 2174 [miles] nc_ocean_border_perimeter <- nc_ocean_border_census |> st_perimeter() nc_ocean_border_perimeter / 2 ## 1930852 [m] units::set_units(nc_ocean_border_perimeter / 2, "miles") ## 1200 [miles] as.numeric((nc_ocean_border_perimeter / 2) / nc_border_perimeter) |> scales::label_percent()() ## [1] "55%"

This approach is great. We don’t need to worry about making sure the coastline map matches the resolution of the state map, since we’re using just one map. Everything aligns perfectly automatically.

BUT there’s one more wrinkle to worry about. This approach gets trickier with states that touch other countries, like Washington, where the western border touches the ocean, the northern border touches Canada, and the eastern and southern borders touch other states. The state data knows about Idaho and Oregon, but it doesn’t know that Canada is there. Washington looks like it has an exterior, ocean-facing western border—Canada is replaced with an ocean.

census_wa <- census_states |> filter(NAME == "Washington") bbox_wa <- census_wa |> st_transform(st_crs("EPSG:4269")) |> st_bbox() wa_ocean_border_census <- census_wa |> st_buffer(dist = 0.001) |> st_difference(census_us_giant) ggplot() + geom_sf(data = census_us_giant, fill = clrs[1], alpha = 0.4) + geom_sf(data = census_wa, linewidth = 0.1, fill = clrs[3]) + geom_sf( data = wa_ocean_border_census, linewidth = 0.4, aes(color = "Coastal border extracted from US shape"), key_glyph = draw_key_path ) + annotate(geom = "text", x = I(0.64), y = I(0.91), label = "↑ Canada (British Columbia) ↑") + annotate(geom = "text", x = I(0.54), y = I(0.18), label = "↓ Oregon ↓") + annotate(geom = "text", x = I(0.87), y = I(0.5), label = "Idaho →") + scale_color_manual(values = c(clrs[8])) + labs(title = "Washington and the fake Canadian ocean", color = NULL) + coord_sf( xlim = bbox_wa[c(1, 3)], ylim = bbox_wa[c(2, 4)], crs = st_crs("EPSG:4269") ) + theme_map() + theme( panel.background = element_rect(fill = clr_ocean), plot.title = element_text(margin = margin(6.5, 0, 6.5, 0)), legend.key.height = unit(0.5, "lines"), legend.key = element_rect(fill = NA, color = NA), legend.position = "top" )

The big US-shaped polygon isn’t big enough for identifying land borders in other countries. In order to identify all the ocean-facing state borders, we need data about Canada and Mexico, and the US Census doesn’t provide that.

So we give up.

Use data from Natural Earth

Just kidding. There’s a better solution.

The Census doesn’t have data on other countries, but other data sources do! The incredible Natural Earth project has a ton of shapefiles for the entire world, both physical (i.e. rivers, coastlines, etc.) and cultural, like borders for countries (what it calls Admin-0) and states/provinces (what it calls Admin-1). It offers three levels of resolution: 1:10 million (high resolution), 1:50 million (medium resolution), and 1:110 million (low resolution).

Plus, like {tigris}, the {rnaturalearth} package makes it really easy to load that data directly into R.

Let’s follow the same process for Hawaiʻi, North Carolina, and Washington using Natural Earth data instead of Census data.

First, we’ll grab high resolution (1:10 million) maps from the Admin-1 (states and provinces) data:

if (!file.exists("ne_data/ne_10m_admin_1_states_provinces_lakes.shp")) { ne_download( type = "admin_1_states_provinces_lakes", scale = 10, destdir = "ne_data", load = FALSE ) } states_provinces <- ne_load( type = "admin_1_states_provinces_lakes", scale = 10, destdir = "ne_data" )

We only need to work with three countries, not the entire world, so we’ll filter it down to just North America before joining everything into a single shape with st_union():

na_giant <- states_provinces |> filter(admin %in% c("United States of America", "Canada", "Mexico")) |> st_union() |> st_transform(crs = st_crs("ESRI:102003")) # Convert to Albers for meters ggplot() + geom_sf(data = na_giant, linewidth = 0, fill = clrs[10]) + labs(title = "One big North America-shaped shape") + coord_sf(crs = st_crs("ESRI:102003")) + theme_map() + theme( panel.background = element_rect(fill = clr_ocean), plot.title = element_text(margin = margin(6.5, 0, 6.5, 0)) )

(That’s so pretty.)

Next, we’ll extract all the US states from the Natural Earth data, since we want this big North American shape to match the states exactly:

ne_states <- states_provinces |> filter(admin == "United States of America") |> st_transform(crs = st_crs("ESRI:102003")) ggplot() + geom_sf(data = na_giant, linewidth = 0, fill = clrs[10]) + geom_sf(data = ne_states, linewidth = 0.05, fill = clrs[8], color = "white") + labs(title = "The United States overlaid on North America") + coord_sf(crs = st_crs("ESRI:102003")) + theme_map() + theme( panel.background = element_rect(fill = clr_ocean), plot.title = element_text(margin = margin(6.5, 0, 6.5, 0)) )

Now we’ll go through the same process as before, using state shapes as cookie cutters from the larger North American shape and extracting non-interior borders.

Hawaiʻi

Here are Hawaiʻi’s extracted borders—because we’re working with high resolution data, we get a lot more detail, including all the Northwestern Hawaiian Islands:

ne_hi <- ne_states |> filter(name == "Hawaii") bbox_ne_hi <- ne_hi |> st_transform(st_crs("EPSG:4269")) |> st_bbox() hi_ocean_border_ne <- ne_hi |> st_buffer(dist = 0.001) |> st_difference(na_giant) ggplot() + geom_sf(data = hi_ocean_border_ne, linewidth = 0.2) + labs(title = "Hawaiʻi’s coastal borders extracted from North America shape") + coord_sf( xlim = bbox_ne_hi[c(1, 3)], ylim = bbox_ne_hi[c(2, 4)], crs = st_crs("EPSG:4269") ) + theme_map()

As before, we can calculate the perimeter of the state, the perimeter of the coastline, and find the proportion. It’s 100%, as expected:

# State border hi_border_perimeter_ne <- ne_hi |> st_perimeter() hi_border_perimeter_ne ## 1458906 [m] units::set_units(hi_border_perimeter_ne, "miles") ## 906.5 [miles] # Coastal border hi_ocean_border_perimeter_ne <- hi_ocean_border_ne |> st_perimeter() hi_ocean_border_perimeter_ne / 2 ## 1458906 [m] units::set_units(hi_ocean_border_perimeter_ne / 2, "miles") ## 906.5 [miles] # Proportion as.numeric((hi_ocean_border_perimeter_ne / 2) / hi_border_perimeter_ne) |> scales::label_percent()() ## [1] "100%"

North Carolina

The same process works for North Carolina

ne_nc <- ne_states |> filter(name == "North Carolina") bbox_ne_nc <- ne_nc |> st_transform(st_crs("EPSG:4269")) |> st_bbox() nc_ocean_border_ne <- ne_nc |> st_buffer(dist = 0.001) |> st_difference(na_giant) ggplot() + geom_sf(data = nc_ocean_border_ne, linewidth = 0.2) + labs(title = "North Carolina’s coastal borders\nextracted from North America shape") + coord_sf(crs = st_crs("EPSG:4269")) + theme_map()

Here’s that extracted coastline with the rest of the state:

ggplot() + geom_sf(data = na_giant, fill = clrs[10], alpha = 0.4) + geom_sf(data = ne_nc, linewidth = 0.1, fill = clrs[6]) + geom_sf( data = nc_ocean_border_ne, linewidth = 0.4, aes(color = "Coastal border extracted from North America shape"), key_glyph = draw_key_path ) + annotate(geom = "text", x = I(0.6), y = I(0.89), label = "↑ Virginia ↑") + annotate(geom = "text", x = I(0.58), y = I(0.3), label = "↓ South Carolina ↓", angle = 313) + annotate(geom = "text", x = I(0.14), y = I(0.5), label = "↓ Georgia") + annotate(geom = "text", x = I(0.27), y = I(0.64), label = "← Tennessee") + scale_color_manual(values = c(clrs[8])) + labs(title = "The actual North Carolina coastline", color = NULL) + coord_sf( xlim = bbox_ne_nc[c(1, 3)], ylim = bbox_ne_nc[c(2, 4)], crs = st_crs("EPSG:4269") ) + theme_map() + theme( panel.background = element_rect(fill = clr_ocean), plot.title = element_text(margin = margin(6.5, 0, 6.5, 0)), legend.key.height = unit(0.5, "lines"), legend.key = element_rect(fill = NA, color = NA), legend.position = "top" )

Now we can calculate proportion of coastline:

# State border nc_border_perimeter_ne <- ne_nc |> st_perimeter() nc_border_perimeter_ne ## 4095093 [m] units::set_units(nc_border_perimeter_ne, "miles") ## 2545 [miles] # Coastal border nc_ocean_border_perimeter_ne <- nc_ocean_border_ne |> st_perimeter() nc_ocean_border_perimeter_ne / 2 ## 2596482 [m] units::set_units(nc_ocean_border_perimeter_ne / 2, "miles") ## 1613 [miles] # Proportion as.numeric((nc_ocean_border_perimeter_ne / 2) / nc_border_perimeter_ne) |> scales::label_percent()() ## [1] "63%"

Sharp-eyed readers will notice something odd, though! Using Census maps at 1:5 million resolution, we found that 55% of North Carolina’s borders were coastline. Using Natural Earth maps at 1:10 million resolution, it becomes 63%. I’ll talk more about that discrepancy later in this post (spoiler: it’s because of the new-to-me Coastline Paradox). The detail in the two maps is different, so the amount of landmass visible along the coastline is different, yielding different perimeters and distances.

Washington

When we were using Census data, we couldn’t use this process with Washington because Canada was invisible and was being treated as an ocean. Now that we have data from Natural Earth, Canada can be accounted for and we’ll get the correct ocean-facing border.

ne_wa <- ne_states |> filter(name == "Washington") bbox_ne_wa <- ne_wa |> st_transform(st_crs("EPSG:4269")) |> st_bbox() wa_ocean_border_ne <- ne_wa |> st_buffer(dist = 0.001) |> st_difference(na_giant) ggplot() + geom_sf(data = wa_ocean_border_ne, linewidth = 0.2) + labs(title = "Washington’s coastal borders\nextracted from North America shape") + coord_sf(crs = st_crs("EPSG:4269")) + theme_map()

Here’s that coastline in context:

ggplot() + geom_sf(data = na_giant, fill = clrs[10], alpha = 0.4) + geom_sf(data = ne_wa, linewidth = 0.1, fill = clrs[3]) + geom_sf( data = wa_ocean_border_ne, linewidth = 0.4, aes(color = "Coastal border extracted from North America shape"), key_glyph = draw_key_path ) + annotate(geom = "text", x = I(0.64), y = I(0.91), label = "↑ Canada (British Columbia) ↑") + annotate(geom = "text", x = I(0.54), y = I(0.18), label = "↓ Oregon ↓") + annotate(geom = "text", x = I(0.87), y = I(0.5), label = "Idaho →") + scale_color_manual(values = c(clrs[8])) + labs(title = "The actual Washington coastline", color = NULL) + coord_sf( xlim = bbox_ne_wa[c(1, 3)], ylim = bbox_ne_wa[c(2, 4)], crs = st_crs("EPSG:4269") ) + theme_map() + theme( panel.background = element_rect(fill = clr_ocean), plot.title = element_text(margin = margin(6.5, 0, 6.5, 0)), legend.key.height = unit(0.5, "lines"), legend.key = element_rect(fill = NA, color = NA), legend.position = "top" )

Since the northern state border is no longer seen as a coastal border, can calculate the correct proportion of coastline:

# State border wa_border_perimeter_ne <- ne_wa |> st_perimeter() wa_border_perimeter_ne ## 3913863 [m] units::set_units(wa_border_perimeter_ne, "miles") ## 2432 [miles] # Coastal border wa_ocean_border_perimeter_ne <- wa_ocean_border_ne |> st_perimeter() wa_ocean_border_perimeter_ne / 2 ## 2591575 [m] units::set_units(wa_ocean_border_perimeter_ne / 2, "miles") ## 1610 [miles] # Proportion as.numeric((wa_ocean_border_perimeter_ne / 2) / wa_border_perimeter_ne) |> scales::label_percent()() ## [1] "66%"

All states

Now that we know that the Natural Earth approach works, let’s apply it to all the states. We’ll nest the geographic data into a list column with a cell for each state, then extract the borders for each state

ne_coastline <- ne_states |> group_by(name) |> nest() |> mutate(ocean_only = map(data, ~ { .x |> st_buffer(dist = 0.001) |> st_difference(na_giant) })) |> unnest(ocean_only) |> ungroup() |> # This special column somehow lost its specialness st_set_geometry("geometry")

ggplot() + geom_sf(data = na_giant, linewidth = 0, fill = clrs[10], alpha = 0.5) + geom_sf( data = ne_coastline, linewidth = 0.2, aes(color = "Coastal border extracted from North America shape"), key_glyph = draw_key_path ) + scale_color_manual(values = c(clrs[8])) + labs(title = "All US coastal borders", color = NULL) + coord_sf(crs = st_crs("ESRI:102003")) + theme_map() + theme( panel.background = element_rect(fill = clr_ocean), plot.title = element_text(margin = margin(6.5, 0, 6.5, 0)), legend.key.height = unit(0.5, "lines"), legend.key = element_rect(fill = NA, color = NA), legend.position = "top" )

That looks great!

Finally, let’s put this in a table. We’ll calculate the perimeter of all the states, then join that data to the coastal data, and then calculate the proportion of coastline.

ne_state_border_lengths <- ne_states |> mutate(border_length = st_perimeter(geometry)) |> st_drop_geometry() |> select(name, border_length) coastline_length_by_state_ne <- ne_coastline |> mutate(coastline_length = st_perimeter(geometry) / 2) |> st_drop_geometry() |> left_join(ne_state_border_lengths, by = join_by(name)) |> mutate(prop_coastline = as.numeric(coastline_length / border_length)) |> mutate(rank = rank(-prop_coastline)) |> mutate( across(c(border_length, coastline_length), list(miles = ~units::set_units(., "miles"))) ) |> select( name, rank, prop_coastline, starts_with("border_length"), starts_with("coastline_length") )

This is fantastic! Hawaiʻi has the highest proportion of coastline, for obvious reasons, and an astounding 95% of Alaska’s borders touch the ocean, even though a huge chunk of the state shares a land border with Canada—there are just so many islands and inlets.

In the contiguous United States, island-y states like Florida, Michigan, and Louisiana have the highest proportion of coastline. Despite their small size, the coastal New England states like New Jersey, Massachusetts, and Rhode Island have a lot of coastline relative to the rest of their borders. Check out Rhode Island at #10—my wife’s offhand observation was pretty accurate! (These columns are sortable.)

Code

Andrew Heiss's blog

Generating universes within universes with a single seed

Introduction

Seeds and reproducible randomness

My (somewhat incorrect) mental model of how seeds work

Making “random” numbers with an equation

Live interactive playground

Cycles and fancier algorithms

Why does it matter if “random” numbers aren’t actually random?

You’re limiting yourself to narrow, known universes

You can seed hack and get any values you want

Real world bad things can happen because of pseudorandom numbers

Can computers even create true randomness?

Moving a mouse around

Lava lamps

Atmospheric noise

How I use true randomness in my own work

“…as an ook cometh of a litel spyr…”

References

Citation

How to automatically convert TikZ images to SVG with Quarto

Quick, easy, ugly diagrams with Mermaid and Graphviz

Nicer diagrams with TikZ

Converting TikZ to both PDF and SVG with Quarto and diagram

Extra settings and complete example

Citation

How to make your data analysis life easier using Positron, Raycast, and Espanso

Positron stuff

How to open a folder as a Positron project with macOS Quick Actions

Projects in Positron

Quickly opening folders in Positron

Workspaces

Project Manager + Raycast

Open a folder in Positron macOS Quick Actions

Making this native?

Citation

How to use Positron’s Connections Pane with DuckDB

DuckDB, {DBI}, and the difficulty of discerning data in a database

DuckDB, {connections}, and the magical Connections Pane

Bonus: Better support for DuckDB in the Connections Pane

The whole game

Citation

Use Positron to run R inside a Docker image through SSH

How to use Docker-based R normally

How to use Docker-based R from Positron

Citation

Open files in external programs with Positron or Visual Studio Code

Quarto and Typora

Running external processes with tasks

Launch tasks more quickly

Do more with tasks

Citation

How to use a histogram as a legend in {ggplot2}

Clean and join data

BLS unemployment data

Census geographic data

Map adjustments

Extract interior state borders

Map with horizontal gradient step legend

Map with histogram legend

Map with automatic histogram legend with {legendry}

Bonus! Use points instead of choropleths

Bonus #2! Use a diverging color scheme + nested legend circles

Citation

How to move Crimea from Russia to Ukraine in maps with R

The Natural Earth Project

Natural Earth’s de facto policy

Natural Earth de jure points of view

Relocating Crimea manually with R and {sf}

Identifying the Crimea POLYGON from a POINT

Extracting the Crimea POLYGON from Russia

Adding the Crimea POLYGON to Ukraine

Updating Russia and Ukraine in the full data

The whole game

Moving Crimea with medium resolution (50m) data

Using the adjusted Natural Earth data as GeoJSON in Observable JS

Broken GeoJSON

GeoJSON and ↻ ↺ winding order ↻ ↺

Clean GeoJSON with correct winding order

Alternative data sources

Converting TikZ to both PDF and SVG with Quarto and `diagram`

Other projections and `.fitExtent()`

Arbitrary areas and `.fitExtent()`

Why not just use `results="asis"`?