我想计算我所在州(北卡罗来纳州(的邮政编码的人口密度。我能够从美国人口普查中提取邮政编码人口和多边形，并使用以下代码绘制北卡罗来纳州的地图：
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
这是按邮政编码绘制的人口地图。为了按邮政编码计算和绘制人口密度图，我需要每个邮政编码多边形的面积(以英里或平方公里为单位(。我正在努力寻找一种方法：(a(从美国人口普查局获得这些数据，(b(在其他地方找到，或者(c(使用多边形几何来计算。
如有任何建议，我们将不胜感激。
另一种方法是在get_acs()中设置keep_geo_vars = TRUE。这将返回每个ZCTA多边形中陆地和水域的面积(平方米(。为了计算人口密度，您可能更喜欢只使用土地面积，而不是每个ZCTA多边形的总面积。
陆地面积变量为ALAND10，水域面积为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...
我找到了一个非常简单的问题答案。请注意，描述每个邮政编码的多边形位于geo数据帧中的变量geometry中。
首先，我们计算面积，默认情况下，面积以m^2为单位返回。
library(sf)
geo$area.m2 <- st_area(geo$geometry)
其次，我们转换为英里^2个单位。
library(units)
geo$area.miles2 <- set_units(geo$area.m2, miles^2)