[R-sig-Geo] Problem with st_union

[Tim Keitt]

When I run the following code:

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)

st_is_valid(states)
st_is_simple(states)

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)) +
  geom_sf(data = states, fill = NA, linetype = "dotted") +
  guides(fill = FALSE)

I find that some of the polygons get dropped from the result. Is this an
issue with the GEOS union code or did I just set it up wrong?

Thanks.

THK

	[[alternative HTML version deleted]]