[jeancochrane]

Background

This PR fixes the geometry ingest by switching from geoarrow to sfarrow. We do this because the version of geoarrow that we rely on (0.1.0) is very outdated and is no longer compatible with R 4.6. This change mirrors an identical change we made in data-architecture recently: ccao-data/data-architecture#1116

The new geometry output is in s3://z-dev-ccao-data-ptaxsim-us-east-1/pin_geometry_raw/part-0.parquet. Once we're ready to construct the 2025 database, we can copy this file into the prod bucket.

Closes #108.

Testing

I had two goals in testing this change:

1. Confirm that the script will resolve the unexpectedly-changed year ranges reported in PTAXSIM 2025 data upload #107 (comment)
2. Confirm that reading and writing with sfarrow produces identical geometry columns

1. Confirm the script resolves unexpectedly-changed year ranges

Save this script to data-raw/pin/test_geometry_ranges.R:

Expand to view the QC script

# Checks that the pin_geometry logic in pin.R (sfarrow read + town fill +
# per-town range collapsing) reproduces the 2024.0.0 geometry ranges and shapes for blocks
# whose PIN10-years went missing in 2025.0.0 (see qc_2025_issues.qmd, issue 7).
#
# Usage: Rscript data-raw/pin/test_geometry_ranges.R [old_db] [block ...]
# Blocks are the first 7 digits of a PIN10. Defaults to 17-09-419.
library(arrow)
library(dplyr)
library(DBI)
library(glue)
library(purrr)
library(sf)
library(tidyr)

args <- commandArgs(trailingOnly = TRUE)
old_db <- if (length(args) >= 1) args[[1]] else "ptaxsim-2024.0.0.db"
blocks <- if (length(args) >= 2) args[-1] else "1709419"
max_year <- 2024

# Same steps as pin.R, restricted to the blocks and to years <= 2024 -----------
remote_bucket_geometry <- "s3://ccao-data-warehouse-us-east-1/spatial/parcel"
full <- map(blocks, \(block) {
  arrow::open_dataset(remote_bucket_geometry) %>%
    filter(
      year >= 2006 & year <= max_year,
      pin10 >= paste0(block, "000"), pin10 <= paste0(block, "999")
    ) %>%
    select(year, x = lon, y = lat, pin10, town_code, geometry) %>%
    sfarrow::read_sf_dataset()
}) %>%
  bind_rows() %>%
  arrange(pin10, year) %>%
  group_by(pin10) %>%
  fill(town_code, .direction = "updown") %>%
  ungroup() %>%
  mutate(town_code = replace_na(town_code, 99))

new_sf <- map(unique(full$town_code), \(town) {
  full %>%
    filter(town_code == town) %>%
    mutate(area = st_area(geometry)) %>%
    group_by(pin10) %>%
    arrange(pin10, year) %>%
    mutate(across(c(area, x, y), lag, .names = "lag_{.col}")) %>%
    mutate(
      diff_area = !(abs(area - lag_area) < units::set_units(0.001, "m^2")),
      diff_cent = !(abs(x - lag_x) < 0.00001 & abs(y - lag_y) < 0.00001),
      pin_group = cumsum(
        (diff_area & diff_cent) |
          (is.na(diff_area) & is.na(diff_cent))
      )
    ) %>%
    group_by(pin10, pin_group) %>%
    mutate(start_year = min(year), end_year = max(year)) %>%
    filter(row_number() == 1) %>%
    ungroup() %>%
    select(pin10, start_year, end_year, geometry)
}) %>%
  bind_rows() %>%
  arrange(pin10, start_year)
new <- st_drop_geometry(new_sf)

