[wrridgeway]

This PR swaps out geoarrow and replaces it with the base arrow package. What this means is that our parquet files in s3 will no longer have spatial metadata where relevant. My solution to this is to make sure that geometry columns and their EPSG identifiers are all recorded for any spatial datasets written as parquets. We can use this data when reading them back in to make sure we are handling spatial data correctly. All of the reading and writing updates are handled in utils.R.

Companion aws-infrastructure PR

Before merging

• delete all _test paths from scripts and raw bucket
• revert .Renviron
• revert spatial schema.yml
• revert non-geoarrow changes to scripts-ccao-data-warehouse-us-east-1/spatial/spatial-environment-road.R
• delete all local geojson and parquet files

After merging

• move data from s3://ccao-data-warehouse-dev-us-east-1/spatial/ into s3://ccao-data-warehouse-us-east-1/spatial/
• replace old raw data with new raw data marked by _test suffix data

Scripts that don't work

Couldn't finish running due to OSM timeouts:

• scripts-ccao-data-raw-us-east-1/spatial/spatial-environment-major_road.R
• scripts-ccao-data-raw-us-east-1/spatial/spatial-environment-secondary_road.R

Are not reproducing old data:

• scripts-ccao-data-warehouse-us-east-1/spatial/spatial-ccao-corner.R
• scripts-ccao-data-warehouse-us-east-1/spatial/spatial-ccao-neighborhood.R

In the case of corners, the new output is not as good as the old output. For neighborhoods, it seems like there's a genuine improvement. I'm still not sure why this is happening, but given the very minor code changes I made to the algorithms that produce these outputs, I am leaning towards this being an issue with updated dependencies (sf, rmapshaper). It's obviously worth getting to the bottom of this, but that work can be done after we've merged in the changes made by this PR.

This does leave one important issue to resolve before merging, though: we will essentially lift and shift the spatial folder in the dev warehouse bucket into our prod warehouse bucket after we merge this, to make sure we aren't using parquet files with attributes that can't be read by the base arrow package without throwing an error. The updated corner lots file is not in a good place, so we'll need to decide how to handle corners and major/secondary osm roads.

1. don't update them: we're not using these features in modeling and we rarely touch warehouse files when using our etl pipeline, leaving them as is won't cause any immediate issues.
2. replace corners with the new file: this will mean the corner lots parquet files are fully integrated into our new etl pipeline, but the data in athena will degrade and we will need to make sure we maintain the old data somewhere for working through what is causing this issue later.
3. convert the old files: rather than rebuilding the files from scratch without geoarrow, we could read them into R with that pakage and then rewrite them out using the new pipeline. this would allow us to preserve the old data while also making it compliant with the new etl pipeline. This will mean what's in athena is not a product of the new pipeline, but that should be fine while we work out what's causing the difference.

Testing

Data generated without any spatial code changes other than swapping arrow for geoarrow

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"

We can see almost all the differences are in tax districts. The source of this divergence is the parcel shapefiles. This is also the source of divergence for combined_municipality_name since it compares municpality names derived from tax districts to those from spatial joins. If we compare tax codes in the spatial parcel files, it appears the new files are significantly more accurate than the older files:

-- Use ptaxsim as our "truth"
with tax as (
	select distinct substr(pin, 1, 10) as pin10,
		year,
		tax_code_num as truth
	from tax.pin
),
diffs as (
	select new.pin10,
		new.year,
		new.tax_code as new,
		old.tax_code as old,
		pins.truth as truth,
		case
			when new.tax_code != pins.truth
			and old.tax_code != pins.truth then 'both'
			-- Assume null tax codes are incorrect
			when (
				(
					new.tax_code != pins.truth
					or new.tax_code is null
				)
				and old.tax_code = pins.truth
			) then 'new'
			when (
				(
					old.tax_code != pins.truth
					or old.tax_code is null
				)
				and new.tax_code = pins.truth
			) then 'old'
		end as which_is_wrong
	from z_dev_wridgeway_spatial.parcel as new
		left join spatial.parcel as old on new.pin10 = old.pin10
		and new.year = old.year
		inner join tax as pins on new.pin10 = pins.pin10
		and new.year = pins.year
	where (new.tax_code != old.tax_code)
		or (
			new.tax_code is null
			and old.tax_code is not null
		)
		or (
			new.tax_code is not null
			and old.tax_code is null
		)
)
select which_is_wrong,
	count() as parcels
from diffs
group by which_is_wrong

WhichIsWrong	Parcels
old	530127
new	416
both	90

This is almost certainly a product of me changing the way the attributes files are constructed in the raw ingest file for the spatial parcel files.

