diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml index bd2aadc..472ced2 100644 --- a/.github/workflows/R-CMD-check.yaml +++ b/.github/workflows/R-CMD-check.yaml @@ -2,7 +2,7 @@ # Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help on: push: - branches: [main, master] + branches: [main] pull_request: schedule: - cron: '0 7 * * *' diff --git a/DESCRIPTION b/DESCRIPTION index c014060..2631df5 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: tongfen Type: Package Title: Make Data Based on Different Geographies Comparable -Version: 0.3.8 +Version: 0.3.9 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. @@ -11,7 +11,7 @@ ByteCompile: yes LazyData: true NeedsCompilation: no Imports: - dplyr (>= 1.0), + dplyr (>= 1.1.0), tidyr (>= 1.0), sf, tibble, @@ -19,9 +19,10 @@ Imports: purrr, stringr, readr, + nanoparquet, + tools, utils, lifecycle -RoxygenNote: 7.3.3 Suggests: knitr, rmarkdown, @@ -34,7 +35,7 @@ Suggests: readxl, scales, microbenchmark, - testthat (>= 3.0.0) + testthat (>= 3.2.0) VignetteBuilder: knitr, rmarkdown URL: https://github.com/mountainMath/tongfen, https://mountainmath.github.io/tongfen/ BugReports: https://github.com/mountainMath/tongfen/issues @@ -42,3 +43,4 @@ Language: en-US RdMacros: lifecycle Depends: R (>= 4.1) +Config/roxygen2/version: 8.1.0 diff --git a/NAMESPACE b/NAMESPACE index 1ba1cc6..16f2796 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -14,9 +14,13 @@ export(meta_for_additive_variables) export(meta_for_ca_census_vectors) export(proportional_reaggregate) export(tongfen_aggregate) +export(tongfen_anomaly_joins) export(tongfen_ca_census_ct) +export(tongfen_detect_anomalies) export(tongfen_estimate) export(tongfen_estimate_ca_census) +export(tongfen_join_correspondence) +export(tongfen_join_regions) export(tongfen_tag_largest_overlap) import(dplyr) import(rlang) diff --git a/NEWS.md b/NEWS.md index e9babd4..2c97b3e 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,57 @@ +# tongfen v.0.3.9 +## Major changes +- new experimental functions to detect and correct for likely geocoding anomalies in timelines on a + common geography, where the same dwellings got assigned to different neighbouring regions in + different years. This shows up as a surprising drop in one region that is offset by a jump in a + neighbouring region. `tongfen_detect_anomalies` lists the candidate regions, + `tongfen_anomaly_joins` determines the regions to join, `tongfen_join_regions` joins them in + data that already is on a common geography and `tongfen_join_correspondence` joins them in a + correspondence for use with `tongfen_aggregate`. See the new "Geocoding anomalies in TongFen + timelines" vignette for details +- StatCan correspondence files are now downloaded as parquet files from a mirror, Statistics Canada + put the original files behind a browser check that blocks programmatic downloads, which broke + `method = "statcan"`. Cached files are checked against the mirror once per session and downloaded + again if they changed. The mirror location can be changed via the `tongfen.statcan_correspondence_url` + option. Previously cached `statcan_correspondence_*.csv` files in the tongfen cache directory are + no longer used and can be removed +## Minor changes +- combining correspondences across three or more datasets no longer runs into a cross join when the + correspondences happen to be ordered so that consecutive ones don't share a geographic identifier. + This gave a dplyr deprecation warning and needlessly large intermediate tables, results are unchanged +- fix `proportional_reaggregate` giving wrong results when the finer level data already has values + for the categories to reaggregate. The values were compared to the parent total across all + categories instead of per category, and a missing value in one child region discarded the + existing values of all its siblings +- fix `tongfen_estimate` underestimating values when the intersection of a source and a target + region is a geometry collection, i.e. several polygons joined by a shared boundary line +- `tongfen_estimate` with `na.rm = FALSE` no longer returns `NA` for target regions that only + share a boundary with a source region with missing values +- `tongfen_estimate` now returns `NA` for target regions that don't overlap the source instead of + erroring out when none of them do, and gives a clear error when `target` already has a column + named like one of the variables to estimate +- fix `tongfen_aggregate` and `aggregate_data_with_meta` scaling averages by their parent variable + more than once when the metadata lists the same variable name for several datasets, as is common + with US census data. `tongfen_aggregate` now only uses the metadata of the dataset being + aggregated, conflicting aggregation rules for the same variable are an error +- averages aggregated with `na.rm = TRUE` are now taken over the regions that have a value. The + parent variable of regions with a missing average still counted toward the total the average + was divided by, pulling the result toward zero. This affects `aggregate_data_with_meta`, + `tongfen_aggregate`, `tongfen_estimate` and the functions built on them +- fix "Average to" variables like percentage changes overwriting each other's base when several of + them share a parent variable. In that case the base columns in the result are named after the + variable (`base_`) instead of the parent +- `estimate_tongfen_correspondence` no longer requires the geometry column to be named `geometry` +- `refresh = TRUE` in `get_tongfen_correspondence_ca_census` and `get_tongfen_ca_census` now also + refreshes the cached StatCan correspondence files +- US Census Bureau relationship files are downloaded to a temporary file first, an interrupted + download no longer leaves a broken file in the cache. They are now cached in the same place as + the StatCan correspondence files, also honouring the `tongfen.cache_path` environment variable + and the `custom_data_path` option +- `tongfen_estimate_ca_census` returns its result visibly +- documentation fixes, among others the `tongfen_aggregate` example now passes a named list of + datasets matching the metadata +- requires dplyr 1.1.0 or newer, which the package already relied on + # tongfen v.0.3.8 ## Breaking changes - `get_tongfen_ca_census` now honours its `base_geo`, `na.rm`, `tolerance`, `crs` and @@ -68,12 +122,12 @@ - squish several edge case bugs # tongfen v.0.3.6 -## Major changs +## Major changes - better downsampling that can also accommodate averages - performance improvements ## Minor changes - better documentation -- allow for datasets vartiables by census year for canadian data +- allow for datasets variables by census year for Canadian data - fix issue where some metadata might get duplicated # tongfen 0.3.2 @@ -85,7 +139,7 @@ ## Major changes - Added `tongfen_estimate_ca_census` function for new CensusMapper endpoint, tying into new {cancensus} functionality. ## Minor changes -- Custom impelementation of `tongfen_etimate` for finer control +- Custom implementation of `tongfen_estimate` for finer control - Fix compatibility issue with changes in {sf} package # tongfen 0.3 diff --git a/R/helpers.R b/R/helpers.R index eb527c3..2e15131 100644 --- a/R/helpers.R +++ b/R/helpers.R @@ -14,6 +14,57 @@ tongfen_cache_dir <- function(){ tempdir() } +tongfen_session <- new.env(parent=emptyenv()) + +# ETag of a remote file, NULL if it can't be determined, e.g. when offline +remote_etag <- function(url){ + headers <- tryCatch(suppressWarnings(curlGetHeaders(url)),error=function(e) NULL) + if (is.null(headers) || !identical(attr(headers,"status"),200L)) return(NULL) + etag <- grep("^etag:",headers,ignore.case=TRUE,value=TRUE) + if (length(etag)==0) return(NULL) + gsub("^etag:\\s*|\"|\\s+$","",etag[length(etag)],ignore.case=TRUE) +} + +# location of the cached US Census Bureau relationship files +us_cache_dir <- function(cache_path=NULL){ + file.path(nullify_blank(cache_path) %||% tongfen_cache_dir(),"us_data") +} + +# Download a remote file to the local path unless the local copy is still current. +# The ETag of the downloaded file is kept next to the cached file and compared to the remote ETag +# the first time the file is requested in a session, the file is only downloaded again if it changed. +# Files that never change don't need to be checked against the remote, with `check_remote=FALSE` +# the cached file is used as is. The download goes to a temporary file first so that an +# interrupted download does not leave a broken file in the cache. +cached_download <- function(url,path,refresh=FALSE,check_remote=TRUE){ + etag_path <- paste0(path,".etag") + cached <- file.exists(path) && !refresh + if (cached && (!check_remote || isTRUE(tongfen_session[[url]]))) return(path) + etag <- if (check_remote) remote_etag(url) + if (cached) { + if (is.null(etag)) { + message(paste0("Could not check ",url," for updates, using cached version.")) + } + local_etag <- if (file.exists(etag_path)) readLines(etag_path,n=1,warn=FALSE) + if (is.null(etag) || identical(etag,local_etag)) { + tongfen_session[[url]] <- TRUE + return(path) + } + } + if (!dir.exists(dirname(path))) dir.create(dirname(path),recursive=TRUE) + tmp <- tempfile(tmpdir=dirname(path)) + on.exit(unlink(tmp)) + utils::download.file(url,tmp,mode="wb",quiet=TRUE) + # S3 ETags of files that were not uploaded in parts are the md5 checksum of the file + if (!is.null(etag) && grepl("^[0-9a-f]{32}$",etag) && !identical(unname(tools::md5sum(tmp)),etag)) { + stop(paste0("Download of ",url," is corrupted, please try again.")) + } + file.copy(tmp,path,overwrite=TRUE) + if (is.null(etag)) unlink(etag_path) else writeLines(etag,etag_path) + tongfen_session[[url]] <- TRUE + path +} + inner_join_tongfen_correspondence <- function(data,correspondence,link){ data %>% inner_join(correspondence %>% @@ -140,6 +191,15 @@ assert <- function (expr, error) { if (! expr) stop(error, call. = FALSE) } +# Pairs of intersecting geometries as row indices into `x` and `y`. The sparse index +# list is turned into a tibble directly, `as.data.frame()` on an empty result drops +# the columns we join on. +intersects_pairs <- function(x, y) { + m <- sf::st_intersects(x, y, sparse = TRUE) + tibble(row.id = rep(seq_along(m), lengths(m)), + col.id = as.integer(unlist(m))) +} + # Dissolve the geometries of `data` by `grouping_var`, the geometric equivalent # of `summarize()`. Groups holding a single geometry - the bulk of the groups @@ -196,18 +256,22 @@ aggregate_correspondences <- function(correspondences){ select(!matches("Tongfen") | matches("TongfenMethod")) } # compute full correspondence, smallest table first to keep intermediate - # join results as small as possible - index_order <- correspondences %>% lapply(nrow) %>% unlist() %>% order() - - correspondence <- correspondences[[index_order[1]]] %>% - clean_correspondence_names() - if (length(correspondences)>1) for (index in index_order[-1]) { - c <- correspondences[[index]] %>% - clean_correspondence_names() - match_columns <- intersect(names(correspondence),names(c)) - match_columns <- match_columns[!grepl("TongfenMethod",match_columns)] - correspondence <- inner_join(correspondence,c,by=match_columns) %>% + # join results as small as possible, but only join tables that share an identifier + # with the tables joined so far, joining unrelated tables gives a cross join + remaining <- correspondences[order(vapply(correspondences,nrow,integer(1)))] %>% + lapply(clean_correspondence_names) + correspondence <- remaining[[1]] + remaining <- remaining[-1] + while (length(remaining)>0) { + match_columns <- lapply(remaining,function(c) { + match_columns <- intersect(names(correspondence),names(c)) + match_columns[!grepl("TongfenMethod",match_columns)] + }) + index <- which(lengths(match_columns)>0)[1] + if (is.na(index)) stop("Correspondences can't be combined, they don't share a common geographic identifier.") + correspondence <- inner_join(correspondence,remaining[[index]],by=match_columns[[index]]) %>% unique() + remaining <- remaining[-index] } method_columns <- names(correspondence)[grepl("TongfenMethod",names(correspondence))] diff --git a/R/tonfen_deprecated.R b/R/tonfen_deprecated.R index 9efe73c..372ed12 100644 --- a/R/tonfen_deprecated.R +++ b/R/tonfen_deprecated.R @@ -1,4 +1,4 @@ -#' Check geographic integrety +#' Check geographic integrity #' #' @description #' \lifecycle{deprecated} diff --git a/R/tongfen.R b/R/tongfen.R index 5640454..f81c484 100644 --- a/R/tongfen.R +++ b/R/tongfen.R @@ -5,8 +5,8 @@ #' #' Generates metadata to be used in tongfen_aggregate. Variables need to be additive like counts. #' -#' @param dataset identifier for the dataset contianing the variable -#' @param variables (named) vecotor with additive variables +#' @param dataset identifier for the dataset containing the variable +#' @param variables (named) vector with additive variables #' @return a tibble to be used in tongfen_aggregate #' @export #' @@ -29,6 +29,21 @@ meta_for_additive_variables <- function(dataset,variables){ +# Averages are aggregated by scaling them by their parent variable, summing up, and dividing +# by the summed up parent again. The parent of regions where the average is missing must not +# count toward the total the sum gets divided by, so each average carries its own weight column. +average_weight_name <- function(variables) { + if (length(variables)==0) return(character(0)) + paste0("...weight_",variables) +} + +average_weight_exprs <- function(to_scale,parent_lookup) { + lapply(setNames(to_scale, average_weight_name(to_scale)), \(col) { + parent <- as.name(unname(parent_lookup[col])) + rlang::expr(replace(as.numeric(!!parent),is.na(!!as.name(col)),NA_real_)) + }) +} + cut_meta <- function(data,meta){ meta <- meta %>% filter(.data$variable %in% names(data)|.data$label %in% names(data)) %>% @@ -52,14 +67,16 @@ pre_scale <- function(data,meta,meta_var="data_var",quiet=FALSE) { if (length(to_scale) > 0) { + weight_exprs <- average_weight_exprs(to_scale,parent_lookup) scale_exprs <- lapply(setNames(to_scale, to_scale), \(col) rlang::expr(!!as.name(col) * !!as.name(unname(parent_lookup[col])))) - data <- data %>% mutate(!!!scale_exprs) + data <- data %>% mutate(!!!weight_exprs) %>% mutate(!!!scale_exprs) } data } +# divides by the weights added in `pre_scale` if they are still around, and by the parent otherwise post_scale <- function(data,meta,meta_var="data_var") { meta_name_lookup <- setNames(meta %>% pull(meta_var),meta$variable) meta$parent_name <- meta_name_lookup[meta$parent] @@ -67,9 +84,14 @@ post_scale <- function(data,meta,meta_var="data_var") { to_scale <- filter(meta,.data$rule %in% c("Median","Average")) %>% pull(meta_var) if (length(to_scale) > 0) { - scale_exprs <- lapply(setNames(to_scale, to_scale), \(col) - rlang::expr(!!as.name(col) / !!as.name(unname(parent_lookup[col])))) - data <- data %>% mutate(!!!scale_exprs) + scale_exprs <- lapply(setNames(to_scale, to_scale), \(col) { + weight <- average_weight_name(col) + if (!(weight %in% names(data))) weight <- unname(parent_lookup[col]) + rlang::expr(!!as.name(col) / !!as.name(weight)) + }) + data <- data %>% + mutate(!!!scale_exprs) %>% + select(-any_of(average_weight_name(to_scale))) } data @@ -102,6 +124,17 @@ post_scale <- function(data,meta,meta_var="data_var") { #'} aggregate_data_with_meta <- function(data,meta,geo=FALSE,na.rm=TRUE,quiet=FALSE){ meta <- meta %>% filter(.data$variable %in% names(data)) + # the same variable can be listed several times if meta spans several datasets, + # each variable must only be aggregated (and scaled) once + duplicates <- duplicated(meta$variable) + if (any(duplicates)) { + rules <- meta %>% select(any_of(c("variable","rule","parent","units"))) %>% unique() + ambiguous <- unique(rules$variable[duplicated(rules$variable)]) + if (length(ambiguous)>0) + stop(paste0("Conflicting aggregation rules in metadata for ",paste0(ambiguous,collapse = ", "), + ", please only pass the metadata for the dataset that is being aggregated.")) + meta <- meta[!duplicates,] + } grouping_var=groups(data) %>% as.character parent_lookup <- setNames(meta$parent,meta$variable) to_scale <- filter(meta,.data$rule %in% c("Median","Average"))$variable @@ -116,16 +149,25 @@ aggregate_data_with_meta <- function(data,meta,geo=FALSE,na.rm=TRUE,quiet=FALSE) message(paste0("Can't TongFen medians, will approximate by treating as averages: ",paste0(median_vars,collapse = ", "))) } + weight_variables <- average_weight_name(to_scale) if (length(to_scale) > 0) { + weight_exprs <- average_weight_exprs(to_scale,parent_lookup) scale_exprs <- lapply(setNames(to_scale, to_scale), \(col) rlang::expr(!!as.name(col) * !!as.name(unname(parent_lookup[col])))) - data <- data %>% mutate(!!!scale_exprs) + data <- data %>% mutate(!!!weight_exprs) %>% mutate(!!!scale_exprs) } + # the base an "Average to" variable gets scaled by depends on the variable, variables + # sharing a parent each need their own base column + scale_from_parents <- unname(parent_lookup[to_scale_from]) + shared_parent <- scale_from_parents %in% scale_from_parents[duplicated(scale_from_parents)] + base_vectors <- setNames(paste0("base_",ifelse(shared_parent,to_scale_from,scale_from_parents)), + to_scale_from) + base_variables <- c() for (x in to_scale_from) { scale_type <- meta %>% filter(.data$variable==x) %>% pull(units) %>% as.character() - base_vector <- paste0("base_",parent_lookup[x]) + base_vector <- unname(base_vectors[x]) base_variables <- c(base_variables,base_vector) if (scale_type=="Percentage ratio (0.0-1.0)") { data <- data %>% mutate(!!base_vector:=!!as.name(parent_lookup[x])/(!!as.name(x)+1)) @@ -144,22 +186,24 @@ aggregate_data_with_meta <- function(data,meta,geo=FALSE,na.rm=TRUE,quiet=FALSE) 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), + summarize_at(c(meta$variable,weight_variables),sum,na.rm=na.rm), by=grouping_var) } else { - data <- data %>% summarize_at(meta$variable,sum,na.rm=na.rm) + data <- data %>% summarize_at(c(meta$variable,weight_variables),sum,na.rm=na.rm) } if (length(to_scale) > 0) { scale_exprs <- lapply(setNames(to_scale, to_scale), \(col) - rlang::expr(!!as.name(col) / !!as.name(unname(parent_lookup[col])))) - data <- data %>% mutate(!!!scale_exprs) + rlang::expr(!!as.name(col) / !!as.name(average_weight_name(col)))) + data <- data %>% + mutate(!!!scale_exprs) %>% + select(-all_of(weight_variables)) } # Optimized: Vectorized division by base vectors if (length(to_scale_from) > 0) { for (x in to_scale_from) { - base_vector <- paste0("base_", parent_lookup[x]) + base_vector <- unname(base_vectors[x]) data[[x]] <- data[[x]] / data[[base_vector]] } } @@ -186,9 +230,12 @@ rename_with_meta <- function(data,meta,ds=NULL){ #' @description #' \lifecycle{maturing} #' -#' Aggregate variables secified in meta for several datasets according to correspondence. +#' Aggregate variables specified in meta for several datasets according to correspondence. #' -#' @param data list of datasets to be aggregated +#' @param data named list of datasets to be aggregated. The names identify the datasets, they are +#' matched against the `geo_dataset` column in `meta` to pick the aggregation rules and labels +#' for each dataset. Without names, or with names not found in `meta`, the rules for all +#' datasets are applied and the variables keep their original names #' @param correspondence correspondence data for gluing up the datasets #' @param meta metadata containing aggregation rules as for example returned by `meta_for_ca_census_vectors` #' @param base_geo identifier for which data element to base the final geography on, @@ -200,17 +247,17 @@ rename_with_meta <- function(data,meta,ds=NULL){ #' @export #' #' @examples -#' # aggregate census tract level 2006 population data on common gepgraphy build through +#' # aggregate census tract level 2006 and 2016 population data on common geography built through #' # correspondence from 2006 and 2016 census tracts in the City of Vancouver. #' \dontrun{ #' regions <- list(CSD="5915022") #' geo1 <- cancensus::get_census("CA06",regions=regions,geo_format='sf',level='CT') #' geo2 <- cancensus::get_census("CA16",regions=regions,geo_format='sf',level='CT') -#' meta <- meta_for_additive_variables("CA06","Population") +#' meta <- meta_for_additive_variables(c("CA06","CA16"),"Population") #' correspondence <- get_tongfen_correspondence_ca_census(geo_datasets=c('CA06','CA16'), #' regions=regions,level='CT') -#' result <- tongfen_aggregate(list(geo1 %>% rename(GeoUIDCA06=GeoUID), -#' geo2 %>% rename(GeoUIDCA16=GeoUID)),correspondence,meta) +#' result <- tongfen_aggregate(list(CA06=geo1 %>% rename(GeoUIDCA06=GeoUID), +#' CA16=geo2 %>% rename(GeoUIDCA16=GeoUID)),correspondence,meta) #'} tongfen_aggregate <- function(data,correspondence,meta=NULL, base_geo = NULL, na.rm = TRUE){ data <- ensure_names(data) @@ -239,7 +286,13 @@ tongfen_aggregate <- function(data,correspondence,meta=NULL, base_geo = NULL, na by=match_column) %>% group_by(.data$TongfenID,.data$TongfenUID) if (!is.null(meta)) { - d <- d %>% aggregate_data_with_meta(meta,na.rm=na.rm) + # only use the aggregation rules for this dataset if meta distinguishes datasets, + # the same variable name can come with different rules in different datasets + ds_meta <- meta + if (ds %in% as.character(meta$geo_dataset)) { + ds_meta <- meta %>% filter(as.character(.data$geo_dataset)==ds) + } + d <- d %>% aggregate_data_with_meta(ds_meta,na.rm=na.rm) } else { if ("sf" %in% class(d)) { d <- summarize_geometry_by_group(d,c("TongfenID","TongfenUID")) @@ -291,7 +344,7 @@ tongfen_aggregate <- function(data,correspondence,meta=NULL, base_geo = NULL, na #' @param geo_match A named string informing on what column names to match data and parent_data #' @param categories Vector of column names to re-aggregate #' @param base Column name to use for proportional weighting when re-aggregating, or named vector with column name for each category. -#' Categries that should be re-aggregated as means should be set to NA and will only be reaggregated if the base data has NA values. +#' Categories that should be re-aggregated as means should be set to NA and will only be reaggregated if the base data has NA values. #' @return dataframe with downsampled variables from parent_data #' @keywords reaggregate proportionally wrt base variable #' @export @@ -397,7 +450,9 @@ proportional_reaggregate <- function(data,parent_data,geo_match,categories,base= values_to="p_value") if (vt %in% c("numeric","integer","integer64")) { d_combined <- full_join(d_base,d_parent,by=c(geo_match,"category"="category")) %>% - mutate(s_value=sum(.data$value),.by=names(geo_match)) %>% + # the parent value gets compared to what the children of the same category + # already add up to, missing child values count as zero + mutate(s_value=sum(.data$value,na.rm=TRUE),.by=c(names(geo_match),"category")) %>% mutate(across(any_of(c("p_value","s_value")),\(x)coalesce(x,0))) %>% mutate(value=case_when(.data$agg_type=="additive" ~ coalesce(.data$value,0) + .data$weight*(.data$p_value-.data$s_value), is.na(.data$value) ~ .data$p_value, @@ -422,7 +477,7 @@ proportional_reaggregate <- function(data,parent_data,geo_match,categories,base= select(-any_of(id)) } -#' Generate togfen correspondence for two geographies +#' Generate tongfen correspondence for two geographies #' #' @description #' \lifecycle{maturing} @@ -465,42 +520,35 @@ estimate_tongfen_single_correspondence <- function(geo1,geo2,geo1_uid,geo2_uid, if (robust) { - if (!st_is_valid(geo1)) geo1 <- geo1 %>% st_make_valid() - if (!st_is_valid(geo2)) geo2 <- geo2 %>% st_make_valid() + if (!isTRUE(all(st_is_valid(geo1)))) geo1 <- geo1 %>% st_make_valid() + if (!isTRUE(all(st_is_valid(geo2)))) geo2 <- geo2 %>% st_make_valid() } - robust_tolerance_buffer <- function(geo,geo_uid,tolerance,max_tries=20) { + # works on the geometries directly, the geometry column can go by any name + robust_tolerance_buffer <- function(geo,tolerance,max_tries=20) { t <- tolerance - d <- geo - d$geometry=st_buffer(geo$geometry,-t) + geometry <- st_geometry(geo) + buffered <- st_buffer(geometry,-t) count=0 - empties <- st_is_empty(d) + empties <- st_is_empty(buffered) while (sum(empties) > 0 & count 0) { stop("Unable to match within given tolerance, some geographies are too fine.") } - d + st_set_geometry(geo,buffered) } - cgeo1 <- geo1 %>% robust_tolerance_buffer(geo_uid = geo1_uid,tolerance = tolerance) - cgeo2 <- geo2 %>% robust_tolerance_buffer(geo_uid = geo2_uid,tolerance = tolerance) - - # Both intersections are necessary (buffered cgeo1 vs geo2, and cgeo2 vs geo1). The - # sparse index list is turned into a tibble directly, `as.data.frame()` on an empty - # result drops the columns we join on. - intersects_pairs <- function(x, y) { - m <- st_intersects(x, y, sparse = TRUE) - tibble(row.id = rep(seq_along(m), lengths(m)), - col.id = as.integer(unlist(m))) - } + cgeo1 <- geo1 %>% robust_tolerance_buffer(tolerance = tolerance) + cgeo2 <- geo2 %>% robust_tolerance_buffer(tolerance = tolerance) - i1 <- intersects_pairs(cgeo1, geo2) %>% + # Both intersections are necessary (buffered cgeo1 vs geo2, and cgeo2 vs geo1). + i1 <-intersects_pairs(cgeo1, geo2) %>% left_join(id1, by = c("row.id" = "id1")) %>% left_join(id2, by = c("col.id" = "id2")) %>% select(-"row.id",-"col.id") @@ -518,7 +566,7 @@ estimate_tongfen_single_correspondence <- function(geo1,geo2,geo1_uid,geo2_uid, correspondence } -#' Generate togfen correspondence for list of geographies +#' Generate tongfen correspondence for list of geographies #' #' @description #' \lifecycle{maturing} @@ -614,7 +662,7 @@ estimate_tongfen_correspondence <- function(data, -#' Check geographic integrety +#' Check geographic integrity #' #' @description #' \lifecycle{maturing} @@ -626,7 +674,7 @@ estimate_tongfen_correspondence <- function(data, #' 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 data a list of geographic 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) #' returned by `estimate_tongfen_correspondence`. diff --git a/R/tongfen_anomalies.R b/R/tongfen_anomalies.R new file mode 100644 index 0000000..e4fcd48 --- /dev/null +++ b/R/tongfen_anomalies.R @@ -0,0 +1,563 @@ +# Geocoding of census data varies over time, the same dwelling units can get assigned to +# different, neighbouring, regions in different years. In a timeline on a common geography +# this shows up as a surprising drop in one region that is offset by a jump in a neighbouring +# region. The functions in this file detect such patterns and join the affected regions, +# extending the idea behind TongFen from changing boundaries to data that got assigned across +# boundaries. See https://doodles.mountainmath.ca/posts/2024-07-26-geocoding-errors-in-aggregate-data/ +# for background, the implementation follows the method described there. + +anomaly_params <- function(rel_scale,abs_scale,p,surprise_cutoff,total_surprise_cutoff, + cutoff_fact,surprise_reduction_const,sum_fact) { + params <- list(rel_scale=rel_scale,abs_scale=abs_scale,p=p, + surprise_cutoff=surprise_cutoff,total_surprise_cutoff=total_surprise_cutoff, + cutoff_fact=cutoff_fact,surprise_reduction_const=surprise_reduction_const, + sum_fact=sum_fact) + valid <- vapply(params,\(x) is.numeric(x) && length(x)==1 && !is.na(x),logical(1)) + assert(all(valid),paste0("Need a single number for ",paste0(names(params)[!valid],collapse=", "),".")) + assert(rel_scale>0 && abs_scale>0 && p>0,"rel_scale, abs_scale and p have to be positive.") + params +} + +# Surprise of a change in a count, on a scale from 0 to 1. Only decreases are surprising, a +# relative decrease of `rel_scale` and an absolute decrease of `abs_scale` each are half way +# to full surprise, and it takes both for a change to be surprising. The relative change is +# taken with respect to `base`, which does not have to be the count the change started from. +anomaly_surprise <- function(change,base,rel_scale,abs_scale) { + decrease_surprise <- function(x,scale) { + r <- x + r[] <- 0 + decrease <- is.finite(x) & x<=0 + r[decrease] <- 1-0.5^(-x[decrease]/scale) + r + } + decrease_surprise(change/base,rel_scale)*decrease_surprise(change,abs_scale) +} + +# One round of looking for regions to join. `V` is a matrix of counts with a row per region +# and a column per year, `edges` a two column matrix with the row indices of neighbouring +# regions. Returns the candidate regions, most surprising first, with the neighbour that +# takes away most of the surprise when joined and if that is enough to join them. +anomaly_round <- function(V,edges,params) { + nt <- ncol(V) + change <- V[,-1,drop=FALSE]-V[,-nt,drop=FALSE] + base <- V[,-nt,drop=FALSE] + count_surprise <- \(S) rowSums(S>params$surprise_cutoff) + total_surprise <- \(S) rowSums(S^params$p)^(1/params$p) + + S <- anomaly_surprise(change,base,params$rel_scale,params$abs_scale) + count <- count_surprise(S) + total <- total_surprise(S) + + index <- which(count>0 & total>params$total_surprise_cutoff) + index <- index[order(-count[index],-total[index],-index)] + result <- tibble(index=index, + surprise_count=as.integer(count[index]), + surprise_total=total[index], + period=max.col(S[index,,drop=FALSE],ties.method="first"), + neighbour=NA_integer_, + surprise_total_joined=NA_real_, + join=FALSE) + + from <- c(edges[,1],edges[,2]) + to <- c(edges[,2],edges[,1]) + keep <- from %in% index + from <- from[keep] + to <- to[keep] + if (length(from)==0) return(result) + + # Joining a region with flat counts lowers the surprise just by growing the denominator + # of the relative change. To keep the surprise comparable the change of the joined region + # is measured against the counts of the candidate region alone. + # A neighbour with unknown change can't make up for a surprising change. + change_neighbour <- change[to,,drop=FALSE] + change_neighbour[is.na(change_neighbour)] <- 0 + S_joined <- anomaly_surprise(change[from,,drop=FALSE]+change_neighbour, + base[from,,drop=FALSE],params$rel_scale,params$abs_scale) + count_joined <- count_surprise(S_joined) + total_joined <- total_surprise(S_joined) + + o <- order(from,count_joined,total_joined,to) + best <- o[!duplicated(from[o])] + best <- best[match(result$index,from[best])] + + result$neighbour <- to[best] + result$surprise_total_joined <- total_joined[best] + join <- result$surprise_total_joined < params$cutoff_fact*result$surprise_total | + result$surprise_total-result$surprise_total_joined > params$surprise_reduction_const | + result$surprise_total_joined < params$cutoff_fact*params$sum_fact* + (result$surprise_total+total[result$neighbour]) + result$join <- !is.na(join) & join + result +} + +# Joins regions until there are no more pairs of regions left to join. Returns the group +# each region ended up in and the round in which it first got joined to another region. +anomaly_joins <- function(V,edges,params) { + membership <- seq_len(nrow(V)) + round_joined <- rep(NA_integer_,nrow(V)) + round <- 0L + repeat { + candidates <- anomaly_round(V,edges,params) + candidates <- candidates[candidates$join,] + # regions can only be part of one join per round, the most surprising ones go first + matched <- logical(nrow(V)) + target <- seq_len(nrow(V)) + for (r in seq_len(nrow(candidates))) { + i <- candidates$index[r] + j <- candidates$neighbour[r] + if (matched[i] || matched[j]) next + matched[c(i,j)] <- TRUE + target[max(i,j)] <- min(i,j) + } + if (!any(matched)) break + + round <- round+1L + round_joined[matched[membership] & is.na(round_joined)] <- round + + # contract the joined regions, the neighbours of a joined region are the neighbours + # of the regions it is made up of + target <- match(target,sort(unique(target))) + membership <- target[membership] + V <- rowsum(V,target) + edges <- cbind(target[edges[,1]],target[edges[,2]]) + edges <- edges[edges[,1]!=edges[,2],,drop=FALSE] + edges <- unique(cbind(pmin(edges[,1],edges[,2]),pmax(edges[,1],edges[,2]))) + } + list(membership=membership,round=round_joined) +} + +# Neighbouring regions as two column matrix of row indices into `ids`, each pair of +# neighbours is listed once. +anomaly_neighbours <- function(data,ids,neighbours) { + if (is.null(neighbours)) { + assert("sf" %in% class(data), + "Need data of class sf to determine neighbouring regions, alternatively specify the neighbours.") + geometry <- sf::st_geometry(data) + pairs <- intersects_pairs(geometry,geometry) + from <- pairs$row.id + to <- pairs$col.id + } else if (is.data.frame(neighbours)) { + assert(ncol(neighbours)>=2, + "Need neighbours to have two columns with the identifiers of neighbouring regions.") + from <- match(as.character(neighbours[[1]]),ids) + to <- match(as.character(neighbours[[2]]),ids) + } else if (is.list(neighbours)) { + # neighbours list like the ones from the spdep package, regions without neighbours hold a 0 + region_ids <- attr(neighbours,"region.id") %||% names(neighbours) %||% ids + assert(length(region_ids)==length(neighbours), + "Need neighbours to have an entry for each region.") + from <- rep(seq_along(neighbours),lengths(neighbours)) + to <- as.integer(unlist(neighbours)) + keep <- !is.na(to) & to>0 + from <- match(as.character(region_ids)[from[keep]],ids) + to <- match(as.character(region_ids)[to[keep]],ids) + } else { + stop("Don't know how to interpret neighbours, need a table with pairs of identifiers or a neighbours list.") + } + keep <- !is.na(from) & !is.na(to) & from!=to + unique(cbind(pmin(from[keep],to[keep]),pmax(from[keep],to[keep]))) +} + +anomaly_input <- function(data,variables,id,neighbours) { + assert(is.character(variables) && length(variables)>=2, + "Need at least two variables making up the timeline to look for anomalies.") + assert(is.character(id) && length(id)==1,"Need id to be the name of the identifier column.") + missing_columns <- setdiff(c(id,variables),names(data)) + assert(length(missing_columns)==0, + paste0("Did not find ",paste0(missing_columns,collapse=", ")," in data.")) + d <- data %>% ungroup() + if ("sf" %in% class(d)) d <- d %>% sf::st_drop_geometry() + not_numeric <- variables[!vapply(variables,\(v) is.numeric(d[[v]]),logical(1))] + assert(length(not_numeric)==0, + paste0("Variables have to be numeric, got ",paste0(not_numeric,collapse=", "),".")) + id_values <- d[[id]] + ids <- as.character(id_values) + assert(!anyNA(ids) && !anyDuplicated(ids), + paste0("Need ",id," to uniquely identify the regions in data.")) + + V <- matrix(as.numeric(unlist(d[variables],use.names=FALSE)),ncol=length(variables)) + edges <- anomaly_neighbours(data,ids,neighbours) + + # Ties are broken by the order of the regions, sort by identifier so that the result + # does not depend on the order the regions come in. + o <- order(ids,method="radix") + position <- order(o) + edges <- cbind(position[edges[,1]],position[edges[,2]]) + list(V=V[o,,drop=FALSE], + edges=cbind(pmin(edges[,1],edges[,2]),pmax(edges[,1],edges[,2])), + ids=ids[o], + id_values=id_values[o]) +} + + +#' Detect likely geocoding anomalies in timelines on a common geography +#' +#' @description +#' \lifecycle{experimental} +#' +#' TongFen is only as good as the geocoding that assigned the underlying data to geographic regions +#' in the first place. Geocoding varies over time, and the same dwelling units, and the people living +#' in them, can get assigned to different neighbouring regions in different years. In a timeline +#' on a common geography this shows up as a surprising drop in one region that is offset by +#' a corresponding jump in a neighbouring region. +#' +#' This function lists the candidate regions with surprising drops in the given count variable, +#' together with the neighbouring region that takes away most of the surprise when both are joined. +#' Use it to check for possible problems and to calibrate the parameters before joining regions with +#' `tongfen_anomaly_joins`. Not all surprising drops are due to geocoding problems, a drop that is +#' not complemented by a neighbouring region is likely real. +#' +#' The surprise of a change between two consecutive years ranges from 0 to 1. Only decreases are surprising, +#' the surprise is the product of the surprise of the relative and of the absolute decrease, so that it takes +#' a decrease that is large in both relative and absolute terms to be surprising. The total surprise of a +#' region is the `p`-norm of the surprises across all changes in the timeline. +#' +#' A candidate region and its neighbour are flagged for joining if joining reduces the total surprise of the +#' candidate region to below `cutoff_fact` times its total surprise, or reduces it by more than +#' `surprise_reduction_const`, or reduces it to below `cutoff_fact * sum_fact` times the sum of +#' the total surprises of both regions. To keep this comparable the surprise after joining is computed +#' from the change of the joined regions relative to the counts of the candidate region alone. +#' +#' @param data data on a common geography, with one row per region, for example as returned by +#' `tongfen_aggregate` or `get_tongfen_ca_census`. Needs to be of class sf unless `neighbours` is specified. +#' @param variables names of the columns holding the timeline of a count variable like population or +#' dwellings, in temporal order. Changes from or to a missing value are not surprising and don't make up for +#' surprising changes in neighbouring regions, and joined regions are missing a value if one of the regions they +#' are made up of is. Replace missing values by zero beforehand if they stand for regions where nothing got counted +#' @param id name of the column that uniquely identifies the regions, default is "TongfenID" +#' @param neighbours optional, neighbouring regions as a table with the identifiers of pairs of neighbouring +#' regions in the first two columns, or as a neighbours list like the ones returned by `spdep::poly2nb`. +#' By default all regions with intersecting geometries are neighbours, which can miss neighbours +#' if the geometries have been simplified and don't share their boundaries any more. +#' @param rel_scale relative decrease that is half way to full surprise, default is `0.25` for a 25\% drop +#' @param abs_scale absolute decrease that is half way to full surprise, default is `200` +#' @param p exponent of the norm used to combine the surprises across the timeline into the total surprise. +#' Large values focus on the most surprising change, 1 adds up the surprises of all changes, default is `4` +#' @param surprise_cutoff changes with larger surprise count as surprising, default is `0.15`. Only +#' regions with at least one surprising change are candidates +#' @param total_surprise_cutoff only regions with larger total surprise are candidates, default is `0.75` +#' @param cutoff_fact join regions if the total surprise after joining is lower than this share of the +#' total surprise of the candidate region, default is `0.6` +#' @param surprise_reduction_const join regions if joining lowers the total surprise by more than this, +#' default is `0.15` +#' @param sum_fact join regions if the total surprise after joining is lower than `cutoff_fact * sum_fact` +#' times the sum of the total surprises of both regions, default is `0.7` +#' @return A tibble with one row for each candidate region, most surprising first, with the identifier of the +#' region, the number of surprising changes `surprise_count`, the total surprise `surprise_total`, +#' the `period` with the most surprising change, the identifier of the `neighbour` that takes away most of +#' the surprise, the total surprise `surprise_total_joined` after joining both and `join` indicating if both +#' regions qualify to get joined. +#' @export +#' +#' @examples +#' # Check 2001 through 2021 dissemination area level population timelines in the +#' # City of Vancouver for possible geocoding problems +#' \dontrun{ +#' datasets <- c("CA01","CA06","CA11","CA16","CA21") +#' meta <- meta_for_additive_variables(datasets,"Population") +#' data <- get_tongfen_ca_census(regions=list(CSD="5915022"),meta=meta,level="DA",base_geo="CA21") +#' +#' anomalies <- tongfen_detect_anomalies(data,paste0("Population_",datasets)) +#' } +tongfen_detect_anomalies <- function(data,variables,id="TongfenID",neighbours=NULL, + rel_scale=0.25,abs_scale=200,p=4, + surprise_cutoff=0.15,total_surprise_cutoff=0.75, + cutoff_fact=0.6,surprise_reduction_const=0.15,sum_fact=0.7) { + params <- anomaly_params(rel_scale,abs_scale,p,surprise_cutoff,total_surprise_cutoff, + cutoff_fact,surprise_reduction_const,sum_fact) + input <- anomaly_input(data,variables,id,neighbours) + candidates <- anomaly_round(input$V,input$edges,params) + periods <- paste0(variables[-length(variables)],"-",variables[-1]) + + tibble(!!id:=input$id_values[candidates$index], + surprise_count=candidates$surprise_count, + surprise_total=candidates$surprise_total, + period=periods[candidates$period], + neighbour=input$id_values[candidates$neighbour], + surprise_total_joined=candidates$surprise_total_joined, + join=candidates$join) +} + + +#' Determine regions to join to correct for likely geocoding anomalies +#' +#' @description +#' \lifecycle{experimental} +#' +#' Looks for regions with surprising drops in the timeline of a count variable that are complemented by +#' a neighbouring region, as explained in `tongfen_detect_anomalies`, and joins them. This gets +#' repeated on the joined regions until there are no more regions left that qualify to get joined. +#' In each round a region only gets joined with one other region, the most surprising regions go first. +#' +#' Joining regions trades geographic detail for timelines that are consistent over time. The parameters +#' control how aggressively regions get joined and are best calibrated on the data at hand, erring on the side of +#' joining too few regions risks keeping geocoding problems, erring on the other side risks removing real +#' changes and needlessly coarsens the geography. +#' +#' The result can be used to join the regions via `tongfen_join_regions`, or to update a correspondence +#' via `tongfen_join_correspondence`. +#' +#' @inheritParams tongfen_detect_anomalies +#' @return A tibble with one row for each region that gets joined with other regions, with the identifier +#' of the region, the identifier of the joined region it becomes part of in the column named like the +#' identifier with suffix `_joined`, by default `TongfenID_joined`, and the `round` in which the region +#' first got joined to another region. The identifier of a joined region is the smallest +#' identifier of the regions it is made up of. +#' @export +#' +#' @examples +#' # Correct 2001 through 2021 dissemination area level population timelines in the +#' # City of Vancouver for likely geocoding problems +#' \dontrun{ +#' datasets <- c("CA01","CA06","CA11","CA16","CA21") +#' meta <- meta_for_additive_variables(datasets,"Population") +#' data <- get_tongfen_ca_census(regions=list(CSD="5915022"),meta=meta,level="DA",base_geo="CA21") +#' +#' joins <- tongfen_anomaly_joins(data,paste0("Population_",datasets)) +#' corrected_data <- tongfen_join_regions(data,joins,meta) +#' } +tongfen_anomaly_joins <- function(data,variables,id="TongfenID",neighbours=NULL, + rel_scale=0.25,abs_scale=200,p=4, + surprise_cutoff=0.15,total_surprise_cutoff=0.75, + cutoff_fact=0.6,surprise_reduction_const=0.15,sum_fact=0.7) { + params <- anomaly_params(rel_scale,abs_scale,p,surprise_cutoff,total_surprise_cutoff, + cutoff_fact,surprise_reduction_const,sum_fact) + input <- anomaly_input(data,variables,id,neighbours) + result <- anomaly_joins(input$V,input$edges,params) + + membership <- result$membership + joined <- tabulate(membership)[membership]>1 + # region with the smallest identifier in each group + o <- order(membership,input$ids,method="radix") + first <- !duplicated(membership[o]) + smallest <- integer(max(membership,0L)) + smallest[membership[o][first]] <- o[first] + + joins <- tibble(!!id:=input$id_values[joined], + !!paste0(id,"_joined"):=input$id_values[smallest[membership[joined]]], + round=result$round[joined]) + joins[order(input$ids[smallest[membership[joined]]],input$ids[joined],method="radix"),] +} + + +# Combine the TongfenUIDs of regions that get joined. A TongfenUID lists the identifiers +# making up a region as ":, :". +merge_tongfen_uids <- function(uids) { + uids <- unique(uids[!is.na(uids)]) + # the result must not depend on the order the regions come in + uids <- uids[order(uids,method="radix")] + parts <- unlist(strsplit(uids," ",fixed=TRUE)) + parts <- parts[nzchar(parts)] + split <- regexpr(":",parts,fixed=TRUE) + # not in the expected format, just string them together + if (length(parts)==0 || any(split<2)) return(paste0(uids,collapse=" ")) + columns <- substr(parts,1,split-1) + values <- strsplit(substring(parts,split+1),",",fixed=TRUE) + vapply(unique(columns),function(column) { + v <- unique(unlist(values[columns==column])) + paste0(column,":",paste0(v[order(v,method="radix")],collapse=",")) + },character(1)) %>% + paste0(collapse=" ") +} + +# metadata for aggregating data that already has been aggregated to a common geography, +# where variables are named by their label +meta_for_joining_regions <- function(data,meta) { + meta <- cut_meta(data,meta) + # the same variable can be part of several datasets, look up parents within the dataset + dataset <- if ("geo_dataset" %in% names(meta)) as.character(meta$geo_dataset) else rep("",nrow(meta)) + key <- paste(dataset,meta$variable,sep="\x1f") + parent <- meta$data_var[match(paste(dataset,meta$parent,sep="\x1f"),key)] + + needs_parent <- meta$rule %in% c("Average","Median","AverageTo") & is.na(parent) + if (any(needs_parent)) + stop(paste0("Can't join regions for ",paste0(unique(meta$data_var[needs_parent]),collapse=", "), + ", the data does not have the parent variables needed to aggregate them. ", + "Use `tongfen_join_correspondence` and `tongfen_aggregate` on the original data instead."), + call.=FALSE) + + meta %>% + mutate(variable=.data$data_var,parent=!!parent) %>% + select(any_of(c("variable","rule","parent","units","type"))) %>% + unique() +} + + +#' Join regions in data on a common geography +#' +#' @description +#' \lifecycle{experimental} +#' +#' Joins regions in data that has already been aggregated to a common geography, for example to correct +#' for likely geocoding anomalies as determined by `tongfen_anomaly_joins`. The data, and the geometries if +#' the data is of class sf, of the regions that get joined are aggregated, all other regions are left as they are. +#' +#' Variables are aggregated according to the metadata, numeric variables that are not part of the metadata +#' are assumed to be additive. Variables that are not additive, like averages, can only be aggregated if their +#' parent variable is part of the data. If that is not the case use `tongfen_join_correspondence` to update the +#' correspondence the data was built from and aggregate the original data again with `tongfen_aggregate`. +#' +#' @param data data on a common geography, with one row per region, for example as returned by +#' `tongfen_aggregate` or `get_tongfen_ca_census` +#' @param joins table with the regions to join as returned by `tongfen_anomaly_joins`, with the identifier +#' of the region and the identifier of the joined region it becomes part of in the column named like the +#' identifier with suffix `_joined` +#' @param meta optional metadata containing aggregation rules as for example returned by `meta_for_ca_census_vectors`, +#' variables are matched by their label. Numeric variables that are not part of the metadata are treated as +#' additive, if `NULL` (the default) that is the case for all numeric variables +#' @param id name of the column that uniquely identifies the regions, default is "TongfenID" +#' @param na.rm logical, determines how NA values should be treated when aggregating variables, +#' default is `TRUE` +#' @return The data with the regions joined. Joined regions take the place and the identifier of the +#' region with the smallest identifier among the regions they are made up of. Variables that are +#' not numeric and not part of the metadata are `NA` for joined regions. +#' @export +#' +#' @examples +#' # Correct 2001 through 2021 dissemination area level population timelines in the +#' # City of Vancouver for likely geocoding problems +#' \dontrun{ +#' datasets <- c("CA01","CA06","CA11","CA16","CA21") +#' meta <- meta_for_additive_variables(datasets,"Population") +#' data <- get_tongfen_ca_census(regions=list(CSD="5915022"),meta=meta,level="DA",base_geo="CA21") +#' +#' joins <- tongfen_anomaly_joins(data,paste0("Population_",datasets)) +#' corrected_data <- tongfen_join_regions(data,joins,meta) +#' } +tongfen_join_regions <- function(data,joins,meta=NULL,id="TongfenID",na.rm=TRUE) { + joined_id <- paste0(id,"_joined") + assert(id %in% names(data),paste0("Did not find ",id," in data.")) + assert(all(c(id,joined_id) %in% names(joins)), + paste0("Need joins to have columns ",id," and ",joined_id,".")) + data <- data %>% ungroup() + ids <- as.character(data[[id]]) + assert(!anyNA(ids) && !anyDuplicated(ids), + paste0("Need ",id," to uniquely identify the regions in data.")) + + match_index <- match(ids,as.character(joins[[id]])) + listed <- !is.na(match_index) + if (!any(listed)) return(data) + + new_id_values <- data[[id]] + new_id_values[listed] <- joins[[joined_id]][match_index[listed]] + # regions that other regions get joined to are part of the join, even if not listed themselves + affected <- listed | ids %in% as.character(new_id_values[listed]) + + is_sf <- "sf" %in% class(data) + geo_column <- if (is_sf) attr(data,"sf_column") else NULL + joined <- data[affected,] + joined[[id]] <- new_id_values[affected] + new_ids <- as.character(joined[[id]]) + + value_columns <- setdiff(names(data),c(id,"TongfenUID",geo_column)) + join_meta <- tibble(variable=character(0),rule=character(0),parent=character(0),type=character(0)) + if (!is.null(meta)) join_meta <- meta_for_joining_regions(joined,meta) + # results of tongfen calls can have count variables that are not part of the metadata, + # like the population, dwelling and household counts for Canadian census data + additive_columns <- setdiff(value_columns,join_meta$variable) + additive_columns <- additive_columns[vapply(additive_columns,\(v) is.numeric(joined[[v]]),logical(1))] + if (length(additive_columns)>0) { + message(paste0(ifelse(is.null(meta),"No metadata given, treating all numeric variables as additive: ", + "Treating numeric variables that are not part of the metadata as additive: "), + paste0(additive_columns,collapse=", "))) + join_meta <- bind_rows(join_meta, + tibble(variable=additive_columns,rule="Additive", + parent=NA_character_,type="Manual")) + } + dropped_columns <- setdiff(value_columns,join_meta$variable) + if (length(dropped_columns)>0) + message(paste0("Don't know how to aggregate ",paste0(dropped_columns,collapse=", "), + ", setting to NA for joined regions.")) + + aggregated <- joined %>% + select(all_of(c(id,intersect(join_meta$variable,names(joined)),geo_column))) %>% + group_by(!!as.name(id)) %>% + aggregate_data_with_meta(join_meta,na.rm=na.rm) %>% + ungroup() + + if ("TongfenUID" %in% names(data)) { + uids <- vapply(split(as.character(joined$TongfenUID),new_ids),merge_tongfen_uids,character(1)) + aggregated$TongfenUID <- unname(uids[as.character(aggregated[[id]])]) + } + + result <- bind_rows(data[!affected,],aggregated) + # joined regions take the place of the region they got their identifier from + result <- result[order(match(as.character(result[[id]]),ids)),names(data)] + if (is_sf) { + geometry_types <- unique(as.character(sf::st_geometry_type(result))) + if (length(geometry_types)>1 && all(geometry_types %in% c("POLYGON","MULTIPOLYGON"))) + result <- sf::st_cast(result,"MULTIPOLYGON") + } + result +} + + +#' Join regions in a correspondence +#' +#' @description +#' \lifecycle{experimental} +#' +#' Updates a correspondence so that the given regions are joined, for example to correct for likely geocoding +#' anomalies as determined by `tongfen_anomaly_joins`. The updated correspondence can be used in `tongfen_aggregate` +#' to aggregate data on the coarser common geography, which works for all variables `tongfen_aggregate` can +#' deal with and for data that was not part of detecting the anomalies. +#' +#' @param correspondence correspondence table with columns the unique geographic identifiers for each of the +#' geographies and the TongfenID and TongfenUID, as for example returned by `estimate_tongfen_correspondence` +#' or `get_tongfen_correspondence_ca_census` +#' @param joins table with the regions to join as returned by `tongfen_anomaly_joins`, with columns +#' `TongfenID` and `TongfenID_joined` +#' @return The correspondence with updated TongfenID and TongfenUID for the regions that got joined. If the +#' correspondence has a TongfenMethod column "anomaly" gets added to the method of the regions that got joined. +#' @export +#' +#' @examples +#' # Correct for likely geocoding problems in dissemination area level population timelines +#' # and use the updated correspondence to aggregate data on the corrected common geography +#' \dontrun{ +#' regions <- list(CSD="5915022") +#' datasets <- c("CA01","CA06","CA11","CA16","CA21") +#' meta <- meta_for_additive_variables(datasets,"Population") +#' data <- get_tongfen_ca_census(regions=regions,meta=meta,level="DA",base_geo="CA21") +#' joins <- tongfen_anomaly_joins(data,paste0("Population_",datasets)) +#' +#' correspondence <- get_tongfen_correspondence_ca_census(geo_datasets=datasets, +#' regions=regions,level="DA") %>% +#' tongfen_join_correspondence(joins) +#' } +tongfen_join_correspondence <- function(correspondence,joins) { + assert("TongfenID" %in% names(correspondence),"Did not find TongfenID in correspondence.") + assert(all(c("TongfenID","TongfenID_joined") %in% names(joins)), + "Need joins to have columns TongfenID and TongfenID_joined.") + ids <- as.character(correspondence$TongfenID) + missing_regions <- setdiff(as.character(joins$TongfenID),ids) + if (length(missing_regions)>0) + warning(paste0("Did not find ",length(missing_regions)," of the regions to join in the correspondence.")) + + match_index <- match(ids,as.character(joins$TongfenID)) + listed <- !is.na(match_index) + if (!any(listed)) return(correspondence) + + ids[listed] <- as.character(joins$TongfenID_joined)[match_index[listed]] + # regions that other regions get joined to are part of the join, even if not listed themselves + affected <- listed | ids %in% ids[listed] + new_ids <- ids[affected] + correspondence$TongfenID[affected] <- new_ids + if ("TongfenUID" %in% names(correspondence)) { + uids <- vapply(split(as.character(correspondence$TongfenUID[affected]),new_ids), + merge_tongfen_uids,character(1)) + correspondence$TongfenUID[affected] <- unname(uids[new_ids]) + } + if ("TongfenMethod" %in% names(correspondence)) { + method <- correspondence$TongfenMethod[affected] + tagged <- grepl("anomaly",method,fixed=TRUE) + method[!tagged] <- paste0(method[!tagged],", anomaly") + correspondence$TongfenMethod[affected] <- method + } + correspondence +} + +#' @import dplyr +#' @importFrom rlang .data +NULL +if(getRversion() >= "2.15.1") utils::globalVariables(c(".")) diff --git a/R/tongfen_ca.R b/R/tongfen_ca.R index 60caeb6..3d3d491 100644 --- a/R/tongfen_ca.R +++ b/R/tongfen_ca.R @@ -1,13 +1,9 @@ -correspondence_ca_census_urls <- list( - "2006"=list("DB"="https://www12.statcan.gc.ca/census-recensement/2011/geo/ref/files-fichiers/2006_92-156_DB_ID_txt.zip", - "DA"="https://www12.statcan.gc.ca/census-recensement/2011/geo/ref/files-fichiers/2006_92-156_DA_AD_txt.zip"), - "2011"=list("DB"="https://www12.statcan.gc.ca/census-recensement/2011/geo/ref/files-fichiers/2011_92-156_DB_ID_txt.zip", - "DA"="https://www12.statcan.gc.ca/census-recensement/2011/geo/ref/files-fichiers/2011_92-156_DA_AD_txt.zip"), - "2016"=list("DB"="https://www12.statcan.gc.ca/census-recensement/2011/geo/ref/files-fichiers/2016/2016_92-156_DB_ID_csv.zip", - "DA"="https://www12.statcan.gc.ca/census-recensement/2011/geo/ref/files-fichiers/2016/2016_92-156_DA_AD_csv.zip"), - "2021"=list("DB"="https://www12.statcan.gc.ca/census-recensement/2021/geo/aip-pia/correspondence-correspondance/files-fichiers/2021_92-156-X_DB_ID.zip", - "DA"="https://www12.statcan.gc.ca/census-recensement/2021/geo/aip-pia/correspondence-correspondance/files-fichiers/2021_92-156-X_DA_AD.zip") -) +# StatCan correspondence files as parquet, built by data-raw/statcan_correspondence.R +correspondence_ca_census_url <- function(year,level){ + base_url <- nullify_blank(getOption("tongfen.statcan_correspondence_url")) %||% + "https://mountainmath.s3.ca-central-1.amazonaws.com/tongfen/statcan_correspondence/v1" + paste0(sub("/+$","",base_url),"/statcan_correspondence_",year,"_",level,".parquet") +} ca_census_base <- c("Population","Dwellings","Households") @@ -30,15 +26,6 @@ datasets_from_vectors <- function(vs){ ds } -GEO_DATASET_LOOKUP <- c( - setNames(rep("CA1996",1),paste0("TX",seq(2000,2000))), - setNames(rep("CA01",5),paste0("TX",seq(2001,2005))), - setNames(rep("CA06",6),paste0("TX",seq(2006,2011))), - setNames(rep("CA11",4),paste0("TX",seq(2012,2015))), - setNames(rep("CA16",5),paste0("TX",seq(2016,2020))), - setNames(rep("CA16",21),paste0("CA",seq(2000,2020),"RMS")) -) - geo_dataset_for_years <- function(years){ require_suggested("cancensus") dataset_list <- cancensus::list_census_datasets() @@ -54,7 +41,6 @@ geo_dataset_for_years <- function(years){ geo_dataset_from_dataset <- function(datasets){ require_suggested("cancensus") - if (TRUE) { # legacy until cancensus updates datasets <- datasets %>% gsub("^CA11[NF]$","CA11",.) %>% gsub("\\d{4}x","",.) dataset_list <- cancensus::list_census_datasets() lapply(datasets, function(ds){ @@ -64,19 +50,9 @@ geo_dataset_from_dataset <- function(datasets){ unique() }) %>% unlist() - } else { - result <- tibble(dataset=datasets,geo_dataset=GEO_DATASET_LOOKUP[datasets]) %>% - mutate(geo_dataset=ifelse(is.na(.data$geo_dataset),.data$dataset %>% - years_from_datasets() %>% - as.character() %>% - substr(3,4) %>% - paste0("CA",.), - .data$geo_dataset)) - result$geo_dataset - } } -#' Generate metadata from Candian census vectors +#' Generate metadata from Canadian census vectors #' #' @description #' \lifecycle{maturing} @@ -92,7 +68,7 @@ geo_dataset_from_dataset <- function(datasets){ #' @examples #' # Build metadata for vectors #' \dontrun{ -#' meta <- meta_for_ca_census_vectors("v_CA16_4836","v_CA16_4838","v_CA16_4899") +#' meta <- meta_for_ca_census_vectors(c("v_CA16_4836","v_CA16_4838","v_CA16_4899")) #'} meta_for_ca_census_vectors <- function(vectors){ require_suggested("cancensus") @@ -104,9 +80,6 @@ meta_for_ca_census_vectors <- function(vectors){ nn[nn==""]=vectors[nn==""] } - if (length(vectors)==0) { - meta <- tibble::tibble(variable=NA,label=NA,dataset=datasets_from_vectors(vectors)) - } meta <- tibble::tibble(variable=vectors,label=nn,dataset=datasets_from_vectors(vectors)) %>% mutate(type="Original", aggregation="0",units=NA) datasets <- meta$dataset %>% @@ -138,13 +111,13 @@ meta_for_ca_census_vectors <- function(vectors){ select(variable="parent","dataset") %>% mutate(type="Extra",aggregation="Additive",rule="Additive") %>% filter(!is.na(.data$variable),!.data$variable %in% meta$variable) %>% - filter(!duplicated(.data$variable,.data$dataset)) %>% + distinct(.data$variable,.data$dataset,.keep_all=TRUE) %>% mutate(label=.data$variable) if (nrow(extras)>0) { meta <- meta %>% bind_rows(extras) %>% - filter(!duplicated(.data$variable,.data$dataset)) + distinct(.data$variable,.data$dataset,.keep_all=TRUE) } meta <- meta %>% @@ -155,13 +128,13 @@ meta_for_ca_census_vectors <- function(vectors){ -#' Generate metadata from Candian census vectors +#' Generate metadata from Canadian census vectors #' #' @description #' \lifecycle{maturing} #' #' Add Population, Dwellings, and Household counts to metadata -#' @param meta ribble with metadata as for example provided by `meta_for_ca_census_vectors` +#' @param meta tibble with metadata as for example provided by `meta_for_ca_census_vectors` #' @return tibble with metadata add_census_ca_base_variables <- function(meta){ new_meta <- meta$geo_dataset %>% @@ -184,6 +157,11 @@ add_census_ca_base_variables <- function(meta){ #' @description #' \lifecycle{maturing} #' +#' The correspondence files are downloaded from a mirror of the Statistics Canada correspondence files +#' and cached in the tongfen cache directory. The cached files are checked against the mirror once per session +#' and get downloaded again if they changed. The location of the mirror can be changed via the +#' `tongfen.statcan_correspondence_url` option. +#' #' @param year census year, only 2006 through 2021 are supported #' @param level geographic level, DA or DB #' @param refresh reload the correspondence files, default is `FALSE` @@ -193,36 +171,9 @@ get_single_correspondence_ca_census_for <- function(year,level=c("DA","DB"),refr year=as.character(year)[1] if (!(level %in% c("DA","DB"))) stop("Level needs to be DA or DB") if (!(year %in% c("2006","2011","2016","2021"))) stop("Year needs to be 2006, 2011, 2016, or 2021") - new_field=paste0(level,"UID",year) - old_field=paste0(level,"UID",as.integer(year)-5) - path=file.path(tongfen_cache_dir(),paste0("statcan_correspondence_",year,"_",level,".csv")) - if (refresh || !file.exists(path)) { - url=correspondence_ca_census_urls[[year]][[level]] - tmp=tempfile() - utils::download.file(url,tmp) - exdir=file.path(tempdir(),paste0("correspondence_",year,"_",level)) - if (dir.exists(exdir)) unlink(exdir,recursive=TRUE) - dir.create(exdir,showWarnings = FALSE) - utils::unzip(tmp,exdir=exdir) - file=dir(exdir,"\\.txt|\\.csv") - if (length(file)==0) { - p<-dir(exdir)[1] - if (dir.exists(file.path(exdir,p))) { - exdir=file.path(exdir,p) - file=dir(exdir,"\\.txt|\\.csv") - } - } - if (level=="DB") headers=c(new_field,old_field,"flag") else headers=c(new_field,old_field,paste0("DBUID",year),"flag") - unwanted <- paste0(level,"UID",year) - d<-readr::read_csv(file.path(exdir,file),col_types = readr::cols(.default = "c"),col_names = headers) %>% - select(all_of(c(new_field,old_field,"flag"))) %>% - unique() %>% - filter(!grepl(unwanted,!!as.name(new_field))) - readr::write_csv(d,path) - unlink(tmp) - unlink(exdir,recursive = TRUE) - } - result <- readr::read_csv(path,col_types = readr::cols(.default = "c")) + path=file.path(tongfen_cache_dir(),paste0("statcan_correspondence_",year,"_",level,".parquet")) + cached_download(correspondence_ca_census_url(year,level),path,refresh=refresh) + result <- tibble::as_tibble(nanoparquet::read_parquet(path)) # manual corrections if (year=="2021" && level=="DB") { @@ -243,7 +194,7 @@ get_single_correspondence_ca_census_for <- function(year,level=c("DA","DB"),refr #' @description #' \lifecycle{maturing} #' -#' Get correspondence file for several Candian censuses on a common geography. Requires sf and cancensus package to be available +#' Get correspondence file for several Canadian censuses on a common geography. Requires sf and cancensus package to be available #' #' @param regions census region list, should be inclusive list of GeoUIDs across censuses #' @param geo_datasets vector of census geography dataset identifiers @@ -253,7 +204,7 @@ 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 tolerance tolerance for `estimate_tongfen_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`) @@ -352,7 +303,7 @@ get_tongfen_correspondence_ca_census <- function(geo_datasets, regions, level="C correspondence_years=all_geo_years[-1] correspondence <- correspondence_years %>% lapply(function(year){ - c <- get_single_correspondence_ca_census_for(year,statcan_level) %>% + c <- get_single_correspondence_ca_census_for(year,statcan_level,refresh=refresh) %>% select(-"flag") previous_year <- all_geo_years[which(all_geo_years==year)-1] ds1 <- all_geo_datasets[all_geo_years==year] @@ -396,15 +347,15 @@ get_tongfen_correspondence_ca_census <- function(geo_datasets, regions, level="C } -#' Togfen data from several Canadian censuses +#' Tongfen data from several Canadian censuses #' #' @description #' \lifecycle{maturing} #' -#' Get data from several Candian censuses on a common geography. Requires sf and cancensus package to be available +#' Get data from several Canadian censuses on a common geography. Requires sf and cancensus package to be available #' #' @param regions census region list, should be inclusive list of GeoUIDs across censuses -#' @param meta metadata for the census veraiables to aggregate, for example as returned +#' @param meta metadata for the census variables to aggregate, for example as returned #' by \code{meta_for_ca_census_vectors}. #' @param level aggregation level to return data on (default is "CT") #' @param method tongfen method, options are "statcan" (the default), "estimate", "identifier". @@ -416,7 +367,7 @@ get_tongfen_correspondence_ca_census <- function(geo_datasets, regions, level="C #' 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, +#' @param tolerance tolerance for `estimate_tongfen_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`) diff --git a/R/tongfen_ca_deprecated.R b/R/tongfen_ca_deprecated.R index 418ff0e..f0cf13b 100644 --- a/R/tongfen_ca_deprecated.R +++ b/R/tongfen_ca_deprecated.R @@ -85,7 +85,7 @@ get_tongfen_census_da <- function(regions,vectors,geo_format=NA,use_cache=TRUE,n #' @return dataframe with variables on common geography #' @export 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()") + lifecycle::deprecate_warn("0.2.0", "get_tongfen_ca_census_ct_from_da()", "get_tongfen_ca_census()") meta <- meta_for_ca_census_vectors(vectors) base_geo <- if (is.na(geo_format)) NULL else meta$geo_dataset %>% unique() %>% sort() %>% first() @@ -106,11 +106,13 @@ get_tongfen_ca_census_ct_from_da <- function(regions,vectors,geo_format=NA,use_c #' \lifecycle{deprecated} #' #' Aggregate variables to common CTs, returns data2 on new tiling matching data1 geography -#' @param data1 cancensus CT level datatset for year1 < year2 to serve as base for common geography -#' @param data2 cancensus CT level datatset for year2 to be aggregated to common geography +#' @param data1 cancensus CT level dataset for year1 < year2 to serve as base for common geography +#' @param data2 cancensus CT level dataset for year2 to be aggregated to common geography #' @param data2_sum_vars vector of variable names to by summed up when aggregating geographies #' @param data2_group_vars optional vector of grouping variables #' @param na.rm optional parameter to remove NA values when summing, default = `TRUE` +#' @return `data2` with the variables in `data2_sum_vars` aggregated to a common geography matching `data1`, +#' identified by the `GeoUID` of `data1` #' @export tongfen_ca_census_ct <- function(data1,data2,data2_sum_vars,data2_group_vars=c(),na.rm=TRUE) { lifecycle::deprecate_warn("0.2.0", "tongfen_ca_census_ct()", "tongfen_aggregate()") @@ -159,7 +161,7 @@ tongfen_ca_census_ct <- function(data1,data2,data2_sum_vars,data2_group_vars=c() #' #' @description #' \lifecycle{deprecated} -#' Joins the StatCan correspodence files for several census years +#' Joins the StatCan correspondence files for several census years #' #' @param years list of census years #' @param level geographic level, DA or DB diff --git a/R/tongfen_ca_estimate.R b/R/tongfen_ca_estimate.R index 6fcaecb..bf24258 100644 --- a/R/tongfen_ca_estimate.R +++ b/R/tongfen_ca_estimate.R @@ -25,11 +25,12 @@ #' one of these for all variables. #' @param na.rm how to deal with NA values, default is \code{FALSE}. #' @param quiet suppress progress messages +#' @return `geometry` with the estimated values for the census variables specified by `meta` #' @export #' #' @examples -#' # Estimate a common geography for 2006 and 2016 dissemination areas in the City of Vancouver -#' # based on the geographic data and check estimation errors +#' # Estimate the 2016 population within 1 km of Toronto City Hall from dissemination area level +#' # census data #' \dontrun{ #' toronto_city_hall <- sf::st_point(c(-79.3839,43.6534)) %>% #' sf::st_sfc(crs=4326) %>% @@ -42,7 +43,7 @@ #' data <- tongfen_estimate_ca_census(toronto_city_hall,meta,level="DA",intersection_level="CT") #' #' print(paste0("Approximately ",scales::comma(data$Population,accuracy=100), -#' " people live within a 1 km radius of Toronto City.")) +#' " people live within a 1 km radius of Toronto City Hall.")) #' #'} tongfen_estimate_ca_census <- function(geometry, meta, level, @@ -85,7 +86,7 @@ tongfen_estimate_ca_census <- function(geometry, meta, level, census_data <- g } - result <- tongfen_estimate(target = geometry, source = census_data, meta = meta,na.rm = na.rm) + tongfen_estimate(target = geometry, source = census_data, meta = meta,na.rm = na.rm) } diff --git a/R/tongfen_estimate.R b/R/tongfen_estimate.R index 8342e84..75f63ce 100644 --- a/R/tongfen_estimate.R +++ b/R/tongfen_estimate.R @@ -13,7 +13,9 @@ #' @param meta metadata for variable aggregation, see `meta_for_additive_variables` and `meta_for_ca_census_vectors` for more information #' on how to construct metadata. #' @param na.rm remove NA values when aggregating, default is FALSE -#' @return `target` with estimated quantities from `source` as specified by `meta` +#' @return `target` with estimated quantities from `source` as specified by `meta`, regions in `target` +#' that don't overlap with `source` have `NA` values. Columns in `target` can't have the same name +#' as the variables to be estimated. #' @export #' #' @examples @@ -46,6 +48,12 @@ tongfen_estimate <- function(target,source,meta,na.rm=FALSE) { cut_meta(source,.) %>% mutate(var_name=paste0("v",row_number())) + clashes <- intersect(meta$data_var,names(target)) + if (length(clashes)>0) { + stop(paste0("Target already has columns named ",paste0(clashes,collapse=", "), + ", please rename or remove these before estimating.")) + } + # rename variables, st_interpolate_aw does not handle column names with special characters safe_rename_vars <- setNames(meta$data_var,meta$var_name) safe_rename_back <- setNames(meta$var_name,meta$data_var) @@ -54,30 +62,33 @@ tongfen_estimate <- function(target,source,meta,na.rm=FALSE) { i = suppressMessages(st_intersection(st_geometry(source), st_geometry(target))) idx = attr(i, "idx") - gc = which(st_is(i, "GEOMETRYCOLLECTION")) - i[gc] = st_collection_extract(i[gc], "POLYGON") + # the area only counts the polygonal parts of intersections, source regions that only + # touch a target region along a boundary don't overlap and must not contribute + area_st <- as.numeric(st_area(i)) + overlaps <- which(area_st > 0) + idx <- idx[overlaps,,drop=FALSE] + area_st <- area_st[overlaps] source <- source %>% rename(!!!safe_rename_vars) - source_area <- unclass(st_area(source)) + source_area <- as.numeric(st_area(source)) x_st <- source[idx[,1],, drop=FALSE] %>% + st_drop_geometry() %>% select(names(safe_rename_vars)) %>% - pre_scale(meta,meta_var = "var_name") %>% - mutate(...area_st = st_area(i) %>% unclass, - ...area_s = source_area[idx[,1]]) %>% - mutate(...factor = .data$...area_st/.data$...area_s) %>% - mutate(...partial = .data$...factor < 0.99) %>% - st_drop_geometry() + pre_scale(meta,meta_var = "var_name") - x_st[meta$var_name] <- lapply(x_st[meta$var_name], `*`, x_st$...factor) + # variables to sum up, including the weights for averages added when pre-scaling + sum_vars <- names(x_st) + x_st[sum_vars] <- lapply(x_st[sum_vars], `*`, area_st/source_area[idx[,1]]) - x_st <- stats::aggregate(x_st, list(idx[,2]), sum, na.rm=na.rm) + # target regions without overlap with the source don't show up here and end up as NA + x_st <- x_st %>% + mutate(!!unique_key:=idx[,2]) %>% + group_by(.data[[unique_key]]) %>% + summarize(across(all_of(sum_vars), \(x) sum(x,na.rm=na.rm)),.groups="drop") result <- target %>% - left_join(x_st %>% - select(-all_of(c("...factor", "...partial", "...area_s", "...area_st"))) %>% - rename(!!unique_key:="Group.1"), - by=unique_key) %>% + left_join(x_st,by=unique_key) %>% select(-all_of(unique_key)) %>% post_scale(meta,meta_var = "var_name") %>% rename(!!!safe_rename_back) @@ -117,17 +128,17 @@ tongfen_estimate <- function(target,source,meta,na.rm=FALSE) { #' @param target custom geography #' @param source input geography #' @param target_id name of the column in `target` table with unique id (character) -#' @return `source` with extra column with name `"target_id"` and column `...overlap_fraction` with -#' the proportion of overlap of the target geometry with the respective `target_id` +#' @return `source` with extra column with the name given by `target_id` and column `...overlap_fraction` with +#' the proportion of the area of the source region that overlaps with the region in `target` with that id #' @export #' #' @examples -#' # Estimate 2006 Populatino in the City of Vancouver dissemination ares on 2016 census geoographies +#' # Tag 2016 dissemination areas in the City of Vancouver by the 2006 census tract they overlap +#' # the most with #' \dontrun{ -#' geo1 <- cancensus::get_census("CA06",regions=list(CSD="5915022"),geo_format='sf',level='DA') +#' geo1 <- cancensus::get_census("CA06",regions=list(CSD="5915022"),geo_format='sf',level='CT') #' geo2 <- cancensus::get_census("CA16",regions=list(CSD="5915022"),geo_format='sf',level='DA') -#' meta <- meta_for_additive_variables("CA06","Population") -#' result <- tongfen_estimate(geo2 %>% rename(Population_2016=Population),geo1,meta) +#' result <- tongfen_tag_largest_overlap(geo2,geo1 %>% select(CT_2006=GeoUID),"CT_2006") #'} tongfen_tag_largest_overlap <- function(source, target, target_id) { target_geo_types <- target %>% sf::st_geometry_type() %>% unique diff --git a/R/tongfen_us.R b/R/tongfen_us.R index 19ecf35..7e9790a 100644 --- a/R/tongfen_us.R +++ b/R/tongfen_us.R @@ -59,14 +59,9 @@ get_us_ct_correspondence_path <- function(state,year){ # 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) - local_path <- file.path(cache_path,basename(path)) - if (!file.exists(local_path)) { - if (!dir.exists(cache_path)) dir.create(cache_path, recursive = TRUE) - utils::download.file(path,local_path,quiet = TRUE) - } + local_path <- file.path(us_cache_dir(cache_path),basename(path)) + cached_download(path,local_path,check_remote=FALSE) 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", @@ -97,13 +92,8 @@ get_us_ct_correspondence_2020 <- function(state,min_area_share=0.01, 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") - local_path <- file.path(cache_path,file) - if (!file.exists(local_path)) { - if (!dir.exists(cache_path)) dir.create(cache_path, recursive = TRUE) - utils::download.file(path,local_path,quiet=TRUE) - } + local_path <- file.path(us_cache_dir(cache_path),basename(path)) + cached_download(path,local_path,check_remote=FALSE) d<-readr::read_csv(local_path,progress=FALSE, col_names=c("STATE00","COUNTY00","TRACT00","GEOID00", "POP00","HU00","PART00","AREA00","AREALAND00", @@ -130,12 +120,8 @@ get_us_ct_correspondence_2010 <- function(state,min_area_share=0.01, 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, recursive = TRUE) - utils::download.file(path,local_path,quiet=TRUE) - } + local_path <- file.path(us_cache_dir(cache_path),basename(path)) + cached_download(path,local_path,check_remote=FALSE) 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), @@ -202,15 +188,21 @@ get_us_ct_correspondence <- function(state, datasets, min_area_share=0.01, # the 2000 to 2010 county subdivision comparability file, covering all states get_us_county_subdivision_correspondence <- function(cache_path=getOption("tongfen.cache_path")){ require_suggested("readxl") - cache_path = file.path(cache_path %||% tempdir(),"us_data") + cache_path = us_cache_dir(cache_path) file <- "Cousub_comparability.xlsx" local_path <- file.path(cache_path,file) if (!file.exists(local_path)) { if (!dir.exists(cache_path)) dir.create(cache_path, recursive = TRUE) - tmp=tempfile(fileext = ".zip") + # download and unpack in a temporary directory, an interrupted download must not leave + # a broken file in the cache + tmp_dir <- tempfile() + dir.create(tmp_dir) + on.exit(unlink(tmp_dir,recursive=TRUE)) + tmp <- file.path(tmp_dir,"cousub_comparabilityxls.zip") path="https://www2.census.gov/geo/docs/maps-data/data/comp/cousub_comparabilityxls.zip" - utils::download.file(path,tmp,quiet=TRUE) - utils::unzip(tmp,exdir = cache_path) + utils::download.file(path,tmp,mode="wb",quiet=TRUE) + utils::unzip(tmp,files=file,exdir=tmp_dir) + file.copy(file.path(tmp_dir,file),local_path,overwrite=TRUE) } readxl::read_xlsx(local_path) } @@ -223,14 +215,10 @@ get_us_county_subdivision_correspondence <- function(cache_path=getOption("tongf # 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, recursive = TRUE) - utils::download.file(path,local_path,quiet=TRUE) - } + local_path <- file.path(us_cache_dir(cache_path),basename(path)) + cached_download(path,local_path,check_remote=FALSE) 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", @@ -300,7 +288,8 @@ get_us_county_subdivision_correspondence_for <- function(state, datasets, min_ar #' 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 +#' `tongfen.cache_path` option. If that is not set the `tongfen.cache_path` environment variable +#' and the `custom_data_path` option are used, falling 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 @@ -351,7 +340,7 @@ sumfile_for_dataset <- function(sumfile, ds){ unname(sumfile[[ds]]) } -#' Get US census data for 2000 and 2010 census on common census tract based geography +#' Get US census data for several censuses on a common geography #' #' @description #' \lifecycle{maturing} @@ -369,14 +358,15 @@ sumfile_for_dataset <- function(sumfile, ds){ #' @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 base_geo dataset to use as base geography, for example `"dec2010"`, has to be one of +#' the datasets in `meta`. Default is `NULL`, which uses the first dataset in `meta`. #' @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 "_"). +#' @return sf object with (wide form) census variables with census year as suffix (separated by underscore "_"). #' @export #' #' @examples @@ -404,7 +394,7 @@ get_tongfen_us_census <- function(regions,meta,level='tract',survey="census", assert(base_geo %in% datasets,paste0("base_geo has to be one of the datasets ",paste0(datasets,collapse=", "))) invalid_datasets <- setdiff(datasets,names(valid_us_census_datasets)) 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(level %in% c('tract','county subdivision'),"Only census tracts and county subdivisions are supported right now.") assert(survey %in% c('census'),"Only census surveys are supported right now.") if (!is.null(sumfile)) { if (is.null(names(sumfile))) { diff --git a/README.md b/README.md index 7a1ad3d..d86e2eb 100644 --- a/README.md +++ b/README.md @@ -29,13 +29,13 @@ library(tongfen) ### Caching correspondence files The `get_tongfen_ca_census` and `get_tongfen_correspondence_ca_census` methods make use of the StatCan correspondence -files when run with `method = "statcan"`. To speed up this process it is useful to permanently cache these files instead of having to download them repeatedly. If caching is desired, set either +files when run with `method = "statcan"`. Statistics Canada no longer allows programmatic downloads of these files, so the package downloads them from a mirror hosting the files in parquet format. To speed up this process it is useful to permanently cache these files instead of having to download them again in every session. If caching is desired, set either * `options("tongfen.cache_path"="")` * `Sys.setenv("tongfen.cache_path"="")` * `options("custom_data_path"="")` -in your `.Rprofile` or `.Renviron` file. +in your `.Rprofile` or `.Renviron` file. Cached files are checked against the mirror once per session and only get downloaded again if they changed, when offline the cached files are used as they are. ## General TongFen @@ -52,7 +52,7 @@ A convenience function to validate geographic TongFen fit via area comparison is Finding a common tiling of several different yet congruent geographies is only one part of the problem TongFen addresses, aggregating up the variables is the other part. The `tongfen` package deals with this using a *metadata* table that specifies how variables should be aggregated. In it's simplest form values are simply added up. The `meta_for_additive_variables` convenience function builds the metadata for additive variables. Metadata for non-additive variables like averages, ratios or percentages needs more care to build, it requires additional information on the **parent variable** that specifies the denominator of the average, ratio or percentage. Other data, like medians, can't be aggregated up, although `tongfen` can provide estimates of medians on aggregated geographies by treating them as averages. ### Packaged data -The package ships with a subset of [voting data from Elections Canada](https://www.elections.ca/content.aspx?section=ele&dir=pas&document=index&lang=e) for the 42nd and 43rd federal elections as well as the polling district geographies for the [42nd](https://open.canada.ca/data/en/dataset/6a78ccfd-6bba-4109-b040-87cb8c71ec35) and [43rd](https://open.canada.ca/data/en/dataset/e70e3263-8584-4f22-94cb-8c15b616cbfc). This facilitates running the example vignette on polling districts without having to download external data. Both are available as open data covered under the [Open Government Licence - Canda](https://open.canada.ca/en/open-government-licence-canada). +The package ships with a subset of [voting data from Elections Canada](https://www.elections.ca/content.aspx?section=ele&dir=pas&document=index&lang=e) for the 42nd and 43rd federal elections as well as the polling district geographies for the [42nd](https://open.canada.ca/data/en/dataset/6a78ccfd-6bba-4109-b040-87cb8c71ec35) and [43rd](https://open.canada.ca/data/en/dataset/e70e3263-8584-4f22-94cb-8c15b616cbfc). This facilitates running the example vignette on polling districts without having to download external data. Both are available as open data covered under the [Open Government Licence - Canada](https://open.canada.ca/en/open-government-licence-canada). ## Data-specific implementations The need for TongFen comes up frequently with certain types of geographies. Census geographies is one such example. In some cases these data sources come with their own correspondence files that go beyond geographic matchup but also join regions to alleviate data integrity problems like geocoding issues. @@ -67,14 +67,15 @@ The package is well-integrated to work with Canadian census data in two essentia ### US census data -* `get_tongfen_us_census` integrates the data acquisition (via the [**tidycensus** package](https://walker-data.com/tidycensus/index.html)) with TongFen, and adds the tongfen `method = "census.gov"` to use the US Census Bureau correspondence files for matching. +* `get_tongfen_us_census` integrates the data acquisition (via the [**tidycensus** package](https://walker-data.com/tidycensus/index.html)) with TongFen, using the US Census Bureau relationship files to build the common geography. +* `get_tongfen_correspondence_us_census` breaks out the correspondence generation from the US Census Bureau relationship files, to tongfen data that comes on census geographies but is obtained by other means. The relationship files are cached in the `us_data` folder of the cache path described above. ## Other implementations The `tongfen` package is open to add extensions for other specialized data sources, as well as extensions of existing ones. ## Fixed target geography estimation -When geographies aren't sufficiently congruent or the target geography is fixed, we won't be able to use the `tongfen` methods to compute the data on a common geography but have to instead rely on estimates. The `tongfen_estimate` makes no assumption on the underlying geographies and returns estimates of the data on the target geography. It uses area-weighted interpolation to achieve this, and can be refined to dasymmetric estimates using the `proportional_reaggregate` function. +When geographies aren't sufficiently congruent or the target geography is fixed, we won't be able to use the `tongfen` methods to compute the data on a common geography but have to instead rely on estimates. The `tongfen_estimate` makes no assumption on the underlying geographies and returns estimates of the data on the target geography. It uses area-weighted interpolation to achieve this, and can be refined to dasymetric estimates using the `proportional_reaggregate` function. This method has the example that it works independent of the nature of the underlying geographies, but comes at the heavy price of only being an estimate. To be useful for research purposes we also need methods to estimate the errors this introduces and the effects this has on subsequent analysis results. @@ -84,8 +85,8 @@ Methods to facilitate this are still under active development. If you wish to cite tongfen: - von Bergmann, J. (2024). tongfen: R package to - Make Data Based on Different Geographies Comparable. v0.3.7. + von Bergmann, J. (2026). tongfen: R package to + Make Data Based on Different Geographies Comparable. v0.3.9. DOI: 10.32614/CRAN.package.tongfen @@ -94,9 +95,9 @@ A BibTeX entry for LaTeX users is @Manual{tongfen, author = {Jens {von Bergmann}}, title = {tongfen: R package to Make Data Based on Different Geographies Comparable}, - year = {2024}, + year = {2026}, doi = {10.32614/CRAN.package.tongfen}, - note = {R package version 0.3.7}, + note = {R package version 0.3.9}, url = {https://mountainmath.github.io/tongfen/}, } ``` diff --git a/cran-comments.md b/cran-comments.md index c6fa36b..298facd 100644 --- a/cran-comments.md +++ b/cran-comments.md @@ -1,63 +1,29 @@ -# 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 +# tongfen v.0.3.9 ## 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 +- new experimental functions `tongfen_detect_anomalies`, `tongfen_anomaly_joins`, `tongfen_join_regions` + and `tongfen_join_correspondence` to detect and correct for likely geocoding anomalies in timelines + on a common geography, together with a new vignette +- StatCan correspondence files are now downloaded as parquet files from a mirror, Statistics Canada + put the original files behind a browser check that blocks programmatic downloads, which broke + `method = "statcan"` ## 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 -- fix `proportional_reaggregate` ignoring all but the first base variable when `base` names a - different variable per category -- fix `estimate_tongfen_correspondence` with `method="identifier"` erroring out when every - geographic identifier matches -- packages in Suggests (`cancensus`, `tidycensus`, `readxl`) are now used conditionally, with an - actionable message when they are not installed -- faster `check_tongfen_areas` and `aggregate_correspondences` - -# tongfen v.0.3.7 -## Major changes -- accommodate factors in proportional_reaggregate -- sizable performance increases -- squish several edge case bugs - -# tongfen v.0.3.6 -## Major changs -- better downsampling that can also accommodate averages -- performance improvements -## Minor changes -- better documentation -- allow for datasets vartiables by census year for canadian data -- fix issue where some metadata might get duplicated - - -# Update v.0.3.3 -- Fix compatibility issue with changes in {sf} package -- More reliable GitHub action CRAN checks - -# Update v.0.3.2 -- Added `tongfen_estimate_ca_census` function for new CensusMapper endpoint, tying into new {cancensus} functionality. -- Custom impelementation of `tongfen_etimate` for finer control -- Fix compatibility issue with changes in {sf} package - -# Submission - v.0.3 +- combining correspondences across three or more datasets no longer runs into a cross join +- fix `proportional_reaggregate` giving wrong results when the finer level data already has values +- fix `tongfen_estimate` mishandling intersections that are geometry collections, regions that only + touch the target, and targets without overlap with the source +- fix averages getting scaled more than once in `tongfen_aggregate` when the metadata lists the same + variable name for several datasets +- fix averages with missing values being pulled toward zero when aggregating with `na.rm = TRUE`, + and "Average to" variables sharing a parent variable overwriting each other's base +- `estimate_tongfen_correspondence` no longer requires the geometry column to be named `geometry` +- `refresh = TRUE` now also refreshes the cached StatCan correspondence files +- US Census Bureau relationship files are downloaded to a temporary file before being moved to the cache +- added missing `\value` documentation for `tongfen_estimate_ca_census` and `tongfen_ca_census_ct` +- fixed typos in the documentation # Test environments * local macOS installation, R 4.6.0 -* GitHub actions (windows-latest, macOS-latest, ubuntu-latest) on release, devel and oldrel +* GitHub actions: macOS-latest (release), windows-latest (release), ubuntu-latest (devel, release, oldrel-1) # R CMD check results 0 errors | 0 warnings | 0 notes diff --git a/data-raw/statcan_correspondence.R b/data-raw/statcan_correspondence.R new file mode 100644 index 0000000..44e64d5 --- /dev/null +++ b/data-raw/statcan_correspondence.R @@ -0,0 +1,51 @@ +## code to prepare the StatCan DA and DB correspondence files hosted on S3 +## +## StatCan put the correspondence files behind a browser check, so they can't be downloaded +## programmatically any more. Download and extract the zip files manually from +## https://www12.statcan.gc.ca/census-recensement/2011/geo/ref/files-fichiers/2006_92-156_DB_ID_txt.zip +## https://www12.statcan.gc.ca/census-recensement/2011/geo/ref/files-fichiers/2006_92-156_DA_AD_txt.zip +## https://www12.statcan.gc.ca/census-recensement/2011/geo/ref/files-fichiers/2011_92-156_DB_ID_txt.zip +## https://www12.statcan.gc.ca/census-recensement/2011/geo/ref/files-fichiers/2011_92-156_DA_AD_txt.zip +## https://www12.statcan.gc.ca/census-recensement/2011/geo/ref/files-fichiers/2016/2016_92-156_DB_ID_csv.zip +## https://www12.statcan.gc.ca/census-recensement/2011/geo/ref/files-fichiers/2016/2016_92-156_DA_AD_csv.zip +## https://www12.statcan.gc.ca/census-recensement/2021/geo/aip-pia/correspondence-correspondance/files-fichiers/2021_92-156-X_DB_ID.zip +## https://www12.statcan.gc.ca/census-recensement/2021/geo/aip-pia/correspondence-correspondance/files-fichiers/2021_92-156-X_DA_AD.zip +## into `input_dir`, then run this script and upload the parquet files in `output_dir` to S3. + +library(dplyr) + +input_dir <- "~/Downloads" +output_dir <- "~/Downloads/tongfen_statcan_correspondence" + +level_codes <- c(DA="DA_AD",DB="DB_ID") + +prepare_statcan_correspondence <- function(year,level){ + new_field <- paste0(level,"UID",year) + old_field <- paste0(level,"UID",year-5) + dirs <- dir(input_dir,paste0("^",year,"_92-156.*",level_codes[[level]]),full.names=TRUE) + dirs <- dirs[dir.exists(dirs)] + file <- dir(dirs,"\\.(txt|csv)$",full.names=TRUE,recursive=TRUE) + stopifnot(length(file)==1) + # DA files are at DB granularity and carry the DBUID in the third column, 2021 files have extra DGUID columns + headers <- if (level=="DB") c(new_field,old_field,"flag") else c(new_field,old_field,paste0("DBUID",year),"flag") + d <- readr::read_csv(file,col_types=readr::cols(.default="c"),col_names=FALSE) %>% + select(all_of(seq_along(headers))) %>% + setNames(headers) %>% + filter(grepl("^\\d+$",.data[[new_field]])) %>% # header row in the 2016 and 2021 files + select(all_of(c(new_field,old_field,"flag"))) %>% + unique() %>% + arrange(.data[[new_field]],.data[[old_field]]) + stopifnot(all(grepl("^\\d+$",d[[old_field]])),all(d$flag %in% as.character(1:4))) + d +} + +dir.create(output_dir,showWarnings=FALSE) +for (year in c(2006,2011,2016,2021)) { + for (level in c("DA","DB")) { + prepare_statcan_correspondence(year,level) %>% + nanoparquet::write_parquet(file.path(output_dir,paste0("statcan_correspondence_",year,"_",level,".parquet")), + compression="zstd", + # nanoparquet does not compress at its default zstd level + options=nanoparquet::parquet_options(compression_level=19)) + } +} diff --git a/docs/404.html b/docs/404.html index d0348fd..008e97f 100644 --- a/docs/404.html +++ b/docs/404.html @@ -30,7 +30,7 @@ tongfen - 0.3.8 + 0.3.9 + + + + + +
+ + + + +
+
+

