I am trying to calculate the area of each census tract in the United States. I used tidycensus to extract a multipolygon for each census tract with its coordinates.
library(tidycensus)
library(purrr)
census_api_key("KEY")
us <- unique(fips_codes$state)[1:51]
geometry <- map_df(us, function(x) {
get_acs(geography = "tract", variables = "B19013_001", geometry = TRUE,
state = x)
})
None of the methods I have been able to find for calculating the area of a multipolygon work with lists, which is the format tidycensus stores the coordinates in. Is there another format that would be easier to use? Is there some way of calculating the area using coordinates in this format?
I have tried gArea() from the rgeos package, area() from the raster package, st_area() from the sf package, and polyarea() from the geometry package.
Calculating area of a multipolygon using sf::st_area() is working fine for me. Here's an example of downloading data from tidycensus and calculating the area of all census tracts in Vermont.
library(tidycensus)
library(sf)
library(dplyr)
vt_tracts <- get_acs(
geography = "tract",
variables = "B19013_001",
geometry = TRUE,
state = "VT"
)
vt_tracts %>%
mutate(area = st_area(.))
#> Simple feature collection with 184 features and 6 fields
#> geometry type: POLYGON
#> dimension: XY
#> bbox: xmin: -73.43774 ymin: 42.72685 xmax: -71.46455 ymax: 45.01666
#> CRS: 4269
#> First 10 features:
#> GEOID NAME variable estimate
#> 1 50001960100 Census Tract 9601, Addison County, Vermont B19013_001 78750
#> 2 50001960200 Census Tract 9602, Addison County, Vermont B19013_001 74958
#> 3 50001960300 Census Tract 9603, Addison County, Vermont B19013_001 59464
#> 4 50001960400 Census Tract 9604, Addison County, Vermont B19013_001 78672
#> 5 50001960500 Census Tract 9605, Addison County, Vermont B19013_001 56066
#> 6 50001960600 Census Tract 9606, Addison County, Vermont B19013_001 58343
#> 7 50001960700 Census Tract 9607, Addison County, Vermont B19013_001 53910
#> 8 50001960800 Census Tract 9608, Addison County, Vermont B19013_001 66713
#> 9 50001960900 Census Tract 9609, Addison County, Vermont B19013_001 65515
#> 10 50001961000 Census Tract 9610, Addison County, Vermont B19013_001 62763
#> moe geometry area
#> 1 8764 POLYGON ((-73.1968 44.26663... 212702063 [m^2]
#> 2 8731 POLYGON ((-73.39963 44.1553... 158329687 [m^2]
#> 3 10068 POLYGON ((-73.27275 44.1768... 6614856 [m^2]
#> 4 4139 POLYGON ((-73.43774 44.0450... 359804506 [m^2]
#> 5 10897 POLYGON ((-73.1495 44.17669... 109378523 [m^2]
#> 6 3601 POLYGON ((-73.06466 44.0408... 530756434 [m^2]
#> 7 7404 POLYGON ((-73.16763 44.0289... 70741373 [m^2]
#> 8 10440 POLYGON ((-73.19061 44.0189... 30961721 [m^2]
#> 9 4949 POLYGON ((-73.41258 43.9827... 478027734 [m^2]
#> 10 7192 POLYGON ((-73.18767 43.8964... 134156214 [m^2]
Although to answer your more narrow question about finding the area of US Census tracts, you could use the tigris pacakge (which is what tidycensus uses under the hood to pull geometry), which in addition to the geometry, returns the land area (ALAND) and water area (AWATER) of each tract (or any other geography requested) directly from the Census files.
tigris::tracts(state = "VT", class = "sf")
#> Simple feature collection with 184 features and 12 fields
#> geometry type: MULTIPOLYGON
#> dimension: XY
#> bbox: xmin: -73.4379 ymin: 42.72693 xmax: -71.46505 ymax: 45.01666
#> CRS: 4269
#> First 10 features:
#> STATEFP COUNTYFP TRACTCE GEOID NAME NAMELSAD MTFCC FUNCSTAT
#> 1 50 007 000600 50007000600 6 Census Tract 6 G5020 S
#> 2 50 007 000800 50007000800 8 Census Tract 8 G5020 S
#> 3 50 007 000900 50007000900 9 Census Tract 9 G5020 S
#> 4 50 007 001100 50007001100 11 Census Tract 11 G5020 S
#> 5 50 001 960400 50001960400 9604 Census Tract 9604 G5020 S
#> 6 50 001 960900 50001960900 9609 Census Tract 9609 G5020 S
#> 7 50 001 960500 50001960500 9605 Census Tract 9605 G5020 S
#> 8 50 001 961000 50001961000 9610 Census Tract 9610 G5020 S
#> 9 50 001 960600 50001960600 9606 Census Tract 9606 G5020 S
#> 10 50 001 960300 50001960300 9603 Census Tract 9603 G5020 S
#> ALAND AWATER INTPTLAT INTPTLON geometry
#> 1 2240206 147061 +44.4833940 -073.1936217 MULTIPOLYGON (((-73.20609 4...
#> 2 1349586 0 +44.4609464 -073.2086195 MULTIPOLYGON (((-73.21524 4...
#> 3 602429 0 +44.4714041 -073.2088310 MULTIPOLYGON (((-73.21393 4...
#> 4 1736310 0 +44.4509355 -073.2197727 MULTIPOLYGON (((-73.23253 4...
#> 5 321773486 38567561 +44.0870615 -073.2671642 MULTIPOLYGON (((-73.43771 4...
#> 6 456817599 21511403 +43.9055842 -073.2904483 MULTIPOLYGON (((-73.41251 4...
#> 7 107525416 1714148 +44.0978095 -073.0485314 MULTIPOLYGON (((-73.1487 44...
#> 8 128778705 5371192 +43.9008582 -073.1035817 MULTIPOLYGON (((-73.18875 4...
#> 9 528821667 1880822 +43.9987500 -072.9398701 MULTIPOLYGON (((-73.06466 4...
#> 10 6413118 203405 +44.1632873 -073.2485362 MULTIPOLYGON (((-73.27275 4...