# Ranges from the 2024.0.0 DB --------------------------------------------------
conn <- dbConnect(RSQLite::SQLite(), old_db)
like <- paste0("substr(pin10, 1, 7) = '", blocks, "'", collapse = " OR ")
old <- dbGetQuery(
  conn,
  glue("
    SELECT pin10, start_year, end_year, geometry AS wkt
    FROM pin_geometry_raw
    WHERE {like}
  ")
) %>%
  as_tibble() %>%
  mutate(across(c(start_year, end_year), as.integer)) %>%
  arrange(pin10, start_year)
dbDisconnect(conn)

# Compare ----------------------------------------------------------------------
cat("PIN10s (old / new):", n_distinct(old$pin10), "/", n_distinct(new$pin10), "\n")
cat("Ranges (old / new):", nrow(old), "/", nrow(new), "\n")
cat("Ranges present in only one side:\n")
print(anti_join(old[1:3], new, by = c("pin10", "start_year", "end_year")) %>%
  mutate(side = "old only") %>%
  bind_rows(anti_join(new, old[1:3], by = c("pin10", "start_year", "end_year")) %>%
    mutate(side = "new only")) %>%
  count(side, start_year, end_year))

stopifnot(
  "New ranges differ from 2024.0.0" = isTRUE(all.equal(
    as.data.frame(old[1:3]), as.data.frame(new),
    check.attributes = FALSE
  ))
)
cat("OK: ranges match the 2024.0.0 DB\n")

# Compare geometries. The old DB stores WKT written by create_db.R with
# st_as_text()'s default precision, so produce the new WKT the same way.
same_geom <- old$wkt == st_as_text(st_geometry(new_sf))
cat("Geometries identical:", sum(same_geom), "/", length(same_geom), "\n")
stopifnot("Geometries differ from 2024.0.0" = all(same_geom))
cat("OK: geometries match the 2024.0.0 DB\n")

Run it on all the blocks we found to have altered ranges:

Rscript data-raw/pin/test_geometry_ranges.R ptaxsim-2024.0.0.db 1709419 1716220 1715103 1709452 1719211 0835201 1722501 1727501

Output:

PIN10s (old / new): 660 / 660
Ranges (old / new): 1790 / 1790
Ranges present in only one side:
# A tibble: 0 × 4
# ℹ 4 variables: side <chr>, start_year <int>, end_year <int>, n <int>
OK: ranges match the 2024.0.0 DB
Geometries identical: 1790 / 1790
OK: geometries match the 2024.0.0 DB

2. Confirm that sfarrow produces identical geometry columns

Save this script to data-raw/pin/qc_geometry_vs_2024.R:

Expand to view the QC script

# QC for the rebuilt pin_geometry_raw parquet file. Confirms that, relative to
# the 2024.0.0 DB:
#   1. Year ranges are identical, except that ranges ending in 2024 may be
#      extended to 2025 and PIN10s/shapes first seen in 2025 are new
#   2. Geometries are identical for every range that overlaps
#
# Usage: Rscript data-raw/pin/qc_geometry_vs_2024.R [new_parquet] [old_db]
library(dplyr)
library(DBI)
library(sf)

# Accept paths via command line, but fall back to the dev bucket and the
# local 2024 db path
args <- commandArgs(trailingOnly = TRUE)
new_path <- if (length(args) >= 1) {
  args[[1]]
} else {
  "s3://z-dev-ccao-data-ptaxsim-us-east-1/pin_geometry_raw/part-0.parquet"
}
old_db <- if (length(args) >= 2) args[[2]] else "ptaxsim-2024.0.0.db"
old_max_year <- 2024

# Old ranges and WKT from the 2024.0.0 DB
conn <- dbConnect(RSQLite::SQLite(), old_db)
old <- dbGetQuery(
  conn,
  "SELECT pin10, start_year, end_year, geometry AS wkt FROM pin_geometry_raw"
) %>%
  as_tibble() %>%
  mutate(across(c(start_year, end_year), as.integer))
dbDisconnect(conn)

# New ranges and WKT. WKT is generated the same way as in create_db.R so the
# strings are comparable. The sfarrow primary geometry warning is harmless
new <- suppressWarnings(sfarrow::st_read_parquet(new_path)) %>%
  mutate(wkt = st_as_text(geometry)) %>%
  st_drop_geometry() %>%
  select(pin10, start_year, end_year, wkt) %>%
  as_tibble() %>%
  mutate(across(c(start_year, end_year), as.integer))

# Cap new ranges at 2024 and drop ranges that start in 2025 or later, since
# those cover only new data
new_capped <- new %>%
  filter(start_year <= old_max_year) %>%
  mutate(end_year = pmin(end_year, old_max_year))

key <- c("pin10", "start_year", "end_year")
cat("Ranges (old / new / new capped at 2024):",
    nrow(old), "/", nrow(new), "/", nrow(new_capped), "\n")
cat("Duplicate keys (old / new):",
    sum(duplicated(old[key])), "/", sum(duplicated(new_capped[key])), "\n")

# 1. Year ranges
old_only <- anti_join(old, new_capped, by = key)
new_only <- anti_join(new_capped, old, by = key)
cat("Ranges only in old DB:", nrow(old_only), "\n")
cat("Ranges only in new file (pre-2025):", nrow(new_only), "\n")
if (nrow(old_only) > 0) print(head(old_only %>% select(all_of(key)), 20))
if (nrow(new_only) > 0) print(head(new_only %>% select(all_of(key)), 20))

# 2. Geometries for overlapping ranges
both <- inner_join(old, new_capped, by = key, suffix = c("_old", "_new"))
diff_geom <- both %>% filter(wkt_old != wkt_new)
cat("Overlapping ranges:", nrow(both), "\n")
cat("Overlapping ranges with different geometry:", nrow(diff_geom), "\n")
if (nrow(diff_geom) > 0) print(head(diff_geom %>% select(all_of(key)), 20))

# 3. Informational: what is new in 2025
cat("PIN10s only in new file:",
    n_distinct(anti_join(new, old, by = "pin10")$pin10), "\n")
cat("New-file ranges starting in 2025:",
    sum(new$start_year > old_max_year), "\n")

stopifnot(
  "Year ranges differ from 2024.0.0" = nrow(old_only) + nrow(new_only) == 0,
  "Geometries differ from 2024.0.0" = nrow(diff_geom) == 0
)
cat("OK: ranges and geometries match the 2024.0.0 DB for all pre-2025 data\n")

Run the script:

Rscript data-raw/pin/qc_geometry_vs_2024.R

Output:

Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.4.0; sf_use_s2() is TRUE
Ranges (old / new / new capped at 2024): 1471760 / 1473815 / 1471760 
Duplicate keys (old / new): 0 / 0 
Ranges only in old DB: 0 
Ranges only in new file (pre-2025): 0 
Overlapping ranges: 1471760 
Overlapping ranges with different geometry: 0 
PIN10s only in new file: 1699 
New-file ranges starting in 2025: 2055 
OK: ranges and geometries match the 2024.0.0 DB for all pre-2025 data

[jeancochrane]

While taking a look at this code, I noticed that the .by param was only added to tidyr::fill() in version 1.3.2, which is currently the latest version in the version list. Switching to the more-verbose group_by() makes this code compatible with a wider range of tidyr versions, which I think is good practice for a package like this that might be installed alongside a variety of different versions of its dependencies.

[jeancochrane]

Note that this read_sf_dataset() call raises a warning, which I think is ultimately benign:

In arrow_to_sf(tbl, geo) :
  Primary geometry column not found, using next available.

I think what's happening here is that geometry_3435 is set as the primary geometry column in the spatial.parcel parquet file metadata, but we're not selecting it in select(). Since it's missing, read_sf_dataset() warns that it can't find the primary geometry column, and it falls back to using the other geometry column that is present instead (i.e. geometry).

Quick script to confirm that geometry_3435 is the primary geometry column:

library(arrow)
ds <- open_dataset("s3://ccao-data-warehouse-us-east-1/spatial/parcel")
geo <- jsonlite::fromJSON(ds$metadata$geo)
print(geo$primary_column)

[1] "geometry_3435"

I thought about two ways of resolving the warning, neither of which seem exactly right to me:

• Suppress it with suppressWarnings(): This would prevent the warning from printing to the console, but might suppress other, more important warnings in the future, so I decided against it.
• Select geometry_3435 in the select() call, and then drop it after setting geometry as the geometry column: This would resolve the warning without suppressing other warnings, but it would roughly double the memory use of the script, since we would need to load both geometry and geometry_3435 into memory. This script already bumps up against the limits of the available memory on our biggest machines, so I don't think doubling memory use is a good idea.

Ultimately I think it's just best to leave the harmless warning as-is, but I'm open to pushback.

[wrridgeway]

You actually saw me fighting with this a bit in my first PR to replace geoarrow in the etl pipeline (scripts that assume geometry_3435 vs geometry will be the default geo column). It's really frustrating that our parcel spatial data is set to use geometry_3435 as the primary geometry column and no other spatial data is. I'd say we should fix this, but there would likely be a lot of subtle downstream effects that are hard to detect. I think it's fine to leave the warning, especially since sfarrow natively moves onto the only other remaining geo column.

[kyrasturgill]

I'm definitely ok with leaving as is!

[wrridgeway]

Thanks Jean! Glad my hunch was right (it seems?) and those diffs were just from the way I had to run this script locally on my laptop.

[wrridgeway]

You actually saw me fighting with this a bit in my first PR to replace geoarrow in the etl pipeline (scripts that assume geometry_3435 vs geometry will be the default geo column). It's really frustrating that our parcel spatial data is set to use geometry_3435 as the primary geometry column and no other spatial data is. I'd say we should fix this, but there would likely be a lot of subtle downstream effects that are hard to detect. I think it's fine to leave the warning, especially since sfarrow natively moves onto the only other remaining geo column.

[kyrasturgill]

Thanks for handling this! Excited to leave geoarrow behind us.