Geocoding anomalies in TongFen timelines

+

Jens von +Bergmann

+ +

2026-10-04

+ + Source: vignettes/tongfen_anomalies.Rmd +
tongfen_anomalies.Rmd
+
+ + + +

TongFen makes data on different geographies comparable by aggregating +it up to a common geography. The result is only as good as the geocoding +that assigned the underlying data to geographic regions in the first +place. Geocoding changes over time, and the same dwelling units, and the +people living in them, can get assigned to different neighbouring +regions in different years. In a timeline on a common geography this +shows up as a surprising drop in one region that is offset by a jump in +a neighbouring region.

+

The fix is in the spirit of TongFen, joining the affected regions +gives a slightly coarser geography on which the data is consistent over +time. The functions in this vignette look for such patterns and join the +regions on demand. The method is explained in more detail in a blog +post, this vignette follows the example from that post.

+
+library(dplyr)
+library(tidyr)
+library(ggplot2)
+library(cancensus)
+library(sf)
+library(tongfen)
+# cancensus::set_api_key("<your cancensus API key>")
+
+

Population timelines for Toronto +

+

As an example we take the population from the 1971 through 2011 +censuses that Statistics Canada tabulated on 2016 dissemination areas, +together with the 2016 population. All data comes on the same geography, +so there is no need to TongFen, but the data for the earlier years is +geocoded from the road network and block face of the time, which does +not always match up with 2016 dissemination areas.

