Setup

library(sf)

## Linking to GEOS 3.8.1, GDAL 3.2.1, PROJ 7.2.1

library(tidyverse)

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──

## ✓ ggplot2 3.3.5     ✓ purrr   0.3.4
## ✓ tibble  3.1.5     ✓ dplyr   1.0.7
## ✓ tidyr   1.1.4     ✓ stringr 1.4.0
## ✓ readr   2.0.2     ✓ forcats 0.5.1

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

library(tigris)

## To enable 
## caching of data, set `options(tigris_use_cache = TRUE)` in your R script or .Rprofile.

library(tidycensus)

## 
## Attaching package: 'tidycensus'

## The following object is masked from 'package:tigris':
## 
##     fips_codes

options(tigris_use_cache = TRUE)
options(tigris_class = "sf")

library(tmap)

## Registered S3 methods overwritten by 'stars':
##   method             from
##   st_bbox.SpatRaster sf  
##   st_crs.SpatRaster  sf

## To install your API key for use in future sessions, run this function with `install = TRUE`.

Tigris

Get the basic data using Tigris.

counties_us = counties(cb=TRUE,resolution = '5m')
tm_shape(counties_us) + tm_borders()

.

Clean up

There are problems because of the outlying areas, which are really part of the US.

Get rid of Alaska, Hawaii and all the islands.

counties_us = counties_us %>% 
  filter(!(STATEFP %in% c("02","15")) & 
             as.numeric(STATEFP) < 60)
tm_shape(counties_us) + tm_borders()

Fix the crs

Make the crs 4326 (WGS84) for compatibility.

counties_us = st_transform(counties_us,4326)
st_crs(counties_us)

## Coordinate Reference System:
##   User input: EPSG:4326 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     USAGE[
##         SCOPE["Horizontal component of 3D system."],
##         AREA["World."],
##         BBOX[-90,-180,90,180]],
##     ID["EPSG",4326]]

Save the file.

save(counties_us,file = "counties_us.Rdata")