Skip to content

Replace geoarrow with sfarrow - #1116

Merged
wrridgeway merged 10 commits into
masterfrom
replace-geoarrow-with-sfarrow
Nov 3, 2026
Merged

wrridgeway merged 10 commits into
masterfrom
replace-geoarrow-with-sfarrow

Conversation

@wrridgeway

@wrridgeway wrridgeway commented Oct 27, 2026 •

Copy link
Copy Markdown
Member

I did significantly less testing for this PR than #1100, but I also made much smaller changes. sfarrow is not having any issues reading our current geoparquet files and handling them the same as geoarrow did as far as I can tell. I did rewrite the parcel geoparquet for 2025 in the dev bucket, as well as the neighborhood geoparquet for 2021. I rebuilt z_de_wridgeway_location.vw_pin10_location_fill and reran my testing script and got exactly the same results as in #1110. Script is pasted below. Note, this PR will not make the changes to the parcel files #1100 did that lead to the tax district changes shown below. Those are persisted from rebuilding raw parcel data in that PR and using a different source for attributes.

vw_pin10_location_fill

The ultimate test to see if using standard, non-spatial parquet files in athena produces the same spatial results as using geo-parquets made with geoarrow is to build our location db and check the new version of vw_pin10_location_fill against the old version.

library(arrow)
library(DBI)
library(dplyr)
library(noctua)
library(purrr)

noctua_options(unload = TRUE)
AWS_ATHENA_CONN_NOCTUA <- dbConnect(noctua::athena(), rstudio_conn_tab = FALSE)