+
+years <- c(1971,seq(1981,2011,5))
+vectors <- c(setNames(paste0("v_CA",years,"x16_1"),years),"2016"="v_CA16_1")
+timeline <- names(vectors)
+
+toronto <- get_census("CA16CT",regions=list(CSD="3520005"),vectors=vectors,
+                      level="DA",geo_format="sf",quiet=TRUE) %>%
+  select(GeoUID,all_of(timeline)) %>%
+  mutate(across(all_of(timeline),\(x) coalesce(x,0)))
+

Dissemination areas without population in a given year come back as +missing values. Changes from or to a missing value are never considered +surprising, so we set them to zero to mark these as areas where nobody +got counted.

+

The area around Crescent Town shows what the problem looks like.

+
+crescent_town <- c("35204370","35204765")
+
+plot_timelines <- function(data) {
+  data %>%
+    st_drop_geometry() %>%
+    pivot_longer(all_of(timeline),names_to="Year",values_to="Population") %>%
+    ggplot(aes(x=Year,y=Population,colour=GeoUID,group=GeoUID)) +
+    geom_line() +
+    geom_point() +
+    scale_y_continuous(labels=scales::comma,limits=c(0,NA))
+}
+
+toronto %>%
+  filter(GeoUID %in% crescent_town) %>%
+  plot_timelines() +
+  labs(title="Population in two neighbouring dissemination areas")
+

