v0.3.8

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.

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)

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

R/helpers.R

@@ -24,116 +24,184 @@ 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<length(hs) && filter(ddd,is.na(.data$TongfenID)) %>% 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) {
   if (! expr) stop(error, call. = FALSE)
 }
 
 
+# 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")
 }