Rewrites to avoid errors while re-running

Neighborhhods

There are a lot of very small differences, primarily along the borders of the neighborhoods. We're talking about essentially zero-width lines. There are other, larger differences - I've chosen a few neighborhoods where these differences appeared and included them below in the check_new and check_old objects. To me, it looks like the new file is doing an objectively better job, for whatever reason. The edits I made to the script didn't change our algorithm, they just handled an error that popped up from trying to negative buffer empty geometries. My only guesses as to why these differences are appearing is that I used different syntax to make absolutely sure we're actually using a 3435 projection during all the buffering or upgrading sf or rmapshaperaffected some of the buffering.

library(dplyr)
library(DBI)
library(noctua)
library(sf)
library(mapview)

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

# Grab old and new neighborhood geometries for 2021
old <- dbGetQuery(conn = AWS_ATHENA_CONN_NOCTUA, "
select town_nbhd, geometry from spatial.neighborhood where year = '2021'
")

new <- dbGetQuery(conn = AWS_ATHENA_CONN_NOCTUA, "
select town_nbhd, geometry from z_dev_wridgeway_spatial.neighborhood where year = '2021'
")

# Convert to sf objects
old_shape <- st_as_sf(old, crs = 4326)
new_shape <- st_as_sf(new, crs = 4326)

old_p <- old_shape %>% st_transform(3435) %>% st_make_valid() %>% arrange(town_nbhd)
new_p <- new_shape %>% st_transform(3435) %>% st_make_valid() %>% arrange(town_nbhd)

go <- st_geometry(old_p)
gn <- st_geometry(new_p)

# Function to compute differences between two sets of geometries
diff_geo <- function(a, b, label) {
  g <- lapply(seq_along(a), function(i) {
    d <- st_difference(a[i], b[i])
    if (length(d) == 0) st_geometrycollection() else d[[1]]
  })
  st_sf(town_nbhd = old_p$town_nbhd, change = label, geometry = st_sfc(g, crs = 3435)) %>%
    filter(!st_is_empty(geometry))
}

diffs <- bind_rows(
  diff_geo(go, gn, "In old only"),
  diff_geo(gn, go, "In new only")
) %>%
  # drop sliver noise (< 1 sq ft)
  filter(as.numeric(st_area(geometry)) > 1)

in_old_only <- diffs %>% filter(change == "In old only")
in_new_only <- diffs %>% filter(change == "In new only")

# Inspect neighborhood parcels where large differences exist
check_new <- dbGetQuery(conn = AWS_ATHENA_CONN_NOCTUA, "
select geometry from z_dev_wridgeway_spatial.parcel where year = '2021' and town_code || nbhd_code IN ('70300', '76012', '16080', '38082')
") %>%
  st_as_sf(crs = 4326)

check_old <- dbGetQuery(conn = AWS_ATHENA_CONN_NOCTUA, "
select geometry from spatial.parcel where year = '2021' and town_code || nbhd_code IN ('70300', '76012', '16080', '38082')
") %>%
  st_as_sf(crs = 4326)

# Plot everything together
mapview(old_shape, col.regions = "blue", alpha.regions = 0.5) + 
  mapview(new_shape, col.regions = "red", alpha.regions = 0.5) +
  mapview(check_new, col.regions = "green", alpha.regions = 0.5) +
  mapview(check_old, col.regions = "purple", alpha.regions = 0.5) +
  mapview(in_old_only, col.regions = "blue", alpha.regions = 0.5) +
  mapview(in_new_only, col.regions = "red", alpha.regions = 0.5)

I thought this might be due to changes in neighborhood #'s between the new and old spatial.parcel tables but that is not the case (though there are some minor differences):

select old.pin10, old.year,
new.town_code || new.nbhd_code as new,
old.town_code || old.nbhd_code as old
from spatial.parcel as old
	full outer join z_dev_wridgeway_spatial.parcel as new on new.pin10 = old.pin10
	and new.year = old.year
where new.town_code || new.nbhd_code != old.town_code || old.nbhd_code
and old.year = '2021'

returns no rows

And we can see from comparing the old and new parcels for these neighborhoods that that's also not the issue.

Issues

Corner Lots

These are not generating as they did in the past. We can see that here:

The script now filters out parcels it wouldn't have in the past during filter(max_degree == 90 | n() >= 3, na.rm = TRUE) because some of these angles are now less than 90. I'm still trying to figure out why this is the case.

[wrridgeway]

This is the only real change to the actual processing here. We project before buffering so sf can use euclidean geometry. This greatly speeds up processing.

[wrridgeway]

Just making sure we built my dev location tables and views using my dev spatial tables.

[wrridgeway]

We typically don't write data to our raw bucket as parquet files, but where we do I used a _test suffix not to overwrite production data.

[wrridgeway]

Fixing a previous oversight.

[jeancochrane]

[Question, non-blocking] How did you determine that the IDOT shapefiles use CRS 3436?

[wrridgeway]

I am moving us away from the old attributes source since our current system of record seems to agree much better with tax district data from other county sources.

[wrridgeway]

Updated with deactivat and cur handling.

[wrridgeway]

Switching over to town code as it's coded in ias.

[wrridgeway]

copilot helped me sort out some issues with rJava being a PitA. We need it to use tabulizer.

[jeancochrane]

[Suggestion, non-blocking] Hmm, this makes me a bit nervous: The absolute path to libjvm.so feels like it could easily change in the future or be different between machines, and the need to set flags to disable a feature of Java 16+ module system feels brittle. I wonder if it's worth trying to switch to the newer tabulapdf package before we persist these changes, to see if it resolves these problems? Part of my hesitation stems from a lack of full understanding of the problems these two settings are intending to address, so if you feel confident in your understanding and think it's the right call, I trust you to move forward.

[wrridgeway]

I'm pretty sure this walk can be reverted, I just needed some flexibility to work with _test data here.

[jeancochrane]

Thank you for your persistence getting this over the line! Basically no major concerns from me, other than one idea about where the corner lots inconsistency might be coming from. Feel free to punt it to a follow-up if you'd rather get the majority of these changes in before moving on to small fixes.

[jeancochrane]

[Thought, non-blocking] You probably noticed this already, but it seems like the package is now developed under the name tabulapdf with a new maintainer: https://docs.ropensci.org/tabulapdf/news/index.html#changes-to-tabulapdf-105

[jeancochrane]

[Question, non-blocking] For the sake of my own clarity: It seems a little strange to drop the 3435 geometry column and then immediately transform the main geometry column (presumably 4326) to 3435, but we're doing this because the st_() functions can only operate on the main geometry column, is that right? Is there any consequence to re-transforming the 4326 column instead of using the 3435 column directly? I notice this is a pattern that is now widespread across the scripts so I want to make sure I understand it.

[jeancochrane]

[Question, non-blocking] How did you determine that the IDOT shapefiles use CRS 3436?

[jeancochrane]

[Question, non-blocking] To clarify my understanding: This transformation step is actually unnecessary because the Cook boundary geometry is already 4326, right? I don't necessarily mind the inconsistency, I just want to make sure I understand why we were able to remove the equivalent transform step in spatial-access.R and export-geojson.R.

[jeancochrane]

[Thought, non-blocking] Good idea! I wonder if this street network might be useful more broadly, such that we might want to save it under a more general location than ccao/corner/? You would have a better sense of its usefulness than I do, so I leave it up to you to decide whether it's worth the trouble.

[jeancochrane]

[Nitpick, optional] Am I right that the raw script is now handling the step of removing way types that are not main roads?

Suggested change

	# Fetch the OSM street network for the county, removing any OSM way types
	# that are not main roads
	# Fetch the OSM street network for the county

[jeancochrane]

[Question, blocking] I wonder if this might be the source of the right-angle issue you mentioned in the PR body? I'm still a geospatial novice, but my understanding is that 4326 operates in degree space, meaning the rectangle will be stretched into a parallelogram on the ground. It seems like the previous code transformed to 4326 after calculating the rectangle, which makes more sense to me (though I could be wrong):

data-architecture/etl/scripts-ccao-data-warehouse-us-east-1/spatial/spatial-ccao-corner.R

Lines 94 to 98 in 75c08a1

	town_mrr <- town_parcels %>%
	st_minimum_rotated_rectangle() %>%
	select(lon, lat, geometry = geometry_3435) %>%
	mutate(id = row_number()) %>%
	st_transform(4326) %>%

[jeancochrane]

[Suggestion, non-blocking] Hmm, this makes me a bit nervous: The absolute path to libjvm.so feels like it could easily change in the future or be different between machines, and the need to set flags to disable a feature of Java 16+ module system feels brittle. I wonder if it's worth trying to switch to the newer tabulapdf package before we persist these changes, to see if it resolves these problems? Part of my hesitation stems from a lack of full understanding of the problems these two settings are intending to address, so if you feel confident in your understanding and think it's the right call, I trust you to move forward.

[wrridgeway]

We decided to take a different approach in #1116. Other changes made in this PR were included in #1123 and #1125.