+

The population jumps back and forth between the two areas, while the +sum of the two is fairly steady from 1981 on. People did not move back +and forth, their homes got geocoded to a different dissemination area in +different years.

+
+
+

Detecting anomalies +

+

tongfen_detect_anomalies lists the regions with +surprising drops. Only decreases are surprising, and a decrease needs to +be large in both relative and absolute terms. For each of these +candidate regions it finds the neighbouring region that takes away most +of the surprise when both are joined, and checks if that reduction is +large enough to justify joining them. Our regions are identified by +their GeoUID instead of the TongfenID the +function looks for by default.

+
+anomalies <- tongfen_detect_anomalies(toronto,timeline,id="GeoUID",total_surprise_cutoff=0.4)
+
+anomalies %>% filter(GeoUID %in% crescent_town)
+#> # A tibble: 2 × 7
+#>   GeoUID   surprise_count surprise_total period  neighbour surprise_total_joined
+#>   <chr>             <int>          <dbl> <chr>   <chr>                     <dbl>
+#> 1 35204370              3          0.933 1991-1… 35204765                  0.616
+#> 2 35204765              2          1.03  1981-1… 35204370                  0.142
+#> # ℹ 1 more variable: join <lgl>
+

Both areas are candidates, and each one is the neighbour that best +explains the surprising drops of the other. Not all candidates find a +neighbour to pair up with. Population does drop for real, for example +when a site gets cleared for redevelopment, and such regions are left +alone.