old <- dbGetQuery(conn = AWS_ATHENA_CONN_NOCTUA, "
select * from location.vw_pin10_location_fill where year = '2025'
") %>%
  arrange(pin10)

new <- dbGetQuery(conn = AWS_ATHENA_CONN_NOCTUA, "
select * from z_dev_wridgeway_location.vw_pin10_location_fill where year = '2025'
") %>%
  arrange(pin10)

diffs <- names(old)[names(old) != "pin10"] %>%
  set_names() %>%
  map(function(x) {
  old <- old %>% select(pin10, "old" = all_of(x))
  new <- new %>% select(pin10, "new" = all_of(x))
  compare <- inner_join(old, new, by = "pin10") %>%
   filter(old != new | (!is.na(old) & is.na(new)) | (!is.na(new) & is.na(old)))

  if (nrow(compare) > 0) {
    return(compare)
  } else {
    return(NULL)
  }

  })

diffs <- diffs[lengths(diffs) != 0]
> names(diffs)
 [1] "tax_municipality_num"                "tax_municipality_name"               "combined_municipality_name"         
 [4] "tax_school_elementary_district_num"  "tax_school_elementary_district_name" "tax_school_secondary_district_num"  
 [7] "tax_school_secondary_district_name"  "tax_school_unified_district_num"     "tax_school_unified_district_name"   
[10] "tax_community_college_district_num"  "tax_community_college_district_name" "tax_fire_protection_district_num"   
[13] "tax_fire_protection_district_name"   "tax_library_district_num"            "tax_library_district_name"          
[16] "tax_park_district_num"               "tax_park_district_name"              "tax_sanitation_district_num"        
[19] "tax_sanitation_district_name"        "tax_special_service_area_num"        "tax_special_service_area_name"      
[22] "tax_tif_district_num"                "tax_tif_district_name"               "tax_data_year"

Comment on lines +9 to +18
# rJava's LD_LIBRARY_PATH from .Renviron isn't picked up by the dynamic
# linker in time, so preload libjvm.so by absolute path first.
dyn.load("/usr/lib/jvm/java-21-openjdk-amd64/lib/server/libjvm.so")
# tabulizer/RJavaTools reflects into JDK internals (e.g. java.util
# iterators), which Java 16+'s module system denies by default (JEP 396).
# Must be set before the JVM starts, i.e. before library(rJava).
options(java.parameters = c(
"--add-opens=java.base/java.util=ALL-UNNAMED",
"--add-opens=java.base/java.lang=ALL-UNNAMED"
))

@wrridgeway wrridgeway Oct 29, 2026 •

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I know you're iffy about this, but I tried tabulapdf and ran into the same exact issue. The good news is that this is one small file that we should never need to rebuilt, not something we have to do on a yearly basis. If you'd like, we could open a new PR to try and figure out some other solution that moves away from tabulizer and tabulapdf. I'd prefer not to get bogged down in fixing this problem in this PR.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Makes sense, I opened #1119 to track the effort to simplify this. I'll dig into it in November if I can find the time.

esri_chicago_geojson <- here(footprint_tmp_dir, "chicago-2008.geojson")
esri_chicago_remote_raw <- file.path(
input_bucket, "esri", "chicago-2008.parquet"
input_bucket, "esri", "chicago-2008.geojson"

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These files were incorrectly specified.

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm am leaving the big changes from #1100 for another PR.

Comment on lines +24 to 29
parcels <- sfarrow::read_sf_dataset(arrow::open_dataset(
paste0("s3://ccao-data-warehouse-us-east-1/spatial/parcel/year=", year)
))
)) %>%
select(-geometry) %>%
rename(geometry = geometry_3435)

@wrridgeway wrridgeway Oct 29, 2026 •

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm guessing it's from updating sf, but this script throws a ton of errors when we're dealing with two geometry columns. Since it immediately sets the parcels object to epsg:3435 for processing on line 30 (gone now) I'm just getting rid of the 4326 column right off the bat here. The raw parcel geoparquets are all saved with geometry_3435 as their active geometry column, so they ingest with a crs of epsg:3435 and renaming geometry_3435 means we don't have to edit all the spatial manipulation below.

@wrridgeway
wrridgeway marked this pull request as ready for review October 29, 2026 17:00
@wrridgeway
wrridgeway requested a review from a team as a code owner October 29, 2026 17:00
Comment thread etl/renv.lock

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We update a number of packages for compatibility with R 4.6.1.

@jeancochrane jeancochrane left a comment

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thank you for your persistence! I like that this change is much simpler than #1100, even though sfarrow is not actively maintained.

library(arrow)
library(ccao)
library(DBI)
library(geoarrow)

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[Question, non-blocking] Do you have any idea why this script doesn't write out files like the rest of the ETL raw scripts do? I assume it's not causing problems but I'm curious.

@wrridgeway wrridgeway Nov 3, 2026 •

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I never even noticed this. There's a comment in the script:

# This is static data and this script does not need to be re-run unless the data
# is no longer availale in the Data Department's raw S3 bucket.

but I don't even know what file in the raw bucket it produces. I'll open an issue to investigate it.

Comment on lines +9 to +18
# rJava's LD_LIBRARY_PATH from .Renviron isn't picked up by the dynamic
# linker in time, so preload libjvm.so by absolute path first.
dyn.load("/usr/lib/jvm/java-21-openjdk-amd64/lib/server/libjvm.so")
# tabulizer/RJavaTools reflects into JDK internals (e.g. java.util
# iterators), which Java 16+'s module system denies by default (JEP 396).
# Must be set before the JVM starts, i.e. before library(rJava).
options(java.parameters = c(
"--add-opens=java.base/java.util=ALL-UNNAMED",
"--add-opens=java.base/java.lang=ALL-UNNAMED"
))

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Makes sense, I opened #1119 to track the effort to simplify this. I'll dig into it in November if I can find the time.

@wrridgeway
wrridgeway merged commit e73d333 into master Nov 3, 2026
11 checks passed
@wrridgeway
wrridgeway deleted the replace-geoarrow-with-sfarrow branch November 3, 2026 15:50
jeancochrane added a commit to ccao-data/ptaxsim that referenced this pull request Nov 6, 2026
## 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
#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`:

<details>
<summary>Expand to view the QC script</summary>

```r
# 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")
```

</details>

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

```bash
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`:

<details>
<summary>Expand to view the QC script</summary>

```r
# 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")
```

</details>

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
```
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants