1

I would like to calculate the population density for ZIP codes in my state (North Carolina). I am able to extract the ZIP code populations and polygons from the US Census and plot the map of North Carolina using the following code:

library(tidycensus)
geo <- get_acs(geography = 'zcta',        # Get zip code-level data
               state = 'NC',              # for NC
               year = 2019,               # for 2019
               table = 'B01003',          # from this US Census table
               geometry = TRUE) %>%       # and return the geospatial polygons for mapping
       dplyr::rename('population' = estimate, 'zipcode' = NAME) %>%
       select(-moe) %>%
       arrange(zipcode) 

p <- tm_shape(geo) +
  tm_polygons('population')
p

This maps population by ZIP code. In order to calculate and map the population density by ZIP code, I need the area (in miles or kilometers squared) of each ZIP code polygon. I am struggling to find a way of (a) getting this data from the US Census site, (b) finding it elsewhere, or (c) using the polygon geometry to calculate it.

Any suggestions will be appreciated.

Brant Inman
  • 53
  • 3

2 Answers2

1

Another approach is to set keep_geo_vars = TRUE in get_acs(). This will return the area (in square meters) of the land and water areas in each ZCTA polygon. For calculating population density, you may prefer to use just the land area rather than the total area of each ZCTA polygon.

The land area variables is ALAND10 and the water area is AWATER10

library(tidycensus)
get_acs(geography = 'zcta',     
        state = 'NC',           
        year = 2019,            
        table = 'B01003',       
        geometry = TRUE,        
        keep_geo_vars = TRUE)

#> Getting data from the 2015-2019 5-year ACS
#> Simple feature collection with 808 features and 9 fields
#> geometry type:  MULTIPOLYGON
#> dimension:      XY
#> bbox:           xmin: -84.32187 ymin: 33.84232 xmax: -75.46062 ymax: 36.58812
#> geographic CRS: NAD83
#> First 10 features:
#>    ZCTA5CE10     AFFGEOID10 GEOID   ALAND10 AWATER10        NAME   variable
#> 1      28906 8600000US28906 28906 864608629 28813485 ZCTA5 28906 B01003_001
#> 2      28721 8600000US28721 28721 285413675    41953 ZCTA5 28721 B01003_001
#> 3      28365 8600000US28365 28365 498948199  2852124 ZCTA5 28365 B01003_001
#> 4      27317 8600000US27317 27317 139042432  9345547 ZCTA5 27317 B01003_001
#> 5      27562 8600000US27562 27562 139182043 11466187 ZCTA5 27562 B01003_001
#> 6      28748 8600000US28748 28748 218992045        0 ZCTA5 28748 B01003_001
#> 7      28025 8600000US28025 28025 282005597   384667 ZCTA5 28025 B01003_001
#> 8      28441 8600000US28441 28441 331231481   282711 ZCTA5 28441 B01003_001
#> 9      27893 8600000US27893 27893 285314738  3173744 ZCTA5 27893 B01003_001
#> 10     28101 8600000US28101 28101   2319755   131290 ZCTA5 28101 B01003_001
#>    estimate  moe                       geometry
#> 1     19701  736 MULTIPOLYGON (((-84.31137 3...
#> 2     10401  668 MULTIPOLYGON (((-82.93706 3...
#> 3     15533 1054 MULTIPOLYGON (((-78.29569 3...
#> 4     16169  875 MULTIPOLYGON (((-79.87638 3...
#> 5      2149  431 MULTIPOLYGON (((-78.99166 3...
#> 6     12606 1020 MULTIPOLYGON (((-82.88801 3...
#> 7     54425 1778 MULTIPOLYGON (((-80.62793 3...
#> 8      3396  588 MULTIPOLYGON (((-78.6127 34...
#> 9     39531 1258 MULTIPOLYGON (((-78.06593 3...
#> 10      970  245 MULTIPOLYGON (((-81.09122 3...
mfherman
  • 392
  • 4
  • 16
0

I was able to find an answer to my question that is very simple. Note that the polygons describing each ZIP code are found in the variable geometry in the geo data frame.

First, we calculate the area, which by default is returned in m^2 units.

library(sf)
geo$area.m2 <- st_area(geo$geometry)

Second, we convert to miles^2 units.

library(units)
geo$area.miles2 <- set_units(geo$area.m2, miles^2)
Brant Inman
  • 53
  • 3