+
+anomalies %>% count(join)
+#> # A tibble: 2 × 2
+#>   join      n
+#>   <lgl> <int>
+#> 1 FALSE   254
+#> 2 TRUE    288
+
+
+

Joining regions +

+

tongfen_anomaly_joins joins the regions that qualify and +looks again, joined regions can have surprising drops that are +complemented by another neighbour. This repeats until there are no more +regions left to join. The result lists the regions that got joined +together with the identifier of the joined region they are now part of +and the round in which they first got joined.

+
+joins <- tongfen_anomaly_joins(toronto,timeline,id="GeoUID",total_surprise_cutoff=0.4)
+
+joins %>% filter(GeoUID %in% crescent_town)
+#> # A tibble: 2 × 3
+#>   GeoUID   GeoUID_joined round
+#>   <chr>    <chr>         <int>
+#> 1 35204370 35204370          1
+#> 2 35204765 35204370          1
+

tongfen_join_regions applies the joins to the data, +aggregating the variables and geometries of the regions that get joined +and leaving all others as they are.

+
+toronto_joined <- tongfen_join_regions(toronto,joins,id="GeoUID")
+
+c(original=nrow(toronto),joined=nrow(toronto_joined))
+#> original   joined 
+#>     3702     3430
+
+toronto_joined %>%
+  filter(GeoUID %in% crescent_town) %>%
+  plot_timelines() +
+  labs(title="Population in the joined region")
+

