Remove the outer boundary of the US using tidycensus, sf, and/or tmap packages

417 views Asked by At

I downloaded the polygons for the contiguous United States from the tidycensus package. However, due to mapping needs, I would like to remove the outline of the US from the polygon, leaving the internal state borders only. Here is my code:

library(tidyverse)
library(tidycensus)
library(tigris)
library(tmap)
library(sf)

`%notin%` <- Negate(`%in%`)

acs_vars <- c(pop_total = "B00001_001",
              age_total = "B01001_001")
us = get_acs(
 geography = "state",
  variables = acs_vars,
  geometry = TRUE,
  output = "wide",
  year= 2018
)

contig_48 = us %>% filter(NAME %notin% c("Alaska", "Hawaii", "Puerto Rico"))

tm_shape(contig_48) + tm_polygons(col = "white")

In other words, I would like to remove the result of running st_union on my contig_48 object. I would like to be left with contig_48 LESS the result of running st_union(contig_48).

I would greatly appreciate any help on this ! Thanks a lot!

1

There are 1 answers

0
dieghernan On BEST ANSWER

So here the strategy would be:

  1. Create a single POLYGON of the contiguos US and cast it to LINESTRING (in fact, this is just to create the country borders using contig_48).
  2. Create a small buffer of the border line to capture well the borders of the state (I did this in EPSG 3857 - Mercator with a buffer of 1 km. since NAD83 is in degrees)
  3. Convert the POLYGONS of the states to LINESTRING and substract the buffered borders.

Find here a reprex:

library(tidyverse)
library(tidycensus)
library(tigris)
library(tmap)
library(sf)
#> Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1; sf_use_s2() is TRUE


acs_vars <- c(pop_total = "B00001_001",
              age_total = "B01001_001")
us = get_acs(
  geography = "state",
  variables = acs_vars,
  geometry = TRUE,
  output = "wide",
  year= 2018
)

contig_48 = us %>% filter(!NAME %in% c("Alaska", "Hawaii", "Puerto Rico"))

# Get internal boundaries
library(dplyr)
whole_pol <- contig_48 %>% mutate(id="country") %>%
  group_by(id) %>% summarise(n=n()) %>%
  st_cast("MULTILINESTRING") %>%
  st_transform(3857) %>%
  # Small buffer to capture the lines
  st_buffer(1000) %>%
  st_transform(st_crs(contig_48))


qtm(whole_pol)

# Now cast to lines all the shapes
contig_48_lines <- st_cast(contig_48, "MULTILINESTRING")

# Difference

inner_boundaries <- st_difference(contig_48_lines, whole_pol)
#> Warning: attribute variables are assumed to be spatially constant throughout all
#> geometries

tm_shape(contig_48) +
  tm_polygons() +
  tm_shape(inner_boundaries) + tm_lines(col = "red")

Created on 2022-03-23 by the reprex package (v2.0.1)