diff --git a/DESCRIPTION b/DESCRIPTION index f9aa3c6..c014060 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: tongfen Type: Package Title: Make Data Based on Different Geographies Comparable -Version: 0.3.7 +Version: 0.3.8 Authors@R: c( person("Jens", "von Bergmann", email = "jens@mountainmath.ca", role = c("aut", "cre"), comment = "creator and maintainer")) Description: Several functions to allow comparisons of data across different geographies, in particular for Canadian census data from different censuses. diff --git a/NAMESPACE b/NAMESPACE index 863178c..1ba1cc6 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -8,6 +8,7 @@ export(get_tongfen_ca_census_ct_from_da) export(get_tongfen_census_ct) export(get_tongfen_census_da) export(get_tongfen_correspondence_ca_census) +export(get_tongfen_correspondence_us_census) export(get_tongfen_us_census) export(meta_for_additive_variables) export(meta_for_ca_census_vectors) diff --git a/NEWS.md b/NEWS.md index 3b64630..64e2c7a 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,52 @@ +# tongfen v.0.3.8 +## Breaking changes +- `get_tongfen_ca_census` now honours its `base_geo`, `na.rm`, `tolerance`, `crs` and + `data_transform` arguments, all of which were silently ignored. Most visibly, the + documented default `base_geo = NULL` now returns data without geographic information, + where previously the geography of the first dataset was returned. Pass `base_geo` to + get an `sf` object back +- removed the `area_mismatch_cutoff` argument from `get_tongfen_ca_census` and + `get_tongfen_correspondence_ca_census`, it never had any effect. Use `check_tongfen_areas` + to inspect area mismatches, keeping in mind that geographies for different years are + simplified independently and differ in how water features are cut out +## Major changes +- correspondence tables are now built via a vectorised connected components pass instead of + a row-by-row union-find, which makes tongfen on large geographies dramatically faster + (dissemination blocks for a large province: minutes down to seconds) +- the "statcan" method no longer downloads census geometries it does not use +- dissolving geometries skips regions that don't need to be merged +- new `get_tongfen_correspondence_us_census` to get correspondence tables for US census + geographies without also fetching the data +- US census tract correspondence tables now reach back to the 1990 census (`dec1990`). The + Census Bureau has retired the 1990 API endpoint, so 1990 data itself has to be brought in + separately, for example from NHGIS, and combined via `tongfen_aggregate` +- US county subdivisions can now be matched across the 2010 and 2020 censuses, previously only + the 2000 and 2010 censuses were available +- US correspondence tables no longer chain regions together over slivers. The Census Bureau + relationship files list every geometric overlap, including boundaries that only shifted + slightly, and matching those up merged unrelated regions into one common geography. The new + `min_area_share` argument controls how much area two regions have to have in common to count + as related, default is `0.01`, and no region is ever dropped. This gives substantially finer + common geographies, for Rhode Island tracts across the 2010 and 2020 censuses 198 instead of + 60, for Vermont 151 instead of 26 +## Minor changes +- `get_tongfen_us_census` gained a `sumfile` argument, passed through to tidycensus, either a + single value for all censuses or a vector named by dataset. Without it 2020 data is read from + the PL 94-171 redistricting file, which carries almost no variables +- `get_tongfen_correspondence_ca_census` gained a `crs` argument for the spatial + intersections, default is `3347` (Statistics Canada Lambert) +- missing geographic identifiers no longer merge unrelated regions into one common geography +- fix crash when tongfen-ing census tracts across non-adjacent censuses +- fix `get_tongfen_census_ct`, `get_tongfen_census_da` and `get_tongfen_ca_census_ct_from_da` + erroring out when called with `geo_format=NA` +- fix 2020 US census tract identifiers getting stripped of their leading zeros when read from + the relationship file, which silently dropped most tracts out of the result. For Rhode Island + 246 of 250 tracts were affected +- US county subdivision data now errors out up front on censuses it can't be matched across, + instead of failing with "Did not find matching geographic identifiers" after downloading + the comparability file and the census data +- faster `check_tongfen_areas` and `aggregate_correspondences` + # tongfen v.0.3.7 ## Major changes - accommodate factors in proportional_reaggregate diff --git a/R/helpers.R b/R/helpers.R index a59078b..2d6e332 100644 --- a/R/helpers.R +++ b/R/helpers.R @@ -24,97 +24,116 @@ inner_join_tongfen_correspondence <- function(data,correspondence,link){ -get_tongfen_correspondence <- function(dd){ - hs <- names(dd)[!grepl("TongfenMethod",names(dd))] - index = 1 - ddd<- dd %>% - mutate(TongfenID=!!as.name(hs[index])) - - while (index% nrow > 0) { - ddd<- ddd %>% - mutate(TongfenID=coalesce(.data$TongfenID,paste0(index,"_",!!as.name(hs[index])))) - index <- index + 1 +# Connected components of the "rows sharing an identifier value" graph. +# +# Each row of `dd` is a node, two rows are adjacent if they agree on the value of +# at least one identifier column. Rather than materialising the (potentially huge) +# edge set, each row points at a parent row and we alternate between hooking every +# row onto the smallest parent in each of its identifier groups and pointer +# jumping to collapse the resulting chains. This computes connected components in +# a logarithmic number of fully vectorised passes and never allocates more than a +# handful of vectors of length `nrow(dd)`. +# +# `codes` is a list of integer vectors, one per identifier column, with 0 marking +# missing values (which must not link rows). +connected_components <- function(codes, n) { + comp <- seq_len(n) + repeat { + changed <- FALSE + # hook each row onto the smallest parent among the rows sharing one of its + # identifier values + for (cd in codes) { + keep <- cd > 0L + if (!any(keep)) next + k <- cd[keep] + current <- comp[keep] + o <- order(k, current, method = "radix") + ko <- k[o] + first <- c(TRUE, ko[-1L] != ko[-length(ko)]) + group_min <- integer(max(k)) + group_min[ko[first]] <- current[o][first] + new <- group_min[k] + if (any(new != current)) { + comp[keep] <- new + changed <- TRUE + } + } + # pointer jumping, halves the depth of every chain per pass + repeat { + jumped <- comp[comp] + if (all(jumped == comp)) break + comp <- jumped + changed <- TRUE + } + if (!changed) break } + comp +} - # Optimized connected components using union-find approach - # Build a mapping of all unique IDs to their component root - # This is much faster than repeated group_by operations - - # Create union-find parent mapping - all_ids <- unique(ddd$TongfenID) - parent <- setNames(all_ids, all_ids) - - # Find root with path compression - find_root <- function(x) { - if (parent[x] == x) return(x) - parent[x] <<- find_root(parent[x]) # Path compression - return(parent[x]) - } +# smallest value of `values` within each group of `comp`, returned as a lookup +# vector indexed by component +group_first_sorted <- function(comp, values, k) { + o <- order(comp, values, method = "radix") + co <- comp[o] + first <- c(TRUE, co[-1L] != co[-length(co)]) + out <- character(k) + out[co[first]] <- values[o][first] + out +} - # Union two components - union_ids <- function(x, y) { - root_x <- find_root(x) - root_y <- find_root(y) - if (root_x != root_y) { - # Always attach to the smaller ID (alphabetically) - if (root_x < root_y) { - parent[root_y] <<- root_x - } else { - parent[root_x] <<- root_y - } - } +get_tongfen_correspondence <- function(dd){ + hs <- names(dd)[!grepl("TongfenMethod",names(dd))] + n <- nrow(dd) + if (n == 0) { + return(dd %>% mutate(TongfenID=character(0),TongfenUID=character(0)) %>% ungroup()) } - # For each identifier column, union all IDs that share the same identifier value - for (nn in hs) { - id_groups <- ddd %>% - select(identifier = !!as.name(nn), "TongfenID") %>% - distinct() %>% - group_by(.data$identifier) %>% - summarise(ids = list(.data$TongfenID), .groups = "drop") - - for (i in seq_len(nrow(id_groups))) { - ids_in_group <- id_groups$ids[[i]] - if (length(ids_in_group) > 1) { - # Union all pairs in this group - for (j in 2:length(ids_in_group)) { - union_ids(ids_in_group[1], ids_in_group[j]) - } - } + values <- lapply(hs, function(nn) as.character(dd[[nn]])) + names(values) <- hs + + # dense integer codes per column, 0 for NA so that missing identifiers never + # link two rows together + codes <- lapply(values, function(v) { + match(v, c(NA_character_, sort(unique(v[!is.na(v)])))) - 1L + }) + + comp <- connected_components(codes, n) + k <- max(comp) + + # label each row by its first available identifier, later columns are prefixed + # by their position so that identifiers from different columns cannot collide + base <- values[[1]] + if (anyNA(base) && length(hs) > 1) { + for (i in seq(2, length(hs))) { + missing <- is.na(base) + if (!any(missing)) break + v <- values[[i]][missing] + base[missing] <- ifelse(is.na(v), NA_character_, paste0(i, "_", v)) } } - - # Apply the final mapping - find root for each ID - final_mapping <- setNames( - vapply(all_ids, find_root, character(1), USE.NAMES = FALSE), - all_ids - ) - - # Map all TongfenIDs to their roots - ddd <- ddd %>% - mutate(TongfenID = final_mapping[.data$TongfenID]) - - # Vectorized UID generation - # Pre-compute all the grouped values at once - uid_parts <- ddd %>% - group_by(.data$TongfenID) %>% - summarise( - across( - all_of(hs), - ~paste0(cur_column(), ":", paste0(sort(unique(.x)), collapse = ",")), - .names = "uid_{.col}" - ), - .groups = "drop" - ) %>% - mutate( - TongfenUID = do.call(paste, c(select(., starts_with("uid_")), sep = " ")) - ) %>% - select("TongfenID", "TongfenUID") - - # Join the UIDs back - ddd %>% - left_join(uid_parts, by = "TongfenID") %>% - ungroup() + if (anyNA(base)) base[is.na(base)] <- paste0("row_", which(is.na(base))) + + # TongfenID is the smallest row label in the component + dd$TongfenID <- group_first_sorted(comp, base, k)[comp] + + # TongfenUID enumerates all identifiers making up the component + uid_parts <- lapply(hs, function(nn) { + v <- values[[nn]] + ok <- !is.na(v) + cc <- comp[ok] + vv <- v[ok] + o <- order(cc, vv, method = "radix") + cc <- cc[o] + vv <- vv[o] + keep <- c(TRUE, cc[-1L] != cc[-length(cc)] | vv[-1L] != vv[-length(vv)]) + collapsed <- vapply(split(vv[keep], cc[keep]), paste0, character(1), collapse = ",") + out <- character(k) + out[as.integer(names(collapsed))] <- paste0(nn, ":", collapsed) + out + }) + dd$TongfenUID <- do.call(paste, c(uid_parts, list(sep = " ")))[comp] + + dd %>% ungroup() } assert <- function (expr, error) { @@ -122,18 +141,67 @@ assert <- function (expr, error) { } +# Dissolve the geometries of `data` by `grouping_var`, the geometric equivalent +# of `summarize()`. Groups holding a single geometry - the bulk of the groups +# when tongfen-ing fine geographies - are passed through directly instead of +# being sent through `st_union()`, which is where most of the time in dissolving +# large geographies goes. +summarize_geometry_by_group <- function(data,grouping_var){ + geo_column <- attr(data,"sf_column") + keys <- data %>% ungroup() %>% sf::st_drop_geometry() %>% select(all_of(grouping_var)) + key <- do.call(paste,c(unname(as.list(keys)),list(sep="\x1f"))) + u <- unique(key) + u <- u[order(u,method="radix")] + index <- match(key,u) + geometry <- data[[geo_column]] + + counts <- tabulate(index,nbins=length(u)) + # first row belonging to each group, groups are in sorted key order + first_row <- integer(length(u)) + first_row[rev(index)] <- rev(seq_along(index)) + + result <- vector("list",length(u)) + singles <- counts==1L + result[singles] <- geometry[first_row[singles]] + multi <- which(counts>1L) + if (length(multi)>0) { + rows <- split(seq_along(index)[index %in% multi],index[index %in% multi]) + result[as.integer(names(rows))] <- lapply(rows,function(i) + suppressMessages(sf::st_union(geometry[i]))[[1]]) + } + + out <- keys[first_row,,drop=FALSE] + out[[geo_column]] <- sf::st_sfc(result,crs=sf::st_crs(geometry)) %>% + sf::st_cast("MULTIPOLYGON") + sf::st_sf(out,sf_column_name=geo_column) +} + + +# row-wise `paste0(unique(...),collapse=", ")` over a set of columns. Values are +# drawn from a small vocabulary, so the collapse is only computed once per +# distinct combination rather than once per row. +collapse_unique_by_row <- function(data,columns){ + sep <- "\x1f" + key <- do.call(paste,c(unname(as.list(data[columns])),list(sep=sep))) + u <- unique(key) + collapsed <- vapply(strsplit(u,sep,fixed=TRUE), + function(x) paste0(unique(x),collapse=", "), + character(1)) + unname(collapsed[match(key,u)]) +} + aggregate_correspondences <- function(correspondences){ clean_correspondence_names <- function(correspondence) { correspondence %>% select(!matches("Tongfen") | matches("TongfenMethod")) } - # compute full correspondence - # order by length to speed up the process - lengths <- correspondences %>% lapply(nrow) %>% unlist %>% rank(ties.method = "first") + # compute full correspondence, smallest table first to keep intermediate + # join results as small as possible + index_order <- correspondences %>% lapply(nrow) %>% unlist() %>% order() - correspondence <- correspondences[[lengths[1]]] %>% + correspondence <- correspondences[[index_order[1]]] %>% clean_correspondence_names() - if (length(correspondences)>1) for (index in lengths[-1]) { + if (length(correspondences)>1) for (index in index_order[-1]) { c <- correspondences[[index]] %>% clean_correspondence_names() match_columns <- intersect(names(correspondence),names(c)) @@ -143,9 +211,9 @@ aggregate_correspondences <- function(correspondences){ } method_columns <- names(correspondence)[grepl("TongfenMethod",names(correspondence))] - correspondence$M <- apply(correspondence[,method_columns],1,function(d)paste0(unique(d),collapse = ", ")) - correspondence %>% select(-method_columns) %>% - rename(TongfenMethod=.data$M) + correspondence$M <- collapse_unique_by_row(correspondence,method_columns) + correspondence %>% select(-all_of(method_columns)) %>% + rename(TongfenMethod="M") } diff --git a/R/tongfen.R b/R/tongfen.R index 79a6e9c..f3ccac3 100644 --- a/R/tongfen.R +++ b/R/tongfen.R @@ -141,11 +141,7 @@ aggregate_data_with_meta <- function(data,meta,geo=FALSE,na.rm=TRUE,quiet=FALSE) meta <- meta %>% bind_rows(tibble(variable=base_variables,type="Base")) if ("sf" %in% class(data)) { - geo_column=attr(data,"sf_column") - data <- left_join(data %>% - select(c(geo_column,grouping_var)) %>% - summarize(!!geo_column:=suppressMessages(sf::st_union(!!as.name(geo_column))) %>% - sf::st_cast("MULTIPOLYGON")), + data <- left_join(summarize_geometry_by_group(data,grouping_var), data %>% sf::st_set_geometry(NULL) %>% summarize_at(meta$variable,sum,na.rm=na.rm), @@ -198,6 +194,8 @@ rename_with_meta <- function(data,meta,ds=NULL){ #' @param base_geo identifier for which data element to base the final geography on, #' uses the first data element if `NULL` (default), #' expects that `base_geo` is an element of `names(data)`. +#' @param na.rm logical, determines how NA values should be treated when aggregating variables, +#' default is `TRUE` #' @return aggregated dataset of class sf if base_geo is not NULL and data is of type sf or tibble otherwise. #' @export #' @@ -214,7 +212,7 @@ rename_with_meta <- function(data,meta,ds=NULL){ #' result <- tongfen_aggregate(list(geo1 %>% rename(GeoUIDCA06=GeoUID), #' geo2 %>% rename(GeoUIDCA16=GeoUID)),correspondence,meta) #'} -tongfen_aggregate <- function(data,correspondence,meta=NULL, base_geo = NULL){ +tongfen_aggregate <- function(data,correspondence,meta=NULL, base_geo = NULL, na.rm = TRUE){ data <- ensure_names(data) nn <- names(data) if (is.null(base_geo)) base_geo <- nn[1] @@ -240,14 +238,10 @@ tongfen_aggregate <- function(data,correspondence,meta=NULL, base_geo = NULL){ by=match_column) %>% group_by(.data$TongfenID,.data$TongfenUID) if (!is.null(meta)) { - d <- d %>% aggregate_data_with_meta(meta) + d <- d %>% aggregate_data_with_meta(meta,na.rm=na.rm) } else { if ("sf" %in% class(d)) { - geo_column=attr(d,"sf_column") - d <- d %>% summarize(!!geo_column:=suppressMessages(sf::st_union(!!as.name(geo_column))) %>% - sf::st_cast("MULTIPOLYGON"), - .groups="drop") - + d <- summarize_geometry_by_group(d,c("TongfenID","TongfenUID")) } else { d <- d %>% summarize(.groups="drop") } @@ -261,7 +255,7 @@ tongfen_aggregate <- function(data,correspondence,meta=NULL, base_geo = NULL){ for (ds in nn[nn!=base_geo]) { aggregated_data <- aggregated_data %>% inner_join(data_new[[ds]] %>% - select(-.data$TongfenUID) %>% + select(-"TongfenUID") %>% rename_with_meta(meta,ds), by="TongfenID") } @@ -351,7 +345,7 @@ proportional_reaggregate <- function(data,parent_data,geo_match,categories,base= na_base <- "...na_base" while (na_base %in% names(data)) { - na_base <- paste0(na_base) + na_base <- paste0("...",na_base) } @@ -615,6 +609,10 @@ estimate_tongfen_correspondence <- function(data, #' Sanity check for areas of estimated tongfen correspondence. This is useful if for example the total extent #' of geo1 and geo2 differ and there are regions at the edges with large difference in overlap. #' +#' The result is a diagnostic, not a pass/fail test. Geographies for different years are +#' simplified independently and differ in how water features are cut out, so a sizable area +#' mismatch does not by itself mean the regions were matched up incorrectly. +#' #' @param data alist of geogrpahic data of class sf #' @param correspondence Correspondence table with columns the unique geographic identifiers for each of the #' geographies and the TongfenID (and optionally TongfenUID and TongfenMethod) @@ -668,14 +666,15 @@ check_tongfen_areas <- function(data,correspondence) { purrr::reduce(full_join,by="TongfenID") method_columns <- names(summary_data)[grepl("TongfenMethod",names(summary_data))] - summary_data$M <- apply(summary_data[,method_columns],1,function(d)paste0(unique(d),collapse = ", ")) + summary_data$M <- collapse_unique_by_row(summary_data,method_columns) + summary_data <- summary_data %>% + select(-all_of(method_columns)) %>% + rename(TongfenMethod="M") + + area_columns <- names(summary_data)[grepl("^area_",names(summary_data))] + areas <- lapply(area_columns,function(cc) as.numeric(summary_data[[cc]])) summary_data %>% - select(-method_columns) %>% - rename(TongfenMethod=.data$M) %>% - mutate(maxa=apply(select(.,matches("^area_")),1,max), - mina=apply(select(.,matches("^area_")),1,min)) %>% - mutate(max_log_ratio=log(.data$maxa/.data$mina)) %>% - select(-.data$maxa,-.data$mina) + mutate(max_log_ratio=log(do.call(pmax,areas)/do.call(pmin,areas))) } diff --git a/R/tongfen_ca.R b/R/tongfen_ca.R index d212812..def0d44 100644 --- a/R/tongfen_ca.R +++ b/R/tongfen_ca.R @@ -198,7 +198,7 @@ get_single_correspondence_ca_census_for <- function(year,level=c("DA","DB"),refr tmp=tempfile() utils::download.file(url,tmp) exdir=file.path(tempdir(),paste0("correspondence_",year,"_",level)) - if (dir.exists(exdir)) file.remove(exdir,recursive=TRUE) + if (dir.exists(exdir)) unlink(exdir,recursive=TRUE) dir.create(exdir,showWarnings = FALSE) utils::unzip(tmp,exdir=exdir) file=dir(exdir,"\\.txt|\\.csv") @@ -250,10 +250,12 @@ get_single_correspondence_ca_census_for <- function(year,level=c("DA","DB"),refr #' this method only works for "DB", "DA" and "CT" levels. #' * "estimate" uses `estimate_tongfen_correspondence` to build up the common geography from scratch based on geographies. #' * "identifier" assumes regions with identical geographic identifier are identical, and builds up the the correspondence for regions with unmatched geographic identifiers. -#' @param tolerance tolerance for `estimate_tongen_correspondence` in metres, default value is 50 metres. -#' @param area_mismatch_cutoff discard areas returned by `estimate_tongfen_correspondence` with area mismatch (log ratio) greater than cutoff. +#' @param tolerance tolerance for `estimate_tongen_correspondence` in metres, default value is 50 metres, +#' only used when method is 'estimate' or 'identifier' #' @param quiet suppress download progress output, default is `FALSE` #' @param refresh optional character, refresh data cache for this call, (default `FALSE`) +#' @param crs CRS to use for the spatial intersections if method is 'identifier' or +#' 'estimate', default is `3347` (Statistics Canada Lambert) #' @return dataframe with the multi-census correspondence file #' @export #' @@ -264,8 +266,8 @@ get_single_correspondence_ca_census_for <- function(year,level=c("DA","DB"),refr #' regions=list(CMA="59933"),level='CT') #'} get_tongfen_correspondence_ca_census <- function(geo_datasets, regions, level="CT", method="statcan", - tolerance = 50, area_mismatch_cutoff = 0.1, - quiet = FALSE, refresh = FALSE) { + tolerance = 50, + quiet = FALSE, refresh = FALSE, crs = 3347) { geo_datasets <- normalize_datasets(geo_datasets) if (method=="statcan") { @@ -282,12 +284,21 @@ get_tongfen_correspondence_ca_census <- function(geo_datasets, regions, level="C use_cache <- !refresh - data <- lapply(geo_datasets,function(g_ds){ - cancensus::get_census(dataset=g_ds, regions=regions, level=level, geo_format='sf', - labels="short", quiet=quiet, use_cache = use_cache) %>% - mutate(!!paste0("GeoUID",g_ds):=.data$GeoUID) - }) %>% - setNames(geo_datasets) + # the "statcan" method only ever looks at geographic identifiers, only the + # geometry based methods need to download the (potentially very large) geometries + geo_format <- if (method=="statcan") NA else 'sf' + + if (method=="statcan" && level=="CT") { + # correspondence is built from DA level links below, the CT level data is not needed + data <- list() + } else { + data <- lapply(geo_datasets,function(g_ds){ + cancensus::get_census(dataset=g_ds, regions=regions, level=level, geo_format=geo_format, + labels="short", quiet=quiet, use_cache = use_cache) %>% + mutate(!!paste0("GeoUID",g_ds):=.data$GeoUID) + }) %>% + setNames(geo_datasets) + } if (method=="statcan") { statcan_level <- level @@ -296,10 +307,12 @@ get_tongfen_correspondence_ca_census <- function(geo_datasets, regions, level="C years<-as.integer(geo_years) all_geo_years=seq(min(years),max(years),5) all_geo_datasets <- geo_dataset_for_years(all_geo_years) - for (g_ds in setdiff(all_geo_datasets,geo_datasets)) { - data[[g_ds]] <- cancensus::get_census(dataset=g_ds, regions=regions, level=level, geo_format='sf', - labels="short", quiet=quiet, use_cache = use_cache) %>% - mutate(!!paste0("GeoUID",g_ds):=.data$GeoUID) + if (level!="CT") { + for (g_ds in setdiff(all_geo_datasets,geo_datasets)) { + data[[g_ds]] <- cancensus::get_census(dataset=g_ds, regions=regions, level=level, geo_format=geo_format, + labels="short", quiet=quiet, use_cache = use_cache) %>% + mutate(!!paste0("GeoUID",g_ds):=.data$GeoUID) + } } prefix=paste0(statcan_level,"UID") @@ -321,7 +334,7 @@ get_tongfen_correspondence_ca_census <- function(geo_datasets, regions, level="C base_column <- paste0(prefix,year) match_column <- paste0("GeoUID",ds) data[[ds]] %>% - st_set_geometry(NULL) %>% + sf::st_drop_geometry() %>% select_at(match_column) %>% mutate(!!base_column:=!!as.name(match_column)) }) %>% @@ -340,7 +353,7 @@ get_tongfen_correspondence_ca_census <- function(geo_datasets, regions, level="C ds2 <- all_geo_datasets[all_geo_years==previous_year] if (!is.null(ds1) && length(ds1)>0) { match_column <- intersect(names(c),names(c_links[[ds1]])) - if (!is.null(match_column)) { + if (length(match_column)>0) { c <- c %>% inner_join(c_links[[ds1]],by=match_column) %>% select(-all_of(match_column)) %>% @@ -349,7 +362,7 @@ get_tongfen_correspondence_ca_census <- function(geo_datasets, regions, level="C } if (!is.null(ds2) && length(ds2)>0) { match_column <- intersect(names(c),names(c_links[[ds2]])) - if (!is.null(match_column)) { + if (length(match_column)>0) { c <- c %>% inner_join(c_links[[ds2]],by=match_column) %>% select(-all_of(match_column)) %>% @@ -369,8 +382,8 @@ get_tongfen_correspondence_ca_census <- function(geo_datasets, regions, level="C correspondence <- estimate_tongfen_correspondence(data, geo_identifiers, method = method, - tolerance=200, - computation_crs=3347) + tolerance=tolerance, + computation_crs=crs) } correspondence @@ -394,16 +407,15 @@ get_tongfen_correspondence_ca_census <- function(geo_datasets, regions, level="C #' * "estimate" uses `estimate_tongfen_correspondence` to build up the common geography from scratch based on geographies. #' * "identifier" assumes regions with identical geographic identifier are identical, and builds up the the correspondence for regions with unmatched geographic identifiers. #' @param base_geo base census year to build up common geography from, `NULL` (the default) to not return -#' any geographi data -#' @param na.rm logical, determines how NA values should be treated when aggregating variables +#' any geographic data +#' @param na.rm logical, determines how NA values should be treated when aggregating variables, +#' default is `FALSE` #' @param tolerance tolerance for `estimate_tongen_correspondence` in metres, default value is 50 metres, #' only used when method is 'estimate' or 'identifier' -#' @param area_mismatch_cutoff discard areas returned by `estimate_tongfen_correspondence` with area mismatch (log ratio) greater than cutoff, -#' only used when method is 'estimate' or 'identifier' #' @param quiet suppress download progress output, default is `FALSE` #' @param refresh optional character, refresh data cache for this call, (default `FALSE`) #' @param crs optional CRS to transform data to, and use for spatial intersections if method is -#' 'identifier' or 'estimate' +#' 'identifier' or 'estimate', defaults to `3347` (Statistics Canada Lambert) for the intersections #' @param data_transform optional transform function to be applied to census data after being returned from cancensus #' @return dataframe with variables on common geography #' @export @@ -423,7 +435,6 @@ get_tongfen_correspondence_ca_census <- function(geo_datasets, regions, level="C get_tongfen_ca_census <- function(regions,meta,level="CT",method="statcan", base_geo=NULL,na.rm=FALSE, tolerance = 50, - area_mismatch_cutoff = 0.1, quiet=FALSE, refresh=FALSE, crs=NULL, @@ -432,6 +443,13 @@ get_tongfen_ca_census <- function(regions,meta,level="CT",method="statcan", geo_datasets <- meta$geo_dataset %>% unique() %>% sort() + if (!is.null(base_geo)) { + base_geo <- normalize_datasets(base_geo) + assert(length(base_geo)==1,"base_geo has to be a single dataset") + assert(base_geo %in% geo_datasets, + paste0("base_geo has to be one of the datasets ",paste0(geo_datasets,collapse=", "))) + } + meta <- meta %>% add_census_ca_base_variables() data <- lapply(geo_datasets,function(g_ds){ @@ -440,13 +458,15 @@ get_tongfen_ca_census <- function(regions,meta,level="CT",method="statcan", .data$type != "Base") %>% pull(.data$variable) %>% as.character() + # only the base geography is returned, no need to download geometries for the others + geo_format <- if (!is.null(base_geo) && g_ds==base_geo) 'sf' else NA c <- cancensus::get_census(dataset=g_ds, regions=regions, vectors=vectors, - level=level, geo_format='sf', + level=level, geo_format=geo_format, labels="short", quiet=quiet, use_cache = use_cache) %>% mutate(!!paste0("GeoUID",g_ds):=.data$GeoUID) - if (!is.null(crs)) c <- c %>% sf::st_transform(crs) - c + if (!is.null(crs) && !is.null(base_geo) && g_ds==base_geo) c <- c %>% sf::st_transform(crs) + c %>% data_transform() }) %>% setNames(geo_datasets) @@ -460,10 +480,11 @@ get_tongfen_ca_census <- function(regions,meta,level="CT",method="statcan", level = level, method = method, tolerance = tolerance, - area_mismatch_cutoff = area_mismatch_cutoff, quiet = quiet, - refresh = refresh) - aggregated_data <- tongfen_aggregate(data,correspondence,meta) + refresh = refresh, + crs = crs %||% 3347) + aggregated_data <- tongfen_aggregate(data,correspondence,meta, + base_geo=base_geo,na.rm=na.rm) } aggregated_data %>% rename_with_meta(meta) diff --git a/R/tongfen_ca_deprecated.R b/R/tongfen_ca_deprecated.R index 5d5e77f..40c4d75 100644 --- a/R/tongfen_ca_deprecated.R +++ b/R/tongfen_ca_deprecated.R @@ -23,7 +23,7 @@ get_tongfen_census_ct <- function(regions, lifecycle::deprecate_warn("0.2.0", "get_tongfen_census_ct()", "get_tongfen_ca_census()") #warning("This method is deprecated, use `get_tongfen_ca_census(regions,vectors,level='CT', method = 'identifier', tolerance = 500, ...)` instead") meta <- meta_for_ca_census_vectors(vectors) - base_geo <- ifelse(is.na(geo_format),NULL,meta$geo_dataset %>% unique %>% sort %>% first ) + base_geo <- if (is.na(geo_format)) NULL else meta$geo_dataset %>% unique() %>% sort() %>% first() get_tongfen_ca_census(regions = regions, meta = meta, @@ -55,7 +55,7 @@ get_tongfen_census_ct <- function(regions, get_tongfen_census_da <- function(regions,vectors,geo_format=NA,use_cache=TRUE,na.rm=TRUE,quiet=TRUE) { lifecycle::deprecate_warn("0.2.0", "get_tongfen_census_da()", "get_tongfen_ca_census()") meta <- meta_for_ca_census_vectors(vectors) - base_geo <- ifelse(is.na(geo_format),NULL,meta$geo_dataset %>% unique %>% sort %>% first ) + base_geo <- if (is.na(geo_format)) NULL else meta$geo_dataset %>% unique() %>% sort() %>% first() get_tongfen_ca_census(regions = regions, meta = meta, @@ -87,7 +87,7 @@ get_tongfen_census_da <- function(regions,vectors,geo_format=NA,use_cache=TRUE,n get_tongfen_ca_census_ct_from_da <- function(regions,vectors,geo_format=NA,use_cache=TRUE,na.rm=TRUE,quiet=TRUE) { lifecycle::deprecate_warn("0.2.0", "get_tongfen_census_da()", "get_tongfen_census_ca()") meta <- meta_for_ca_census_vectors(vectors) - base_geo <- ifelse(is.na(geo_format),NULL,meta$geo_dataset %>% unique %>% sort %>% first ) + base_geo <- if (is.na(geo_format)) NULL else meta$geo_dataset %>% unique() %>% sort() %>% first() get_tongfen_ca_census(regions = regions, meta = meta, diff --git a/R/tongfen_ca_estimate.R b/R/tongfen_ca_estimate.R index 90cd31c..badd6b5 100644 --- a/R/tongfen_ca_estimate.R +++ b/R/tongfen_ca_estimate.R @@ -67,7 +67,8 @@ tongfen_estimate_ca_census <- function(geometry, meta, level, # So maybe a function like get_tongfen_correspondence_from_seed census_data <- get_tongfen_ca_census(regions = regions, meta = meta, - level = level, na.rm = na.rm, quiet = quiet) %>% + level = level, base_geo = datasets, + na.rm = na.rm, quiet = quiet) %>% sf::st_transform(st_crs(geometry)) if (!is.null(downsample_level)){ diff --git a/R/tongfen_estimate.R b/R/tongfen_estimate.R index fa5842b..eb6f749 100644 --- a/R/tongfen_estimate.R +++ b/R/tongfen_estimate.R @@ -56,38 +56,26 @@ tongfen_estimate <- function(target,source,meta,na.rm=FALSE) { gc = which(st_is(i, "GEOMETRYCOLLECTION")) i[gc] = st_collection_extract(i[gc], "POLYGON") - two_d = which(st_dimension(i) == 2) - i[two_d] = st_cast(i[two_d], "MULTIPOLYGON") source <- source %>% rename(!!!safe_rename_vars) + source_area <- unclass(st_area(source)) x_st <- source[idx[,1],, drop=FALSE] %>% select(names(safe_rename_vars)) %>% pre_scale(meta,meta_var = "var_name") %>% mutate(...area_st = st_area(i) %>% unclass, - ...area_s = unclass(st_area(.))) %>% + ...area_s = source_area[idx[,1]]) %>% mutate(...factor = .data$...area_st/.data$...area_s) %>% mutate(...partial = .data$...factor < 0.99) %>% st_drop_geometry() - - for (var in meta$var_name) { - # for (ci in naive_CI) { - # c <- ci/100 - # x_st <- x_st %>% - # mutate(!!paste0(var,"_lower_",ci) := !!as.name(var) * ifelse(.data$...partial,(1-c) * .data$...factor, 1), - # !!paste0(var,"_upper_",ci) := !!as.name(var) * ifelse(.data$...partial,(1-c) * .data$...factor + c, 1)) - # - # } - x_st <- x_st %>% - mutate(!!var:=!!as.name(var) * .data$...factor) - } + x_st[meta$var_name] <- lapply(x_st[meta$var_name], `*`, x_st$...factor) x_st <- stats::aggregate(x_st, list(idx[,2]), sum, na.rm=na.rm) result <- target %>% left_join(x_st %>% - select(-.data$...factor,-.data$...partial,-.data$...area_s,-.data$...area_st) %>% + select(-all_of(c("...factor", "...partial", "...area_s", "...area_st"))) %>% rename(!!unique_key:="Group.1"), by=unique_key) %>% select(-all_of(unique_key)) %>% diff --git a/R/tongfen_us.R b/R/tongfen_us.R index 4b27553..a54673a 100644 --- a/R/tongfen_us.R +++ b/R/tongfen_us.R @@ -1,31 +1,63 @@ fips_code_for_state <- function(s){ tidycensus::fips_codes %>% filter(.data$state==s | .data$state_code==s) %>% - select(.data$state,.data$state_code) %>% + select("state","state_code") %>% unique() } +# census vintages the correspondence layer can bridge, in chronological order, together +# with the name of the GEOID column identifying geographies of that vintage +us_geoid_columns <- c(dec1990 = "GEOID90", + dec2000 = "GEOID00", + dec2010 = "GEOID10", + dec2020 = "GEOID20") + +# Relationship files list every geometric overlap between two censuses, including slivers +# along boundaries that only shifted slightly. Chaining those merges unrelated regions into +# one common geography, so they get cut. The share is taken over the larger of the two sides, +# a region carved out of a bigger one is only a small share of the old one but most of the +# new one. Regions never get dropped: if all parts of a region are slivers its largest part +# is kept, so every region still ends up in some common geography. +area_share <- function(part, total) ifelse(total > 0, part/total, 0) + +cut_correspondence_slivers <- function(d, share, min_area_share){ + keep <- share >= min_area_share + for (column in names(d)) { + g <- d[[column]] + largest <- share == unname(tapply(share, g, max)[g]) + keep <- keep | (!(g %in% g[keep]) & largest) + } + d[keep,,drop=FALSE] %>% unique() +} + +# `year` is the later of the two censuses the relationship file links get_us_ct_correspondence_path <- function(state,year){ - if (year=="2010") { - states <- fips_code_for_state(state) - if (nrow(states)!= 1) { - stop(paste0("Could not determine state: ",state)) - } + states <- fips_code_for_state(state) + if (nrow(states)!= 1) { + stop(paste0("Could not determine state: ",state)) + } + if (year=="2000") { + path <- paste0("https://www2.census.gov/geo/relfiles/tract/", + tolower(states$state),"/", + tolower(states$state), + states$state_code,"pop.txt") + } else if (year=="2010") { path <- paste0("https://www2.census.gov/geo/docs/maps-data/data/rel/trf_txt/", tolower(states$state), states$state_code,"trf.txt") } else if (year=="2020") { - states <- fips_code_for_state(state) - if (nrow(states)!= 1) { - stop(paste0("Could not determine state: ",state)) - } path <- paste0("https://www2.census.gov/geo/docs/maps-data/data/rel2020/t10t20/TAB2010_TAB2020_ST", states$state_code,".zip") + } else { + stop(paste0("No census tract relationship file available for ",year)) } + path } -get_us_ct_correspondence_2020 <- function(state,cache_path=getOption("tongfen.cache_path")) { - states <- fips_code_for_state(state) +# the 2010 to 2020 tract relationship file is block based, tract areas get summed up from +# the blocks they are made up of +get_us_ct_correspondence_2020 <- function(state,min_area_share=0.01, + cache_path=getOption("tongfen.cache_path")) { cache_path = file.path(cache_path %||% tempdir(),"us_data") path <- get_us_ct_correspondence_path(state,2020) @@ -34,14 +66,35 @@ get_us_ct_correspondence_2020 <- function(state,cache_path=getOption("tongfen.ca if (!dir.exists(cache_path)) dir.create(cache_path) utils::download.file(path,local_path,quiet = TRUE) } - readr::read_delim(local_path,delim="|", col_types = "cccccnncccnncnn") %>% + blocks <- readr::read_delim(local_path,delim="|",progress=FALSE, + col_types=readr::cols_only( + STATE_2010="c",COUNTY_2010="c",TRACT_2010="c",BLK_2010="c", + AREALAND_2010="d",AREAWATER_2010="d", + STATE_2020="c",COUNTY_2020="c",TRACT_2020="c",BLK_2020="c", + AREALAND_2020="d",AREAWATER_2020="d", + AREALAND_INT="d",AREAWATER_INT="d")) %>% mutate(GEOID10=paste0(.data$STATE_2010,.data$COUNTY_2010,.data$TRACT_2010), - GEOID20=paste0(.data$STATE_2020,.data$COUNTY_2020,.data$TRACT_2020)) %>% - select(.data$GEOID10,.data$GEOID20)%>% - unique + GEOID20=paste0(.data$STATE_2020,.data$COUNTY_2020,.data$TRACT_2020), + area10=.data$AREALAND_2010+.data$AREAWATER_2010, + area20=.data$AREALAND_2020+.data$AREAWATER_2020, + area_part=.data$AREALAND_INT+.data$AREAWATER_INT) + tracts10 <- blocks %>% select("GEOID10","BLK_2010","area10") %>% unique() %>% + group_by(.data$GEOID10) %>% summarize(area10=sum(.data$area10),.groups="drop") + tracts20 <- blocks %>% select("GEOID20","BLK_2020","area20") %>% unique() %>% + group_by(.data$GEOID20) %>% summarize(area20=sum(.data$area20),.groups="drop") + d <- blocks %>% + group_by(.data$GEOID10,.data$GEOID20) %>% + summarize(area_part=sum(.data$area_part),.groups="drop") %>% + left_join(tracts10,by="GEOID10") %>% + left_join(tracts20,by="GEOID20") + cut_correspondence_slivers(d %>% select("GEOID10","GEOID20"), + pmax(area_share(d$area_part,d$area10), + area_share(d$area_part,d$area20)), + min_area_share) } -get_us_ct_correspondence_2010 <- function(state,cache_path=getOption("tongfen.cache_path")){ +get_us_ct_correspondence_2010 <- function(state,min_area_share=0.01, + cache_path=getOption("tongfen.cache_path")){ path <- get_us_ct_correspondence_path(state,"2010") file <- basename(path) cache_path = file.path(cache_path %||% tempdir(),"us_data") @@ -50,7 +103,7 @@ get_us_ct_correspondence_2010 <- function(state,cache_path=getOption("tongfen.ca if (!dir.exists(cache_path)) dir.create(cache_path) utils::download.file(path,local_path,quiet=TRUE) } - d<-readr::read_csv(local_path, + d<-readr::read_csv(local_path,progress=FALSE, col_names=c("STATE00","COUNTY00","TRACT00","GEOID00", "POP00","HU00","PART00","AREA00","AREALAND00", "STATE10","COUNTY10","TRACT10","GEOID10", @@ -58,32 +111,94 @@ get_us_ct_correspondence_2010 <- function(state,cache_path=getOption("tongfen.ca "AREAPT","AREALANDPT","AREAPCT00PT", "AREALANDPCT00PT","AREAPCT10PT","AREALANDPCT10PT", "POP10PT","POPPCT00","POPPCT10","HU10PT","HUPCT00","HUPCT10"), - col_types = "cccciiccccccciicccnnnnnnnnnnnn") + col_types = "cccciiccccccciicccnnnnnnnnnnnn") %>% + group_by(.data$GEOID00,.data$GEOID10) %>% + summarize(area_part=sum(.data$AREAPT), + area00=max(as.numeric(.data$AREA00)), + area10=max(as.numeric(.data$AREA10)), + .groups="drop") + cut_correspondence_slivers(d %>% select("GEOID00","GEOID10"), + pmax(area_share(d$area_part,d$area00), + area_share(d$area_part,d$area10)), + min_area_share) } -get_us_ct_correspondence <- function(state, datasets, - cache_path=getOption("tongfen.cache_path")){ - c <- NULL - if (setdiff(datasets,c("dec2000","dec2010","dec2020")) %>% length() > 0) { - stop("Invalid census years, can only match censuses 2000 through 2020") +# the 1990 to 2000 relationship files are fixed width, the "pop" variant is the complete one, +# listing every tract rather than only the ones that changed. It only carries the land area of +# each part, tract areas get summed up from those +get_us_ct_correspondence_2000 <- function(state,min_area_share=0.01, + cache_path=getOption("tongfen.cache_path")){ + path <- get_us_ct_correspondence_path(state,"2000") + cache_path = file.path(cache_path %||% tempdir(),"us_data") + local_path <- file.path(cache_path,basename(path)) + if (!file.exists(local_path)) { + if (!dir.exists(cache_path)) dir.create(cache_path) + utils::download.file(path,local_path,quiet=TRUE) } - if ("dec2000" %in% datasets) { - c<-get_us_ct_correspondence_2010(state,cache_path=cache_path) %>% - select(matches("^GEOID\\d{2}$")) - if ("dec2020" %in% datasets) { - c2 <- get_us_ct_correspondence_2020(state,cache_path=cache_path) %>% - select(matches("^GEOID\\d{2}$")) - c <- full_join(c,c2,by="GEOID10") - } - if (!("dec2010" %in% datasets)) c <- c %>% select(-.data$GEOIOD10) - c <- c %>% unique - } else if ("dec2020" %in% datasets) { - c<-get_us_ct_correspondence_2020(state,cache_path=cache_path) %>% - select(matches("^GEOID\\d{2}$")) - } else stop("Invalid census years, can only match censuses 2000 through 2020") - c + d <- readr::read_fwf(local_path, + readr::fwf_cols(STATE90=c(1,2),COUNTY90=c(3,5),TRACT90BASE=c(6,9), + TRACT90SUF=c(10,11),PART90=c(12,12),POP90TRACT=c(13,21), + PCT90=c(22,25),STATE00=c(26,27),COUNTY00=c(28,30), + TRACT00BASE=c(31,34),TRACT00SUF=c(35,36),PART00=c(37,37), + POP00TRACT=c(38,46),PCT00=c(47,50),POPPART=c(51,59), + AREALAND=c(60,73),STAB=c(74,75),COUNTYNAME=c(76,135)), + col_types=readr::cols(.default="c"),progress=FALSE) %>% + mutate(GEOID90=paste0(.data$STATE90,.data$COUNTY90,.data$TRACT90BASE, + coalesce(.data$TRACT90SUF,"00")), + GEOID00=paste0(.data$STATE00,.data$COUNTY00,.data$TRACT00BASE, + coalesce(.data$TRACT00SUF,"00"))) %>% + group_by(.data$GEOID90,.data$GEOID00) %>% + summarize(area_part=sum(as.numeric(.data$AREALAND)),.groups="drop") %>% + group_by(.data$GEOID90) %>% mutate(area90=sum(.data$area_part)) %>% + group_by(.data$GEOID00) %>% mutate(area00=sum(.data$area_part)) %>% + ungroup() + cut_correspondence_slivers(d %>% select("GEOID90","GEOID00"), + pmax(area_share(d$area_part,d$area90), + area_share(d$area_part,d$area00)), + min_area_share) +} + +# stitch relationship files for consecutive censuses into one table spanning all +# requested censuses, dropping the vintages that only served as stepping stones +join_us_correspondence <- function(links, datasets){ + c <- links[[1]] + for (l in links[-1]) c <- full_join(c,l,by=intersect(names(c),names(l))) + c %>% + select(all_of(unname(us_geoid_columns[datasets]))) %>% + unique() } +# the censuses that have to be traversed to get from the earliest to the latest requested +# one, the Census Bureau only publishes relationship files between consecutive censuses +us_correspondence_span <- function(datasets, available){ + invalid_datasets <- setdiff(datasets,available) + if (length(invalid_datasets) > 0) { + stop(paste0("Invalid census years ",paste0(invalid_datasets,collapse=", "), + ", can only match censuses ",paste0(available,collapse=", "))) + } + datasets <- intersect(available,datasets) + if (length(datasets) < 2) { + stop("Need at least two censuses to build a correspondence table.") + } + available[seq(match(datasets[1],available),match(utils::tail(datasets,1),available))] +} + +get_us_ct_correspondence <- function(state, datasets, min_area_share=0.01, + cache_path=getOption("tongfen.cache_path")){ + available <- names(us_geoid_columns) + span <- us_correspondence_span(datasets,available) + links <- utils::head(span,-1) %>% + lapply(function(year){ + f <- switch(year, + dec1990 = get_us_ct_correspondence_2000, + dec2000 = get_us_ct_correspondence_2010, + dec2010 = get_us_ct_correspondence_2020) + f(state,min_area_share=min_area_share,cache_path=cache_path) + }) + join_us_correspondence(links,intersect(available,datasets)) +} + +# the 2000 to 2010 county subdivision comparability file, covering all states get_us_county_subdivision_correspondence <- function(cache_path=getOption("tongfen.cache_path")){ cache_path = file.path(cache_path %||% tempdir(),"us_data") file <- "Cousub_comparability.xlsx" @@ -98,6 +213,123 @@ get_us_county_subdivision_correspondence <- function(cache_path=getOption("tongf readxl::read_xlsx(local_path) } +# The 2010 to 2020 county subdivision relationship file, covering all states. Unlike the +# 2000 to 2010 comparability file this is a geometric overlay, most rows are slivers along +# boundaries that shifted slightly rather than actual relationships. Keeping them chains +# unrelated subdivisions into one common geography, so they get cut. A subdivision carved +# out of a larger one is only a small share of the old one but most of the new one, hence +# the share is taken over the larger of the two. +get_us_county_subdivision_correspondence_2020 <- function(min_area_share=0.01, + cache_path=getOption("tongfen.cache_path")){ + cache_path = file.path(cache_path %||% tempdir(),"us_data") + path <- paste0("https://www2.census.gov/geo/docs/maps-data/data/rel2020/cousub/", + "tab20_cousub20_cousub10_natl.txt") + local_path <- file.path(cache_path,basename(path)) + if (!file.exists(local_path)) { + if (!dir.exists(cache_path)) dir.create(cache_path) + utils::download.file(path,local_path,quiet=TRUE) + } + d <- readr::read_delim(local_path,delim="|",progress=FALSE, + col_types=readr::cols_only(GEOID_COUSUB_10="c",GEOID_COUSUB_20="c", + AREALAND_COUSUB_10="d",AREAWATER_COUSUB_10="d", + AREALAND_COUSUB_20="d",AREAWATER_COUSUB_20="d", + AREALAND_PART="d",AREAWATER_PART="d")) %>% + group_by(GEOID10=.data$GEOID_COUSUB_10,GEOID20=.data$GEOID_COUSUB_20) %>% + summarize(area_part=sum(.data$AREALAND_PART+.data$AREAWATER_PART), + area10=max(.data$AREALAND_COUSUB_10+.data$AREAWATER_COUSUB_10), + area20=max(.data$AREALAND_COUSUB_20+.data$AREAWATER_COUSUB_20), + .groups="drop") + cut_correspondence_slivers(d %>% select("GEOID10","GEOID20"), + pmax(area_share(d$area_part,d$area10), + area_share(d$area_part,d$area20)), + min_area_share) +} + +get_us_county_subdivision_correspondence_for <- function(state, datasets, min_area_share=0.01, + cache_path=getOption("tongfen.cache_path")){ + available <- setdiff(names(us_geoid_columns),"dec1990") + span <- us_correspondence_span(datasets,available) + fips <- fips_code_for_state(state)$state_code + links <- utils::head(span,-1) %>% + lapply(function(year){ + link <- switch(year, + dec2000 = get_us_county_subdivision_correspondence(cache_path=cache_path) %>% + select("GEOID00","GEOID10"), + dec2010 = get_us_county_subdivision_correspondence_2020( + min_area_share=min_area_share,cache_path=cache_path)) + # both files are national, county subdivisions don't cross state lines + link %>% filter(substr(.data$GEOID10,1,2)==fips) %>% unique() + }) + join_us_correspondence(links,intersect(available,datasets)) +} + + +#' Get correspondence table for US census geographies +#' +#' @description +#' \lifecycle{maturing} +#' +#' Builds a correspondence table matching US census geographies across censuses, based on the +#' relationship files published by the US Census Bureau. Censuses that aren't requested but sit +#' in between two that are get traversed on the way, the Census Bureau only publishes +#' relationship files between consecutive censuses. +#' +#' The relationship files are geometric overlays that list every sliver along boundaries that +#' only shifted slightly. Those get cut via `min_area_share`, keeping them would chain +#' unrelated regions into one common geography. +#' +#' The correspondence layer reaches back one census further than +#' \code{\link{get_tongfen_us_census}}. The 1990 census is available as `dec1990` here, but the +#' Census Bureau has retired the 1990 API endpoint, so 1990 data has to be brought in by other +#' means, for example from NHGIS via the ipumsr package, and handed to +#' \code{\link{tongfen_aggregate}} together with this correspondence table. +#' +#' @param datasets vector of censuses to match up, valid values are `dec1990`, `dec2000`, +#' `dec2010` and `dec2020` for census tracts, `dec2000` through `dec2020` for county +#' subdivisions. At least two censuses are needed. +#' @param regions list with regions to query the correspondence for. At this stage, the only +#' valid list is a vector of states, i.e. `regions = list(state=c("CA","OR"))` +#' @param level aggregation level, at this stage the only valid levels are 'tract' and +#' 'county subdivision'. +#' @param min_area_share minimum share of area two geographies have to have in common to count +#' as related, default is `0.01`. The Census Bureau relationship files list every geometric +#' overlap, lowering this pulls in slivers along boundaries that only shifted slightly and +#' chains unrelated regions into one common geography. Raising it gives finer common +#' geographies at the risk of separating regions that did change. No region is ever dropped, +#' if all of its parts are slivers its largest part is kept. +#' @param cache_path optional path to cache the relationship files in, defaults to the +#' `tongfen.cache_path` option and falls back to a temporary directory +#' @return tibble with one row per census geography, a GEOID column for each requested census, +#' and the common geography identified by `TongfenID` and `TongfenUID`. +#' @export +#' +#' @examples +#' # Match up census tracts for the 1990 and 2000 censuses in Rhode Island +#' \dontrun{ +#' correspondence <- get_tongfen_correspondence_us_census(datasets = c("dec1990","dec2000"), +#' regions = list(state="RI")) +#'} +get_tongfen_correspondence_us_census <- function(datasets, regions, level='tract', + min_area_share=0.01, + cache_path=getOption("tongfen.cache_path")){ + assert(level %in% c('tract','county subdivision'), + "Only census tracts and county subdivisions are supported right now.") + + regions$state %>% + lapply(function(state){ + if (level=='tract') { + get_us_ct_correspondence(state,datasets,min_area_share=min_area_share, + cache_path=cache_path) + } else { + get_us_county_subdivision_correspondence_for(state,datasets, + min_area_share=min_area_share, + cache_path=cache_path) + } + }) %>% + bind_rows() %>% + get_tongfen_correspondence() +} + valid_us_census_datasets <- c( dec2000 = "US decentennial census 2000", @@ -105,6 +337,18 @@ valid_us_census_datasets <- c( dec2020 = "US decentennial census 2020" ) +# Censuses that published several summary files need to be told which one to read. tidycensus +# picks a default per census year, for 2020 that is the PL 94-171 redistricting file which only +# carries a handful of variables, most 2020 variables live in the DHC file. `sumfile` is either a +# single value for all censuses or a vector named by dataset, `NULL` leaves the choice to +# tidycensus. +sumfile_for_dataset <- function(sumfile, ds){ + if (is.null(sumfile)) return(NULL) + if (is.null(names(sumfile))) return(unname(sumfile)) + if (!(ds %in% names(sumfile))) return(NULL) + unname(sumfile[[ds]]) +} + #' Get US census data for 2000 and 2010 census on common census tract based geography #' #' @description @@ -113,12 +357,23 @@ valid_us_census_datasets <- c( #' This wraps data acquisition via the tidycensus package and tongfen on a common geography into #' a single convenience function. #' +#' Data is only available for the 2000, 2010 and 2020 censuses, the Census Bureau has retired the +#' 1990 API endpoint. To tongfen 1990 data, obtain it elsewhere and combine it with a +#' correspondence table from \code{\link{get_tongfen_correspondence_us_census}} via +#' \code{\link{tongfen_aggregate}}. +#' #' @param regions list with regions to query the data for. At this stage, the only #' valid list is a vector of states, i.e. `regions = list(state=c("CA","OR"))`` #' @param meta metadata for variables to retrieve #' @param level aggregation level to return the data on. At this stage, the only valid levels are 'tract' and 'county subdivision'. #' @param survey survey to get data for, supported options is "census" #' @param base_geo census year to use as base geography, default is `2010`. +#' @param min_area_share minimum share of area two geographies have to have in common to count +#' as related, default is `0.01`, see \code{\link{get_tongfen_correspondence_us_census}}. +#' @param sumfile summary file to read the variables from, either a single value used for all +#' censuses or a vector named by dataset, for example `c(dec2010="sf1", dec2020="dhc")`. Default +#' is `NULL`, which leaves the choice to tidycensus. Note that tidycensus defaults the 2020 +#' census to the PL 94-171 redistricting file, most 2020 variables need `sumfile="dhc"`. #' @return sf object with (wide form) census variables with census year as suffix (separated by underdcore "_"). #' @export #' @@ -139,7 +394,7 @@ valid_us_census_datasets <- c( #' #'} get_tongfen_us_census <- function(regions,meta,level='tract',survey="census", - base_geo = NULL){ + base_geo = NULL, min_area_share = 0.01, sumfile = NULL){ datasets <- meta$dataset %>% unique if (is.null(base_geo)) base_geo=datasets[1] @@ -148,20 +403,22 @@ get_tongfen_us_census <- function(regions,meta,level='tract',survey="census", assert(length(invalid_datasets)==0, paste0("Invalid datasets :",paste0(invalid_datasets,collapse = ", "))) assert(level %in% c('tract','county subdivision'),"Only census tracts and counties are supported right now.") assert(survey %in% c('census'),"Only census surveys are supported right now.") - - regions$state %>% lapply(function(state){ - if (level=='tract') { - correspondence <- get_us_ct_correspondence(state,datasets) - } else if (level=="county subdivision") { - fips <- fips_code_for_state(state)$state_code - correspondence <- get_us_county_subdivision_correspondence() %>% - filter(.data$STATEFP10==fips) %>% - select(.data$GEOID00,.data$GEOID10) + if (!is.null(sumfile)) { + if (is.null(names(sumfile))) { + assert(length(sumfile)==1, + "sumfile has to be a single value or a vector named by dataset") } else { - stop("Ooops, should have caught this earler.") + invalid_sumfiles <- setdiff(names(sumfile),datasets) + assert(length(invalid_sumfiles)==0, + paste0("Invalid datasets in sumfile: ",paste0(invalid_sumfiles,collapse=", "))) } - correspondence <- correspondence %>% - get_tongfen_correspondence() + } + + regions$state %>% lapply(function(state){ + correspondence <- get_tongfen_correspondence_us_census(datasets = datasets, + regions = list(state=state), + level = level, + min_area_share = min_area_share) data <- datasets %>% lapply(function(ds){ @@ -170,6 +427,7 @@ get_tongfen_us_census <- function(regions,meta,level='tract',survey="census", short_year <- substr(as.character(year),3,4) tidycensus::get_decennial(geography=level, state=state, county=regions$county, variables = m$variable, year = year, + sumfile = sumfile_for_dataset(sumfile,ds), geometry = base_geo==ds, output="wide") %>% rename(!!paste0("GEOID",short_year):=.data$GEOID) }) %>% diff --git a/cran-comments.md b/cran-comments.md index 4459ea3..568f0c3 100644 --- a/cran-comments.md +++ b/cran-comments.md @@ -1,3 +1,27 @@ +# tongfen v.0.3.8 +## Breaking changes +- `get_tongfen_ca_census` now honours its `base_geo`, `na.rm`, `tolerance`, `crs` and + `data_transform` arguments, all of which were silently ignored +- removed the `area_mismatch_cutoff` argument from `get_tongfen_ca_census` and + `get_tongfen_correspondence_ca_census`, it never had any effect +## Major changes +- correspondence tables are now built via a vectorised connected components pass instead of + a row-by-row union-find, making tongfen on large geographies dramatically faster +- the "statcan" method no longer downloads census geometries it does not use +- dissolving geometries skips regions that don't need to be merged +- new `get_tongfen_correspondence_us_census`, US tract correspondence tables now reach back to + the 1990 census and county subdivisions forward to the 2020 census +- US correspondence tables no longer chain regions together over slivers, and no longer strip + leading zeros off 2020 census tract identifiers +## Minor changes +- `get_tongfen_us_census` gained a `sumfile` argument, passed through to tidycensus +- `get_tongfen_correspondence_ca_census` gained a `crs` argument for the spatial intersections +- missing geographic identifiers no longer merge unrelated regions into one common geography +- fix crash when tongfen-ing census tracts across non-adjacent censuses +- US county subdivision data errors out up front on censuses it can't be matched across +- several fixes to the deprecated `get_tongfen_census_*` functions +- faster `check_tongfen_areas` and `aggregate_correspondences` + # tongfen v.0.3.7 ## Major changes - accommodate factors in proportional_reaggregate @@ -26,7 +50,12 @@ # Submission - v.0.3 # Test environments -* local R installation, R 4.0.2 -* GitHub action release +* local macOS installation, R 4.6.0 +* GitHub actions (windows-latest, macOS-latest, ubuntu-latest) on release, devel and oldrel + +# R CMD check results +0 errors | 0 warnings | 0 notes + +There are no reverse dependencies. diff --git a/docs/404.html b/docs/404.html index c868523..d0348fd 100644 --- a/docs/404.html +++ b/docs/404.html @@ -12,8 +12,8 @@ - - + + @@ -30,7 +30,7 @@ tongfen - 0.3.7 + 0.3.8 + + + + + +
+
+
+