+

The map shows the regions that got joined around Crescent Town.

+
+bbox <- toronto %>% filter(GeoUID %in% crescent_town) %>% st_buffer(1500) %>% st_bbox()
+
+ggplot(toronto_joined %>% mutate(joined=GeoUID %in% joins$GeoUID_joined)) +
+  geom_sf(aes(fill=joined),linewidth=0.1) +
+  geom_sf(data=toronto,fill=NA,linewidth=0.1,linetype="dotted") +
+  scale_fill_manual(values=c("TRUE"="steelblue","FALSE"="whitesmoke"),guide="none") +
+  coord_sf(datum=NA,xlim=bbox[c("xmin","xmax")],ylim=bbox[c("ymin","ymax")]) +
+  labs(title="Joined regions around Crescent Town",
+       caption="Joined regions in blue, original dissemination areas dotted")
+

+
+
+

Tuning +

+

Joining regions trades geographic detail for consistency over time, +and how to best make that trade depends on the data and the application. +The parameters are documented in tongfen_detect_anomalies, +the most important ones are

+
    +
  • +rel_scale and abs_scale, the relative and +absolute decrease at which a change is half way to being fully +surprising. The defaults of a 25% drop and a drop of 200 are tuned to +population counts in regions of the size of dissemination areas.
    • +
  • +total_surprise_cutoff, how surprising the timeline of a +region needs to be to become a candidate. The default of 0.75 is +conservative, above we used 0.4 to also pick up less pronounced +cases.
    • +
  • +cutoff_fact, surprise_reduction_const and +sum_fact determine how much of the surprise a neighbour +needs to take away for the regions to get joined.
    • +
+
+c(0.4,0.6,0.75) %>%
+  lapply(\(cutoff) tibble(total_surprise_cutoff=cutoff,
+                          regions_joined=tongfen_anomaly_joins(toronto,timeline,id="GeoUID",
+                                                               total_surprise_cutoff=cutoff) %>%
+                            nrow())) %>%
+  bind_rows()
+#> # A tibble: 3 × 2
+#>   total_surprise_cutoff regions_joined
+#>                   <dbl>          <int>
+#> 1                  0.4             477
+#> 2                  0.6             273
+#> 3                  0.75            152
+

Neighbours are by default determined by intersecting the geometries +of the regions. This can miss neighbours if the geometries have been +simplified, in that case the neighbours argument takes a +table with the identifiers of neighbouring regions or a neighbours list +from the spdep package.

+
+
+

Anomalies in TongFen data +

+

The functions work the same way on data on a common geography built +by TongFen, where the regions are identified by their +TongfenID. As an example we look at the dissemination area +level population in the City of Vancouver for the 2001 through 2021 +censuses.

+
+regions <- list(CSD="5915022")
+datasets <- c("CA01","CA06","CA11","CA16","CA21")
+meta <- meta_for_additive_variables(datasets,"Population")
+
+vancouver <- get_tongfen_ca_census(regions=regions,meta=meta,level="DA",base_geo="CA21",quiet=TRUE)
+
+joins <- tongfen_anomaly_joins(vancouver,paste0("Population_",datasets),total_surprise_cutoff=0.4)
+joins
+#> # A tibble: 4 × 3
+#>   TongfenID TongfenID_joined round
+#>   <chr>     <chr>            <int>
+#> 1 59150762  59150762             1
+#> 2 59153181  59150762             1
+#> 3 59150765  59150765             1
+#> 4 59150770  59150765             1
+

Passing the metadata to tongfen_join_regions makes sure +the variables get aggregated the right way, numeric variables that are +not part of the metadata are assumed to be additive.

+
+vancouver_joined <- tongfen_join_regions(vancouver,joins,meta)
+

The TongfenUID of the joined regions lists all the +dissemination areas they are made up of.

+
+vancouver_joined %>%
+  st_drop_geometry() %>%
+  filter(TongfenID %in% joins$TongfenID_joined) %>%
+  select(TongfenID,TongfenUID,starts_with("Population"))
+#> # A tibble: 2 × 7
+#>   TongfenID TongfenUID           Population_CA21 Population_CA01 Population_CA06
+#>   <chr>     <chr>                          <int>           <dbl>           <dbl>
+#> 1 59150762  GeoUIDCA01:59150762…            1412            1148            1266
+#> 2 59150765  GeoUIDCA01:59150765…            4082            1895            2499
+#> # ℹ 2 more variables: Population_CA11 <dbl>, Population_CA16 <dbl>
+

Variables that are not additive, like averages, can only be +aggregated this way if the variable they are averaged over is part of +the data. The alternative that always works is to join the regions in +the correspondence the common geography was built from, and use the +joined correspondence to aggregate the original data. This also is the +way to use the joins for data other than the one that was used to detect +the anomalies, for example to get average rents on the corrected +geography.

