Aggregate ACS tract estimates to drive-time areas by area weighting

Combines the American Community Survey (ACS) estimates of the census tracts that overlap each drive-time area (isochrone) into an estimate and a margin of error for every site, drive time, and variable. Counts, such as total population (B01003_001), are added up with each tract weighted by the share of its area inside the drive-time area (its coverage weight), which assumes that whatever a variable counts is spread evenly over each tract's area. The assumption is made for each variable separately, and it is a stronger one for a subgroup, such as the people below the poverty level, than for the population as a whole. Medians and per-person values, such as median household income (B19013_001) and per capita income (B19301_001), are averaged over the overlapping tracts with weights proportional to the area each tract shares with the drive-time area. This average stands in for the median or per-person value of the drive-time area, but its weights follow area and ignore how many people live in each tract. It can therefore be far from that value when the overlapping tracts differ in population density.

Usage

cacs_intersect_weight(
  iso_sf,
  acs_sf,
  bg_pop_sf = NULL,
  weight_method = c("area", "population"),
  min_weight = 1e-06,
  verbose = TRUE,
  keep_tract_audit = FALSE,
  cache_dir = NULL
)

Arguments

iso_sf: An sf object of drive-time areas as returned by cacs_isochrone(), with one row for each site and drive time. Data from other sources must have the same columns, which cacs_validate_iso() checks, and each area must contain the shorter ones (ring_topology = "cumulative"; see cacs_rings_to_cumulative()).
acs_sf: An sf object of ACS estimates for census tracts as returned by cacs_acs_prefetch(), with one row for each tract and variable. Data from other sources must have the same columns, which cacs_acs_validate() checks.
bg_pop_sf: An sf object of block-group population estimates, or NULL (the default). It has no effect on the result (see weight_method).
weight_method: A string giving the weighting method: "area" (the default) or "population", which is not implemented yet; using it gives an error before any other work is done.
min_weight: A single number in [0, 1); the default is 1e-6. Tracts whose coverage weight is at or below it are dropped (see step 4 of Details).
verbose: A logical value. If TRUE (the default), progress messages are shown: a summary line when the step finishes and, with five or more site and drive-time pairs, also a line as each pair finishes. See the "Progress messages" section of cacs_run() for how to turn them off or show the lines for any number of pairs.
keep_tract_audit: A logical value, FALSE (the default) or TRUE. If TRUE, the result has a cacs_tract_audit attribute: a table with one row for each tract used in each drive-time area. Its columns are site_id, drive_time_min, GEOID, the tract's coverage weight (area_wt), and the areas of the overlap (int_area_m2) and of the tract (tract_area_m2) in square meters.
cache_dir: A path to the cache folder, or NULL (the default) to use cacs_cache_dir(). The result is saved in, and looked for in, its intersect subfolder. A folder outside the temporary folder of the R session is tidied as described in cacs_cache_dir().

Value

A tibble with one row for each site, drive time, and ACS variable, sorted by site_id, drive_time_min, and variable, followed by one row for each site and drive-time pair with no tract left (see Details). It has the columns that cacs_run() returns in its long form, described in its Value section, except est_total, var_total_raw, est_mean, and var_mean_raw, which cacs_propagate_moe() adds. The main columns are:

estimate, moe: The estimate and its margin of error, at the 90 percent level (see Details).
weight_sum, n_tracts: The sum of the coverage weights of the tracts combined, also on the rows for medians and per-person values, and the number of those tracts, including any with a missing estimate.
estimand_family, weight_basis: The kind of quantity and the weights used: "spatial_total" (a count) with "coverage"; "median_proxy" (a median) or "area_weighted_scalar_proxy" (a per-person value) with "area_mean"; or "metadata_only" (a code that is not combined) with "none", and cacs_propagate_moe() gives an error for such rows.
failure_origin: "none", or "intersection" on the row of a pair with no tract left.

The result has these attributes:

cacs_aggregation_carriers: A table of the weighted sums and averages of the tract estimates, with their variances, for each site, drive time, and variable. cacs_propagate_moe() and cacs_derive_rates() read it, and cacs_derive_rates() removes it.
cacs_aggregation_provenance: A list recording weight_method, min_weight, acs_year, the time the result was computed, the versions of catchmentACS, sf, GEOS, and PROJ, and counts of variables and of missing estimates. n_sites_input is the number of site and drive-time pairs, and n_sites_with_data and n_sites_empty are the numbers of rows of the result with and without tracts.
skipped_geoids: The GEOID of each row of acs_sf skipped in step 3 of Details, so a skipped tract appears once for each of its variables; an empty character vector when nothing is skipped.

It also has cacs_schema_version, the version label ("1.0") of the column layout, and, with keep_tract_audit = TRUE, cacs_tract_audit (see that argument).

Details

Like the ACS margins they are combined from, the margins of error are half-widths of 90 percent confidence intervals, and they treat the weights as fixed and the tract estimates as independent (see cacs_propagate_moe()). Rates, such as the poverty rate, are computed later by cacs_derive_rates(), each as the ratio of two of these weighted counts.

The function works in these steps.

Checks the inputs. iso_sf must be in EPSG:4326 (longitude and latitude on WGS 84) and acs_sf in EPSG:4269 (NAD83), as returned by cacs_isochrone() and cacs_acs_prefetch(), and two rows of acs_sf for the same tract and variable give an error. Both must lie within a box around the contiguous United States and the District of Columbia, so data for Alaska, Hawaii, or Puerto Rico give an error. The codes that the Census Bureau's data API puts in place of some estimates and margins of error (-222222222, -333333333, -555555555, -666666666, -888888888, and -999999999) are set to NA, with a warning that counts them; a margin-of-error code next to a missing estimate is not counted. Other negative values are used as they are. A warning is also given when the bounding box of the drive-time areas extends beyond that of the tracts (see below).
Transforms both inputs to EPSG:5070 (NAD83 / Conus Albers), an equal-area projection, and measures all areas there, in square meters. Invalid geometries are repaired, with a warning. For iso_sf, the change from WGS 84 to NAD83 is the one that PROJ chooses. Without datum grid files, PROJ treats the two as the same; with them, it can move the areas by a meter or two. Results can therefore differ between computers in the fourth or fifth significant digit.
Skips, with a warning, the tracts whose area is zero or not finite, such as a water tract with an empty boundary, and stops with an error if every tract is like that. The skipped tracts are listed in the skipped_geoids attribute.
For each site and drive time, finds the tracts that overlap the area. A tract's coverage weight is the area of its overlap divided by the tract's area. Tracts with a coverage weight at or below min_weight, including tracts that only touch the edge of the area, are dropped. Each area is used whole, so the 10-minute results include the tracts of the 5-minute area.
Combines the remaining tracts for each variable. A count is the sum of the tract estimates multiplied by their coverage weights, and a median or per-person value is the average of the tract estimates weighted by the area of each overlap. The margin of error is $\sqrt{\sum_i (w_i M_i)^2}$, where $M_i$ is the margin of error of tract $i$ and $w_i$ its weight. For a count the weight is the coverage weight, and for a median or per-person value it is the area share (see the "Coverage weights and area shares" section). A missing estimate in any of the tracts, including a code set to NA in step 1, makes the variable's estimate and margin of error NA; a missing margin of error makes only the margin of error NA. A rate that uses the variable is NA in both cases (see cacs_derive_rates()). cacs_acs_prefetch() sets negative margins of error to NA.

Whether a variable is a count, a median, or a per-person value is decided by its ACS code. Tables B19013 and B25077 are medians, and table B19301 is a per-person value. Every other code of the form B, five digits, an underscore, and three digits (such as B17001_002) is treated as a count, so a median or per-person value from another table is added up like a count. Other codes, such as those of tables whose names begin with C or S or end with a letter (such as B17001A), are not combined and get NA values.

Only the tracts in acs_sf are used: any part of a drive-time area outside them (for example, across a state line) adds nothing, so counts come out too small, and rates and medians come from the remaining tracts. The only warning about this, in step 1, compares the bounding box of all the areas in iso_sf with the bounding box of all the tracts in acs_sf: it is given when the first extends more than 0.05 degrees (about 5 km) beyond the second on some side. The tracts of a neighboring state can be included in acs_sf, for example by combining two cacs_acs_prefetch() results with dplyr::bind_rows(); the combined data no longer record the ACS year, so acs_year in the result is NA.

A site and drive time can have no tract left after step 4: no tract overlaps the area, the area is empty or appears twice in iso_sf, or every tract is at or below min_weight. Such a pair gives one row with variable = NA, NA values, n_tracts = 0, and failure_origin = "intersection". No warning is given. The same row is given when the overlap calculation for a pair fails, so one failing pair does not stop the others.

Results are cached (see cacs_set_cache()) in the folder given by cache_dir or cacs_cache_dir(): a call that matches an earlier one returns the saved result without repeating steps 2 to 5 or their warnings. The drive-time areas are matched by the cache_key that cacs_isochrone() attaches to its result together with their contents: the geometries, the coordinate reference system, and the columns site_id, drive_time_min, provider, profile, osm_snapshot_date, and ring_topology. A subset or an edited copy of a cacs_isochrone() result is therefore computed again, while the same areas in another row order use the saved result. options(catchmentACS.cache_intersect = FALSE) turns this cache off.

Coverage weights and area shares

For one site, drive time, and variable, let $a_i$ be the area of the overlap between tract $i$ and the drive-time area and $A_i$ the area of the tract. The tracts are those kept in step 4 of Details that have a row for the variable in acs_sf; a tract without such a row is left out, without a warning. A count then does not include that tract, and n_tracts on the rows of that variable is smaller than on the rows of the variables that have a row for the tract (cacs_acs_validate() describes a check). Counts are combined with the coverage weights (weight_basis = "coverage") $$c_i = \frac{a_i}{A_i},$$ and medians and per-person values with the area shares (weight_basis = "area_mean") $$s_i = \frac{a_i}{\sum_k a_k},$$ which sum to one.

Totals kept for later steps

The cacs_aggregation_carriers attribute has the columns site_id, drive_time_min, variable, estimand_family, weight_sum, and n_tracts, with the same values as in the result. It has four more: est_total, the sum $\sum_i c_i X_i$ of the tract estimates $X_i$; est_mean, the average $\sum_i s_i X_i$; and var_total_raw and var_mean_raw, their variances (the squares of their standard errors). The sum and the average are computed for every variable; estimate and moe come from the sum for counts and from the average for medians and per-person values.

Examples

# Turn the cache off while this example runs (see ?cacs_set_cache).
old <- options(catchmentACS.cache_enabled = FALSE)

# Example data bundled with the package: the drive-time areas are circles
# with a radius of 1 km per minute, and the ACS data are made up.
library(sf)
iso <- readRDS(system.file("extdata", "legacy_2025_isochrones.rds",
                           package = "catchmentACS"))
acs <- readRDS(system.file("extdata", "sample_alabama_subset.rds",
                           package = "catchmentACS"))
iso_07 <- iso[iso$site_id == "AL_SITE_07" & iso$drive_time_min == 10, ]
agg <- cacs_intersect_weight(iso_sf = iso_07, acs_sf = acs,
                             keep_tract_audit = TRUE, verbose = FALSE)

# The coverage weight and the area share of each tract in the area (the
# three tracts have the same area, so here the area shares are the
# coverage weights divided by their sum)
tracts <- attr(agg, "cacs_tract_audit")
tracts$area_share <- tracts$int_area_m2 / sum(tracts$int_area_m2)
tracts[, c("GEOID", "area_wt", "area_share")]
#> # A tibble: 3 × 3
#>   GEOID       area_wt area_share
#>   <chr>         <dbl>      <dbl>
#> 1 01073000401  0.0145     0.0141
#> 2 01075000401  1          0.972 
#> 3 01077000401  0.0145     0.0141

# A count, a median, and a per-person value, with the weights each one uses
agg[agg$variable %in% c("B01003_001", "B19013_001", "B19301_001"),
    c("variable", "estimate", "moe", "weight_basis", "n_tracts")]
#> # A tibble: 3 × 5
#>   variable   estimate    moe weight_basis n_tracts
#>   <chr>         <dbl>  <dbl> <chr>           <int>
#> 1 B01003_001    4643.   612. coverage            3
#> 2 B19013_001   68037. 16160. area_mean           3
#> 3 B19301_001   16725.  2356. area_mean           3

options(old)