Get correspondence table for US census geographies

+ Source: R/tongfen_us.R +
get_tongfen_correspondence_us_census.Rd
+
+ +
+

[Maturing]

+

Builds a correspondence table matching US census geographies across censuses, based on the +relationship files published by the US Census Bureau. Censuses that aren't requested but sit +in between two that are get traversed on the way, the Census Bureau only publishes +relationship files between consecutive censuses.

+

The relationship files are geometric overlays that list every sliver along boundaries that +only shifted slightly. Those get cut via `min_area_share`, keeping them would chain +unrelated regions into one common geography.

+

The correspondence layer reaches back one census further than +get_tongfen_us_census. The 1990 census is available as `dec1990` here, but the +Census Bureau has retired the 1990 API endpoint, so 1990 data has to be brought in by other +means, for example from NHGIS via the ipumsr package, and handed to +tongfen_aggregate together with this correspondence table.

+
+ +
+

Usage

+ 
get_tongfen_correspondence_us_census(
+  datasets,
+  regions,
+  level = "tract",
+  min_area_share = 0.01,
+  cache_path = getOption("tongfen.cache_path")
+)
+
+ +
+

Arguments

+ + +
datasets
+

vector of censuses to match up, valid values are `dec1990`, `dec2000`, +`dec2010` and `dec2020` for census tracts, `dec2000` through `dec2020` for county +subdivisions. At least two censuses are needed.

