1

I have written the following code that downloads state boundaries and population data from the US Census and then recursively merges the least populated states with their least populated neighbors. The idea is to come up with 20 districts of more equal population than the 50 states. As implemented below, some polygons disappear. Anyone know how to fix that problem?

library(sf)
library(tidyverse)
library(tidycensus)

key <- "your key goes here (see tidycensus docs)"

census_api_key(key)

states <- get_estimates("state", "population", geometry = TRUE) %>%
  st_transform(2163) %>%
  spread(variable, value) %>%
  st_simplify(TRUE, 1e3) %>%
  st_buffer(1e3) %>%
  rename(id = GEOID, name = NAME, population = POP, density = DENSITY)

merge_states <- function(x) {
  i <- which.min(x[["population"]])
  j <- st_intersects(x[i, ], x[-i, ])[[1]]
  k <- j[which.min(x[["population"]][j])]
  if (length(k) < 1) {
    l <- seq.int(1, nrow(x))[-i]
    m <- which.min(st_distance(x[i, ], x[l, ]))
    k <- l[m]
  }
  x$geometry[k] <- st_union(x$geometry[k], x$geometry[i])
  x$population[k] <- x$population[k] + x$population[i]
  x <- x[-i, ]
  return(x)
}

zones <- states

while (nrow(zones) > 20) zones <- merge_states(zones)

#plot(zones[, "id"], key.pos = NULL)

ggplot(zones) + geom_sf(aes(fill = id)) + guides(fill = FALSE)

Output of code

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.