+
+correspondence <- get_tongfen_correspondence_ca_census(geo_datasets=datasets,regions=regions,
+                                                       level="DA",quiet=TRUE) %>%
+  tongfen_join_correspondence(joins)
+
+rent_meta <- meta_for_ca_census_vectors(c(rent_2006="v_CA06_2050",rent_2016="v_CA16_4901"))
+
+rent_data <- c("CA06","CA16") %>%
+  lapply(\(ds) get_census(ds,regions=regions,level="DA",labels="short",quiet=TRUE,
+                          vectors=rent_meta %>% filter(geo_dataset==ds) %>% pull(variable),
+                          geo_format=if (ds=="CA16") "sf" else NA) %>%
+           rename(!!paste0("GeoUID",ds):="GeoUID")) %>%
+  setNames(c("CA06","CA16"))
+
+rents <- tongfen_aggregate(rent_data,correspondence,rent_meta,base_geo="CA16")
+
+rents %>%
+  st_drop_geometry() %>%
+  filter(TongfenID %in% joins$TongfenID_joined) %>%
+  select(TongfenID,rent_2006,rent_2016)
+#> # A tibble: 2 × 3
+#>   TongfenID rent_2006 rent_2016
+#>   <chr>         <dbl>     <dbl>
+#> 1 59150762       452.      674.
+#> 2 59150765       474.      913.
+
+
+
+ + + + +
+ + + + + + + diff --git a/docs/articles/tongfen_anomalies.md b/docs/articles/tongfen_anomalies.md new file mode 100644 index 0000000..4e1c12f --- /dev/null +++ b/docs/articles/tongfen_anomalies.md @@ -0,0 +1,311 @@ +# Geocoding anomalies in TongFen timelines + +TongFen makes data on different geographies comparable by aggregating it +up to a common geography. The result is only as good as the geocoding +that assigned the underlying data to geographic regions in the first +place. Geocoding changes over time, and the same dwelling units, and the +people living in them, can get assigned to different neighbouring +regions in different years. In a timeline on a common geography this +shows up as a surprising drop in one region that is offset by a jump in +a neighbouring region. + +The fix is in the spirit of TongFen, joining the affected regions gives +a slightly coarser geography on which the data is consistent over time. +The functions in this vignette look for such patterns and join the +regions on demand. The method is explained in more detail in a [blog +post](https://doodles.mountainmath.ca/posts/2024-07-26-geocoding-errors-in-aggregate-data/), +this vignette follows the example from that post. + +``` r + +library(dplyr) +library(tidyr) +library(ggplot2) +library(cancensus) +library(sf) +library(tongfen) +# cancensus::set_api_key("") +``` + +## Population timelines for Toronto + +As an example we take the population from the 1971 through 2011 censuses +that Statistics Canada tabulated on 2016 dissemination areas, together +with the 2016 population. All data comes on the same geography, so there +is no need to TongFen, but the data for the earlier years is geocoded +from the road network and block face of the time, which does not always +match up with 2016 dissemination areas. + +``` r + +years <- c(1971,seq(1981,2011,5)) +vectors <- c(setNames(paste0("v_CA",years,"x16_1"),years),"2016"="v_CA16_1") +timeline <- names(vectors) + +toronto <- get_census("CA16CT",regions=list(CSD="3520005"),vectors=vectors, + level="DA",geo_format="sf",quiet=TRUE) %>% + select(GeoUID,all_of(timeline)) %>% + mutate(across(all_of(timeline),\(x) coalesce(x,0))) +``` + +Dissemination areas without population in a given year come back as +missing values. Changes from or to a missing value are never considered +surprising, so we set them to zero to mark these as areas where nobody +got counted. + +The area around Crescent Town shows what the problem looks like. + +``` r + +crescent_town <- c("35204370","35204765") + +plot_timelines <- function(data) { + data %>% + st_drop_geometry() %>% + pivot_longer(all_of(timeline),names_to="Year",values_to="Population") %>% + ggplot(aes(x=Year,y=Population,colour=GeoUID,group=GeoUID)) + + geom_line() + + geom_point() + + scale_y_continuous(labels=scales::comma,limits=c(0,NA)) +} + +toronto %>% + filter(GeoUID %in% crescent_town) %>% + plot_timelines() + + labs(title="Population in two neighbouring dissemination areas") +``` + +![](tongfen_anomalies_files/figure-html/unnamed-chunk-3-1.png) + +The population jumps back and forth between the two areas, while the sum +of the two is fairly steady from 1981 on. People did not move back and +forth, their homes got geocoded to a different dissemination area in +different years. + +## Detecting anomalies + +`tongfen_detect_anomalies` lists the regions with surprising drops. Only +decreases are surprising, and a decrease needs to be large in both +relative and absolute terms. For each of these candidate regions it +finds the neighbouring region that takes away most of the surprise when +both are joined, and checks if that reduction is large enough to justify +joining them. Our regions are identified by their `GeoUID` instead of +the `TongfenID` the function looks for by default. + +``` r + +anomalies <- tongfen_detect_anomalies(toronto,timeline,id="GeoUID",total_surprise_cutoff=0.4) + +anomalies %>% filter(GeoUID %in% crescent_town) +#> # A tibble: 2 × 7 +#> GeoUID surprise_count surprise_total period neighbour surprise_total_joined +#> +#> 1 35204370 3 0.933 1991-1… 35204765 0.616 +#> 2 35204765 2 1.03 1981-1… 35204370 0.142 +#> # ℹ 1 more variable: join +``` + +Both areas are candidates, and each one is the neighbour that best +explains the surprising drops of the other. Not all candidates find a +neighbour to pair up with. Population does drop for real, for example +when a site gets cleared for redevelopment, and such regions are left +alone. + +``` r + +anomalies %>% count(join) +#> # A tibble: 2 × 2 +#> join n +#> +#> 1 FALSE 254 +#> 2 TRUE 288 +``` + +## Joining regions + +`tongfen_anomaly_joins` joins the regions that qualify and looks again, +joined regions can have surprising drops that are complemented by +another neighbour. This repeats until there are no more regions left to +join. The result lists the regions that got joined together with the +identifier of the joined region they are now part of and the round in +which they first got joined. + +``` r + +joins <- tongfen_anomaly_joins(toronto,timeline,id="GeoUID",total_surprise_cutoff=0.4) + +joins %>% filter(GeoUID %in% crescent_town) +#> # A tibble: 2 × 3 +#> GeoUID GeoUID_joined round +#> +#> 1 35204370 35204370 1 +#> 2 35204765 35204370 1 +``` + +`tongfen_join_regions` applies the joins to the data, aggregating the +variables and geometries of the regions that get joined and leaving all +others as they are. + +``` r + +toronto_joined <- tongfen_join_regions(toronto,joins,id="GeoUID") + +c(original=nrow(toronto),joined=nrow(toronto_joined)) +#> original joined +#> 3702 3430 +``` + +``` r + +toronto_joined %>% + filter(GeoUID %in% crescent_town) %>% + plot_timelines() + + labs(title="Population in the joined region") +``` + +![](tongfen_anomalies_files/figure-html/unnamed-chunk-8-1.png) + +The map shows the regions that got joined around Crescent Town. + +``` r + +bbox <- toronto %>% filter(GeoUID %in% crescent_town) %>% st_buffer(1500) %>% st_bbox() + +ggplot(toronto_joined %>% mutate(joined=GeoUID %in% joins$GeoUID_joined)) + + geom_sf(aes(fill=joined),linewidth=0.1) + + geom_sf(data=toronto,fill=NA,linewidth=0.1,linetype="dotted") + + scale_fill_manual(values=c("TRUE"="steelblue","FALSE"="whitesmoke"),guide="none") + + coord_sf(datum=NA,xlim=bbox[c("xmin","xmax")],ylim=bbox[c("ymin","ymax")]) + + labs(title="Joined regions around Crescent Town", + caption="Joined regions in blue, original dissemination areas dotted") +``` + +![](tongfen_anomalies_files/figure-html/unnamed-chunk-9-1.png) + +## Tuning + +Joining regions trades geographic detail for consistency over time, and +how to best make that trade depends on the data and the application. The +parameters are documented in `tongfen_detect_anomalies`, the most +important ones are + +- `rel_scale` and `abs_scale`, the relative and absolute decrease at + which a change is half way to being fully surprising. The defaults of + a 25% drop and a drop of 200 are tuned to population counts in regions + of the size of dissemination areas. +- `total_surprise_cutoff`, how surprising the timeline of a region needs + to be to become a candidate. The default of 0.75 is conservative, + above we used 0.4 to also pick up less pronounced cases. +- `cutoff_fact`, `surprise_reduction_const` and `sum_fact` determine how + much of the surprise a neighbour needs to take away for the regions to + get joined. + +``` r + +c(0.4,0.6,0.75) %>% + lapply(\(cutoff) tibble(total_surprise_cutoff=cutoff, + regions_joined=tongfen_anomaly_joins(toronto,timeline,id="GeoUID", + total_surprise_cutoff=cutoff) %>% + nrow())) %>% + bind_rows() +#> # A tibble: 3 × 2 +#> total_surprise_cutoff regions_joined +#> +#> 1 0.4 477 +#> 2 0.6 273 +#> 3 0.75 152 +``` + +Neighbours are by default determined by intersecting the geometries of +the regions. This can miss neighbours if the geometries have been +simplified, in that case the `neighbours` argument takes a table with +the identifiers of neighbouring regions or a neighbours list from the +**spdep** package. + +## Anomalies in TongFen data + +The functions work the same way on data on a common geography built by +TongFen, where the regions are identified by their `TongfenID`. As an +example we look at the dissemination area level population in the City +of Vancouver for the 2001 through 2021 censuses. + +``` r + +regions <- list(CSD="5915022") +datasets <- c("CA01","CA06","CA11","CA16","CA21") +meta <- meta_for_additive_variables(datasets,"Population") + +vancouver <- get_tongfen_ca_census(regions=regions,meta=meta,level="DA",base_geo="CA21",quiet=TRUE) + +joins <- tongfen_anomaly_joins(vancouver,paste0("Population_",datasets),total_surprise_cutoff=0.4) +joins +#> # A tibble: 4 × 3 +#> TongfenID TongfenID_joined round +#> +#> 1 59150762 59150762 1 +#> 2 59153181 59150762 1 +#> 3 59150765 59150765 1 +#> 4 59150770 59150765 1 +``` + +Passing the metadata to `tongfen_join_regions` makes sure the variables +get aggregated the right way, numeric variables that are not part of the +metadata are assumed to be additive. + +``` r + +vancouver_joined <- tongfen_join_regions(vancouver,joins,meta) +``` + +The `TongfenUID` of the joined regions lists all the dissemination areas +they are made up of. + +``` r + +vancouver_joined %>% + st_drop_geometry() %>% + filter(TongfenID %in% joins$TongfenID_joined) %>% + select(TongfenID,TongfenUID,starts_with("Population")) +#> # A tibble: 2 × 7 +#> TongfenID TongfenUID Population_CA21 Population_CA01 Population_CA06 +#> +#> 1 59150762 GeoUIDCA01:59150762… 1412 1148 1266 +#> 2 59150765 GeoUIDCA01:59150765… 4082 1895 2499 +#> # ℹ 2 more variables: Population_CA11 , Population_CA16 +``` + +Variables that are not additive, like averages, can only be aggregated +this way if the variable they are averaged over is part of the data. The +alternative that always works is to join the regions in the +correspondence the common geography was built from, and use the joined +correspondence to aggregate the original data. This also is the way to +use the joins for data other than the one that was used to detect the +anomalies, for example to get average rents on the corrected geography. + +``` r + +correspondence <- get_tongfen_correspondence_ca_census(geo_datasets=datasets,regions=regions, + level="DA",quiet=TRUE) %>% + tongfen_join_correspondence(joins) + +rent_meta <- meta_for_ca_census_vectors(c(rent_2006="v_CA06_2050",rent_2016="v_CA16_4901")) + +rent_data <- c("CA06","CA16") %>% + lapply(\(ds) get_census(ds,regions=regions,level="DA",labels="short",quiet=TRUE, + vectors=rent_meta %>% filter(geo_dataset==ds) %>% pull(variable), + geo_format=if (ds=="CA16") "sf" else NA) %>% + rename(!!paste0("GeoUID",ds):="GeoUID")) %>% + setNames(c("CA06","CA16")) + +rents <- tongfen_aggregate(rent_data,correspondence,rent_meta,base_geo="CA16") + +rents %>% + st_drop_geometry() %>% + filter(TongfenID %in% joins$TongfenID_joined) %>% + select(TongfenID,rent_2006,rent_2016) +#> # A tibble: 2 × 3 +#> TongfenID rent_2006 rent_2016 +#> +#> 1 59150762 452. 674. +#> 2 59150765 474. 913. +``` diff --git a/docs/articles/tongfen_anomalies_files/figure-html/unnamed-chunk-3-1.png b/docs/articles/tongfen_anomalies_files/figure-html/unnamed-chunk-3-1.png new file mode 100644 index 0000000..f4533c8 Binary files /dev/null and b/docs/articles/tongfen_anomalies_files/figure-html/unnamed-chunk-3-1.png differ diff --git a/docs/articles/tongfen_anomalies_files/figure-html/unnamed-chunk-8-1.png b/docs/articles/tongfen_anomalies_files/figure-html/unnamed-chunk-8-1.png new file mode 100644 index 0000000..12a93f2 Binary files /dev/null and b/docs/articles/tongfen_anomalies_files/figure-html/unnamed-chunk-8-1.png differ diff --git a/docs/articles/tongfen_anomalies_files/figure-html/unnamed-chunk-9-1.png b/docs/articles/tongfen_anomalies_files/figure-html/unnamed-chunk-9-1.png new file mode 100644 index 0000000..eea4128 Binary files /dev/null and b/docs/articles/tongfen_anomalies_files/figure-html/unnamed-chunk-9-1.png differ diff --git a/docs/articles/tongfen_ca.html b/docs/articles/tongfen_ca.html index 8f7bdea..a7a225e 100644 --- a/docs/articles/tongfen_ca.html +++ b/docs/articles/tongfen_ca.html @@ -25,7 +25,7 @@ tongfen - 0.3.8 + 0.3.9 + + + + + +
+
+
+

Determine regions to join to correct for likely geocoding anomalies

+ Source: R/tongfen_anomalies.R +
tongfen_anomaly_joins.Rd
+
+ +
+

[Experimental]

+

Looks for regions with surprising drops in the timeline of a count variable that are complemented by +a neighbouring region, as explained in `tongfen_detect_anomalies`, and joins them. This gets +repeated on the joined regions until there are no more regions left that qualify to get joined. +In each round a region only gets joined with one other region, the most surprising regions go first.

+

Joining regions trades geographic detail for timelines that are consistent over time. The parameters +control how aggressively regions get joined and are best calibrated on the data at hand, erring on the side of +joining too few regions risks keeping geocoding problems, erring on the other side risks removing real +changes and needlessly coarsens the geography.

+

The result can be used to join the regions via `tongfen_join_regions`, or to update a correspondence +via `tongfen_join_correspondence`.

+
+ +
+

Usage

+ 
tongfen_anomaly_joins(
+  data,
+  variables,
+  id = "TongfenID",
+  neighbours = NULL,
+  rel_scale = 0.25,
+  abs_scale = 200,
+  p = 4,
+  surprise_cutoff = 0.15,
+  total_surprise_cutoff = 0.75,
+  cutoff_fact = 0.6,
+  surprise_reduction_const = 0.15,
+  sum_fact = 0.7
+)
+
+ +
+

Arguments

+ + +
data
+

data on a common geography, with one row per region, for example as returned by +`tongfen_aggregate` or `get_tongfen_ca_census`. Needs to be of class sf unless `neighbours` is specified.

+ + +
variables
+

names of the columns holding the timeline of a count variable like population or +dwellings, in temporal order. Changes from or to a missing value are not surprising and don't make up for +surprising changes in neighbouring regions, and joined regions are missing a value if one of the regions they +are made up of is. Replace missing values by zero beforehand if they stand for regions where nothing got counted

+ + +
id
+

name of the column that uniquely identifies the regions, default is "TongfenID"

+ + +
neighbours
+

optional, neighbouring regions as a table with the identifiers of pairs of neighbouring +regions in the first two columns, or as a neighbours list like the ones returned by `spdep::poly2nb`. +By default all regions with intersecting geometries are neighbours, which can miss neighbours +if the geometries have been simplified and don't share their boundaries any more.

+ + +
rel_scale
+

relative decrease that is half way to full surprise, default is `0.25` for a 25% drop

+ + +
abs_scale
+

absolute decrease that is half way to full surprise, default is `200`

+ + +
p
+

exponent of the norm used to combine the surprises across the timeline into the total surprise. +Large values focus on the most surprising change, 1 adds up the surprises of all changes, default is `4`

+ + +
surprise_cutoff
+

changes with larger surprise count as surprising, default is `0.15`. Only +regions with at least one surprising change are candidates

+ + +
total_surprise_cutoff
+

only regions with larger total surprise are candidates, default is `0.75`

+ + +
cutoff_fact
+

join regions if the total surprise after joining is lower than this share of the +total surprise of the candidate region, default is `0.6`

+ + +
surprise_reduction_const
+

join regions if joining lowers the total surprise by more than this, +default is `0.15`

+ + +
sum_fact
+

join regions if the total surprise after joining is lower than `cutoff_fact * sum_fact` +times the sum of the total surprises of both regions, default is `0.7`

+ +
+
+

Value

+

A tibble with one row for each region that gets joined with other regions, with the identifier +of the region, the identifier of the joined region it becomes part of in the column named like the +identifier with suffix `_joined`, by default `TongfenID_joined`, and the `round` in which the region +first got joined to another region. The identifier of a joined region is the smallest +identifier of the regions it is made up of.

+
+ +
+

Examples

+ 
# Correct 2001 through 2021 dissemination area level population timelines in the
+# City of Vancouver for likely geocoding problems
+if (FALSE) { # \dontrun{
+datasets <- c("CA01","CA06","CA11","CA16","CA21")
+meta <- meta_for_additive_variables(datasets,"Population")
+data <- get_tongfen_ca_census(regions=list(CSD="5915022"),meta=meta,level="DA",base_geo="CA21")
+
+joins <- tongfen_anomaly_joins(data,paste0("Population_",datasets))
+corrected_data <- tongfen_join_regions(data,joins,meta)
+} # }
+
+
+
+ + +
+ + + + + + + diff --git a/docs/reference/tongfen_anomaly_joins.md b/docs/reference/tongfen_anomaly_joins.md new file mode 100644 index 0000000..5e39412 --- /dev/null +++ b/docs/reference/tongfen_anomaly_joins.md @@ -0,0 +1,139 @@ +# Determine regions to join to correct for likely geocoding anomalies + +**\[experimental\]** + +Looks for regions with surprising drops in the timeline of a count +variable that are complemented by a neighbouring region, as explained in +\`tongfen_detect_anomalies\`, and joins them. This gets repeated on the +joined regions until there are no more regions left that qualify to get +joined. In each round a region only gets joined with one other region, +the most surprising regions go first. + +Joining regions trades geographic detail for timelines that are +consistent over time. The parameters control how aggressively regions +get joined and are best calibrated on the data at hand, erring on the +side of joining too few regions risks keeping geocoding problems, erring +on the other side risks removing real changes and needlessly coarsens +the geography. + +The result can be used to join the regions via \`tongfen_join_regions\`, +or to update a correspondence via \`tongfen_join_correspondence\`. + +## Usage + +``` r +tongfen_anomaly_joins( + data, + variables, + id = "TongfenID", + neighbours = NULL, + rel_scale = 0.25, + abs_scale = 200, + p = 4, + surprise_cutoff = 0.15, + total_surprise_cutoff = 0.75, + cutoff_fact = 0.6, + surprise_reduction_const = 0.15, + sum_fact = 0.7 +) +``` + +## Arguments + +- data: + + data on a common geography, with one row per region, for example as + returned by \`tongfen_aggregate\` or \`get_tongfen_ca_census\`. Needs + to be of class sf unless \`neighbours\` is specified. + +- variables: + + names of the columns holding the timeline of a count variable like + population or dwellings, in temporal order. Changes from or to a + missing value are not surprising and don't make up for surprising + changes in neighbouring regions, and joined regions are missing a + value if one of the regions they are made up of is. Replace missing + values by zero beforehand if they stand for regions where nothing got + counted + +- id: + + name of the column that uniquely identifies the regions, default is + "TongfenID" + +- neighbours: + + optional, neighbouring regions as a table with the identifiers of + pairs of neighbouring regions in the first two columns, or as a + neighbours list like the ones returned by \`spdep::poly2nb\`. By + default all regions with intersecting geometries are neighbours, which + can miss neighbours if the geometries have been simplified and don't + share their boundaries any more. + +- rel_scale: + + relative decrease that is half way to full surprise, default is + \`0.25\` for a 25% drop + +- abs_scale: + + absolute decrease that is half way to full surprise, default is + \`200\` + +- p: + + exponent of the norm used to combine the surprises across the timeline + into the total surprise. Large values focus on the most surprising + change, 1 adds up the surprises of all changes, default is \`4\` + +- surprise_cutoff: + + changes with larger surprise count as surprising, default is \`0.15\`. + Only regions with at least one surprising change are candidates + +- total_surprise_cutoff: + + only regions with larger total surprise are candidates, default is + \`0.75\` + +- cutoff_fact: + + join regions if the total surprise after joining is lower than this + share of the total surprise of the candidate region, default is + \`0.6\` + +- surprise_reduction_const: + + join regions if joining lowers the total surprise by more than this, + default is \`0.15\` + +- sum_fact: + + join regions if the total surprise after joining is lower than + \`cutoff_fact \* sum_fact\` times the sum of the total surprises of + both regions, default is \`0.7\` + +## Value + +A tibble with one row for each region that gets joined with other +regions, with the identifier of the region, the identifier of the joined +region it becomes part of in the column named like the identifier with +suffix \`\_joined\`, by default \`TongfenID_joined\`, and the \`round\` +in which the region first got joined to another region. The identifier +of a joined region is the smallest identifier of the regions it is made +up of. + +## Examples + +``` r +# Correct 2001 through 2021 dissemination area level population timelines in the +# City of Vancouver for likely geocoding problems +if (FALSE) { # \dontrun{ +datasets <- c("CA01","CA06","CA11","CA16","CA21") +meta <- meta_for_additive_variables(datasets,"Population") +data <- get_tongfen_ca_census(regions=list(CSD="5915022"),meta=meta,level="DA",base_geo="CA21") + +joins <- tongfen_anomaly_joins(data,paste0("Population_",datasets)) +corrected_data <- tongfen_join_regions(data,joins,meta) +} # } +``` diff --git a/docs/reference/tongfen_ca_census_ct.html b/docs/reference/tongfen_ca_census_ct.html index cc9ebcc..090d3f2 100644 --- a/docs/reference/tongfen_ca_census_ct.html +++ b/docs/reference/tongfen_ca_census_ct.html @@ -9,7 +9,7 @@ tongfen - 0.3.8 + 0.3.9 + + + + + +
+
+
+

Detect likely geocoding anomalies in timelines on a common geography

+ Source: R/tongfen_anomalies.R +
tongfen_detect_anomalies.Rd
+
+ +
+

[Experimental]

+

TongFen is only as good as the geocoding that assigned the underlying data to geographic regions +in the first place. Geocoding varies over time, and the same dwelling units, and the people living +in them, can get assigned to different neighbouring regions in different years. In a timeline +on a common geography this shows up as a surprising drop in one region that is offset by +a corresponding jump in a neighbouring region.

+

This function lists the candidate regions with surprising drops in the given count variable, +together with the neighbouring region that takes away most of the surprise when both are joined. +Use it to check for possible problems and to calibrate the parameters before joining regions with +`tongfen_anomaly_joins`. Not all surprising drops are due to geocoding problems, a drop that is +not complemented by a neighbouring region is likely real.

+

The surprise of a change between two consecutive years ranges from 0 to 1. Only decreases are surprising, +the surprise is the product of the surprise of the relative and of the absolute decrease, so that it takes +a decrease that is large in both relative and absolute terms to be surprising. The total surprise of a +region is the `p`-norm of the surprises across all changes in the timeline.

+

A candidate region and its neighbour are flagged for joining if joining reduces the total surprise of the +candidate region to below `cutoff_fact` times its total surprise, or reduces it by more than +`surprise_reduction_const`, or reduces it to below `cutoff_fact * sum_fact` times the sum of +the total surprises of both regions. To keep this comparable the surprise after joining is computed +from the change of the joined regions relative to the counts of the candidate region alone.

+
+ +
+

Usage

+ 
tongfen_detect_anomalies(
+  data,
+  variables,
+  id = "TongfenID",
+  neighbours = NULL,
+  rel_scale = 0.25,
+  abs_scale = 200,
+  p = 4,
+  surprise_cutoff = 0.15,
+  total_surprise_cutoff = 0.75,
+  cutoff_fact = 0.6,
+  surprise_reduction_const = 0.15,
+  sum_fact = 0.7
+)
+
+ +
+

