Repository navigation
PTAXSIM 2025 data upload - #107
kyrasturgill wants to merge 38 commits into
Conversation
| rename_with(~ rep("agency", length(.x)), any_of(c( | ||
| "mf_num" | ||
| ))) %>% | ||
| rename_with(~ rep("agency_name", length(.x)), any_of(c( | ||
| "auth_name" | ||
| ))) %>% | ||
| rename_with(~ rep("home_rule_ind", length(.x)), any_of(c( | ||
| "home_rule" | ||
| ))) %>% | ||
| rename_with(~ rep("percent_burden", length(.x)), any_of(c( | ||
| "b" | ||
| ))) %>% |
There was a problem hiding this comment.
These are the field name changes that need to be accounted for in the ingest
| agency_fund <- agency_fund %>% | ||
| semi_join(agency, by = c("year", "agency_num")) |
There was a problem hiding this comment.
There were a few former agencies (now funds) included in the agency fund:
- 020090003, 408000 (ROAD AND BRIDGE HANOVER)
-
- 408000 (GENERAL ASSISTANCE PROVISO)
- 020290003, 408000 (ROAD AND BRIDGE WHEELING)
- 030210001, 003001 (CITY OF CHICAGO LIBRARY FUND)
- 050200001, 003005 (CHICAGO PARK DIST AQUARIUM & MUSEUM BOND)
This code ensures any defunct agencies that are not in the agency table for a certain year are removed. This also requires bumping finalizing and writing agency_fund to after agency table is created and finalized.
There was a problem hiding this comment.
[Suggestion, optional] Good call! I wonder if we should use a warning() call to alert the script runner to any agency funds that are being filtered out, so we have a chance to investigate unexpected mismatches in future runs? If we do this, we should probably also create a list of keys for known mismatches and exclude those from the warning, so that we don't pollute the warning with mismatches that we've already investigated.
There was a problem hiding this comment.
Added warning! Lmk if it is what you were thinking on your next review (disclaimer an LLM wrote it almost in entirety).
| agency_num varchar(9) NOT NULL, | ||
| tax_code_num varchar(5) NOT NULL, | ||
| tax_code_rate double CHECK(tax_code_rate >= 0) NOT NULL, | ||
| tax_code_rate double CHECK(tax_code_rate >= 0) , |
There was a problem hiding this comment.
For some reason this field is blank for a few tax codes. Ultimately this table is no longer used for any functions in 2025, and it's not the primary place to find a tax code rate, which is why I didn't build a workaround to fill these.
There was a problem hiding this comment.
I'm now leaning towards following same protocol as described above, leaving in the NOT NULL requirement and including this as one of the data quality issues in our request to the Clerk's office.
| AND pct_burden <= 1) | ||
| AND (year >= 2024 OR pct_burden IS NOT NULL)), | ||
| total_levy bigint CHECK(total_levy >= 0) NOT NULL, | ||
| total_levy bigint CHECK(total_levy >= 0) , |
There was a problem hiding this comment.
This field no longer exists in the 2025 report which renders it NULL for 2025. The alternative here would be removing this field altogether.
There was a problem hiding this comment.
[Thought, non-blocking] Hmm, without more context it's hard for me to tell whether it would be better to remove this constraint or remove the field. I lean toward removing the field if it seems like the Clerk won't be publishing it in the future, but I trust your judgment about what the most useful choice will be for us and our users.
There was a problem hiding this comment.
I think you're correct that removing the field is better call!
There was a problem hiding this comment.
I believe this should be resolved now with the renaming the field!
jeancochrane
left a comment
There was a problem hiding this comment.
Thank you for this excellent work! I'm going to mess around with the data and the tests more extensively on Monday, but I wanted to get a first pass at review over to you ASAP in the meantime.
| mutate( | ||
| all_zero = rowSums(across(where(is.numeric)) != 0, na.rm = TRUE) == 0 | ||
| ) %>% | ||
| add_count(year, agency_num, fund_num, name = "n_dupe_key") %>% | ||
| filter(!(n_dupe_key > 1 & all_zero)) %>% |
There was a problem hiding this comment.
[Question, non-blocking] Do you think it's possible that we might ever see a fund that has zero-value duplicate rows, but where the numeric columns are all 0 for all rows, not just the dupes? In that case, this code would remove all the rows, leaving no representation of the fund in the data.
This may be fine behavior, or it may be a situation that you're confident can never occur, so I don't feel strongly about the correct fix; same as above, I leave it up to you to decide if you want to handle this edge case.
There was a problem hiding this comment.
Ah good point. I think it's unlikely, but I can't be certain. I wonder if we remove this workaround and add this to our list of requested fixes in the data quality?
| ) | ||
| !is.na(levy_plus_loss), | ||
| # Remove rows added in 2025 for ALL FUND TYPES | ||
| fund_num != "000000" |
There was a problem hiding this comment.
[Question, non-blocking] Is there any chance that 00000 might be a valid fund number in the future? If so, we could consider adding a year = '2025' clause to this filter condition to protect ourselves from accidentally stripping out a valid fund in future years, but I trust your judgement regarding the likelihood of future valid 00000 funds.
There was a problem hiding this comment.
I highly doubt it and think my preference would be to remove it for any year. Though they are inconsistent with their data I feel like I cannot foresee them using '00000' for a fund.
| agency_fund <- agency_fund %>% | ||
| semi_join(agency, by = c("year", "agency_num")) |
There was a problem hiding this comment.
[Suggestion, optional] Good call! I wonder if we should use a warning() call to alert the script runner to any agency funds that are being filtered out, so we have a chance to investigate unexpected mismatches in future runs? If we do this, we should probably also create a list of keys for known mismatches and exclude those from the warning, so that we don't pollute the warning with mismatches that we've already investigated.
| AND pct_burden <= 1) | ||
| AND (year >= 2024 OR pct_burden IS NOT NULL)), | ||
| total_levy bigint CHECK(total_levy >= 0) NOT NULL, | ||
| total_levy bigint CHECK(total_levy >= 0) , |
There was a problem hiding this comment.
[Thought, non-blocking] Hmm, without more context it's hard for me to tell whether it would be better to remove this constraint or remove the field. I lean toward removing the field if it seems like the Clerk won't be publishing it in the future, but I trust your judgment about what the most useful choice will be for us and our users.
|
@kyrasturgill I did some QC today and identified a few small issues, with suggested fixes! To make it easier to rerun the code, I encapsulated my findings in a Quarto doc. Here's the rendered output: qc_2025_issues.html Expand to see the Quarto code---
title: "PTAXSIM 2025.0.0 data QC: issues to review"
date: today
format:
html:
toc: true
toc-depth: 2
code-fold: true
code-summary: "Show code"
embed-resources: true
df-print: kable
execute:
warning: false
message: false
params:
old_db: "ptaxsim-2024.0.0.db"
new_db: "ptaxsim-2025.0.0.db"
new_year: 2025
---
This document summarizes issues found while QC'ing the 2025 data year added
to PTAXSIM in the `2025.0.0` database. Note that I relied heavily on gen AI
to help me draft this doc, but I reviewed everything for correctness and
copyedited the prose, so everything here represents something I actually think.
Comparisons use the published `2025.0.0` and `2024.0.0` databases, which I
assume you have already downloaded from
`https://ccao-data-public-us-east-1.s3.amazonaws.com/ptaxsim/ptaxsim-{version}.db.bz2`.
Decompress them into the repo root before rendering. Raw Clerk files are read
from `data-raw/`.
```{r setup}
suppressPackageStartupMessages({
library(DBI)
library(dplyr)
library(ggplot2)
library(glue)
library(purrr)
library(readxl)
library(snakecase)
library(stringr)
library(tidyr)
})
options(scipen = 99)
new_year <- params$new_year
prev_year <- new_year - 1
conn <- dbConnect(RSQLite::SQLite(), params$new_db, bigint = "numeric")
invisible(dbExecute(conn, glue("ATTACH DATABASE '{params$old_db}' AS old")))
# Helper function to run a query against the DB
q <- function(sql, .envir = parent.frame()) {
dbGetQuery(conn, glue(sql, .envir = .envir))
}
# Helper functions for data formatting
fmt_dollar <- scales::label_dollar(accuracy = 1)
fmt_num <- scales::label_comma(accuracy = 1)
fmt_pct <- scales::label_percent(accuracy = 0.1)
```
# Summary
| \# | Issue | Tables | Severity |
|----|--------------------------------------------------------------|--------------------------------|--------------------------------------|
| 1 | `total_levy`, null in 2025, may use a different column | `agency` | Bug, with a proposed fix |
| 2 | Non-numeric tax code `RR` | `tax_code` | Warning |
| 3 | 5,048 PINs with property class `RR` | `pin` | Warning |
| 4 | Vacant/exempt class changed from `0` to `000` | `pin` | Warning |
| 5 | Large year-over-year changes | `agency`, `agency_fund`, `pin` | Question |
| 6 | 2024 veteran exemption values changed for 7 PINs | `pin` | Unexpected change to prior-year data |
| 7 | 491 pre-2025 PIN10-years lost geometry | `pin_geometry_raw` | Unexpected change to prior-year data |
# 1. `agency.total_levy` is NULL for all 2025 agencies, but may use a different column
In the 2025 Agency Rate Report, the Overview column "Grand Total Levy" appears
to have been replaced by a new column "Grand Total Levy+Loss". `agency.R` turns
this new column into `total_levy_loss` and then drops it
(`select(-cty_overall_eav, -total_levy_loss)`). As a result, `total_levy` is
missing for every 2025 agency. To account for these null values, we relaxed the
`NOT NULL` constraint on `agency.total_levy` in `create_db.sql`.
```{r total-levy-nulls}
q("
SELECT year, COUNT(*) AS n_agencies,
SUM(total_levy IS NULL) AS n_null_total_levy
FROM main.agency
WHERE year >= {new_year - 4}
GROUP BY year
")
```
However, it appears that "Grand Total Levy+Loss" may be the same column as
"Grand Total Levy", just renamed. In every year, the Overview grand total
equals the sum of the Detail sheet's fund-level levy+loss, rather than the plain
fund levy. The subtotal headers support this too: the 2024 file already labels
the non-capped subtotal `Total Non-Cap Levy+Loss`, even though its grand total
is still labeled `Grand Total Levy`.
```{r total-levy-raw, cache=TRUE}
first_present <- function(df, candidates) intersect(candidates, names(df))[1]
agency_files <- list.files("data-raw/agency", "\\.xlsx?$", full.names = TRUE)
raw_levy <- map_dfr(agency_files, function(file) {
ov <- read_xlsx(file, sheet = 1) %>% set_names(to_snake_case(names(.)))
gt_col <- first_present(ov, c("grand_total_levy", "grand_total_levy_loss"))
# 2013 has no grand total columns
if (is.na(gt_col)) {
return(NULL)
}
dt <- read_xlsx(file, sheet = 2) %>% set_names(to_snake_case(names(.)))
levy_col <- first_present(dt, c("levy", "fund_levy"))
lpl_col <- first_present(
dt, c("levy_and_loss", "fund_levy_plus_loss", "levy_loss")
)
dt_id <- first_present(dt, c("agency", "agency_num"))
ov_id <- first_present(ov, c("agency", "agency_num", "mf_num"))
year_col <- first_present(ov, c("tax_year", "year"))
funds <- dt %>%
transmute(
agency_num = str_pad(as.character(.data[[dt_id]]), 9, "left", "0"),
fund = as.character(fund),
levy = as.numeric(.data[[levy_col]]),
levy_plus_loss = as.numeric(.data[[lpl_col]])
) %>%
# 2025 adds a "TOTAL OF ALL FUNDS" row per agency
filter(!fund %in% c("0", "000000")) %>%
group_by(agency_num) %>%
summarize(
fund_levy = sum(levy, na.rm = TRUE),
fund_levy_plus_loss = sum(levy_plus_loss, na.rm = TRUE)
)
ov %>%
transmute(
year = as.integer(.data[[year_col]]),
overview_column = gt_col,
agency_num = str_pad(as.character(.data[[ov_id]]), 9, "left", "0"),
grand_total = as.numeric(.data[[gt_col]])
) %>%
inner_join(funds, by = "agency_num")
})
raw_levy %>%
group_by(year, overview_column) %>%
summarize(
n_agencies = n(),
pct_eq_fund_levy_plus_loss = mean(abs(grand_total - fund_levy_plus_loss) <= 2),
pct_eq_fund_levy = mean(abs(grand_total - fund_levy) <= 2),
.groups = "drop"
) %>%
mutate(across(starts_with("pct_"), fmt_pct))
```
**Proposed fix:** in the Overview rename block of `data-raw/agency/agency.R`,
map the new header onto the old one, before the
`starts_with("grand_total_")` select:
```r
rename_with(~"grand_total_levy", any_of("grand_total_levy_loss")) %>%
```
Then remove `-total_levy_loss` from the `select()` call near the end of the
agency section, and restore `NOT NULL` on `agency.total_levy` in
`data-raw/create_db.sql`. After the rebuild, 2025 `total_levy` should equal
the per-agency sum of `agency_fund.levy_plus_loss`, plus or minus a few dollars
(`agency.R` truncates the Clerk's fractional fund values to integers, so the
sum can fall a few dollars short).
# 2. Non-numeric tax code `RR` in `tax_code`
The 2025 Clerk tax code file (`2025-tax-code-agency-rates.xlsx`) contains a
railroad tax code `RR`, while no earlier year's file has one. No PIN is assigned
to tax code `RR`, so it may not affect our users, but it could cause problems
for any code that assumes tax codes only contain numeric characters. This tax
code also has a maximum rate that seems very high to me (927%) but I might just
be misunderstanding how tax rates work.
```{r rr-tax-code}
q("
SELECT year, tax_code_num, COUNT(*) AS n_agency_rows,
MAX(tax_code_rate) AS tax_code_rate,
(SELECT COUNT(*) FROM main.pin p
WHERE p.year = t.year AND p.tax_code_num = t.tax_code_num) AS n_pins
FROM main.tax_code t
WHERE tax_code_num NOT GLOB '[0-9][0-9][0-9][0-9][0-9]'
GROUP BY year, tax_code_num
")
```
```{r rr-tax-code-mean}
q("
SELECT year,
AVG(tax_code_rate) AS mean_rate_all,
AVG(CASE WHEN tax_code_num GLOB '[0-9][0-9][0-9][0-9][0-9]'
THEN tax_code_rate END) AS mean_rate_excl_non_numeric
FROM main.tax_code
WHERE year IN ({prev_year}, {new_year})
GROUP BY year
")
```
**Suggested fix:** Investigate this new `RR` tax code and document it in the
2025 changelog.
# 3. 5,048 PINs with property class `RR`
In 2025, 5,048 PINs have class `RR`, presumably railroad property. Only 25 of
these PINs were in the 2024 `pin` table, where they had ordinary classes; the
rest are new. These PINs have normal numeric tax codes (not `RR`) and real
tax bills. Every other class in PTAXSIM is numeric, so `RR` could break
downstream code that assumes numeric classes, such as `substr(class, 1, 1)`
to get the major class.
```{r rr-class}
q("
SELECT year, COUNT(*) AS n_pins,
COUNT(DISTINCT tax_code_num) AS n_tax_codes,
SUM(av_clerk) AS total_av_clerk,
ROUND(SUM(tax_bill_total)) AS total_tax_bill
FROM main.pin
WHERE class = 'RR'
GROUP BY year
")
```
```{r rr-class-prev}
q("
SELECT p.class AS class_{prev_year}, COUNT(*) AS n_pins
FROM main.pin p
JOIN main.pin n ON p.pin = n.pin AND n.year = {new_year}
WHERE p.year = {prev_year} AND n.class = 'RR'
GROUP BY 1
ORDER BY 2 DESC
")
```
**Suggested fix:** Investigate `RR` class PINs to make sure they are correct,
particularly ones that changed class in 2025, and document the new class in
the 2025 changelog.
# 4. Vacant/exempt class changed from `0` to `000`
From 2006 through 2024, vacant/exempt PINs have class `'0'`. In 2025 the same
PINs have class `'000'`. Code that filters on `class == "0"` will silently drop
these PINs for 2025.
```{r class-zero}
q("
SELECT year, class, COUNT(*) AS n_pins
FROM main.pin
WHERE class IN ('0', '00', '000') AND year >= {new_year - 3}
GROUP BY year, class
ORDER BY year
")
```
**Suggested fix:** Decide whether to coerce `000` to `0` in 2025, or document
the discrepancy in the changelog.
# 5. Large year-over-year changes
These aren't necessarily errors, but they are large enough changes from
2024 to 2025 that we should confirm them against the source data.
## 5.1. `agency.curr_new_prop` dropped 64%
New property totals are down by quite a lot from 2024. This may be correct,
since a historical view reveals that new property fluctuates a lot
year-to-year. (Get a load of that 2007 spike!) Still, if the 2025 numbers are
correct, then new property is at its lowest level since 2016.
```{r curr-new-prop, fig.height=3}
new_prop <- q("
SELECT year, SUM(curr_new_prop) AS curr_new_prop
FROM main.agency
GROUP BY year
")
ggplot(new_prop, aes(x = as.integer(year), y = curr_new_prop)) +
geom_col() +
scale_y_continuous(labels = scales::label_dollar(scale = 1e-9, suffix = "B")) +
labs(x = "Year", y = "Sum of curr_new_prop") +
theme_minimal()
```
## 5.2. MWRD bond fund `prelim_rate`
Not exactly a year-over-year change, but feels similar: One MWRD bond fund
has a preliminary rate of about 468, but its final rate is about 0.05.
Is this a source data issue, or do I just not understand how preliminary
rates work?
```{r prelim-rate}
q("
SELECT f.year, f.agency_num, i.agency_name, f.fund_num, fi.fund_name,
f.prelim_rate, f.final_rate, f.levy, f.final_levy
FROM main.agency_fund f
LEFT JOIN main.agency_info i ON f.agency_num = i.agency_num
LEFT JOIN main.agency_fund_info fi
ON f.agency_num = fi.agency_num AND f.fund_num = fi.fund_num
WHERE f.year IN ({prev_year}, {new_year})
ORDER BY f.prelim_rate DESC
LIMIT 5
")
```
## 5.3. Exemption totals
Several exemption totals moved sharply from 2024 to 2025. For example,
returning veterans fell 92% and WWII rose 113%. It's possible that some of
this may come from the way that `pin.R` splits the tax roll's total veteran
disability exemption into types using Athena values, so we should confirm that
these changes are correct.
```{r exemptions}
exe_cols <- q("PRAGMA main.table_info(pin)") %>%
filter(str_starts(name, "exe_")) %>%
pull(name)
exe_sums <- q("
SELECT year, {paste0('SUM(', exe_cols, ') AS ', exe_cols, collapse = ', ')}
FROM main.pin
WHERE year >= {new_year - 3}
GROUP BY year
")
exe_sums %>%
pivot_longer(-year, names_to = "exemption") %>%
pivot_wider(names_from = year) %>%
mutate(
pct_change = .data[[as.character(new_year)]] /
.data[[as.character(prev_year)]] - 1
) %>%
arrange(desc(abs(pct_change))) %>%
mutate(
across(matches("^\\d{4}$"), fmt_num),
pct_change = fmt_pct(pct_change)
)
```
# 6. 2024 veteran exemption values changed in the 2025.0.0 DB
`pin.R` was set to process only 2025 (`start_year <- 2025`), so the 2024 PIN
data should be identical in both releases. However, 7 PINs have veteran
disability exemption amounts that moved between buckets, for example from
`exe_vet_dis_lt50` to `exe_vet_dis_ge70`. Did we rerun 2024 at some point?
```{r vet-changes}
vet_cols <- c(
"exe_vet_dis_lt50", "exe_vet_dis_50_69", "exe_vet_dis_ge70",
"exe_vet_dis_100", "exe_wwii"
)
vet_where <- paste0("o.", vet_cols, " != n.", vet_cols, collapse = " OR ")
vet_select <- paste0(
"o.", vet_cols, " AS ", vet_cols, "_old, n.", vet_cols, " AS ", vet_cols, "_new",
collapse = ", "
)
vet_changes <- q("
SELECT o.pin, {vet_select}
FROM old.pin o
JOIN main.pin n ON o.year = n.year AND o.pin = n.pin
WHERE o.year = {prev_year} AND ({vet_where})
")
# Show only the columns that actually changed, in long format
vet_changes %>%
pivot_longer(
-pin,
names_to = c("exemption", ".value"),
names_pattern = "(.*)_(old|new)$"
) %>%
filter(old != new)
```
# 7. Pre-2025 PIN-years lost geometry coverage
`pin_geometry_raw` stores each parcel shape with a `(start_year, end_year)`
validity range. `pin.R` rebuilds this whole table on every run from the
warehouse's yearly parcel data (`spatial.parcel`, 2006 onward), collapsing
consecutive years with the same shape into one range. That means any change
to past years of upstream parcel data also changes past years of PTAXSIM
geometry, even when only one new year of data is meant to be added.
Extending `end_year` from 2024 to 2025 for unchanged parcels is expected.
However, 491 PIN10-years from before 2025 that had a geometry in `2024.0.0` have
none in `2025.0.0`. In every case the shape itself is unchanged. What changed
is that the PIN's first range now starts later, by one year in almost every
case, and in a few cases the range also ends earlier. Almost all of the lost
PIN10-years (487) are from tax year 2009 in downtown Chicago. 469 of them are
on a single block, `17-09-419`, and nearly all are class 592 commercial
property.
I spot-checked a few of the missing PIN10-years in the `spatial.parcel` table,
and it appears that they all have rows for the missing years, but those rows are
missing key attributes like `town_code` and `nbhd_code` (and they have
`has_attributes` set to `FALSE`). We should check with Billy to figure out
exactly what's going on here, but in the meantime, I'm wondering if perhaps the
null town/neighborhood codes are somehow causing these rows to get dropped from
the `pin_geometry` generation script.
```{r geometry, cache=TRUE}
invisible(dbExecute(conn, "
CREATE TEMP TABLE IF NOT EXISTS geo_missing AS
SELECT o.pin10, o.start_year, o.end_year
FROM old.pin_geometry_raw o
LEFT JOIN main.pin_geometry_raw n
ON o.pin10 = n.pin10 AND o.start_year = n.start_year
WHERE n.pin10 IS NULL
"))
invisible(dbExecute(conn, glue("
CREATE TEMP TABLE IF NOT EXISTS geo_pin_years AS
SELECT DISTINCT SUBSTR(pin, 1, 10) AS pin10, year
FROM main.pin
WHERE year < {new_year}
AND SUBSTR(pin, 1, 10) IN (SELECT pin10 FROM temp.geo_missing)
")))
# PIN10-years from the pin table that have a geometry in the old DB but not
# the new DB
invisible(dbExecute(conn, "
CREATE TEMP TABLE IF NOT EXISTS geo_lost AS
SELECT p.pin10, p.year
FROM temp.geo_pin_years p
WHERE EXISTS (
SELECT 1 FROM old.pin_geometry_raw g
WHERE g.pin10 = p.pin10 AND p.year BETWEEN g.start_year AND g.end_year
) AND NOT EXISTS (
SELECT 1 FROM main.pin_geometry_raw g
WHERE g.pin10 = p.pin10 AND p.year BETWEEN g.start_year AND g.end_year
)
"))
geo_lost <- q("SELECT * FROM temp.geo_lost ORDER BY year, pin10")
# For each lost PIN10-year, compare the old range that covered it to the
# PIN's first range in the new DB
geo_ranges <- q("
SELECT l.pin10, l.year AS lost_year,
o.start_year AS old_start_year, o.end_year AS old_end_year,
n.start_year AS new_start_year, n.end_year AS new_end_year,
o.geometry = n.geometry AS same_geometry
FROM temp.geo_lost l
JOIN old.pin_geometry_raw o
ON o.pin10 = l.pin10 AND l.year BETWEEN o.start_year AND o.end_year
JOIN main.pin_geometry_raw n
ON n.pin10 = l.pin10
AND n.start_year = (
SELECT MIN(start_year) FROM main.pin_geometry_raw g
WHERE g.pin10 = l.pin10
)
")
geo_classes <- q("
SELECT l.year, SUBSTR(l.pin10, 1, 7) AS block, p.class
FROM temp.geo_lost l
JOIN main.pin p
ON p.year = l.year AND SUBSTR(p.pin, 1, 10) = l.pin10
")
```
Lost PIN10-years by year and block:
```{r geometry-by-block}
geo_lost %>%
mutate(block = str_sub(pin10, 1, 7)) %>%
count(year, block, name = "n_pin10s") %>%
arrange(desc(n_pin10s))
```
Classes of the affected PINs:
```{r geometry-classes}
count(geo_classes, year, class, name = "n_pins")
```
How the affected ranges changed. The geometry is identical in every case:
```{r geometry-ranges}
geo_ranges %>%
count(
lost_year, old_start_year, old_end_year, new_start_year, new_end_year,
same_geometry = as.logical(same_geometry),
name = "n_pin10s"
) %>%
arrange(desc(n_pin10s))
```
All affected PIN10-years:
```{r geometry-table}
DT::datatable(geo_ranges, rownames = FALSE, options = list(pageLength = 10))
```
```{r teardown, include=FALSE}
dbDisconnect(conn)
``` |
|
@jeancochrane thank you for the additional QC! Here are my responses to your main finds and next steps:
|
|
Thanks for the detailed response Kyra, that all sounds great to me! It seems like I will be most helpful investigating the lost geometry coverage (issue 8) so Billy and I will coordinate on figuring that out when we redo the spatial ingest. |
|
OK upon further investigation of RR PINs, it actually seems that last year's provided tax roll excluded these while prior years classed them as 0. Because of that, I'm leaning towards adding another exception to reclassify RR as 0. |
| # A PIN that has not existed since 2018 shows up in the 2025 tax roll with | ||
| # NA values for av_mailed, av_certified, av_board which causes a test failure | ||
| filter(pin != "28224230088002") |
There was a problem hiding this comment.
I realized when running create_db.R again that this PIN was still causing issues - I can't recall the temporary workaround I did before to get it to run but ultimately I think it just needs to be removed based on the fact that we have no record of it in our parcel universe since 2018.
## 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 ```
This PR handles all of the data uploads for 2025 PTAXSIM data, as well as several minor changes to the ingest code to account for some new data oddities introduced in 2025.
Data Ingest Changes
agency.Ragency_fundtableData file changes