Obstetrical care in Iowa

The third article in Kalen McCain’s excellent series on maternal healthcare in southeastern Iowa, Barriers aren’t limited to birthdays, features a map from numbered page 14 of ACCESS TO OBSTETRICAL CARE IN IOWA: A REPORT to the IOWA STATE LEGISLATURE – CALENDAR YEAR 2021 [PDF].

This got me wondering what it would take to reproduce this figure and perhaps supplement it with others.

First, I hacked together some Python code to parse out the data from the PDF. Based on Automate the Boring Stuff with Python, I’ve been using PyPDF2, but as noted in those very docs, it sounds like we should be using pypdf going forward. (Some nice project history here.)

Maybe there’s a better way — just go to OCR? — but I found I could identify rows of data and position within the row based on the small set of values for the OB Unit status.

You can see the results here: https://github.com/fadend/iowa_ob_unit_2023_report/blob/main/iowa_ob_unit_2023_report/generated/ob_unit.csv.

Then to make a plot, I ended up mostly relying on some lecture notes from Luke Tierney at University of Iowa. (I’m looking forward to looking at the other notes from his course.)

library(ggplot2)
library(sf)
library(tidyr)

# Copied from Tierney's notes
map_data_sf <- function(map, region = ".", exact = FALSE, crs = 4326, ...) {
    maps::map(map, region, exact = exact, plot = FALSE, fill = TRUE, ...) |>
        sf::st_as_sf(crs = crs) |>
        separate(ID, c("region", "subregion"), sep = ",") |>
        sf::st_transform(crs)
}

ob_unit <- read.csv(url("https://raw.githubusercontent.com/fadend/iowa_ob_unit_2023_report/refs/heads/main/iowa_ob_unit_2023_report/generated/ob_unit.csv"), stringsAsFactors=FALSE)

ob_unit[is.na(ob_unit)] <- 0

# "O'Brien" county is "obrien" in the map data.
ob_unit$subregion <- gsub("'", "", tolower(ob_unit$county))

iowa_sf <- map_data_sf("county", "iowa")
iowa_ob_unit_sf <- merge(iowa_sf, ob_unit, on="subregion")

ggplot(iowa_ob_unit_sf) + geom_sf(aes(fill=ob_unit_status)) + geom_sf_text(aes(label=county), size=2) + coord_sf() + theme_void() + scale_fill_manual(name="ob_unit_status", values=setNames(c("#495d7e", "#d0a35c", "#86d1bf"), c("Closed", "Never", "OB Unit"))) + guides(fill=guide_legend(title="OB Unit status"))

This gives us something pretty close to the original:

And we can easily make other plots based on the data, like the percent of births per county attended by an obgyn:

ggplot(iowa_ob_unit_sf) + geom_sf(aes(fill= obgyn_percent)) + geom_sf_text(aes(label=county), size=2) + coord_sf() + theme_void() + guides(fill=guide_legend(title="% births with obgyn"))