+ + +
regions
+

list with regions to query the correspondence for. At this stage, the only +valid list is a vector of states, i.e. `regions = list(state=c("CA","OR"))`

+ + +
level
+

aggregation level, at this stage the only valid levels are 'tract' and +'county subdivision'.

+ + +
min_area_share
+

minimum share of area two geographies have to have in common to count +as related, default is `0.01`. The Census Bureau relationship files list every geometric +overlap, lowering this pulls in slivers along boundaries that only shifted slightly and +chains unrelated regions into one common geography. Raising it gives finer common +geographies at the risk of separating regions that did change. No region is ever dropped, +if all of its parts are slivers its largest part is kept.

+ + +
cache_path
+

optional path to cache the relationship files in, defaults to the +`tongfen.cache_path` option and falls back to a temporary directory

+ +
+
+

Value

+

tibble with one row per census geography, a GEOID column for each requested census, +and the common geography identified by `TongfenID` and `TongfenUID`.

+
+ +
+

Examples

+ 
# Match up census tracts for the 1990 and 2000 censuses in Rhode Island
+if (FALSE) { # \dontrun{
+correspondence <- get_tongfen_correspondence_us_census(datasets = c("dec1990","dec2000"),
+                                                       regions = list(state="RI"))
+} # }
+
+
+
+ + +
+ + + + + + + diff --git a/docs/reference/get_tongfen_correspondence_us_census.md b/docs/reference/get_tongfen_correspondence_us_census.md new file mode 100644 index 0000000..aeaeede --- /dev/null +++ b/docs/reference/get_tongfen_correspondence_us_census.md @@ -0,0 +1,86 @@ +# Get correspondence table for US census geographies + +**\[maturing\]** + +Builds a correspondence table matching US census geographies across +censuses, based on the relationship files published by the US Census +Bureau. Censuses that aren't requested but sit in between two that are +get traversed on the way, the Census Bureau only publishes relationship +files between consecutive censuses. + +The relationship files are geometric overlays that list every sliver +along boundaries that only shifted slightly. Those get cut via +\`min_area_share\`, keeping them would chain unrelated regions into one +common geography. + +The correspondence layer reaches back one census further than +[`get_tongfen_us_census`](https://mountainmath.github.io/tongfen/reference/get_tongfen_us_census.md). +The 1990 census is available as \`dec1990\` here, but the Census Bureau +has retired the 1990 API endpoint, so 1990 data has to be brought in by +other means, for example from NHGIS via the ipumsr package, and handed +to +[`tongfen_aggregate`](https://mountainmath.github.io/tongfen/reference/tongfen_aggregate.md) +together with this correspondence table. + +## Usage + +``` r +get_tongfen_correspondence_us_census( + datasets, + regions, + level = "tract", + min_area_share = 0.01, + cache_path = getOption("tongfen.cache_path") +) +``` + +## Arguments + +- datasets: + + vector of censuses to match up, valid values are \`dec1990\`, + \`dec2000\`, \`dec2010\` and \`dec2020\` for census tracts, + \`dec2000\` through \`dec2020\` for county subdivisions. At least two + censuses are needed. + +- regions: + + list with regions to query the correspondence for. At this stage, the + only valid list is a vector of states, i.e. \`regions = + list(state=c("CA","OR"))\` + +- level: + + aggregation level, at this stage the only valid levels are 'tract' and + 'county subdivision'. + +- min_area_share: + + minimum share of area two geographies have to have in common to count + as related, default is \`0.01\`. The Census Bureau relationship files + list every geometric overlap, lowering this pulls in slivers along + boundaries that only shifted slightly and chains unrelated regions + into one common geography. Raising it gives finer common geographies + at the risk of separating regions that did change. No region is ever + dropped, if all of its parts are slivers its largest part is kept. + +- cache_path: + + optional path to cache the relationship files in, defaults to the + \`tongfen.cache_path\` option and falls back to a temporary directory + +## Value + +tibble with one row per census geography, a GEOID column for each +requested census, and the common geography identified by \`TongfenID\` +and \`TongfenUID\`. + +## Examples + +``` r +# Match up census tracts for the 1990 and 2000 censuses in Rhode Island +if (FALSE) { # \dontrun{ +correspondence <- get_tongfen_correspondence_us_census(datasets = c("dec1990","dec2000"), + regions = list(state="RI")) +} # } +``` diff --git a/docs/reference/get_tongfen_us_census.html b/docs/reference/get_tongfen_us_census.html index 8186ad9..57f67e0 100644 --- a/docs/reference/get_tongfen_us_census.html +++ b/docs/reference/get_tongfen_us_census.html @@ -1,9 +1,17 @@ -Get US census data for 2000 and 2010 census on common census tract based geography — get_tongfen_us_census • tongfenGet US census data for 2000 and 2010 census on common census tract based geography — get_tongfen_us_census • tongfen +a single convenience function. +Data is only available for the 2000, 2010 and 2020 censuses, the Census Bureau has retired the +1990 API endpoint. To tongfen 1990 data, obtain it elsewhere and combine it with a +correspondence table from get_tongfen_correspondence_us_census via +tongfen_aggregate."> Skip to contents @@ -11,7 +19,7 @@ tongfen - 0.3.7 + 0.3.8