Search code examples
rgisr-sftmaptidycensus

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


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!


Solution

  • 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)