Arguments

+ + +
data
+

data on a common geography, with one row per region, for example as returned by +`tongfen_aggregate` or `get_tongfen_ca_census`. Needs to be of class sf unless `neighbours` is specified.

+ + +
variables
+

names of the columns holding the timeline of a count variable like population or +dwellings, in temporal order. Changes from or to a missing value are not surprising and don't make up for +surprising changes in neighbouring regions, and joined regions are missing a value if one of the regions they +are made up of is. Replace missing values by zero beforehand if they stand for regions where nothing got counted

+ + +
id
+

name of the column that uniquely identifies the regions, default is "TongfenID"

+ + +
neighbours
+

optional, neighbouring regions as a table with the identifiers of pairs of neighbouring +regions in the first two columns, or as a neighbours list like the ones returned by `spdep::poly2nb`. +By default all regions with intersecting geometries are neighbours, which can miss neighbours +if the geometries have been simplified and don't share their boundaries any more.

+ + +
rel_scale
+

relative decrease that is half way to full surprise, default is `0.25` for a 25% drop

+ + +
abs_scale
+

absolute decrease that is half way to full surprise, default is `200`

+ + +
p
+

exponent of the norm used to combine the surprises across the timeline into the total surprise. +Large values focus on the most surprising change, 1 adds up the surprises of all changes, default is `4`

+ + +
surprise_cutoff
+

changes with larger surprise count as surprising, default is `0.15`. Only +regions with at least one surprising change are candidates

+ + +
total_surprise_cutoff
+

only regions with larger total surprise are candidates, default is `0.75`

+ + +
cutoff_fact
+

join regions if the total surprise after joining is lower than this share of the +total surprise of the candidate region, default is `0.6`

+ + +
surprise_reduction_const
+

join regions if joining lowers the total surprise by more than this, +default is `0.15`

+ + +
sum_fact
+

join regions if the total surprise after joining is lower than `cutoff_fact * sum_fact` +times the sum of the total surprises of both regions, default is `0.7`

+ +
+
+

Value

+

A tibble with one row for each candidate region, most surprising first, with the identifier of the +region, the number of surprising changes `surprise_count`, the total surprise `surprise_total`, +the `period` with the most surprising change, the identifier of the `neighbour` that takes away most of +the surprise, the total surprise `surprise_total_joined` after joining both and `join` indicating if both +regions qualify to get joined.

+
+ +
+

Examples

+ 
# Check 2001 through 2021 dissemination area level population timelines in the
+# City of Vancouver for possible geocoding problems
+if (FALSE) { # \dontrun{
+datasets <- c("CA01","CA06","CA11","CA16","CA21")
+meta <- meta_for_additive_variables(datasets,"Population")
+data <- get_tongfen_ca_census(regions=list(CSD="5915022"),meta=meta,level="DA",base_geo="CA21")
+
+anomalies <- tongfen_detect_anomalies(data,paste0("Population_",datasets))
+} # }
+
+
+
+ + +
+ + + + + + + diff --git a/docs/reference/tongfen_detect_anomalies.md b/docs/reference/tongfen_detect_anomalies.md new file mode 100644 index 0000000..6ba66b3 --- /dev/null +++ b/docs/reference/tongfen_detect_anomalies.md @@ -0,0 +1,152 @@ +# Detect likely geocoding anomalies in timelines on a common geography + +**\[experimental\]** + +TongFen is only as good as the geocoding that assigned the underlying +data to geographic regions in the first place. Geocoding varies over +time, and the same dwelling units, and the people living in them, can +get assigned to different neighbouring regions in different years. In a +timeline on a common geography this shows up as a surprising drop in one +region that is offset by a corresponding jump in a neighbouring region. + +This function lists the candidate regions with surprising drops in the +given count variable, together with the neighbouring region that takes +away most of the surprise when both are joined. Use it to check for +possible problems and to calibrate the parameters before joining regions +with \`tongfen_anomaly_joins\`. Not all surprising drops are due to +geocoding problems, a drop that is not complemented by a neighbouring +region is likely real. + +The surprise of a change between two consecutive years ranges from 0 +to 1. Only decreases are surprising, the surprise is the product of the +surprise of the relative and of the absolute decrease, so that it takes +a decrease that is large in both relative and absolute terms to be +surprising. The total surprise of a region is the \`p\`-norm of the +surprises across all changes in the timeline. + +A candidate region and its neighbour are flagged for joining if joining +reduces the total surprise of the candidate region to below +\`cutoff_fact\` times its total surprise, or reduces it by more than +\`surprise_reduction_const\`, or reduces it to below \`cutoff_fact \* +sum_fact\` times the sum of the total surprises of both regions. To keep +this comparable the surprise after joining is computed from the change +of the joined regions relative to the counts of the candidate region +alone. + +## Usage + +``` r +tongfen_detect_anomalies( + data, + variables, + id = "TongfenID", + neighbours = NULL, + rel_scale = 0.25, + abs_scale = 200, + p = 4, + surprise_cutoff = 0.15, + total_surprise_cutoff = 0.75, + cutoff_fact = 0.6, + surprise_reduction_const = 0.15, + sum_fact = 0.7 +) +``` + +## Arguments + +- data: + + data on a common geography, with one row per region, for example as + returned by \`tongfen_aggregate\` or \`get_tongfen_ca_census\`. Needs + to be of class sf unless \`neighbours\` is specified. + +- variables: + + names of the columns holding the timeline of a count variable like + population or dwellings, in temporal order. Changes from or to a + missing value are not surprising and don't make up for surprising + changes in neighbouring regions, and joined regions are missing a + value if one of the regions they are made up of is. Replace missing + values by zero beforehand if they stand for regions where nothing got + counted + +- id: + + name of the column that uniquely identifies the regions, default is + "TongfenID" + +- neighbours: + + optional, neighbouring regions as a table with the identifiers of + pairs of neighbouring regions in the first two columns, or as a + neighbours list like the ones returned by \`spdep::poly2nb\`. By + default all regions with intersecting geometries are neighbours, which + can miss neighbours if the geometries have been simplified and don't + share their boundaries any more. + +- rel_scale: + + relative decrease that is half way to full surprise, default is + \`0.25\` for a 25% drop + +- abs_scale: + + absolute decrease that is half way to full surprise, default is + \`200\` + +- p: + + exponent of the norm used to combine the surprises across the timeline + into the total surprise. Large values focus on the most surprising + change, 1 adds up the surprises of all changes, default is \`4\` + +- surprise_cutoff: + + changes with larger surprise count as surprising, default is \`0.15\`. + Only regions with at least one surprising change are candidates + +- total_surprise_cutoff: + + only regions with larger total surprise are candidates, default is + \`0.75\` + +- cutoff_fact: + + join regions if the total surprise after joining is lower than this + share of the total surprise of the candidate region, default is + \`0.6\` + +- surprise_reduction_const: + + join regions if joining lowers the total surprise by more than this, + default is \`0.15\` + +- sum_fact: + + join regions if the total surprise after joining is lower than + \`cutoff_fact \* sum_fact\` times the sum of the total surprises of + both regions, default is \`0.7\` + +## Value + +A tibble with one row for each candidate region, most surprising first, +with the identifier of the region, the number of surprising changes +\`surprise_count\`, the total surprise \`surprise_total\`, the +\`period\` with the most surprising change, the identifier of the +\`neighbour\` that takes away most of the surprise, the total surprise +\`surprise_total_joined\` after joining both and \`join\` indicating if +both regions qualify to get joined. + +## Examples + +``` r +# Check 2001 through 2021 dissemination area level population timelines in the +# City of Vancouver for possible geocoding problems +if (FALSE) { # \dontrun{ +datasets <- c("CA01","CA06","CA11","CA16","CA21") +meta <- meta_for_additive_variables(datasets,"Population") +data <- get_tongfen_ca_census(regions=list(CSD="5915022"),meta=meta,level="DA",base_geo="CA21") + +anomalies <- tongfen_detect_anomalies(data,paste0("Population_",datasets)) +} # } +``` diff --git a/docs/reference/tongfen_estimate.html b/docs/reference/tongfen_estimate.html index e9eeac4..603f34c 100644 --- a/docs/reference/tongfen_estimate.html +++ b/docs/reference/tongfen_estimate.html @@ -13,7 +13,7 @@ tongfen - 0.3.8 + 0.3.9 + + + + + +
+
+
+

Join regions in a correspondence

+ Source: R/tongfen_anomalies.R +
tongfen_join_correspondence.Rd
+
+ +
+

[Experimental]

+

Updates a correspondence so that the given regions are joined, for example to correct for likely geocoding +anomalies as determined by `tongfen_anomaly_joins`. The updated correspondence can be used in `tongfen_aggregate` +to aggregate data on the coarser common geography, which works for all variables `tongfen_aggregate` can +deal with and for data that was not part of detecting the anomalies.

+
+ +
+

Usage

+ 
tongfen_join_correspondence(correspondence, joins)
+
+ +
+

Arguments

+ + +
correspondence
+

correspondence table with columns the unique geographic identifiers for each of the +geographies and the TongfenID and TongfenUID, as for example returned by `estimate_tongfen_correspondence` +or `get_tongfen_correspondence_ca_census`

+ + +
joins
+

table with the regions to join as returned by `tongfen_anomaly_joins`, with columns +`TongfenID` and `TongfenID_joined`

+ +
+
+

Value

+

The correspondence with updated TongfenID and TongfenUID for the regions that got joined. If the +correspondence has a TongfenMethod column "anomaly" gets added to the method of the regions that got joined.

+
+ +
+

Examples

+ 
# Correct for likely geocoding problems in dissemination area level population timelines
+# and use the updated correspondence to aggregate data on the corrected common geography
+if (FALSE) { # \dontrun{
+regions <- list(CSD="5915022")
+datasets <- c("CA01","CA06","CA11","CA16","CA21")
+meta <- meta_for_additive_variables(datasets,"Population")
+data <- get_tongfen_ca_census(regions=regions,meta=meta,level="DA",base_geo="CA21")
+joins <- tongfen_anomaly_joins(data,paste0("Population_",datasets))
+
+correspondence <- get_tongfen_correspondence_ca_census(geo_datasets=datasets,
+                                                       regions=regions,level="DA") %>%
+  tongfen_join_correspondence(joins)
+} # }
+
+
+
+ + +
+ + + + + + + diff --git a/docs/reference/tongfen_join_correspondence.md b/docs/reference/tongfen_join_correspondence.md new file mode 100644 index 0000000..772d411 --- /dev/null +++ b/docs/reference/tongfen_join_correspondence.md @@ -0,0 +1,55 @@ +# Join regions in a correspondence + +**\[experimental\]** + +Updates a correspondence so that the given regions are joined, for +example to correct for likely geocoding anomalies as determined by +\`tongfen_anomaly_joins\`. The updated correspondence can be used in +\`tongfen_aggregate\` to aggregate data on the coarser common geography, +which works for all variables \`tongfen_aggregate\` can deal with and +for data that was not part of detecting the anomalies. + +## Usage + +``` r +tongfen_join_correspondence(correspondence, joins) +``` + +## Arguments + +- correspondence: + + correspondence table with columns the unique geographic identifiers + for each of the geographies and the TongfenID and TongfenUID, as for + example returned by \`estimate_tongfen_correspondence\` or + \`get_tongfen_correspondence_ca_census\` + +- joins: + + table with the regions to join as returned by + \`tongfen_anomaly_joins\`, with columns \`TongfenID\` and + \`TongfenID_joined\` + +## Value + +The correspondence with updated TongfenID and TongfenUID for the regions +that got joined. If the correspondence has a TongfenMethod column +"anomaly" gets added to the method of the regions that got joined. + +## Examples + +``` r +# Correct for likely geocoding problems in dissemination area level population timelines +# and use the updated correspondence to aggregate data on the corrected common geography +if (FALSE) { # \dontrun{ +regions <- list(CSD="5915022") +datasets <- c("CA01","CA06","CA11","CA16","CA21") +meta <- meta_for_additive_variables(datasets,"Population") +data <- get_tongfen_ca_census(regions=regions,meta=meta,level="DA",base_geo="CA21") +joins <- tongfen_anomaly_joins(data,paste0("Population_",datasets)) + +correspondence <- get_tongfen_correspondence_ca_census(geo_datasets=datasets, + regions=regions,level="DA") %>% + tongfen_join_correspondence(joins) +} # } +``` diff --git a/docs/reference/tongfen_join_regions.html b/docs/reference/tongfen_join_regions.html new file mode 100644 index 0000000..4aff240 --- /dev/null +++ b/docs/reference/tongfen_join_regions.html @@ -0,0 +1,144 @@ + +Join regions in data on a common geography — tongfen_join_regions • tongfen + Skip to contents + + +
+
+
+

Join regions in data on a common geography

+ Source: R/tongfen_anomalies.R +
tongfen_join_regions.Rd
+
+ +
+

[Experimental]

+

Joins regions in data that has already been aggregated to a common geography, for example to correct +for likely geocoding anomalies as determined by `tongfen_anomaly_joins`. The data, and the geometries if +the data is of class sf, of the regions that get joined are aggregated, all other regions are left as they are.

+

Variables are aggregated according to the metadata, numeric variables that are not part of the metadata +are assumed to be additive. Variables that are not additive, like averages, can only be aggregated if their +parent variable is part of the data. If that is not the case use `tongfen_join_correspondence` to update the +correspondence the data was built from and aggregate the original data again with `tongfen_aggregate`.

+
+ +
+

Usage

+ 
tongfen_join_regions(data, joins, meta = NULL, id = "TongfenID", na.rm = TRUE)
+
+ +
+

Arguments

+ + +
data
+

data on a common geography, with one row per region, for example as returned by +`tongfen_aggregate` or `get_tongfen_ca_census`

+ + +
joins
+

table with the regions to join as returned by `tongfen_anomaly_joins`, with the identifier +of the region and the identifier of the joined region it becomes part of in the column named like the +identifier with suffix `_joined`

+ + +
meta
+

optional metadata containing aggregation rules as for example returned by `meta_for_ca_census_vectors`, +variables are matched by their label. Numeric variables that are not part of the metadata are treated as +additive, if `NULL` (the default) that is the case for all numeric variables

+ + +
id
+

name of the column that uniquely identifies the regions, default is "TongfenID"

+ + +
na.rm
+

logical, determines how NA values should be treated when aggregating variables, +default is `TRUE`

+ +
+
+

Value

+

The data with the regions joined. Joined regions take the place and the identifier of the +region with the smallest identifier among the regions they are made up of. Variables that are +not numeric and not part of the metadata are `NA` for joined regions.

+
+ +
+

Examples

+ 
# Correct 2001 through 2021 dissemination area level population timelines in the
+# City of Vancouver for likely geocoding problems
+if (FALSE) { # \dontrun{
+datasets <- c("CA01","CA06","CA11","CA16","CA21")
+meta <- meta_for_additive_variables(datasets,"Population")
+data <- get_tongfen_ca_census(regions=list(CSD="5915022"),meta=meta,level="DA",base_geo="CA21")
+
+joins <- tongfen_anomaly_joins(data,paste0("Population_",datasets))
+corrected_data <- tongfen_join_regions(data,joins,meta)
+} # }
+
+
+
+ + +
+ + + + + + + diff --git a/docs/reference/tongfen_join_regions.md b/docs/reference/tongfen_join_regions.md new file mode 100644 index 0000000..40d63d0 --- /dev/null +++ b/docs/reference/tongfen_join_regions.md @@ -0,0 +1,77 @@ +# Join regions in data on a common geography + +**\[experimental\]** + +Joins regions in data that has already been aggregated to a common +geography, for example to correct for likely geocoding anomalies as +determined by \`tongfen_anomaly_joins\`. The data, and the geometries if +the data is of class sf, of the regions that get joined are aggregated, +all other regions are left as they are. + +Variables are aggregated according to the metadata, numeric variables +that are not part of the metadata are assumed to be additive. Variables +that are not additive, like averages, can only be aggregated if their +parent variable is part of the data. If that is not the case use +\`tongfen_join_correspondence\` to update the correspondence the data +was built from and aggregate the original data again with +\`tongfen_aggregate\`. + +## Usage + +``` r +tongfen_join_regions(data, joins, meta = NULL, id = "TongfenID", na.rm = TRUE) +``` + +## Arguments + +- data: + + data on a common geography, with one row per region, for example as + returned by \`tongfen_aggregate\` or \`get_tongfen_ca_census\` + +- joins: + + table with the regions to join as returned by + \`tongfen_anomaly_joins\`, with the identifier of the region and the + identifier of the joined region it becomes part of in the column named + like the identifier with suffix \`\_joined\` + +- meta: + + optional metadata containing aggregation rules as for example returned + by \`meta_for_ca_census_vectors\`, variables are matched by their + label. Numeric variables that are not part of the metadata are treated + as additive, if \`NULL\` (the default) that is the case for all + numeric variables + +- id: + + name of the column that uniquely identifies the regions, default is + "TongfenID" + +- na.rm: + + logical, determines how NA values should be treated when aggregating + variables, default is \`TRUE\` + +## Value + +The data with the regions joined. Joined regions take the place and the +identifier of the region with the smallest identifier among the regions +they are made up of. Variables that are not numeric and not part of the +metadata are \`NA\` for joined regions. + +## Examples + +``` r +# Correct 2001 through 2021 dissemination area level population timelines in the +# City of Vancouver for likely geocoding problems +if (FALSE) { # \dontrun{ +datasets <- c("CA01","CA06","CA11","CA16","CA21") +meta <- meta_for_additive_variables(datasets,"Population") +data <- get_tongfen_ca_census(regions=list(CSD="5915022"),meta=meta,level="DA",base_geo="CA21") + +joins <- tongfen_anomaly_joins(data,paste0("Population_",datasets)) +corrected_data <- tongfen_join_regions(data,joins,meta) +} # } +``` diff --git a/docs/reference/tongfen_tag_largest_overlap.html b/docs/reference/tongfen_tag_largest_overlap.html index c580c9c..4d4ff1a 100644 --- a/docs/reference/tongfen_tag_largest_overlap.html +++ b/docs/reference/tongfen_tag_largest_overlap.html @@ -9,7 +9,7 @@ tongfen - 0.3.8 + 0.3.9