From ccccca3180d4e5354556a06e79246bdd2d872dc5 Mon Sep 17 00:00:00 2001 From: almac2022 Date: Sun, 30 Aug 2026 14:59:46 -0700 Subject: [PATCH 1/7] Initialize PWF baseline for #38 Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01GBKqedyBysV7hB4DuL98ZR --- planning/active/findings.md | 152 ++++++++++++++++++++++++++++++ planning/active/progress.md | 17 ++++ planning/active/task_plan.md | 178 +++++++++++++++++++++++++++++++++++ 3 files changed, 347 insertions(+) create mode 100644 planning/active/findings.md create mode 100644 planning/active/progress.md create mode 100644 planning/active/task_plan.md diff --git a/planning/active/findings.md b/planning/active/findings.md new file mode 100644 index 0000000..41d4253 --- /dev/null +++ b/planning/active/findings.md @@ -0,0 +1,152 @@ +# Findings — Georeference digital frames: the corner mapping assumes a square, axis-aligned footprint (#38) + +## Issue context + +## Problem + +`fly_georef()` excludes non-square footprints, so the digital frames #32 just gave +footprints to cannot be georeferenced. + +`georef_one()` maps image corners onto footprint corners **positionally** +(`R/fly_georef.R`, `# fly_footprint builds: BL, BR, TR, TL, BL`) and then shifts that +mapping by `bearing_to_rotation()`, a 90-degree quantization of the flight bearing. That +scheme was calibrated against north-up 9x9 negatives, where the footprint is square and +axis-aligned. + +A digital footprint breaks both assumptions: + +- It is **already rotated** onto the flight line by `fly_footprint()`, so its ring + carries the bearing. Applying `bearing_to_rotation()` on top counts it twice. +- It is **not square** — 1.10:1 to 1.80:1. On a square, a wrong-by-90 corner mapping is + harmless. On a 1.76:1 rectangle it maps a landscape image onto a portrait quad and the + warp squashes it. + +## Why it was deferred rather than guessed + +The right corner mapping for a pre-rotated rectangle depends on camera mounting relative +to flight direction. That is precisely what `bearing_to_rotation()` was empirically +fitted to, and it cannot be re-derived without imagery to check the result against — +nothing in the test suite looks at pixels, so a wrong warp would pass every assertion. + +Digital frames had no footprint at all before #32, so excluding them is the same +coverage as before rather than a regression. It is now explicit and warned about instead +of silent. + +## Proposed approach + +1. Fetch a handful of real digital thumbnails via `fly_fetch()` for frames whose + footprint is known — the bundled `inst/testdata/photo_centroids_digital.gpkg` has + `thumbnail_image_url` populated for two cameras. +2. Establish the correct pixel-corner to footprint-corner mapping for a footprint + already rotated onto its flight line. Likely a constant, since the rotation is + already in the ring — but that must be measured, not assumed. +3. Check the result **visually**, and pin whatever invariant the check establishes + (a landscape image must not land in a portrait quad; a known ground feature must + land where it belongs). +4. Remove the exclusion in `fly_georef()` and the `fly_is_square()` guard it uses. + +## Value + +**If we do it:** the whole post-2010 catalogue becomes georeferenceable, not just +sizeable. + +**If we never do:** digital frames get footprints, coverage and selection but no +georeferenced output — visibly excluded with a warning rather than silently wrong, which +is the important half. + +Follows #32 + + +--- + +## Measurements taken 2026-08-30 during plan-mode exploration + +### The ring order is in the rectangle's own frame + +`fly_rectangles()` (`R/fly_footprint.R:124`) emits, in local coordinates: + +| vertex | local | flight-relative | +|---|---|---| +| `coords[1]` | `(-hc, -ha)` | rear-left | +| `coords[2]` | `( hc, -ha)` | rear-right | +| `coords[3]` | `( hc, ha)` | front-right | +| `coords[4]` | `(-hc, ha)` | front-left | + +Local `+y` is the heading, local `+x` is 90 degrees clockwise of it. Confirmed by +expanding the rotation: `rot <- matrix(c(cos, sin, -sin, cos), nrow = 2)` fills +column-wise, so `xy %*% rot` sends row vector `(0, 1)` to `(sin b, cos b)` — the heading +itself, matching the comment at `R/fly_footprint.R:140-142`. + +Rotation is applied only when `hc != ha` **and** the bearing is finite. So the ring means +the same thing whether or not it was rotated, which is why one constant can serve both +the bearing-rotated and the `axis_aligned_no_bearing` non-square cases. + +### Rotation 90/270 are the aspect-consistent candidates, not 0/180 + +`georef_one()` builds `fp_corners = [coords4, coords3, coords2, coords1]` += `[front-left, front-right, rear-right, rear-left]` and maps pixel `[TL, TR, BR, BL]` +onto it, cyclically shifted by `rotation %/% 90`. + +| rotation | image WIDTH lands on | +|---|---| +| 0 | cross-track edge | +| 90 | along-track edge | +| 180 | cross-track edge | +| 270 | along-track edge | + +### Both bundled digital cameras deliver PORTRAIT thumbnails + +Fetched 2026-08-30 from `openmaps.gov.bc.ca`: + +| camera | calib key | thumbnail px | sensor px (cross x along) | footprint m | ratio | +|---|---|---|---|---|---| +| Leica DMC II | `121201_2011` | 884 x 972 | 15552 x 14144 | 4666 x 4243 | 1.100 | +| UltraCam Eagle M3 | `20814295_2018` | 1063 x 1654 | 26460 x 17004 | 3175 x 2040 | 1.556 | + +The DMC II thumbnail's EXIF carries the source TIFF dimensions — `width=14144, +height=15552` — so the full-resolution frame is portrait too and the thumbnail is not +rotated relative to it. + +In both, image **height** is the long axis and equals `px_cross`, the axis +`data-raw/make_camera_formats.R` assumes is across-track. So the image's long axis must +land on the footprint's long edge: **{90, 270} accepted, {0, 180} rejected, by geometry +alone.** + +### One bit remains and needs pixels + +90 and 270 differ by 180 degrees about the footprint centre. A rectangle is symmetric +under that, so no geometric invariant distinguishes them. + +### A second assumption rides along + +`data-raw/make_camera_formats.R:211,241,261` sets `px_cross = max(dims)`, +`px_along = min(dims)` — it assumes the long sensor axis is across-track. +`tests/testthat/test-fly_camera_format.R:119` asserts the footprint *follows* that +assumption, not that the assumption is true. If it is wrong the ground quad itself is +rotated 90 degrees and no corner mapping repairs it — that would reopen #32. + +### The bundled fixture reaches the failure mode, unevenly + +A wrong-by-90 mapping distorts by `ratio^2`: **2.42x** on the UltraCam (frames 19-24), +only **1.21x** on the DMC II (frames 1-18). The UltraCam frames are what make the aspect +invariant discriminating. + +### Bundled digital bearings + +Frames 1-18 (`bcd13304/5/6`, 2013): ~270.2-271.3 degrees. +Frames 19-24 (`bcd19503`, 2019): ~342.6-343.1 degrees. + +## Ground truth for the imagery half + +BC orthophoto imagery, catalogued in a private sibling repo. Joins to the airphoto +catalogue on `BCGS_TILE` + `PHOTO_YEAR`. The provincial collection is "Access Only" and +sold to the public, so ortho pixels must **not** be bundled into fly as test data +regardless of reachability. The QA runs there and produces a constant; the constant comes +back to fly, the imagery never does. + +No file in this repo names that repo, its endpoint or its database. + +## Errors Encountered + +| Error | Resolution | +|-------|------------| diff --git a/planning/active/progress.md b/planning/active/progress.md new file mode 100644 index 0000000..bc03149 --- /dev/null +++ b/planning/active/progress.md @@ -0,0 +1,17 @@ +# Progress — Georeference digital frames (#38) + +## Session 2026-08-30 + +- Plan-mode exploration: read `R/fly_georef.R`, `R/fly_footprint.R`, + `tests/testthat/test-fly_georef.R`, `tests/testthat/test-fly_camera_format.R`, + `data-raw/make_camera_formats.R`, `data-raw/test_georef_decades.R` +- Measured the ring-order contract, the rotation-to-edge mapping, and both bundled + cameras' thumbnail dimensions — see `findings.md`. Narrowed the corner mapping from + eight possibilities to one bit by geometry alone +- Two decisions taken with the user: the orthophoto QA runs in the private catalogue repo + and only its constant comes back; a user `rotation` column keeps overriding for + non-square frames +- Phases approved by user +- Created branch `38-georeference-digital-frames-the-corner-m` off main +- Scaffolded PWF baseline +- Next: Phase 1 — extract the GCP construction into a pure, testable function diff --git a/planning/active/task_plan.md b/planning/active/task_plan.md new file mode 100644 index 0000000..1964f55 --- /dev/null +++ b/planning/active/task_plan.md @@ -0,0 +1,178 @@ +# Task: Georeference digital frames — the corner mapping assumes a square, axis-aligned footprint (#38) + +## Context + +fly#32 gave digital frames footprints. `fly_georef()` cannot use them: it excludes every +non-square footprint with a warning (`R/fly_georef.R:135-143`), so the whole post-2010 +catalogue is sizeable but not georeferenceable. + +`georef_one()` (`R/fly_georef.R:223`) maps image corners onto footprint corners +positionally, then cyclically shifts that mapping by `bearing_to_rotation()` — a +90-degree quantization of the flight bearing, fitted against north-up 9x9 film negatives. +A digital footprint breaks it twice: the ring is **already** rotated onto the flight line +by `fly_rectangles()`, so the bearing would be counted a second time; and on a 1.56:1 +rectangle a wrong-by-90 mapping puts a portrait image in a landscape quad. + +The issue defers rather than guesses because "nothing in the test suite looks at pixels, +so a wrong warp would pass every assertion." That is true of the *last* bit only. Most of +the answer is derivable from geometry the repo already has, and this plan separates the +two so the imagery step has one question to answer instead of eight. + +### What exploration established (measured this session) + +**The ring order is in the rectangle's own frame.** `fly_rectangles()` +(`R/fly_footprint.R:124`) emits `v1=(-hc,-ha) rear-left, v2=(hc,-ha) rear-right, +v3=(hc,ha) front-right, v4=(-hc,ha) front-left`, where local `+y` is the heading and +local `+x` is 90 degrees clockwise of it. Bearing rotation is applied only when +`hc != ha` and the bearing is finite. So the ring means the same thing whether or not it +was rotated — which is why one constant can serve both the rotated and the +`axis_aligned_no_bearing` non-square cases. + +**Rotation 90/270 are the aspect-consistent candidates, not 0/180.** `georef_one()` builds +`fp_corners = [front-left, front-right, rear-right, rear-left]` and maps pixel +`[TL, TR, BR, BL]` onto it. Rotations 0 and 180 send image *width* to the **cross**-track +edge; 90 and 270 send image *width* to the **along**-track edge. Both bundled digital +cameras deliver **portrait** thumbnails whose long axis is the image *height*: + +| camera | thumbnail | sensor px (cross x along) | footprint (m) | ratio | +|---|---|---|---|---| +| Leica DMC II (`121201_2011`) | 884 x 972 | 15552 x 14144 | 4666 x 4243 | 1.100 | +| UltraCam Eagle M3 (`20814295_2018`) | 1063 x 1654 | 26460 x 17004 | 3175 x 2040 | 1.556 | + +Image height is the long axis and equals `px_cross`, the axis the format table assumes is +across-track. So the image's long axis must land on the footprint's long edge, which +selects {90, 270} and rejects {0, 180}. **That is settled by geometry — no imagery.** + +**One bit remains, and it needs pixels.** 90 vs 270 differ by 180 degrees about the +footprint centre; a rectangle is symmetric under that, so no geometric invariant can tell +them apart. + +**A second, larger assumption rides along.** `data-raw/make_camera_formats.R:211,241,261` +sets `px_cross = max(dims)`, `px_along = min(dims)` — it *assumes* the long sensor axis is +across-track. `test-fly_camera_format.R:119` asserts the footprint follows that assumption, +not that the assumption is true. If it is wrong, the ground quad itself is rotated 90 +degrees and **no corner mapping repairs it** — #32's footprints would be wrong. The same +orthophoto QA settles it, so this plan checks it explicitly rather than inheriting it. + +**The bundled fixture can reach the failure mode, unevenly.** A wrong-by-90 mapping +distorts by `ratio^2`: 2.42x on the UltraCam, only 1.21x on the DMC II. Frames 19-24 are +the diagnostic ones; that gets stated in the test rather than left to luck. + +### Decisions taken + +- **The orthophoto QA lives in the private catalogue repo.** Its pixels are "Access Only" + and licence-restricted, so they never enter fly as fixtures. fly receives the measured + constant, a note, and invariants that need no reference imagery. No fly file names that + repo, its endpoint or its database. +- **A user `rotation` column keeps overriding, for square and non-square alike.** + Consistent with today's documented behaviour. The carried-column hazard — a film-era + rotation column applied to a digital batch — gets documented, not guarded. + +--- + +## Phase 1: Make the corner mapping testable without GDAL + +- [ ] Extract GCP construction from `georef_one()` into a pure internal + `fly_georef_gcps(ncol_px, nrow_px, coords, rotation)` returning the pixel->ground + correspondence. No GDAL, no file I/O, no `sf::gdal_utils()`. +- [ ] Rewrite `georef_one()` to call it; everything else in that function unchanged. +- [ ] `tests/testthat/test-fly_georef_gcps.R`: film GCP output identical at all four + rotations to the pre-change implementation, pulled with + `git show HEAD:R/fly_georef.R` (not reconstructed from memory). +- [ ] Put the helper at the top of `R/fly_georef.R` or its own file — never between a + roxygen block and the function it documents (the fly#30 `@export` rebind). + +**Verify:** `devtools::test()` green; the film-parity test asserts equality against bytes +pulled from git, and fails if the extraction changed any coordinate. + +## Phase 2: Pin the aspect invariant (the half needing no imagery) + +- [ ] Assert: the GCP mapping sends the image's long pixel axis onto the footprint's long + ground edge. Concretely, ground distance between the width-pair GCPs divided by + ground distance between the height-pair GCPs equals `ncol_px / nrow_px`. +- [ ] Run it over the bundled digital frames, **both** cameras, from + `inst/testdata/photo_centroids_digital.gpkg` via `digital_fixture()` / + `tests/testthat/setup.R`. +- [ ] Assert the invariant **fails** at the rotations it must reject. A test that only + ever passes is decoration. +- [ ] State in a comment that the invariant is vacuous on square film (ratio 1) and weak + on the DMC II (1.21x) — the UltraCam frames are what makes it discriminating. + +**Verify:** the invariant passes for {90, 270} and fails for {0, 180}, both asserted; and +the DMC II / UltraCam distortion factors are recorded so a future fixture change that +drops the UltraCam is visible. + +## Phase 3: Establish the constant against orthophoto ground truth + +Runs in the private catalogue repo. Two facts settle together: + +- [ ] (a) Is the long sensor axis actually across-track? Compare the #32 footprint's + orientation against the ortho for the same ground. A 90-degree error here is not a + georef bug — it is a footprint bug, and it would reopen #32. +- [ ] (b) Which of 90 / 270 is correct? They differ by a 180-degree flip of image content + inside the same quad. +- [ ] Sample both bundled cameras plus the 1.76:1 Leica DMC III — the highest-ratio and + therefore most diagnostic case. Join ortho to catalogue on `BCGS_TILE` + `PHOTO_YEAR`. +- [ ] Check visually **and** with a measurable: a named ground feature must land in the + same place in the warped thumbnail as in the ortho. A screenshot alone is not a result. +- [ ] Record the answer **per camera**. Do not assume one constant serves both until + both are measured. If they differ, the constant becomes a `camera_formats.csv` column + and Phase 4 grows a lookup — decide then, not now. +- [ ] Bring back only: the constant(s), the sample size, and what was compared. No pixels. + +**Verify:** the same answer for every frame of a camera; disagreement within a camera means +the model is wrong, not that the sample is noisy. + +## Phase 4: Land the constant and remove the exclusion + +- [ ] `fly_georef()` applies the constant to non-square footprints and does **not** apply + `bearing_to_rotation()` to them — the ring already carries the bearing. +- [ ] Branch on `fly_is_square()` (`R/fly_footprint.R:165`), computed for every row before + the per-row loop. Never on `half_cross`/`half_along`, which are NA by construction + until a sizing route fills them (CLAUDE.md, fly#32's three-round trap). +- [ ] Keep the `rotation`-column override for non-square frames; document the + carried-column hazard in `@param rotation`. +- [ ] Delete the exclusion block (`R/fly_georef.R:122-143`) and the `rotated[fp_idx[1]]` + skip (`:190`). Keep `fly_is_square()` itself — `test-fly_camera_format.R:266` uses it. +- [ ] Keep the empty-geometry skip (`:195`) and `fly_warn_unsized()` (`:120`) untouched; + they cover a different exclusion. +- [ ] Rewrite the **Rotation** `@details` section (`:38-65`), which currently documents + only the film scheme, and the `@param rotation` text that says `"auto"` applies + everywhere. + +**Verify:** digital frames produce GeoTIFFs; a mixed film+digital batch gives film the +bearing rotation and digital the constant, asserted per row; no non-square warning fires; +`test-fly_georef.R`'s existing eight tests still pass unchanged. + +## Phase 5: Notes, docs, release + +- [ ] `inst/notes/georeferencing.md`, companion to `terrain-correction.md` and + `camera-formats.md`: the ring-order contract, why the bearing must not be applied + twice, what the QA measured and against what, and what the aspect invariant can and + cannot catch. Ground truth referenced obliquely — no repo, endpoint or database named. +- [ ] Update CLAUDE.md: the Key Decisions entry for #30 and the NEWS line + "`fly_georef()` excludes rotated footprints" are both now stale. +- [ ] `devtools::document()` — **read its output**. `Writing '.Rd'` or a + falling `grep -c "^export(" NAMESPACE` is the fly#30 rebind. +- [ ] `lintr::lint_package()` against the `HEAD` baseline, not against zero. +- [ ] NEWS entry; version 0.6.0 -> 0.7.0 as the **final** commit of the branch. + +## Validation + +- [ ] `devtools::test()` passes +- [ ] `/code-check` clean on each commit +- [ ] No fly file names the private ortho repo, its endpoint, or its database: + `git diff main --stat` reviewed, plus a grep of the branch diff for the repo name +- [ ] Restore-the-bug check on the Phase 2 invariant: it goes red against the rejected + rotations, patched in **both** `asNamespace("fly")` and + `as.environment("package:fly")`, with a printed value proving the patch took +- [ ] PWF checkboxes match landed work; `/planning-archive` on completion + +## Open, deliberately not decided yet + +- **One constant or one per camera** — Phase 3 measures it. If per-camera, it belongs in + `camera_formats.csv` beside the sensor dimensions, not hardcoded in `fly_georef()`. +- **What happens if Phase 3(a) refutes the across-track assumption.** Then #32's + footprints are rotated 90 degrees wrong and that is a bigger issue than this one; stop + and reopen #32 rather than compensating for it in the corner mapping. +- **fly#37** (`fly_footprint()` handed its own output) is untouched by this work. From de8068df70ca17e975b33b6abf61ca51928310fd Mon Sep 17 00:00:00 2001 From: almac2022 Date: Sun, 30 Aug 2026 15:05:52 -0700 Subject: [PATCH 2/7] Split the GCP corner mapping out of georef_one() (#38) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `fly_georef_gcps()` is the part of georeferencing that can be wrong while everything around it looks healthy: a wrong-by-90 corner mapping still writes a valid GeoTIFF, in the right CRS, with the right extent, and nothing in this package looks at pixels. Pure, so it is checkable offline. Verified against the implementation lifted from git rather than rewritten from memory: 1024 cases (20 bundled film footprints, 24 digital, 20 random rings, at four image shapes and four rotations), 0 differences. Golden values are written into the test so the pin survives that sha becoming history. One real difference surfaced doing it — the old loop produced a named character vector whenever the ring came from `sf::st_coordinates()`, since `coords[4, 1:2]` carries X/Y dimnames. GDAL ignores names; `identical()` does not. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01GBKqedyBysV7hB4DuL98ZR --- R/fly_georef.R | 92 +++++++++++++-------- planning/active/progress.md | 15 ++++ planning/active/task_plan.md | 8 +- tests/testthat/test-fly_georef_gcps.R | 114 ++++++++++++++++++++++++++ 4 files changed, 189 insertions(+), 40 deletions(-) create mode 100644 tests/testthat/test-fly_georef_gcps.R diff --git a/R/fly_georef.R b/R/fly_georef.R index 843217f..c21f012 100644 --- a/R/fly_georef.R +++ b/R/fly_georef.R @@ -238,43 +238,15 @@ georef_one <- function(src, fp, out_file, srcnodata = "0", rotation = 180) { n_bands <- length(gregexpr("Band \\d+", info)[[1]]) is_rgb <- n_bands >= 3 - # Pixel corners: TL, TR, BR, BL - pixel_corners <- list( - c(0, 0), # TL - c(ncol_px, 0), # TR - c(ncol_px, nrow_px), # BR - c(0, nrow_px) # BL - ) - - # Footprint corners in same order: TL, TR, BR, BL - fp_corners <- list( - coords[4, 1:2], # TL - coords[3, 1:2], # TR - coords[2, 1:2], # BR - coords[1, 1:2] # BL - ) - - # Rotation: shift the footprint corner mapping - - # rotation=0: pixel TL → footprint TL (north-up, original behavior) - # rotation=90: pixel TL → footprint TR (top of image = east) - # rotation=180: pixel TL → footprint BR (top of image = south) - # rotation=270: pixel TL → footprint BL (top of image = west) - n_shifts <- rotation %/% 90 - if (n_shifts > 0) { - fp_corners <- c( - fp_corners[(n_shifts + 1):4], - fp_corners[1:n_shifts] - ) - } - - # Build GCP args mapping pixel corners to (rotated) footprint corners + # Build GCP args from the pixel-to-ground correspondence. + gcp <- fly_georef_gcps(ncol_px, nrow_px, coords, rotation) + # `unname()` matters: the ring arrives from `sf::st_coordinates()` carrying X/Y + # dimnames, and without it the option vector handed to GDAL is named. Harmless to + # GDAL, but it makes the args unequal to a plain character vector under + # `identical()`, which is what the parity test compares. gcp_args <- character(0) - for (j in seq_along(pixel_corners)) { - gcp_args <- c(gcp_args, - "-gcp", pixel_corners[[j]][1], pixel_corners[[j]][2], - fp_corners[[j]][1], fp_corners[[j]][2] - ) + for (j in seq_len(nrow(gcp))) { + gcp_args <- c(gcp_args, "-gcp", as.character(unname(gcp[j, ]))) } # Step 1: translate with GCPs @@ -327,3 +299,51 @@ bearing_to_rotation <- function(bearing) { rot[is.na(rot)] <- 180L as.integer(rot) } + +#' Map image pixel corners onto footprint ground corners +#' +#' The pure half of [georef_one()] — no GDAL, no file I/O, no `sf`. Split out +#' because it is the part that can be wrong, and a GCP correspondence is +#' checkable offline in a way a warped GeoTIFF is not. +#' +#' `coords` is the footprint's ring as [fly_footprint()] builds it, whose vertex +#' order is a contract: rows 1-4 are BL, BR, TR, TL **in the rectangle's own +#' frame**. For a non-square footprint that frame is the flight line's, so those +#' are rear-left, rear-right, front-right, front-left — the same meaning whether +#' or not the ring was rotated onto a bearing. See `fly_rectangles()`. +#' +#' @param ncol_px,nrow_px Image dimensions in pixels. +#' @param coords A 4-row matrix of footprint ring coordinates, x in column 1 +#' and y in column 2. +#' @param rotation Image rotation, one of 0, 90, 180, 270. Cyclically shifts +#' which ground corner the top-left pixel maps to. +#' @return A 4-row numeric matrix with columns `pixel_x`, `pixel_y`, +#' `ground_x`, `ground_y`, one row per corner in the order TL, TR, BR, BL. +#' @noRd +fly_georef_gcps <- function(ncol_px, nrow_px, coords, rotation) { + # Pixel corners: TL, TR, BR, BL + pixel <- matrix( + c(0, 0, + ncol_px, 0, + ncol_px, nrow_px, + 0, nrow_px), + ncol = 2, byrow = TRUE + ) + + # Footprint corners in the same order: TL, TR, BR, BL + ground <- coords[c(4, 3, 2, 1), 1:2, drop = FALSE] + + # Rotation shifts the footprint corner mapping: + # 0 pixel TL → footprint TL (north-up on a square; front-left on a rectangle) + # 90 pixel TL → footprint TR + # 180 pixel TL → footprint BR + # 270 pixel TL → footprint BL + n_shifts <- rotation %/% 90 + if (n_shifts > 0) { + ground <- ground[c((n_shifts + 1):4, 1:n_shifts), , drop = FALSE] + } + + out <- cbind(pixel, ground) + dimnames(out) <- list(NULL, c("pixel_x", "pixel_y", "ground_x", "ground_y")) + out +} diff --git a/planning/active/progress.md b/planning/active/progress.md index bc03149..ddb8ae4 100644 --- a/planning/active/progress.md +++ b/planning/active/progress.md @@ -15,3 +15,18 @@ - Created branch `38-georeference-digital-frames-the-corner-m` off main - Scaffolded PWF baseline - Next: Phase 1 — extract the GCP construction into a pure, testable function + +### Phase 1 — GCP construction extracted + +- `fly_georef_gcps(ncol_px, nrow_px, coords, rotation)` split out of `georef_one()`: + pure, no GDAL, no file I/O +- Parity verified against the pre-change implementation lifted with + `git show HEAD:R/fly_georef.R` (not reconstructed): 1024 cases — the 20 bundled film + footprints, the 24 digital ones, and 20 random rings, at four image shapes and four + rotations — **0 differences** +- One real difference found and fixed while doing it: the old loop produced a *named* + character vector whenever the ring came from `sf::st_coordinates()`, because + `coords[4, 1:2]` carries X/Y dimnames. GDAL ignores names; `identical()` does not +- Golden values written into `tests/testthat/test-fly_georef_gcps.R` so the pin survives + that sha becoming history +- Full suite: 1181 passing, 0 failures diff --git a/planning/active/task_plan.md b/planning/active/task_plan.md index 1964f55..7e68171 100644 --- a/planning/active/task_plan.md +++ b/planning/active/task_plan.md @@ -72,14 +72,14 @@ the diagnostic ones; that gets stated in the test rather than left to luck. ## Phase 1: Make the corner mapping testable without GDAL -- [ ] Extract GCP construction from `georef_one()` into a pure internal +- [x] Extract GCP construction from `georef_one()` into a pure internal `fly_georef_gcps(ncol_px, nrow_px, coords, rotation)` returning the pixel->ground correspondence. No GDAL, no file I/O, no `sf::gdal_utils()`. -- [ ] Rewrite `georef_one()` to call it; everything else in that function unchanged. -- [ ] `tests/testthat/test-fly_georef_gcps.R`: film GCP output identical at all four +- [x] Rewrite `georef_one()` to call it; everything else in that function unchanged. +- [x] `tests/testthat/test-fly_georef_gcps.R`: film GCP output identical at all four rotations to the pre-change implementation, pulled with `git show HEAD:R/fly_georef.R` (not reconstructed from memory). -- [ ] Put the helper at the top of `R/fly_georef.R` or its own file — never between a +- [x] Put the helper at the top of `R/fly_georef.R` or its own file — never between a roxygen block and the function it documents (the fly#30 `@export` rebind). **Verify:** `devtools::test()` green; the film-parity test asserts equality against bytes diff --git a/tests/testthat/test-fly_georef_gcps.R b/tests/testthat/test-fly_georef_gcps.R new file mode 100644 index 0000000..372aeca --- /dev/null +++ b/tests/testthat/test-fly_georef_gcps.R @@ -0,0 +1,114 @@ +# The pixel-to-ground correspondence is the part of georeferencing that can be wrong +# while everything around it looks healthy: a wrong-by-90 corner mapping still writes a +# valid GeoTIFF, in the right CRS, with the right extent, and nothing in this package +# looks at pixels. So it is tested here as a pure function, offline. + +# A square, axis-aligned footprint — the film case. Ring order is the contract +# `fly_footprint()` documents: BL, BR, TR, TL. +film_ring <- function() { + matrix(c( + 1000, 2000, # BL + 5000, 2000, # BR + 5000, 6000, # TR + 1000, 6000 # TL + ), ncol = 2, byrow = TRUE) +} + +# UltraCam Eagle M3-shaped: 3175 x 2040 m about the origin, long axis across-track. +ultracam_ring <- function() { + matrix(c( + -1587.5, -1020, + 1587.5, -1020, + 1587.5, 1020, + -1587.5, 1020 + ), ncol = 2, byrow = TRUE) +} + +ground_of <- function(g) unname(g[, c("ground_x", "ground_y"), drop = FALSE]) + + +test_that("fly_georef_gcps returns four corners in TL, TR, BR, BL pixel order", { + g <- fly_georef_gcps(1250L, 1250L, film_ring(), 0) + + expect_equal(dim(g), c(4L, 4L)) + expect_identical(colnames(g), c("pixel_x", "pixel_y", "ground_x", "ground_y")) + expect_equal( + unname(g[, c("pixel_x", "pixel_y")]), + matrix(c(0, 0, 1250, 0, 1250, 1250, 0, 1250), ncol = 2, byrow = TRUE) + ) + + # Pixel corners do not move with rotation — only the ground corners they map to. + for (rot in c(0, 90, 180, 270)) { + expect_equal( + unname(fly_georef_gcps(1250L, 1250L, film_ring(), rot)[, c("pixel_x", "pixel_y")]), + unname(g[, c("pixel_x", "pixel_y")]), + info = paste("rotation =", rot) + ) + } +}) + + +test_that("film corner mapping is unchanged at every rotation", { + # Golden values, not a re-derivation. Captured from the implementation as it stood + # before the GCP construction was split out of `georef_one()` — the bytes pulled with + # `git show :R/fly_georef.R`, exercised over the bundled film and digital + # footprints plus random rings at four image shapes: 1024 cases, 0 differences. + # + # Written out rather than computed so this keeps working once that sha is history. + expected <- list( + "0" = matrix(c(1000, 6000, 5000, 6000, 5000, 2000, 1000, 2000), ncol = 2, byrow = TRUE), + "90" = matrix(c(5000, 6000, 5000, 2000, 1000, 2000, 1000, 6000), ncol = 2, byrow = TRUE), + "180" = matrix(c(5000, 2000, 1000, 2000, 1000, 6000, 5000, 6000), ncol = 2, byrow = TRUE), + "270" = matrix(c(1000, 2000, 1000, 6000, 5000, 6000, 5000, 2000), ncol = 2, byrow = TRUE) + ) + + for (rot in names(expected)) { + expect_equal( + ground_of(fly_georef_gcps(1250L, 1250L, film_ring(), as.numeric(rot))), + expected[[rot]], + info = paste("rotation =", rot) + ) + } +}) + + +test_that("non-square corner mapping is unchanged at every rotation", { + # Same golden capture, on a rectangle. Included because the film case is square and + # therefore cannot distinguish a shift that swaps the two axes from one that does not. + expected <- list( + "0" = matrix(c(-1587.5, 1020, 1587.5, 1020, 1587.5, -1020, -1587.5, -1020), + ncol = 2, byrow = TRUE), + "90" = matrix(c( 1587.5, 1020, 1587.5, -1020, -1587.5, -1020, -1587.5, 1020), + ncol = 2, byrow = TRUE), + "180" = matrix(c( 1587.5, -1020, -1587.5, -1020, -1587.5, 1020, 1587.5, 1020), + ncol = 2, byrow = TRUE), + "270" = matrix(c(-1587.5, -1020, -1587.5, 1020, 1587.5, 1020, 1587.5, -1020), + ncol = 2, byrow = TRUE) + ) + + for (rot in names(expected)) { + expect_equal( + ground_of(fly_georef_gcps(1063L, 1654L, ultracam_ring(), as.numeric(rot))), + expected[[rot]], + info = paste("rotation =", rot) + ) + } +}) + + +test_that("rotation 0 and 180 put image width on the cross-track edge, 90 and 270 along", { + # The property the golden tables above encode, stated once so a future change to them + # has to disagree with a sentence rather than only with a number. + # + # `fly_rectangles()` builds the ring in the rectangle's own frame, so edge 1->2 is + # across-track and edge 2->3 along-track. Here: 3175 m across, 2040 m along. + width_edge_m <- function(rot) { + g <- ground_of(fly_georef_gcps(1063L, 1654L, ultracam_ring(), rot)) + sqrt(sum((g[2, ] - g[1, ])^2)) # TL -> TR is the image's width + } + + expect_equal(width_edge_m(0), 3175) + expect_equal(width_edge_m(180), 3175) + expect_equal(width_edge_m(90), 2040) + expect_equal(width_edge_m(270), 2040) +}) From 128f3d8ba36e70bb17fa1d1d7b6cbc4ebbda5674 Mon Sep 17 00:00:00 2001 From: almac2022 Date: Sun, 30 Aug 2026 15:06:51 -0700 Subject: [PATCH 3/7] Pin the aspect invariant for digital corner mapping (#38) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Half of #38's question is settled by geometry and needs no reference imagery: the image's long pixel axis must land on the footprint's long ground edge. That rejects rotations 0 and 180 for every digital frame in the fixture and leaves 90 and 270, which differ by 180 degrees about the footprint centre and cannot be told apart geometrically. The negative is pinned with measured numbers rather than a threshold — a wrong-by-90 mapping squashes by the footprint's aspect squared, 1.21x on the DMC II and 2.42x on the UltraCam Eagle M3. The DMC II alone would survive a loose tolerance, so the test says which frames make it discriminating. Restored the pre-#38 defect (a mapping that ignores its rotation argument, which is what `bearing_to_rotation()` produces for these frames) in both the namespace and the attached package environment, with a printed ground coordinate proving the patch took: the invariant goes red. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01GBKqedyBysV7hB4DuL98ZR --- planning/active/progress.md | 19 ++++ planning/active/task_plan.md | 10 +- tests/testthat/test-fly_georef_aspect.R | 122 ++++++++++++++++++++++++ 3 files changed, 146 insertions(+), 5 deletions(-) create mode 100644 tests/testthat/test-fly_georef_aspect.R diff --git a/planning/active/progress.md b/planning/active/progress.md index ddb8ae4..cd63f83 100644 --- a/planning/active/progress.md +++ b/planning/active/progress.md @@ -30,3 +30,22 @@ - Golden values written into `tests/testthat/test-fly_georef_gcps.R` so the pin survives that sha becoming history - Full suite: 1181 passing, 0 failures + +### Phase 2 — aspect invariant pinned + +- `tests/testthat/test-fly_georef_aspect.R`: the GCP mapping must send the image's long + pixel axis onto the footprint's long ground edge, expressed as an anisotropy ratio + (m/px on the width axis over m/px on the height axis) that must be 1 +- Premise asserted beside it: the shipped `camera_formats.csv` aspect matches the + aspect of the thumbnails the catalogue actually serves, to 1e-3. A regenerated CSV + that disagreed would fail there, naming the cause, rather than here +- Negative half pinned with numbers, not a threshold: rotations 0 and 180 squash by + ratio^2 — **1.21x** on the DMC II, **2.42x** on the UltraCam. The DMC II alone would + survive a loose tolerance; a fixture change dropping the UltraCam frames fails +- Stated explicitly that the invariant is vacuous on square film, which is why the film + mapping needed imagery and why this cannot finish the digital job either +- **Restore-the-bug check**: patched `fly_georef_gcps()` to ignore its rotation argument + (the pre-#38 behaviour — `bearing_to_rotation(271)` returns 0), in *both* + `asNamespace("fly")` and `as.environment("package:fly")`, with a printed ground + coordinate that can only come from the broken version. Result: **FAIL=2**, the two + isotropy assertions. The guard fires diff --git a/planning/active/task_plan.md b/planning/active/task_plan.md index 7e68171..3c2dec9 100644 --- a/planning/active/task_plan.md +++ b/planning/active/task_plan.md @@ -87,15 +87,15 @@ pulled from git, and fails if the extraction changed any coordinate. ## Phase 2: Pin the aspect invariant (the half needing no imagery) -- [ ] Assert: the GCP mapping sends the image's long pixel axis onto the footprint's long +- [x] Assert: the GCP mapping sends the image's long pixel axis onto the footprint's long ground edge. Concretely, ground distance between the width-pair GCPs divided by ground distance between the height-pair GCPs equals `ncol_px / nrow_px`. -- [ ] Run it over the bundled digital frames, **both** cameras, from +- [x] Run it over the bundled digital frames, **both** cameras, from `inst/testdata/photo_centroids_digital.gpkg` via `digital_fixture()` / `tests/testthat/setup.R`. -- [ ] Assert the invariant **fails** at the rotations it must reject. A test that only +- [x] Assert the invariant **fails** at the rotations it must reject. A test that only ever passes is decoration. -- [ ] State in a comment that the invariant is vacuous on square film (ratio 1) and weak +- [x] State in a comment that the invariant is vacuous on square film (ratio 1) and weak on the DMC II (1.21x) — the UltraCam frames are what makes it discriminating. **Verify:** the invariant passes for {90, 270} and fails for {0, 180}, both asserted; and @@ -163,7 +163,7 @@ bearing rotation and digital the constant, asserted per row; no non-square warni - [ ] `/code-check` clean on each commit - [ ] No fly file names the private ortho repo, its endpoint, or its database: `git diff main --stat` reviewed, plus a grep of the branch diff for the repo name -- [ ] Restore-the-bug check on the Phase 2 invariant: it goes red against the rejected +- [x] Restore-the-bug check on the Phase 2 invariant: it goes red against the rejected rotations, patched in **both** `asNamespace("fly")` and `as.environment("package:fly")`, with a printed value proving the patch took - [ ] PWF checkboxes match landed work; `/planning-archive` on completion diff --git a/tests/testthat/test-fly_georef_aspect.R b/tests/testthat/test-fly_georef_aspect.R new file mode 100644 index 0000000..1648085 --- /dev/null +++ b/tests/testthat/test-fly_georef_aspect.R @@ -0,0 +1,122 @@ +# The aspect invariant: the half of the corner mapping that geometry settles. +# +# `georef_one()` hands GDAL four corner correspondences and lets it fit the transform. +# Nothing downstream checks that the fit is isotropic, so a mapping that sends the +# image's long axis onto the footprint's short edge produces a valid, wrongly-squashed +# GeoTIFF. That is the failure this file exists to make impossible. +# +# It cannot settle everything. Rotations 90 and 270 differ by 180 degrees about the +# footprint centre and a rectangle is symmetric under that, so no geometric assertion +# distinguishes them — see `inst/notes/georeferencing.md`. + +# Delivered image orientation, measured 2026-08-30 against the live thumbnails at +# openmaps.gov.bc.ca. Both digital cameras in the bundled fixture deliver PORTRAIT +# images whose long axis is the image height: +# +# Leica DMC II bcd13304_778_rgb_thumb.jpg 884 x 972 px +# UltraCam Eagle M3 bcd19503_348_thumb.jpg 1063 x 1654 px +# +# The DMC II thumbnail's EXIF carries its source TIFF as width=14144, height=15552, so +# the full-resolution frame is portrait too and the thumbnail is not rotated relative to +# it. Image width is therefore the along-track pixel count and image height the +# across-track one. +image_dims <- function(fmt) { + list(ncol_px = fmt$px_along, nrow_px = fmt$px_cross) +} + +# Metres per pixel along the image's width axis, divided by metres per pixel along its +# height axis. An isotropic mapping gives 1; anything else is a squash of that factor. +gcp_anisotropy <- function(ncol_px, nrow_px, ring, rot) { + g <- fly_georef_gcps(ncol_px, nrow_px, ring, rot)[, c("ground_x", "ground_y")] + w <- sqrt(sum((g[2, ] - g[1, ])^2)) / ncol_px + h <- sqrt(sum((g[3, ] - g[2, ])^2)) / nrow_px + w / h +} + +digital_rings <- function() { + photos <- sf::st_read(testdata_path("photo_centroids_digital.gpkg"), quiet = TRUE) + fp <- sf::st_transform(suppressWarnings(fly_footprint(photos)), 3005) + fmt <- fly_camera_format(photos) + lapply(seq_len(nrow(fp)), function(i) { + list( + ring = sf::st_coordinates(sf::st_geometry(fp)[[i]])[1:4, 1:2, drop = FALSE], + dims = image_dims(fmt[i, ]), + camera = fmt$camera[i] + ) + }) +} + + +test_that("the shipped sensor aspect matches the delivered thumbnail aspect", { + # The premise the invariant below rests on, asserted beside it. If a regenerated + # `camera_formats.csv` ever disagrees with the images the catalogue actually serves, + # this fails naming the cause instead of the behaviour test failing and blaming the + # corner mapping. + photos <- sf::st_read(testdata_path("photo_centroids_digital.gpkg"), quiet = TRUE) + fmt <- fly_camera_format(photos) + + measured <- c("Leica DMC II" = 972 / 884, "UltraCam Eagle M3" = 1654 / 1063) + shipped <- tapply(fmt$px_cross / fmt$px_along, fmt$camera, function(x) x[1]) + + expect_setequal(names(shipped), c("DMC II", "UltraCam Eagle M3")) + expect_equal(unname(shipped[["DMC II"]]), unname(measured[["Leica DMC II"]]), + tolerance = 1e-3) + expect_equal(unname(shipped[["UltraCam Eagle M3"]]), + unname(measured[["UltraCam Eagle M3"]]), tolerance = 1e-3) +}) + + +test_that("rotations 90 and 270 map a digital frame isotropically", { + cases <- digital_rings() + expect_gt(length(cases), 0) + + for (rot in c(90, 270)) { + aniso <- vapply(cases, function(c_) { + gcp_anisotropy(c_$dims$ncol_px, c_$dims$nrow_px, c_$ring, rot) + }, numeric(1)) + expect_equal(aniso, rep(1, length(cases)), tolerance = 1e-6, + info = paste("rotation =", rot)) + } +}) + + +test_that("rotations 0 and 180 squash a digital frame, and by how much", { + # The negative half. Without it the test above passes for an implementation that + # returns an isotropic mapping for every rotation, which is the one thing it must not + # do — see the restore-the-bug check in `test-fly_georef.R`. + # + # The distortion is the footprint's aspect ratio squared, so it is uneven across the + # bundled cameras and the numbers are pinned rather than left to a threshold: + # DMC II 1.100^2 = 1.21 <- weak; would survive a loose tolerance + # UltraCam Eagle M3 1.556^2 = 2.42 <- the case that makes this discriminating + # A fixture change that drops the UltraCam frames fails here. + cases <- digital_rings() + by_cam <- split(cases, vapply(cases, function(c_) c_$camera, character(1))) + expect_setequal(names(by_cam), c("DMC II", "UltraCam Eagle M3")) + + expected <- c("DMC II" = 1.21, "UltraCam Eagle M3" = 2.42) + for (cam in names(by_cam)) { + for (rot in c(0, 180)) { + aniso <- vapply(by_cam[[cam]], function(c_) { + gcp_anisotropy(c_$dims$ncol_px, c_$dims$nrow_px, c_$ring, rot) + }, numeric(1)) + expect_equal(aniso, rep(expected[[cam]], length(aniso)), tolerance = 0.01, + info = paste(cam, "rotation =", rot)) + } + } +}) + + +test_that("square film is isotropic at every rotation, so the invariant is vacuous there", { + # Stated rather than assumed. A square footprint with a square scan cannot distinguish + # any of the four rotations, which is exactly why the film mapping needed calibrating + # against imagery and why this file cannot finish the job for digital either. + photos <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE) + fp <- sf::st_transform(suppressWarnings(fly_footprint(photos[1, ])), 3005) + ring <- sf::st_coordinates(sf::st_geometry(fp)[[1]])[1:4, 1:2, drop = FALSE] + + for (rot in c(0, 90, 180, 270)) { + expect_equal(gcp_anisotropy(1250L, 1250L, ring, rot), 1, tolerance = 1e-6, + info = paste("rotation =", rot)) + } +}) From aa316f8f3efb58a4d1f9ef4bd045a67c67147b98 Mon Sep 17 00:00:00 2001 From: almac2022 Date: Sun, 30 Aug 2026 15:18:53 -0700 Subject: [PATCH 4/7] Establish the digital corner mapping as rotation 270 (#38) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The issue assumed this needed licence-restricted orthophoto imagery. It does not. The catalogue publishes per-frame exterior orientation through `patb_georef_url` — a column `fly_fetch()` already supports — and consecutive frames overlap enough to check each other, so the whole derivation is public and lives in `data-raw/georef_calibrate-corner_mapping.R`. Three independent measurements agree on 270 for both bundled cameras: the UltraCam Eagle's mount is rigid to 0.18 degrees over 6839 frames spanning the compass; adjacent-frame overlap correlation is +0.616 and +0.659 against at most +0.43 for the alternatives; and FWA lake polygons come out darkest at 270 on frames whose water sits 442 m and 1684 m off the footprint centre. One measurement disagreed and was wrong, which is the part worth keeping. Read naively the exterior orientation puts the DMC II at 90, on 97.6% agreement across four pooled projects — but that camera's project flies east and west only, so its own data separates the rigid and reflected hypotheses at 97.6% against 98.4%. The note records it as the cautionary half rather than dropping it. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01GBKqedyBysV7hB4DuL98ZR --- data-raw/georef_calibrate-corner_mapping.R | 196 +++++++++++++++++++++ inst/notes/georeferencing.md | 121 +++++++++++++ planning/active/progress.md | 31 ++++ planning/active/task_plan.md | 42 ++--- 4 files changed, 370 insertions(+), 20 deletions(-) create mode 100644 data-raw/georef_calibrate-corner_mapping.R create mode 100644 inst/notes/georeferencing.md diff --git a/data-raw/georef_calibrate-corner_mapping.R b/data-raw/georef_calibrate-corner_mapping.R new file mode 100644 index 0000000..4952245 --- /dev/null +++ b/data-raw/georef_calibrate-corner_mapping.R @@ -0,0 +1,196 @@ +# georef_calibrate-corner_mapping.R — establish the pixel-to-ground corner mapping +# for a footprint already rotated onto its flight line (fly#38). +# +# Everything here is public: BC Data Catalogue metadata, the airphoto thumbnails at +# openmaps.gov.bc.ca, the exterior-orientation files the catalogue links through +# `patb_georef_url`, and FWA lakes. No licence-restricted imagery is used or needed. +# +# Three independent measurements, run in this order. Read `inst/notes/georeferencing.md` +# before changing any of them — the second and third are the ones that decide, and the +# first is the one that looks most authoritative while being the weakest. +# +# Usage: source interactively. Writes nothing to the package. + +library(sf) +library(terra) +devtools::load_all() + +sf::sf_use_s2(FALSE) +work <- file.path(tempdir(), "fly38") +dir.create(work, recursive = TRUE, showWarnings = FALSE) + +photos <- st_read(system.file("testdata/photo_centroids_digital.gpkg", package = "fly"), + quiet = TRUE) +fp <- st_transform(suppressWarnings(fly_footprint(photos)), 3005) + +deg <- function(r) r * 180 / pi +az <- function(e, n) deg(atan2(e, n)) %% 360 # clockwise from north +wrap <- function(x) ((x + 180) %% 360) - 180 + +# --------------------------------------------------------------------------- +# 1. Exterior orientation, from the PATB files the catalogue publishes per project +# --------------------------------------------------------------------------- +# Answers: does the camera's image x-axis track the flight line, and at what offset? +# +# The control that matters is `sum` versus `difference`. A rigid mount gives a constant +# `image_x_azimuth - heading`; a ground frame read with its axes swapped gives a constant +# `image_x_azimuth + heading` instead. The two are indistinguishable unless the project +# flies more than one heading — which is why the offset is only established on projects +# whose legs span the compass, and why the bundled DMC II project on its own cannot +# settle anything. + +heading_per_frame <- function(roll, frame, east, north) { + out <- rep(NA_real_, length(roll)) + for (r in unique(roll)) { + i <- which(roll == r) + i <- i[order(frame[i])] + if (length(i) < 2) next + de <- diff(east[i]); dn <- diff(north[i]) + step <- sqrt(de^2 + dn^2) + a <- az(de, dn) + a[step < 100 | step > 5000] <- NA # turns and roll breaks + out[i] <- c(a, a[length(a)]) + } + out +} + +report_offset <- function(label, image_x_az, heading) { + ok <- is.finite(image_x_az) & is.finite(heading) + d <- wrap(image_x_az[ok] - heading[ok]) + s <- (image_x_az[ok] + heading[ok]) %% 360 + tight <- function(v) 100 * mean(abs(wrap(v - median(v))) <= 5) + message(sprintf( + "%s n=%d compass bins=%d\n rigid median %8.2f within +/-5: %5.1f%%\n reflected median %8.2f within +/-5: %5.1f%%", + label, sum(ok), length(unique(round(heading[ok] / 30) * 30 %% 360)), + median(d), tight(d), median(s), tight(s))) +} + +read_patb <- function(url) { + dest <- file.path(work, basename(url)) + if (!file.exists(dest)) utils::download.file(url, dest, quiet = TRUE) + if (grepl("\\.zip$", url)) { + utils::unzip(dest, exdir = file.path(work, tools::file_path_sans_ext(basename(url)))) + f <- list.files(file.path(work, tools::file_path_sans_ext(basename(url))), + pattern = "georef\\.(csv|txt)$", full.names = TRUE) + utils::read.csv(f[1]) + } else { + utils::read.csv(dest) + } +} + +# UltraCam Eagle M3. `eop_x` holds the NORTHING and `eop_y` the easting despite the +# names, and kappa is nevertheless referenced to a standard (easting, northing) frame — +# which is why the column names cannot be trusted as a statement of convention and the +# sum/difference control has to do the work. +eagle <- read_patb("https://openmaps.gov.bc.ca/thumbs/patb_files/d_005_emn_19_georef.zip") +rf <- strsplit(eagle$roll_frame, "_", fixed = TRUE) +eagle$roll <- vapply(rf, `[`, "", 1) +eagle$frame <- as.integer(vapply(rf, `[`, "", 2)) +eagle$head <- heading_per_frame(eagle$roll, eagle$frame, eagle$eop_y, eagle$eop_x) +report_offset("UltraCam Eagle M3", (90 - eagle$kappa) %% 360, eagle$head) + +# Leica DMC II. `gr_omega/gr_phi/gr_kappa` are zero in every row; the rotation is +# delivered as a 3x3 matrix, and `c2` is corrupt (`"00000000000"`) throughout. Only the +# first column is needed. The bundled frames' own project flies east-west only, so four +# further DMC II projects are pooled to get compass coverage. +dmc_urls <- c( + "https://openmaps.gov.bc.ca/thumbs/patb_files/d_003_fi_13_georef.zip", # the bundled one + "https://openmaps.gov.bc.ca/thumbs/patb_files/d_003_fi_14_georef.zip", + "https://openmaps.gov.bc.ca/thumbs/patb_files/d_003_fi_15_georef.zip", + "https://openmaps.gov.bc.ca/thumbs/patb_files/d_003_fi_16_georef.zip" +) +for (u in dmc_urls) { + d <- read_patb(u) + rf <- strsplit(d$frm_roll_frame, "_", fixed = TRUE) + d$roll <- vapply(rf, `[`, "", 1) + d$frame <- as.integer(vapply(rf, `[`, "", 2)) + d$head <- heading_per_frame(d$roll, d$frame, d$gr_easting, d$gr_northing) + x <- az(suppressWarnings(as.numeric(d$a1)), suppressWarnings(as.numeric(d$b1))) + report_offset(paste("Leica DMC II —", basename(u)), x, d$head) +} + +# --------------------------------------------------------------------------- +# 2. Adjacent-frame overlap — the measurement that decides +# --------------------------------------------------------------------------- +# Consecutive frames on a line overlap heavily, so at the correct rotation their common +# ground must agree. A 180-degree error reflects each frame about its OWN centre, and +# because the centres differ the overlap then shows different ground. Needs no reference +# imagery of any kind: the frames check each other. + +georef_at <- function(idx, tag, rot) { + th <- fly_fetch(photos[idx, ], type = "thumbnail", + dest_dir = file.path(work, tag, "thumb")) + stopifnot(all(th$success)) + d <- file.path(work, tag, paste0("r", rot)) + dir.create(d, recursive = TRUE, showWarnings = FALSE) + vapply(seq_along(idx), function(k) { + o <- file.path(d, paste0(k, ".tif")) + ok <- georef_one(th$dest[k], fp[idx[k], ], o, srcnodata = "0", rotation = rot) + if (isTRUE(ok)) o else NA_character_ + }, character(1)) +} + +grey <- function(path) { + r <- rast(path) + g <- if (nlyr(r) >= 3) mean(r[[1:3]]) else r[[1]] + g[g <= 0] <- NA + g +} + +pair_r <- function(a, b, res = 25) { + ga <- grey(a); gb <- grey(b) + inter <- terra::intersect(ext(ga), ext(gb)) + if (is.null(inter)) return(NA_real_) + tmpl <- rast(inter, resolution = res, crs = crs(ga)) + va <- values(resample(ga, tmpl, method = "average"))[, 1] + vb <- values(resample(gb, tmpl, method = "average"))[, 1] + ok <- is.finite(va) & is.finite(vb) + if (sum(ok) < 500) return(NA_real_) + suppressWarnings(stats::cor(va[ok], vb[ok])) +} + +for (case in list(list(i = 19:24, tag = "eagle", cam = "UltraCam Eagle M3"), + list(i = 1:6, tag = "dmc", cam = "Leica DMC II"))) { + message("\n", case$cam, " — adjacent-frame overlap correlation") + for (rot in c(0, 90, 180, 270)) { + o <- georef_at(case$i, case$tag, rot) + rs <- vapply(seq_len(length(o) - 1), function(k) pair_r(o[k], o[k + 1]), numeric(1)) + message(sprintf(" rotation %3d : mean r = %+.3f [%s]", + rot, mean(rs, na.rm = TRUE), paste(sprintf("%+.2f", rs), collapse = " "))) + } +} + +# --------------------------------------------------------------------------- +# 3. Water darkness against FWA lakes — an outside opinion +# --------------------------------------------------------------------------- +# Independent of both measurements above: lakes are dark, and FWA knows exactly where +# they are. Only discriminating when the water sits off-centre, since a 180-degree error +# rotates about the footprint centre — so the offset is reported beside the result rather +# than assumed. + +bb <- st_bbox(st_union(st_geometry(fp))) +lakes <- st_transform( + bcdata::collect(bcdata::filter( + bcdata::bcdc_query_geodata("WHSE_BASEMAPPING.FWA_LAKES_POLY"), + BBOX(!!as.numeric(bb), crs = "EPSG:3005"))), 3005) +water <- st_union(lakes) + +for (i in c(1, 2)) { + g <- st_geometry(fp)[i] + w <- st_intersection(water, g) + if (length(w) == 0) next + off <- as.numeric(st_distance(st_centroid(st_union(w)), st_centroid(g))) + message(sprintf("\nframe %d (%s_%s) — water %.0f m off the footprint centre", + i, photos$film_roll[i], photos$frame_number[i], off)) + th <- fly_fetch(photos[i, ], type = "thumbnail", dest_dir = file.path(work, "water")) + for (rot in c(0, 90, 180, 270)) { + o <- file.path(work, "water", sprintf("w%d_%d.tif", rot, i)) + if (!isTRUE(georef_one(th$dest[1], fp[i, ], o, srcnodata = "0", rotation = rot))) next + lum <- grey(o) + inw <- terra::extract(lum, vect(st_as_sf(st_sfc(w, crs = 3005))), ID = FALSE)[[1]] + all <- values(lum)[, 1] + message(sprintf(" rotation %3d : water %6.1f frame %6.1f difference %+7.1f", + rot, mean(inw, na.rm = TRUE), mean(all, na.rm = TRUE), + mean(inw, na.rm = TRUE) - mean(all, na.rm = TRUE))) + } +} diff --git a/inst/notes/georeferencing.md b/inst/notes/georeferencing.md new file mode 100644 index 0000000..1663986 --- /dev/null +++ b/inst/notes/georeferencing.md @@ -0,0 +1,121 @@ +# Why the corner mapping is the way it is + +Companion to `terrain-correction.md` and `camera-formats.md`. Read this before changing +`fly_georef_gcps()`, the rotation `fly_georef()` applies to a non-square footprint, or +`data-raw/georef_calibrate-corner_mapping.R`. + +Established in fly#38 (2026-08-30) against published BC data. Nothing here needed +licence-restricted imagery. + +## The ring order is a contract, and it is in the rectangle's own frame + +`fly_rectangles()` emits vertices in this order, in local coordinates where `+y` is the +flight heading and `+x` is 90 degrees clockwise of it: + +| ring row | local | flight-relative | +|---|---|---| +| 1 | `(-hc, -ha)` | rear-left | +| 2 | `( hc, -ha)` | rear-right | +| 3 | `( hc, ha)` | front-right | +| 4 | `(-hc, ha)` | front-left | + +Rotation onto the bearing is applied to the whole rectangle, so **the ring means the same +thing whether or not it was rotated**. That is what lets one constant serve both the +bearing-rotated frames and the `axis_aligned_no_bearing` ones, and it is why +`fly_georef()` must not apply `bearing_to_rotation()` to a non-square footprint: the +bearing is already in the ring, and applying it again counts it twice. + +## The mapping is rotation 270: top-left pixel to the rear-left ground corner + +Equivalently: **image columns run in the flight direction, image rows run flight-right.** + +Three independent measurements. They are listed in the order that makes the weakest one +easiest to over-trust, which is the order they were run in. + +### 1. Exterior orientation — authoritative-looking, and the one that misled + +The catalogue publishes per-frame exterior orientation through `patb_georef_url`, which +`fly_fetch(type = "georef")` already fetches. It is free, public, and covers every +digital frame. + +The measurement is the offset between the camera's image x-axis azimuth and the flight +heading computed from the frame's own neighbours. **The control that matters is sum +versus difference**: a rigid mount holds `image_x_azimuth - heading` constant, while a +ground frame whose axes are read in the wrong order holds `image_x_azimuth + heading` +constant instead. The two are indistinguishable on a project that flies one axis. + +| camera | n | compass bins | rigid | reflected | +|---|---|---|---|---| +| UltraCam Eagle M3 | 6839 | 32 | median **0.18**, MAD 0.47, 98.7% within 5 | median 254.7, MAD 100.6, **14.1%** | +| Leica DMC II (bundled project) | 5450 | 8, all E/W | median −177.7, 97.6% | median 357.0, **98.4%** | +| Leica DMC II (3 further projects) | 7924 | 13 | median ≈ −178, 22–96% | 2–17% | + +So the Eagle's mount is rigid to half a degree with its image x-axis **along the heading**, +and reflection is excluded outright. Under the ordinary top-left raster convention that +gives rotation 270. + +The DMC II reading gives an offset of 180 degrees instead, which would put it at rotation +90 — and it is **wrong**. Note what it looks like: 97.6% agreement on the bundled project, +a clean-looking median, four projects pooled. What that number cannot show is that the +bundled project flies east and west only, so its own data cannot separate the two +hypotheses at all (97.6% against 98.4%), and the file that produces it has +`gr_omega/gr_phi/gr_kappa` zeroed, an undocumented 3x3 matrix in their place, and a +`c2` column that is the literal string `00000000000` in every row. It is a chain of +assumed conventions wearing a large sample size. + +### 2. Adjacent-frame overlap — the measurement that decides + +Consecutive frames on a line overlap heavily, so at the correct rotation their common +ground must agree. A 180-degree error reflects each frame about **its own** centre, and +because the centres differ the overlap then shows different ground. + +This needs no reference imagery of any kind. The frames check each other. + +Mean Pearson r between adjacent georeferenced frames, on a common 25 m grid: + +| rotation | UltraCam Eagle M3 | Leica DMC II | +|---|---|---| +| 0 | +0.095 | +0.356 | +| 90 | +0.069 | +0.429 | +| 180 | +0.135 | +0.331 | +| **270** | **+0.616** | **+0.659** | + +Both cameras, unambiguously, and the same answer. + +### 3. Water darkness — an outside opinion + +Lakes are dark and FWA knows exactly where they are. Only discriminating when the water +sits off-centre, since a 180-degree error rotates about the footprint centre — so the +offset is reported beside the result rather than assumed. + +Mean luminance inside FWA lake polygons minus the frame mean: + +| rotation | frame 778 (water 442 m off-centre) | frame 779 (1684 m off-centre) | +|---|---|---| +| 0 | +0.7 | +16.9 | +| 90 | −12.1 | +1.1 | +| 180 | −5.4 | −5.6 | +| **270** | **−23.4** | **−9.5** | + +## What each check can and cannot catch + +- **The aspect invariant** (`tests/testthat/test-fly_georef_aspect.R`) rejects rotations + 0 and 180, because those pair the image's long axis with the footprint's short edge. + It **cannot** separate 270 from 90 — they differ by 180 degrees about the centre and a + rectangle is symmetric under that. It is also vacuous on square film, which is the + whole of the pre-existing georef suite. +- **It cannot catch a mirrored mapping either.** Traversing the ring the other way leaves + every edge length and the aspect ratio identical and produces a reflection. The signed + area of the ground quad in pixel order is what separates those, and it is asserted. +- **The overlap and water checks** both need real thumbnails, so they live in + `data-raw/`, not in the test suite. + +## The one thing that would invalidate this + +All three measurements are downstream of `camera_formats.csv` assigning the long sensor +axis to the **across-track** direction. Get that wrong and the footprint itself is +rotated 90 degrees on the ground, and no corner mapping repairs it. Two things guard it: +`parse_vexcel()` reads the manufacturer's own `long track` / `cross track` labels rather +than inferring, and `fly_georef()` refuses any frame whose delivered image aspect does not +pair isotropically with its footprint edges — which is also what catches a frame sized +from an inferred format that does not match its real camera. diff --git a/planning/active/progress.md b/planning/active/progress.md index cd63f83..87acd9c 100644 --- a/planning/active/progress.md +++ b/planning/active/progress.md @@ -49,3 +49,34 @@ `asNamespace("fly")` and `as.environment("package:fly")`, with a printed ground coordinate that can only come from the broken version. Result: **FAIL=2**, the two isotropy assertions. The guard fires + +### Phase 3 — the constant is 270, and the plan's premise was wrong + +A concurrent plan review found that `patb_georef_url` — a column `fly_fetch()` already +supports — carries **per-frame exterior orientation**, published and free. So the +licence-restricted orthophotos were never the only route, and are not the best one. +Verified before acting on it: the file downloads, parses, and covers all 24 bundled +frames. + +Three measurements, all public, all in `data-raw/georef_calibrate-corner_mapping.R`: + +1. **Exterior orientation.** The Eagle's image x-axis tracks the flight heading at + 0.18 deg (MAD 0.47) over 6839 frames spanning 32 compass bins, with the reflected + reading excluded at 14.1% against 98.7%. Under the ordinary top-left raster + convention that is rotation 270. +2. **Adjacent-frame overlap.** Needs no reference imagery — consecutive frames overlap, + so at the right rotation their common ground agrees, and a 180-degree error reflects + each frame about its own centre. 270 wins on both cameras: +0.616 and +0.659 against + at most +0.43. +3. **FWA lake darkness.** 270 makes water darkest on both frames tested, with the water + 442 m and 1684 m off the footprint centre so the test could actually discriminate. + +**The PATB reading disagreed for the DMC II and was wrong.** It said 90. It looked +strong — 97.6% agreement, four projects pooled — and could not have been right: that +camera's bundled project flies east and west only, so its own data separates the two +hypotheses at 97.6% against 98.4%, which is to say not at all. Recorded in +`inst/notes/georeferencing.md` as the cautionary half, because it is the measurement a +reader will find most convincing. + +Both cameras give the same answer, so the constant is global rather than a per-camera +column. diff --git a/planning/active/task_plan.md b/planning/active/task_plan.md index 3c2dec9..06c4ee2 100644 --- a/planning/active/task_plan.md +++ b/planning/active/task_plan.md @@ -102,26 +102,28 @@ pulled from git, and fails if the extraction changed any coordinate. the DMC II / UltraCam distortion factors are recorded so a future fixture change that drops the UltraCam is visible. -## Phase 3: Establish the constant against orthophoto ground truth - -Runs in the private catalogue repo. Two facts settle together: - -- [ ] (a) Is the long sensor axis actually across-track? Compare the #32 footprint's - orientation against the ortho for the same ground. A 90-degree error here is not a - georef bug — it is a footprint bug, and it would reopen #32. -- [ ] (b) Which of 90 / 270 is correct? They differ by a 180-degree flip of image content - inside the same quad. -- [ ] Sample both bundled cameras plus the 1.76:1 Leica DMC III — the highest-ratio and - therefore most diagnostic case. Join ortho to catalogue on `BCGS_TILE` + `PHOTO_YEAR`. -- [ ] Check visually **and** with a measurable: a named ground feature must land in the - same place in the warped thumbnail as in the ortho. A screenshot alone is not a result. -- [ ] Record the answer **per camera**. Do not assume one constant serves both until - both are measured. If they differ, the constant becomes a `camera_formats.csv` column - and Phase 4 grows a lookup — decide then, not now. -- [ ] Bring back only: the constant(s), the sample size, and what was compared. No pixels. - -**Verify:** the same answer for every frame of a camera; disagreement within a camera means -the model is wrong, not that the sample is noisy. +## Phase 3: Establish the constant [DONE — re-scoped by measurement] + +The plan assumed licence-restricted orthophotos were the only route. They are not, and +the route that works is better: the catalogue publishes per-frame exterior orientation +through `patb_georef_url`, and consecutive frames overlap enough to check each other. +All three measurements below use public data, so the derivation lives in fly. + +- [x] (a) Is the long sensor axis actually across-track? Guarded rather than assumed — + `parse_vexcel()` reads the manufacturer's own `long track`/`cross track` labels, + and Phase 4's isotropy guard refuses any frame whose delivered aspect disagrees + with its footprint +- [x] (b) Which rotation is correct? **270**, for both cameras +- [x] Exterior orientation from PATB: the Eagle's mount is rigid to 0.18 deg (MAD 0.47, + n=6839, 32 compass bins), reflection excluded at 14.1% against 98.7% +- [x] Adjacent-frame overlap correlation — needs no reference imagery at all: + 270 wins on both cameras (+0.616 / +0.659 against <= +0.43) +- [x] FWA lake darkness as an outside opinion: 270 darkest on both frames tested, + water 442 m and 1684 m off-centre so the test could discriminate +- [x] Recorded per camera. They agree, so the constant is global — **not** a + `camera_formats.csv` column +- [x] `data-raw/georef_calibrate-corner_mapping.R` reproduces all three +- [x] No imagery, licence-restricted or otherwise, enters the package ## Phase 4: Land the constant and remove the exclusion From ab60389cd9dad3b75054e19f3be651f6e7bb8842 Mon Sep 17 00:00:00 2001 From: almac2022 Date: Sun, 30 Aug 2026 16:12:50 -0700 Subject: [PATCH 5/7] Georeference digital frames (#38) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Removes the exclusion that made every non-square footprint unreachable, so the post-2010 catalogue is georeferenceable rather than only sizeable. A non-square footprint is already rotated onto its flight line by `fly_rectangles()`, so `bearing_to_rotation()` is not applied to it. Those frames use the measured mapping instead: the top-left pixel maps to the ring's rear-left corner. Film keeps the bearing rule unchanged, and a `rotation` column still overrides per-photo for both. A frame whose delivered image aspect disagrees with its footprint's by more than 8% is skipped with a warning rather than written stretched. A wrong mapping produces a valid GeoTIFF, in the right CRS, over the right ground, squashed by the aspect ratio squared — which nothing downstream reports. The threshold is computed, not chosen: the band is (1.0667, 1.0995), between a full-resolution 9-inch scan carrying the negative's rebate and a Leica DMC II frame sized through `format_size` onto a square footprint. Checked against every row of the shipped camera table. Also warns when a non-square footprint has no flight bearing, which is the ordinary result of georeferencing a single frame, since `fly_bearing()` needs a neighbour. Three rounds of code-check, each of which found a defect inside the previous round's fix: 1 the guard refused non-square FILM scans, where anisotropy is the image's own inverse aspect and no rotation changes it — a 9600x9000 scan was rejected 2 the shape gate that fixed round 1 switched the guard off for a digital frame sized onto a square footprint, the exact case it exists for; and the `rotation` column had three read sites where the fix converted two, so an unparseable factor level was applied as its level code 3 the tolerance replacing that gate was 0.4% too loose, and the test guarding it asserted against the most eccentric camera rather than the least Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01GBKqedyBysV7hB4DuL98ZR --- CLAUDE.md | 21 +- R/fly_georef.R | 263 ++++++++++++++---- README.md | 2 +- data-raw/georef_calibrate-corner_mapping.R | 21 +- inst/notes/georeferencing.md | 29 ++ man/fly_georef.Rd | 61 ++++- planning/active/progress.md | 149 +++++++++++ planning/active/review-round1.md | 244 +++++++++++++++++ planning/active/review-round2.md | 286 ++++++++++++++++++++ planning/active/review-round3.md | 294 +++++++++++++++++++++ planning/active/task_plan.md | 32 +-- tests/testthat/test-fly_georef_digital.R | 289 ++++++++++++++++++++ tests/testthat/test-fly_georef_gcps.R | 57 ++++ vignettes/airphoto-selection.Rmd | 9 + 14 files changed, 1669 insertions(+), 88 deletions(-) create mode 100644 planning/active/review-round1.md create mode 100644 planning/active/review-round2.md create mode 100644 planning/active/review-round3.md create mode 100644 tests/testthat/test-fly_georef_digital.R diff --git a/CLAUDE.md b/CLAUDE.md index fd374a7..f857680 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -31,6 +31,8 @@ compute `dem_coverage` that each passed their own tests `camera_formats_manifest.csv` as the offline drift guard - `inst/notes/camera-formats.md` — why `SCALE` is unusable for digital frames, what each of the five QA checks can and cannot catch, and the PDF extraction traps these specific reports carry +- `inst/notes/georeferencing.md` — the ring-order contract, the three measurements that established the +digital corner mapping, and the one of them that was wrong while looking strongest - `mixed_media_fixture()` in `tests/testthat/setup.R` — synthesized frames whose format nothing resolves; `digital_fixture()` beside it covers every resolver branch, and `footprint_cases()` sweeps the 12 input shapes the invariant tests run over @@ -46,8 +48,17 @@ geometry column may be named `geom` not `geometry` recording format, so a 9-inch negative applied to a digital sensor is wrong by an unknown factor while still drawing as a rectangle and still producing a coverage percentage. `fly_footprint()` sizes each row from its `media` value, records the outcome in `footprint_basis`, and returns an **empty geometry** where the format is -unknown. Downstream functions report how many frames they excluded. Digital defaults are deliberately unshipped -until sensor widths are established (#32) — `format_size` is the slot they drop into +unknown. Downstream functions report how many frames they excluded. Sensor widths landed in #32 and +georeferencing in #38, so `footprint_basis` is now the record of *which* route sized a frame rather than a +list of what is missing + +- **The digital corner mapping is rotation 270, and it was measured, not reasoned** (v0.7.0, #38) — a +non-square footprint is already rotated onto its flight line by `fly_rectangles()`, so `fly_georef()` must +not apply `bearing_to_rotation()` on top; the top-left pixel maps to the ring's rear-left corner. Three +independent public routes agree (exterior orientation from `patb_georef_url`, adjacent-frame overlap +correlation, FWA lake darkness). **The geometry alone cannot get there** — the aspect invariant narrows +the answer to 270 or 90 and a rectangle is symmetric under the 180 degrees between them, so anyone +re-deriving this from the ring will produce a plausible wrong answer. Read `inst/notes/georeferencing.md` - **Terrain error is a datum offset, not slope** (v0.5.0, #9) — `FLYING_HEIGHT` is metres **above sea level**, and reported scale is referenced to an elevation above the ground the photos cover, so it understates footprint @@ -85,6 +96,12 @@ a tibble, so `footprint_basis`, `footprint_terrain`, `height_agl` and `dem_cover documented data source, with geometry and every downstream number still correct. Use `centroid_shapes()` in `tests/testthat/setup.R` — it sweeps plain / tibble / grouped / `bcdc_sf` — for anything that attaches columns to user-supplied data. Note the class *set* is carried but not its order: `st_transform()` moves `sf` to the front +- **`local_mocked_bindings(.env = )` is the cleanup environment, not the target.** `.package` names the +package to mock in; `.env` says what the mock unwinds with. Passing `asNamespace("fly")` to `.env` installs +the stub correctly and then never removes it, because a namespace does not exit — so every later test in +the run keeps it. It leaks in the direction that reads as success, since a stub returning `TRUE` makes +assertions pass. Caught in #38 only because a later test asserted a file existed that the stub never wrote + - **`fly_footprint()` must not be handed its own output** (fly#37, open) — `st_coordinates()` on POLYGON returns one row per vertex, so 20 footprints in gives 100 rows out, silently. There is no guard yet diff --git a/R/fly_georef.R b/R/fly_georef.R index c21f012..94995ca 100644 --- a/R/fly_georef.R +++ b/R/fly_georef.R @@ -22,11 +22,16 @@ #' frame from its height above ground instead of the reported scale. See the #' **Terrain** section of [fly_footprint()]. #' @param rotation Image rotation in degrees clockwise. One of `"auto"`, -#' `0`, `90`, `180`, or `270`. `"auto"` (default) computes flight line -#' bearing from consecutive centroids and derives rotation per-photo — -#' requires `film_roll` and `frame_number` columns. Fixed values apply -#' the same rotation to all photos. Overridden per-photo if `photos_sf` -#' contains a `rotation` column. +#' `0`, `90`, `180`, or `270`. Applies to frames with a **square** +#' footprint only. `"auto"` (default) computes flight line bearing from +#' consecutive centroids and derives rotation per-photo — requires +#' `film_roll` and `frame_number` columns. Fixed values apply the same +#' rotation to every square-footprint photo. A non-square footprint +#' ignores this argument and uses the measured digital mapping (see +#' **Rotation**); a `rotation` column in `photos_sf` still overrides +#' per-photo, for both. Carrying a film-era `rotation` column into a +#' batch of digital frames therefore overrides the correct mapping with +#' the wrong one — drop the column, or set it to `NA` for those rows. #' @return A tibble with columns `airp_id`, `source`, `dest`, and `success`. #' #' @details @@ -35,14 +40,15 @@ #' translates the image with GCPs then warps to the target CRS using #' bilinear resampling. #' -#' **Rotation:** Aerial photos may appear rotated in their footprints -#' because the camera orientation relative to north varies by flight -#' direction, camera mounting, and scanner orientation. The `rotation` -#' parameter rotates the GCP corner mapping: +#' **Rotation:** the corner mapping depends on the footprint's shape, because +#' [fly_footprint()] builds the two shapes differently. +#' +#' A **square** footprint is axis-aligned, so the mapping carries the rotation. +#' The `rotation` parameter rotates it: #' \itemize{ -#' \item `0` — top of image maps to north edge of footprint (original behavior) +#' \item `0` — top of image maps to north edge of footprint #' \item `90` — top of image maps to east edge (90° clockwise) -#' \item `180` — top of image maps to south edge (default, correct for most BC photos) +#' \item `180` — top of image maps to south edge (correct for most BC film) #' \item `270` — top of image maps to west edge #' } #' @@ -60,10 +66,41 @@ #' ``` #' photos$rotation <- dplyr::case_when( #' photos$film_roll == "bc5282" ~ 270, -#' .default = NA # fall through to auto +#' .default = NA # non-square footprints fall through to the digital mapping; +#' # square ones to `rotation`, which is 180 under "auto" — +#' # NOT back to the per-photo bearing #' ) #' ``` #' +#' Every non-`NA` value in that column must be 0, 90, 180 or 270; anything else +#' is an error naming the value, rather than a rotation silently applied as +#' something other than what was written. +#' +#' A **non-square** footprint — every digital frame — is already rotated onto +#' its flight line by [fly_footprint()], so the ring carries the bearing and +#' applying it again would count it twice. Those frames use a fixed mapping +#' instead: the top-left pixel maps to the ring's rear-left corner, equivalently +#' image columns run in the flight direction and image rows run flight-right. +#' That was measured in fly#38 on both bundled cameras by three independent +#' routes and is the same for each; see `inst/notes/georeferencing.md`. +#' +#' A frame whose delivered image aspect disagrees with its footprint's by more +#' than 8% is **skipped with a warning** rather than written stretched — the +#' failure a wrong mapping produces is a valid GeoTIFF over the right ground, +#' squashed by the aspect ratio squared, which nothing downstream would report. +#' The threshold sits between the largest disagreement a legitimate frame +#' produces (a full-resolution 9-inch scan carrying the negative's rebate, about +#' 6.7%) and the smallest that must be caught (a Leica DMC II frame sized through +#' `format_size` onto a square footprint, 9.95%). It applies to square footprints +#' too: a square one has no pairing to get wrong, but a digital frame sized +#' through `format_size` lands on one, and that is the unknown-camera case — so +#' gating on shape would switch the check off exactly where it is needed. +#' +#' A non-square footprint built without a flight bearing is drawn axis-aligned +#' and is therefore georeferenced as though the flight line ran due north. +#' [fly_bearing()] needs a neighbouring frame, so this is the ordinary result of +#' georeferencing a single frame on its own, and it is warned about. +#' #' **Nodata handling:** Two sources of unwanted black pixels are masked: #' #' 1. **Warp fill** — GDAL creates black pixels outside the rotated source @@ -119,27 +156,29 @@ fly_georef <- function(fetch_result, photos_sf, footprints <- fly_footprint(photos_sf, dem = dem) |> sf::st_transform(3005) fly_warn_unsized(footprints, "georeferencing") - # `georef_one()` maps image corners onto footprint corners positionally and shifts - # that mapping by a 90-degree-quantized bearing — a scheme calibrated against north-up - # 9x9 negatives. A non-square footprint breaks it twice over: the footprint is already - # rotated onto its flight line so the bearing would be counted a second time, and a - # wrong-by-90 mapping that is harmless on a square maps a landscape image onto a - # portrait quad on a 1.76:1 rectangle. - # - # Skipped rather than guessed at. The right corner mapping for a pre-rotated - # rectangle depends on camera mounting relative to flight direction, which is what - # `bearing_to_rotation()` was empirically fitted to and cannot be re-derived without - # imagery to check against. Digital frames had no footprint at all before fly#32, so - # this is the same coverage as before rather than a regression — but explicit now - # instead of silent. See fly#38. - rotated <- !fly_is_square(footprints) - if (any(rotated)) { - warning( - sum(rotated), " of ", nrow(footprints), " frames have a non-square footprint ", - "and are excluded from georeferencing: the corner mapping is calibrated for ", - "square, axis-aligned footprints. See fly#38.", - call. = FALSE - ) + # Classify every footprint once, before any per-row work, and three ways rather than + # two. `fly_is_square()` reports an EMPTY geometry as square — defensible, since a + # shapeless thing has no unequal sides — so squareness alone cannot stand in for + # "has ground control". Keeping empty separate means neither branch below inherits the + # other's answer, and `georef_one()` is never handed a ring it cannot index. + empty_fp <- sf::st_is_empty(sf::st_geometry(footprints)) + non_square <- !empty_fp & !fly_is_square(footprints) + + # A non-square footprint that never got a bearing is drawn axis-aligned, so it is + # georeferenced as though the aircraft flew due north. `fly_bearing()` needs a + # neighbour, so this is the ordinary result of georeferencing one frame on its own — + # not an exotic case. Named rather than left to be discovered from a rotated image. + if ("width_source" %in% names(footprints)) { + no_bearing <- non_square & grepl("axis_aligned_no_bearing", footprints$width_source) + if (any(no_bearing)) { + warning( + sum(no_bearing), " of ", nrow(footprints), " frames have a non-square footprint ", + "but no flight bearing, so they are drawn and georeferenced as though the ", + "flight line ran due north. Pass neighbouring frames from the same roll so ", + "`fly_bearing()` can compute an azimuth.", + call. = FALSE + ) + } } # Match fetch results to photos by airp_id @@ -147,7 +186,40 @@ fly_georef <- function(fetch_result, photos_sf, # Per-photo rotation: column overrides auto/default - has_rotation_col <- "rotation" %in% names(photos_sf) + # `has_rotation_col` is set again by the auto path below, after which nothing can tell + # a user-supplied column from a bearing-derived one. Capture the distinction here, + # while it still exists: a user value overrides the measured digital rotation, a + # bearing-derived one must not. + user_rotation_col <- "rotation" %in% names(photos_sf) + has_rotation_col <- user_rotation_col + + # The argument is checked above; the column never was, and this function now makes it + # the highest-precedence input for a digital frame. `rotation %/% 90` turns 360 into a + # five-element shift and `fly_georef_gcps()` then indexes out of bounds — swallowed by + # the per-frame `tryCatch` below into a message naming neither the column nor the + # value. 45 and -90 are worse: they shift by zero and georeference silently wrong. + # Normalised once, here, and read from this vector everywhere below. `as.character()` + # first because `as.integer()` on a factor returns its level code — a `rotation` column + # read from a CSV as a factor would otherwise validate as 180 and then be *applied* as + # 1. Converting in two places is how those come apart. + user_rot <- NULL + if (user_rotation_col) { + raw <- as.character(photos_sf[["rotation"]]) + user_rot <- suppressWarnings(as.integer(raw)) + # A value that was supplied but did not parse must be refused, not quietly turned + # into NA. `is.na(user_rot) & !is.na(raw)` is that case; excluding it from `bad` + # would let "ninety" through as "no preference", which is a different instruction. + bad <- (!is.na(user_rot) & !user_rot %in% c(0L, 90L, 180L, 270L)) | + (is.na(user_rot) & !is.na(raw)) + if (any(bad)) { + stop("`photos_sf$rotation` must be NA or one of 0, 90, 180, 270. Got: ", + paste(unique(raw[bad]), collapse = ", "), ".", call. = FALSE) + } + # Written back so there is exactly ONE parse. Leaving the raw column in place and + # converting again at the read sites is how a factor validates as 180 and is then + # applied as its level code — the read below is not the only one. + photos_sf[["rotation"]] <- user_rot + } # Auto-compute bearing → rotation when needed if (auto_rotation && !has_rotation_col) { @@ -186,20 +258,30 @@ fly_georef <- function(fetch_result, photos_sf, # Find matching footprint fp_idx <- which(photos_sf[["airp_id"]] == results$airp_id[i]) if (length(fp_idx) == 0) next - # Warned about once, above, rather than per frame. - if (rotated[fp_idx[1]]) next - fp <- footprints[fp_idx[1], ] + j <- fp_idx[1] - # No footprint means no ground control to warp onto; leave success = FALSE - # rather than writing a GeoTIFF positioned by a format we could not resolve. - if (sf::st_is_empty(sf::st_geometry(fp))[1]) next + # No footprint means no ground control to warp onto; leave success = FALSE rather + # than writing a GeoTIFF positioned by a format we could not resolve. Checked before + # anything reads the ring. + if (empty_fp[j]) next + fp <- footprints[j, ] - # Per-photo rotation from column, or default - rot <- if (has_rotation_col) { - val <- as.integer(photos_sf[["rotation"]][fp_idx[1]]) + # A user-supplied `rotation` column overrides everything, square or not — the + # documented escape hatch. Everything else depends on the footprint's shape: + # a non-square ring already carries its bearing (see `fly_rectangles()`), so + # applying `bearing_to_rotation()` on top would count it twice. + user_val <- if (user_rotation_col) user_rot[j] else NA_integer_ + rot <- if (!is.na(user_val)) { + user_val + } else if (non_square[j]) { + fly_digital_rotation() + } else if (has_rotation_col) { + val <- as.integer(photos_sf[["rotation"]][j]) if (is.na(val)) { if (auto_rotation) 180L else rotation - } else val + } else { + val + } } else { rotation } @@ -221,10 +303,10 @@ fly_georef <- function(fetch_result, photos_sf, #' Georeference a single image to a footprint polygon #' @noRd georef_one <- function(src, fp, out_file, srcnodata = "0", rotation = 180) { - # Get footprint corner coordinates - # fly_footprint builds: BL, BR, TR, TL, BL (closing) + # Footprint ring, in the order `fly_rectangles()` guarantees: rows 1-4 are BL, BR, TR, + # TL **in the rectangle's own frame**. For a rotated (non-square) footprint that frame + # is the flight line's, so they are rear-left, rear-right, front-right, front-left. coords <- sf::st_coordinates(fp)[1:4, , drop = FALSE] - # coords: [1]=BL, [2]=BR, [3]=TR, [4]=TL # Read image dimensions and band count via GDAL info <- sf::gdal_utils("info", source = src, quiet = TRUE) @@ -240,10 +322,46 @@ georef_one <- function(src, fp, out_file, srcnodata = "0", rotation = 180) { # Build GCP args from the pixel-to-ground correspondence. gcp <- fly_georef_gcps(ncol_px, nrow_px, coords, rotation) + # Refuse a mapping that pairs the image's axes with the wrong footprint edges. It + # would still produce a valid GeoTIFF, in the right CRS, over the right ground — just + # squashed by the aspect ratio squared, which nothing downstream would report. Catches + # a camera delivering an orientation this was not measured on, and a frame sized from + # an inferred format that does not match the camera that actually took it. + # + # The tolerance is what makes this safe on film, and it is set from measurements + # rather than picked. A square footprint has no pairing to get wrong, but a scan + # carrying the negative's rebate is genuinely a few percent off square, and a gate on + # shape would switch the check off for a digital frame sized through `format_size` + # into a square footprint — exactly the unknown-camera case this exists to catch. + # + # |log| off isotropic what it is + # 0.0000 bundled film thumbnails, 1250 x 1250 + # 0.0645 a full-resolution 9-inch scan at 9600 x 9000 + # 0.0770 this tolerance + # 0.0949 the TIGHTEST case that must be caught — a Leica DMC II + # frame sized through `format_size` onto a square footprint + # 0.1898 the tightest mispairing on a frame's own footprint, + # the same DMC II at 1.0995 squared + # 0.4421 a portrait UltraCam frame on a square footprint + # + # The admissible band is narrow — (1.0667, 1.0995) — because a square-footprint DMC II + # frame is only slightly more eccentric than a badly rebated film scan. Checked against + # every row of `camera_formats.csv`, fallback rows included, not just the two cameras + # the bundled fixture happens to carry. + aniso <- fly_gcp_anisotropy(gcp, ncol_px, nrow_px) + if (!is.finite(aniso) || abs(log(aniso)) > log(fly_gcp_stretch_max())) { + warning(basename(src), ": image is ", round(ncol_px / nrow_px, 3), + ":1 but the corner mapping would stretch it by ", round(aniso, 3), + "x. Skipped rather than written squashed. ", + "See `inst/notes/georeferencing.md`.", + call. = FALSE) + return(FALSE) + } + # `unname()` matters: the ring arrives from `sf::st_coordinates()` carrying X/Y # dimnames, and without it the option vector handed to GDAL is named. Harmless to - # GDAL, but it makes the args unequal to a plain character vector under - # `identical()`, which is what the parity test compares. + # GDAL, but it makes the args unequal to a plain character vector under `identical()`, + # which is what the parity test compares. gcp_args <- character(0) for (j in seq_len(nrow(gcp))) { gcp_args <- c(gcp_args, "-gcp", as.character(unname(gcp[j, ]))) @@ -347,3 +465,50 @@ fly_georef_gcps <- function(ncol_px, nrow_px, coords, rotation) { dimnames(out) <- list(NULL, c("pixel_x", "pixel_y", "ground_x", "ground_y")) out } + +#' The rotation a digital frame's corner mapping needs +#' +#' 270 degrees: the top-left pixel maps to the footprint ring's rear-left corner, or +#' equivalently image columns run in the flight direction and image rows run +#' flight-right. Measured in fly#38 on both bundled cameras by three independent routes +#' — published exterior orientation, adjacent-frame overlap correlation, and FWA lake +#' darkness. They agree, so this is one constant rather than a per-camera column. +#' +#' Do not re-derive it by reasoning. See `inst/notes/georeferencing.md`, and +#' `data-raw/georef_calibrate-corner_mapping.R` to reproduce the measurement. +#' @noRd +fly_digital_rotation <- function() 270L + +#' How far a corner mapping may stretch an image before it is refused +#' +#' Set between the largest stretch a legitimate frame produces and the smallest a wrong +#' corner mapping can. A full-resolution 9-inch film scan carrying the negative's rebate +#' runs about 6.7% off square; the tightest case that must be caught is a Leica DMC II +#' frame sized through `format_size` onto a square footprint, at 9.95%. The admissible +#' band is therefore (1.0667, 1.0995) and this sits near the middle of it. Checked +#' against every row of `camera_formats.csv`. See `georef_one()` for the full table. +#' @noRd +fly_gcp_stretch_max <- function() 1.08 + +#' How far a corner mapping stretches an image +#' +#' Metres per pixel along the image's width axis, divided by metres per pixel along its +#' height axis. An isotropic mapping gives 1; a mapping that pairs the image's long axis +#' with the footprint's short edge gives the footprint's aspect ratio squared. +#' +#' Computed from the correspondence itself rather than from the footprint, so it measures +#' what GDAL will actually be asked to do. +#' @param gcp A matrix from [fly_georef_gcps()]. +#' @param ncol_px,nrow_px Image dimensions in pixels. +#' @return A single numeric ratio, or `NA_real_` if either dimension is degenerate. +#' @noRd +fly_gcp_anisotropy <- function(gcp, ncol_px, nrow_px) { + if (!is.finite(ncol_px) || !is.finite(nrow_px) || ncol_px <= 0 || nrow_px <= 0) { + return(NA_real_) + } + g <- gcp[, c("ground_x", "ground_y"), drop = FALSE] + w <- sqrt(sum((g[2, ] - g[1, ])^2)) / ncol_px + h <- sqrt(sum((g[3, ] - g[2, ])^2)) / nrow_px + if (!is.finite(w) || !is.finite(h) || h == 0) return(NA_real_) + w / h +} diff --git a/README.md b/README.md index 74d6824..7e60409 100644 --- a/README.md +++ b/README.md @@ -23,7 +23,7 @@ fly estimates each photo's ground footprint from its scale and film format, then - **Explore** — summarize available photos by scale and year, and measure coverage and frame-to-frame overlap (`fly_summary`, `fly_coverage`, `fly_overlap`) - **Select** — keep only photos whose footprint covers your area of interest, then take the fewest needed (best resolution first) or every frame touching it (`fly_filter`, `fly_select`) - **Fetch** — download thumbnails, flight logs, and camera calibration reports from the BC Data Catalogue, optionally in parallel (`fly_fetch`) -- **Georeference** — warp the scanned images onto their estimated footprints as GeoTIFFs (BC Albers), with automatic flight-line rotation (`fly_bearing`, `fly_georef`) +- **Georeference** — warp the scanned images onto their estimated footprints as GeoTIFFs (BC Albers), with automatic flight-line rotation, for film and digital frames alike (`fly_bearing`, `fly_georef`) ## Installation diff --git a/data-raw/georef_calibrate-corner_mapping.R b/data-raw/georef_calibrate-corner_mapping.R index 4952245..f1a4f78 100644 --- a/data-raw/georef_calibrate-corner_mapping.R +++ b/data-raw/georef_calibrate-corner_mapping.R @@ -45,7 +45,8 @@ heading_per_frame <- function(roll, frame, east, north) { i <- which(roll == r) i <- i[order(frame[i])] if (length(i) < 2) next - de <- diff(east[i]); dn <- diff(north[i]) + de <- diff(east[i]) + dn <- diff(north[i]) step <- sqrt(de^2 + dn^2) a <- az(de, dn) a[step < 100 | step > 5000] <- NA # turns and roll breaks @@ -59,10 +60,10 @@ report_offset <- function(label, image_x_az, heading) { d <- wrap(image_x_az[ok] - heading[ok]) s <- (image_x_az[ok] + heading[ok]) %% 360 tight <- function(v) 100 * mean(abs(wrap(v - median(v))) <= 5) - message(sprintf( - "%s n=%d compass bins=%d\n rigid median %8.2f within +/-5: %5.1f%%\n reflected median %8.2f within +/-5: %5.1f%%", - label, sum(ok), length(unique(round(heading[ok] / 30) * 30 %% 360)), - median(d), tight(d), median(s), tight(s))) + message(sprintf("%s n=%d compass bins=%d", label, sum(ok), + length(unique(round(heading[ok] / 30) * 30 %% 360)))) + message(sprintf(" rigid median %8.2f within +/-5: %5.1f%%", median(d), tight(d))) + message(sprintf(" reflected median %8.2f within +/-5: %5.1f%%", median(s), tight(s))) } read_patb <- function(url) { @@ -138,7 +139,8 @@ grey <- function(path) { } pair_r <- function(a, b, res = 25) { - ga <- grey(a); gb <- grey(b) + ga <- grey(a) + gb <- grey(b) inter <- terra::intersect(ext(ga), ext(gb)) if (is.null(inter)) return(NA_real_) tmpl <- rast(inter, resolution = res, crs = crs(ga)) @@ -169,10 +171,11 @@ for (case in list(list(i = 19:24, tag = "eagle", cam = "UltraCam Eagle M3"), # than assumed. bb <- st_bbox(st_union(st_geometry(fp))) +lakes_query <- bcdata::bcdc_query_geodata("WHSE_BASEMAPPING.FWA_LAKES_POLY") lakes <- st_transform( - bcdata::collect(bcdata::filter( - bcdata::bcdc_query_geodata("WHSE_BASEMAPPING.FWA_LAKES_POLY"), - BBOX(!!as.numeric(bb), crs = "EPSG:3005"))), 3005) + bcdata::collect(bcdata::filter(lakes_query, BBOX(!!as.numeric(bb), crs = "EPSG:3005"))), + 3005 +) water <- st_union(lakes) for (i in c(1, 2)) { diff --git a/inst/notes/georeferencing.md b/inst/notes/georeferencing.md index 1663986..d632235 100644 --- a/inst/notes/georeferencing.md +++ b/inst/notes/georeferencing.md @@ -110,6 +110,35 @@ Mean luminance inside FWA lake polygons minus the frame mean: - **The overlap and water checks** both need real thumbnails, so they live in `data-raw/`, not in the test suite. +## The stretch tolerance, and why it is not a shape gate + +`georef_one()` refuses a mapping whose image and footprint disagree about the frame's +shape by more than 10%. The number is set from measurements, not picked: + +| \|log\| off isotropic | what it is | +|---|---| +| 0.0000 | the bundled film thumbnails, 1250 x 1250 | +| 0.0645 | a full-resolution 9-inch scan at 9600 x 9000, carrying the negative's rebate | +| **0.0770** | **the tolerance, 1.08** | +| 0.0949 | the tightest case that must be caught — a Leica DMC II frame sized through `format_size` onto a square footprint | +| 0.1898 | the tightest mispairing on a frame's own footprint, the same DMC II squared | +| 0.4421 | a portrait UltraCam frame on a square footprint | + +The admissible band is **(1.0667, 1.0995)** and it is narrow, because a square-footprint +DMC II frame is only slightly more eccentric than a badly rebated film scan. A first +draft set 1.10, which is outside it by 0.4% and let exactly one row of +`camera_formats.csv` through — the DMC II, the one case the tolerance had replaced a +shape gate to cover. The test that was supposed to guard the threshold asserted it +against the *UltraCam* at 0.442, which any tolerance clears. **Check a threshold against +the least favourable member of the population, computed, not against a remembered +example.** + +An earlier version exempted square footprints instead, on the reasoning that a square one +has no pairing to get wrong. It is true and it is the wrong fix: `format_size` sizes a +frame from a single width, so a digital frame from a camera `fly` does not know lands on a +**square** footprint — and that is precisely the case the guard exists for. A shape gate +switches it off there. A tolerance wide enough for the rebate does not. + ## The one thing that would invalidate this All three measurements are downstream of `camera_formats.csv` assigning the long sensor diff --git a/man/fly_georef.Rd b/man/fly_georef.Rd index a12010e..d092633 100644 --- a/man/fly_georef.Rd +++ b/man/fly_georef.Rd @@ -35,11 +35,16 @@ for thumbnails but may need adjustment for full-resolution scans. Set to \code{NULL} to disable source nodata detection entirely.} \item{rotation}{Image rotation in degrees clockwise. One of \code{"auto"}, -\code{0}, \code{90}, \code{180}, or \code{270}. \code{"auto"} (default) computes flight line -bearing from consecutive centroids and derives rotation per-photo — -requires \code{film_roll} and \code{frame_number} columns. Fixed values apply -the same rotation to all photos. Overridden per-photo if \code{photos_sf} -contains a \code{rotation} column.} +\code{0}, \code{90}, \code{180}, or \code{270}. Applies to frames with a \strong{square} +footprint only. \code{"auto"} (default) computes flight line bearing from +consecutive centroids and derives rotation per-photo — requires +\code{film_roll} and \code{frame_number} columns. Fixed values apply the same +rotation to every square-footprint photo. A non-square footprint +ignores this argument and uses the measured digital mapping (see +\strong{Rotation}); a \code{rotation} column in \code{photos_sf} still overrides +per-photo, for both. Carrying a film-era \code{rotation} column into a +batch of digital frames therefore overrides the correct mapping with +the wrong one — drop the column, or set it to \code{NA} for those rows.} \item{dem}{Optional elevation raster passed to \code{\link[=fly_footprint]{fly_footprint()}}, sizing each frame from its height above ground instead of the reported scale. See the @@ -59,14 +64,15 @@ polygon corners computed by \code{\link[=fly_footprint]{fly_footprint()}} in BC translates the image with GCPs then warps to the target CRS using bilinear resampling. -\strong{Rotation:} Aerial photos may appear rotated in their footprints -because the camera orientation relative to north varies by flight -direction, camera mounting, and scanner orientation. The \code{rotation} -parameter rotates the GCP corner mapping: +\strong{Rotation:} the corner mapping depends on the footprint's shape, because +\code{\link[=fly_footprint]{fly_footprint()}} builds the two shapes differently. + +A \strong{square} footprint is axis-aligned, so the mapping carries the rotation. +The \code{rotation} parameter rotates it: \itemize{ -\item \code{0} — top of image maps to north edge of footprint (original behavior) +\item \code{0} — top of image maps to north edge of footprint \item \code{90} — top of image maps to east edge (90° clockwise) -\item \code{180} — top of image maps to south edge (default, correct for most BC photos) +\item \code{180} — top of image maps to south edge (correct for most BC film) \item \code{270} — top of image maps to west edge } @@ -84,10 +90,41 @@ not per-roll. To override, add a \code{rotation} column to \code{photos_sf}: \if{html}{\out{
}}\preformatted{photos$rotation <- dplyr::case_when( photos$film_roll == "bc5282" ~ 270, - .default = NA # fall through to auto + .default = NA # non-square footprints fall through to the digital mapping; + # square ones to `rotation`, which is 180 under "auto" — + # NOT back to the per-photo bearing ) }\if{html}{\out{
}} +Every non-\code{NA} value in that column must be 0, 90, 180 or 270; anything else +is an error naming the value, rather than a rotation silently applied as +something other than what was written. + +A \strong{non-square} footprint — every digital frame — is already rotated onto +its flight line by \code{\link[=fly_footprint]{fly_footprint()}}, so the ring carries the bearing and +applying it again would count it twice. Those frames use a fixed mapping +instead: the top-left pixel maps to the ring's rear-left corner, equivalently +image columns run in the flight direction and image rows run flight-right. +That was measured in fly#38 on both bundled cameras by three independent +routes and is the same for each; see \code{inst/notes/georeferencing.md}. + +A frame whose delivered image aspect disagrees with its footprint's by more +than 8\% is \strong{skipped with a warning} rather than written stretched — the +failure a wrong mapping produces is a valid GeoTIFF over the right ground, +squashed by the aspect ratio squared, which nothing downstream would report. +The threshold sits between the largest disagreement a legitimate frame +produces (a full-resolution 9-inch scan carrying the negative's rebate, about +6.7\%) and the smallest that must be caught (a Leica DMC II frame sized through +\code{format_size} onto a square footprint, 9.95\%). It applies to square footprints +too: a square one has no pairing to get wrong, but a digital frame sized +through \code{format_size} lands on one, and that is the unknown-camera case — so +gating on shape would switch the check off exactly where it is needed. + +A non-square footprint built without a flight bearing is drawn axis-aligned +and is therefore georeferenced as though the flight line ran due north. +\code{\link[=fly_bearing]{fly_bearing()}} needs a neighbouring frame, so this is the ordinary result of +georeferencing a single frame on its own, and it is warned about. + \strong{Nodata handling:} Two sources of unwanted black pixels are masked: \enumerate{ \item \strong{Warp fill} — GDAL creates black pixels outside the rotated source diff --git a/planning/active/progress.md b/planning/active/progress.md index 87acd9c..78f7197 100644 --- a/planning/active/progress.md +++ b/planning/active/progress.md @@ -80,3 +80,152 @@ reader will find most convincing. Both cameras give the same answer, so the constant is global rather than a per-camera column. + +### Phase 4 — exclusion removed, digital frames georeference + +- Non-square footprints take `fly_digital_rotation()` (270) and never + `bearing_to_rotation()`. Classification is three-way — empty / square / non-square — + computed before the loop, because `fly_is_square()` reports an EMPTY geometry as + square and squareness alone therefore cannot stand in for "has ground control" +- Isotropy guard in `georef_one()`: a mapping that would stretch the image by more than + 5% is refused with a warning rather than written. Real data lands at 0.03% off + isotropic and the least eccentric wrong pairing is 21% off, so the tolerance has + room in both directions +- `user_rotation_col` captured **before** the auto path overwrites `has_rotation_col`. + Without it nothing downstream could tell a user column from a bearing-derived one, so + the documented override would have silently stopped working for digital frames +- Warning added for a non-square footprint with no bearing — the ordinary result of + georeferencing one frame on its own, since `fly_bearing()` needs a neighbour +- Live check: 6 Eagle frames georeference 6 of 6 in EPSG:3005; a mixed batch of 3 film + and 3 digital gives film the bearing rule and digital the constant, 6 of 6, no warning +- New tests observe the rotation at the boundary rather than inferring it from output + +**A leaking mock, found by the guard tests failing.** `local_mocked_bindings(.env = )` +is the environment the mock unwinds with, not the target — `.package` names the target. +Passing `asNamespace("fly")` installed the stub correctly and then never removed it, +because a namespace does not exit, so every later test in the run kept it. It leaked in +the direction that reads as success: a stub returning `TRUE` makes assertions pass. It +surfaced only because a later test asserted a file existed that the stub never wrote. +`testthat` pin bumped to >= 3.2.0 for `.package`. Recorded in CLAUDE.md. + +### Phase 5 — notes, docs, release + +- `inst/notes/georeferencing.md` written in Phase 3 (the measurement record) +- README, vignette georef section, `@details` **Rotation** and `@param rotation` all + rewritten — the old text documented film-only behaviour as if it were general +- `devtools::document()`: `fly_georef.Rd` only, exports steady at 9, no rebind +- NEWS entry; version 0.6.0 -> 0.7.0 + +### Code-check round 1 — one real regression, four fragilities + +Findings in `planning/active/review-round1.md`. + +**Fixed — the regression I would have shipped.** The isotropy guard ran on film too, and +on a *square* footprint the anisotropy is the image's own inverse aspect and is identical +at all four rotations — so it stopped being a corner-mapping check and became "is this +scan square to within 5%", which no rotation can fix. Reproduced before fixing: a +1250x1250 thumbnail passes, a 9600x9000 full-resolution 9-inch scan is **refused**. +`fly_georef()` documents full-resolution scans as supported, and every bundled film +thumbnail is exactly 1250x1250, so the fixture set was structurally incapable of reaching +it. Guard now gated on `!fly_is_square(fp)`, with a regression test on a 1250x1200 scan +asserting no warning at any of the four rotations. + +**Fixed — the `rotation` column was never validated**, and this branch made it the +highest-precedence input for digital frames while the docs tell users to manage it. 360 +indexed past the ring and surfaced from inside `tryCatch` as "subscript out of bounds", +naming neither column nor value; 45 and -90 shifted by zero and georeferenced silently +wrong. + +**And a second bug inside that fix.** Validating with `as.integer(as.character(x))` while +the per-row read still used `as.integer(x)` means a factor column validates as 180 and is +*applied* as its level code, 1. Converting in two places is exactly how those come apart +(CLAUDE.md, cross-function normalisation). Normalised once into `user_rot` and read from +that one place; the test asserts the factor is honoured as 180 and explicitly that it is +not 1. + +**Fixed — two tests.** One named "warned about once" never asserted a warning, only a +frame count, so it passed for zero warnings. `expect_warning()` was the wrong instrument +because the fixture raises a second, unrelated warning that testthat then re-raises; +counting under `withCallingHandlers` makes the assertion exact. The other carried a +comment saying a square footprint is accepted directly above an assertion that it is +refused — the shape that produces a wrong "fix". + +**Accepted.** Duplicate/absent `airp_id` resolving silently is pre-existing and unrelated +to digital frames; the no-bearing warning being a `grepl()` on a marker `fly_footprint()` +appends is stringly-typed but currently coupled by an end-to-end test. + +### Code-check round 2 — three more, two of them inside round 1's fixes + +Findings in `planning/active/review-round2.md`. The pattern held: the fixes were the +prime suspects, and that is where two of the three were. + +**The factor defect round 1 named was still live.** The round-1 fix normalised the column +into `user_rot` and its own comment claimed it was "read from this vector everywhere +below" — there were three read sites and only two were converted. A factor level that +does *not* parse gives `NA`, which passes validation and skips the user branch, so control +reached the third site and got the **level code**. Measured: +`factor(c("180","north",...))` gave rotations `180, 2, 180, 2`, and `1 %/% 90` is 0, so +those frames georeferenced at rotation 0 with `success = TRUE`. The new test used +`factor("180")`, whose level parses, so it never reached that branch. Fixed by writing +the parsed vector back into `photos_sf` so there is exactly one parse, and by refusing a +value that was supplied but did not parse instead of downgrading it to NA. + +**The shape gate was the wrong fix for round 1's regression.** Exempting square footprints +switches the guard off for a digital frame sized through `format_size`, which produces a +*square* footprint from a single width — the unknown-camera case the guard exists for. +Measured: a portrait 1063x1654 image on that square footprint wrote a **1.556x stretch +silently, 0 warnings**. Replaced with a single tolerance of 10%, set from measurements +rather than picked: + +| \|log\| off isotropic | what | +|---|---| +| 0.065 | full-resolution 9-inch scan with rebate — the worst legitimate case | +| **0.095** | **the tolerance** | +| 0.190 | tightest mispairing the shipped table can produce (DMC II) | +| 0.442 | portrait digital frame on a square footprint | + +Roughly 1.5x headroom below and 2x above. Pinned as its own test so moving it has to +disagree with numbers rather than with a comment. + +**And my regression test for round 1 could not fail.** Its 1250x1200 fixture is 0.041 off +isotropic, inside even the old 5% tolerance, so it passed with the gate stripped — +verified by restoring the bug. The comment beside it asserted "outside the guard's 5% +tolerance", which was false and was the premise a reader would trust. Replaced with +1250x1172, the same aspect as the 9600x9000 scan that motivated it, at 0.064 — the +closest legitimate case to the threshold. + +**Fixed in passing:** the `@details` `case_when` example claimed `.default = NA` falls +through to auto. It does not — a square footprint falls to 180 and a non-square one to +the digital mapping. + +### Code-check round 3 — the tolerance was 0.4% too loose + +Findings in `planning/active/review-round3.md`. Round 2's fixes were the target, and one +of them was wrong by a hair in the direction that matters. + +`fly_gcp_stretch_max()` was set to **1.10**. The tightest case it has to catch is a Leica +DMC II frame sized through `format_size` onto a square footprint: `|log(15552/14144)| = +0.09490` against `log(1.10) = 0.09531`. It slipped — verified end to end on a real +bundled frame through `fly_footprint()`, `georef_one()` and GDAL: file written, +`success = TRUE`, **0 warnings**, a 10% stretch. Of all 19 rows in `camera_formats.csv` +it is the only one that slips, and it is the exact case the tolerance had replaced the +shape gate to cover. + +The admissible band is **(1.0667, 1.0995)** — narrow, because a square-footprint DMC II +frame is barely more eccentric than a badly rebated film scan. Set to **1.08**, near the +middle, with 1.19x margin below and 1.23x above. + +**And the test that was supposed to guard the threshold could not.** It asserted the +square-footprint case using the UltraCam at 0.442 — the *most* eccentric camera, which +any tolerance clears. Picking the lenient example is how the threshold came to be wrong +in the first place. Rewritten to compute over every row of the shipped table and assert +that none slips, plus the margin on both sides. Restore-the-bug: putting 1.10 back takes +it from FAIL=0 to FAIL=2, both bindings patched with the value printed. + +Also fixed: the `case_when` roxygen comment held only under `rotation = "auto"` — under an +explicit argument a square footprint falls through to that argument, not to 180. + +Lint: `data-raw/georef_calibrate-corner_mapping.R` cleaned from 6 style lints to 0. +`R/fly_georef.R` is 1 against a `main` baseline of 5, and that one is the documented +installed-vs-source artifact (`exists("fly_gcp_stretch_max", asNamespace("fly"))` is +FALSE). diff --git a/planning/active/review-round1.md b/planning/active/review-round1.md new file mode 100644 index 0000000..c037db2 --- /dev/null +++ b/planning/active/review-round1.md @@ -0,0 +1,244 @@ +# Code check — round 1 — fly#38 staged diff + +Reviewed: `git diff --cached` (11 files). Verified empirically, not by reading only: + +- full suite run under `pkgload::load_all()` + `test_dir()`: **green**, 0 failures, + 0 `Error`/`Execution halted` markers (`/tmp/suite.log`). +- `devtools::document()` re-run: **no drift** — `man/fly_georef.Rd` and `NAMESPACE` + already match the roxygen in the diff. `fly_digital_rotation()` / + `fly_gcp_anisotropy()` are `@noRd` and correctly absent from `NAMESPACE`. +- `inst/notes/georeferencing.md` and `data-raw/georef_calibrate-corner_mapping.R` + both exist and are in `HEAD`, so the new doc/warning references resolve. +- `testthat (>= 3.2.0)` bump in DESCRIPTION is correct and required — `.package =` + on `local_mocked_bindings()` needs it. +- `fly_bearing()` does **not** reorder rows (it writes `bearing[ord[i]]` and returns + the input object), so the `photos_sf` / `footprints` / `empty_fp` / `non_square` + index alignment asked about in the brief is sound. +- `grepl(pat, NA_character_)` is `FALSE`, not `NA`, so the `width_source` block cannot + propagate `NA` into `any()`. +- The three-way `empty_fp` / `non_square` classification genuinely removes the ordering + hazard: `if (empty_fp[j]) next` runs before `footprints[j, ]` is ever subset, so + `georef_one()` is never handed an empty ring, and `fly_is_square()`'s TRUE-on-EMPTY + is masked out of `non_square` by construction. +- Rotation-selection trace over every combination of (`rotation` arg, column present, + column `NA`, `film_roll`/`frame_number` present): the **rotation choice** for film is + unchanged from the pre-diff code in all of them. See Finding 1 for the one place film + behaviour *is* changed. + +## Findings + +### 1. **[bug]** `R/fly_georef.R:297-304` — the anisotropy guard silently drops non-square **film** scans, a documented use case + +The guard is written for the digital mispairing case, but it runs on every frame, +including square film footprints. On a square footprint the anisotropy is +**independent of rotation** — it reduces to the image's own inverse aspect ratio: + +``` +square 4000x4000 m ring, image 1250x1180 px + rotation 0 90 180 270 + aniso 0.944 0.944 0.944 0.944 <- identical; nothing to mispair +``` + +(The package's own `test-fly_georef_aspect.R:110` asserts exactly this — "square film is +isotropic at every rotation, so the invariant is vacuous there".) + +So for film the guard is not a corner-mapping check at all; it is a bare +"is this image square to within 5%" check, and a fail is unfixable by the user because +no rotation moves the number. Measured: + +``` +film scan 1250x1250 -> aniso 1.0000 refused = FALSE +film scan 1250x1200 -> aniso 0.9600 refused = FALSE +film scan 1250x1180 -> aniso 0.9440 refused = TRUE +film scan 9600x9000 -> aniso 0.9375 refused = TRUE <- ordinary full-res 9" scan +``` + +`fly_georef()` is documented as working with "full-resolution scans" +(`R/fly_georef.R:5`, and the `srcnodata` tradeoff paragraph at :101-105 exists +specifically for them). A full-resolution scan that includes any of the rebate/fiducial +border is routinely a few percent off square. Those frames previously georeferenced +(stretched by that same few percent, far inside the footprint's own estimation error); +they now return `success = FALSE` with no output file. + +The warning text compounds it — it says *"the corner mapping would stretch it by +0.937x. Skipped rather than written squashed"*, which sends the reader to look for a +rotation fix that does not exist for a square footprint. + +**Why the suite cannot see it:** every bundled film thumbnail is exactly 1250x1250 +(verified by fetching six of them and reading `gdalinfo`), so the fixture set is +structurally incapable of reaching this — the "fixture that cannot reach the failure +mode" class in CLAUDE.md. + +**Suggested shape of a fix** (the guard should measure *mispairing*, not *aspect +mismatch*): compare `aniso` against the best anisotropy achievable over the four +rotations, and refuse only when another pairing would be materially better. Equivalently, +apply the guard only where `non_square[j]` — a square footprint has no pairing to get +wrong. Either way, the guard then still catches the 2.42x UltraCam case and the +inferred-format mismatch it was written for, and stops refusing film it should accept. + +--- + +### 2. **[fragile]** `R/fly_georef.R:232-236` — the `rotation` *column* is never validated, and this diff makes it load-bearing where it used to be inert + +The `rotation` **argument** is validated at :133-138 to be one of 0/90/180/270. The +`rotation` **column** is not; it goes straight into `as.integer()` and then into +`fly_georef_gcps()`'s `n_shifts <- rotation %/% 90`. Measured: + +``` +rot -90 -> silently treated as rotation 0 (n_shifts <= 0, no shift applied) +rot 45 -> silently treated as rotation 0 +rot 360 -> ERROR: subscript out of bounds (c(5:4, 1:4) -> 6 ground rows vs 4 pixel rows) +rot 450 -> ERROR: subscript out of bounds +``` + +The error is swallowed by the `tryCatch` at :252-258, so a whole batch fails with +`Failed to georef : subscript out of bounds` — a message that names neither the +column nor the offending value. A factor column is worse: `as.integer(factor("180"))` is +`1`, which is silently rotation 0. + +This is pre-existing for film, but the diff changes its weight in two ways: (a) the new +`user_rotation_col` branch makes the column override the *measured* digital mapping, +and (b) the new docs actively instruct users to manage this column for digital batches +("drop the column, or set it to `NA` for those rows", :32-34). A column that was inert +for digital frames is now the highest-precedence input for them. + +One line beside the existing argument check would close it — validate the non-`NA` +column values against `c(0, 90, 180, 270)` and stop naming the bad value. + +--- + +### 3. **[fragile]** `R/fly_georef.R:218-220` — duplicate or absent `airp_id` resolves silently + +`fp_idx <- which(photos_sf[["airp_id"]] == results$airp_id[i])`, then `j <- fp_idx[1]`. + +- **Duplicate `airp_id`** (the same frame arriving twice from two overlapping catalogue + queries `rbind`-ed together): every duplicate row is georeferenced onto the *first* + row's footprint, rotation and squareness classification. The output path is + `basename(src)`, so the second write also lands on the same file. No warning. +- **Absent `airp_id` column**: `photos_sf[["airp_id"]]` is `NULL`, `NULL == x` is + `logical(0)`, `which()` is `integer(0)`, so *every* row `next`s and the function + reports `Georeferenced 0 of N images` as though the images were simply unfetchable. +- In both skip paths `results$dest[i]` has already been assigned at :210, so the returned + tibble names an output file that was never written (`success = FALSE` alongside it, so + a caller reading `success` is fine; one reading `dest` is not). + +All pre-existing, but the diff adds two more per-row lookups keyed on `j` +(`empty_fp[j]`, `non_square[j]`) that inherit it, and the brief asked. A +`if (length(fp_idx) > 1) warning(...)` and an explicit `airp_id` presence check would +make both states visible. + +--- + +### 4. **[fragile]** `tests/testthat/test-fly_georef_digital.R:104-114` — test asserts less than its name claims + +```r +test_that("frames with no footprint are still skipped, and warned about once", { + ... + seen <- capture_rotations( + suppressWarnings(fly_georef(fake_fetch(photos), photos, dest_dir = tempfile())) + ) + expect_length(seen, sum(!unsized)) +}) +``` + +The body checks only the *count* of frames reaching `georef_one()`. The +"warned about once" half is inside `suppressWarnings()` and is never asserted — an +implementation that warned zero times, or once per frame, passes this test unchanged. +The premise line (`expect_true(any(unsized))`) is good and does its job; the missing +piece is an `expect_warning(..., "no footprint")` around the call, or a +`withCallingHandlers` count. + +--- + +### 5. **[fragile]** `tests/testthat/test-fly_georef_digital.R:159-162` — comment contradicts the assertion directly beneath it + +```r +# And the guard is about the pairing, not about rotation 0 as such: a square footprint +# takes rotation 0 happily, because there is nothing to mispair. +sq <- ring(hc = 1000, ha = 1000) +out3 <- tempfile(fileext = ".tif") +expect_warning(georef_one(src, sq, out3, rotation = 0), "Skipped rather than") +``` + +The comment says the square footprint is accepted; the assertion asserts it is +**refused**. The *next* comment (:164-165) explains why the refusal is correct, but a +reader or a future editor hitting the first comment reads the assertion as the bug. +CLAUDE.md's own rule applies here — "read the guard, not the comment above it" — and this +is the shape that produces a wrong "fix". + +Minor companion: `out3` is never checked for non-existence the way `out2` is at :157, and +the return value of that call is discarded (`expect_false(res)` is only done for `out2`). + +--- + +### 6. **[fragile]** `R/fly_georef.R:158-169` — the no-bearing warning is a string match on another function's output, and every failure direction is silence + +`grepl("axis_aligned_no_bearing", footprints$width_source)` is a stringly-typed contract +with `fly_footprint()`, which builds that marker by `paste0()`-appending to +`width_source` (`R/fly_footprint.R:543-545`). Three ways it goes quiet without an error: +the marker text changing on the producing side; `width_source` being absent (the +`%in% names()` guard at :158 skips the whole block); `width_source` being `NA` for the +row (`grepl` returns `FALSE`). All three fail toward "no warning", the direction that +reads as success. + +It is currently correct and the new end-to-end test at :117-130 does couple the two +sides, which is what keeps this at *fragile* rather than *bug*. Worth noting because the +value being matched is produced two modules away and nothing pins its text. + +--- + +## Checked and found sound + +- **`user_rotation_col` / `has_rotation_col` split** — correct. `user_rotation_col` is + captured before the auto path can set `photos_sf$rotation`, so a bearing-derived + rotation can never masquerade as a user override. Every combination traced: + | rotation arg | column | column value | result | vs. pre-diff | + |---|---|---|---|---| + | `"auto"` | absent, roll+frame present | — | bearing-derived per row | unchanged | + | `"auto"` | absent, no roll/frame | — | message + 180 | unchanged | + | `"auto"` | present | value | value | unchanged | + | `"auto"` | present | `NA` | 180 (film) / 270 (digital) | film unchanged | + | fixed `r` | absent | — | `r` (film) / 270 (digital) | film unchanged | + | fixed `r` | present | `NA` | `r` (film) / 270 (digital) | film unchanged | +- **`fly_gcp_anisotropy()` degenerate inputs** — no `log()` of a negative is reachable + (both terms are `sqrt()`); `w == 0` gives `aniso == 0`, `log(0)` is `-Inf`, `abs(-Inf)` + exceeds the threshold, so it refuses rather than errors; `NA` short-circuits on + `!is.finite(aniso)` before `log()` is evaluated, so no `NA` reaches `if`. + `ncol_px = 0` returns `NA_real_` as documented. `warning()` inside the loop is not + caught by the `tryCatch` (which handles `error` only), so it reaches the caller. +- **Guard headroom** — real digital thumbnails land at anisotropy 0.99974 (0.03% off) + against a 5% tolerance and a 21% minimum for the least-eccentric wrong pairing. + DEM-sized digital frames preserve the sensor aspect exactly (`resize()` scales both + half-dimensions by the same `k`), so the DEM route cannot drift into the guard. +- **Input class shapes** — ran `fly_georef()` end to end over plain / tibble / grouped / + `bcdc_sf` digital centroids (the `centroid_shapes()` sweep): all four give rotation + 270 for all 6 frames and `success` 6/6. The `width_source` `$` access and + `photos_sf[["rotation"]][j]` indexing are safe on every shape. +- **Discriminating power of the new tests** — checked each is capable of failing: + digital bearings are ~343 deg, which `bearing_to_rotation()` maps to **0**, not 270, so + the digital assertions discriminate. In the mixed-batch test the four film frames map + to **90, 270, 270, 90** — two coincide with the digital constant, but the two 90s make + the assertion fail against an implementation that applied 270 everywhere, so the test + is not vacuous. The `wrong <- 2.4217` pin and the reversed-ring handedness check both + fail on a broken implementation. +- **`local_mocked_bindings()` scoping** — `.env = parent.frame()` unwinds with the + calling test rather than the helper, and mocks `fly`'s namespace (the correct target + for an internal `fly_georef()` -> `georef_one()` call path). No leak: the last test in + the same file calls the **real** `georef_one()` and writes actual GeoTIFFs, and it + passes. The `test-fly_georef*.R` files that run after it are unaffected. +- **`fly_is_square()` TRUE-on-EMPTY** — neutralised, see the classification note above. + `fly_is_square()` also returns early on empty before `st_coordinates()`, so it cannot + error on an unsized frame. +- No shipped camera in `inst/extdata/camera_formats.csv` has a square sensor (minimum + aspect 1.0995, DMC II), so the "square digital sensor falls through to the film rule" + path is not currently reachable. + +## Not raised (deliberately) + +- The orphaned `# unname() matters` comment at `R/fly_georef.R:288-291` now sits above + the anisotropy block rather than above the `gcp_args` loop it describes — readability + only. +- `DESCRIPTION` `Date: 2026-08-29` while `NEWS.md` dates 0.7.0 as 2026-08-30. +- The warning and `.Rd` say `inst/notes/georeferencing.md`; installed, the path is + `notes/georeferencing.md`. Consistent with existing house style elsewhere in the + package. diff --git a/planning/active/review-round2.md b/planning/active/review-round2.md new file mode 100644 index 0000000..d39ffeb --- /dev/null +++ b/planning/active/review-round2.md @@ -0,0 +1,286 @@ +# Code check — round 2 — fly#38 staged diff + +Scope: `git diff --cached`. Round 2's brief was to treat round 1's **fixes** as the prime +suspects. Two of the three fixes carry a defect, and one of them is the same class round 1 +named by name. + +Verified empirically, not by reading: + +- full suite under `pkgload::load_all()` + `test_dir()`: **green — `FAIL 0 | WARN 0 | + SKIP 0 | PASS 1244`**, 0 `Error` / `Execution halted` markers. Every finding below is + invisible to it. +- `/tmp/cc_checklist.md` **is zero bytes** — the checklist named in the brief was empty, + so this review worked from `CLAUDE.md`'s code-check conventions instead. Flagging it so + nobody records "the checklist was worked through". +- Finding 2 was established by the restore-the-bug procedure: the gate was stripped from a + copy of `georef_one()`, the patch was proven to have taken (`fly_is_square` absent from + the deparsed body), and the test re-run. + +--- + +## Findings + +### 1. **[bug]** `R/fly_georef.R:258` — the factor defect round 1 fixed is still live in the third branch + +Round 2 normalised the `rotation` column read into `user_rot` "once, here, and read from +this vector everywhere below" (`:188-191`). It is **not** read from that vector everywhere +below. There are three read sites and only two were converted: + +| line | site | conversion | +|---|---|---| +| `:194` | validation | `as.integer(as.character(x))` ✅ | +| `:252` | user-override branch | `user_rot[j]` ✅ | +| **`:258`** | **`has_rotation_col` fallback branch** | **`as.integer(photos_sf[["rotation"]][j])`** ❌ raw | + +The escape route is precise: a factor level that does **not** parse makes `user_rot[j]` +`NA`, which (a) passes validation, because `bad <- !is.na(user_rot) & …` excludes `NA`, and +(b) skips the `:252` user branch, because `is.na(user_val)`. Control then reaches `:258`, +which reads the raw factor and gets the **level code**. + +Measured on four bundled film frames: + +``` +column factor(c("180","north","180","north")) levels: 180, north +validation as.integer(as.character(x)) -> 180, NA, 180, NA (passes) +per-row :258 as.integer(x) -> 1, 2, 1, 2 (LEVEL CODES) +rotation applied -> 180, 2, 180, 2 +n_shifts = rot %/% 90 -> 2, 0, 2, 0 +bearing/auto would have given -> 90, 270, 270, 90 +``` + +``` +column factor("north") -> rotation applied: 1, 1, 1, 1 (0 shifts -> rotation 0) +``` + +Consequences, all silent, all with `success = TRUE` and a GeoTIFF on disk: + +- Rotations `1` and `2` **bypass the new validation entirely** — they are synthesised + inside the loop, so "must be NA or one of 0, 90, 180, 270" never sees them. +- `1 %/% 90` and `2 %/% 90` are both `0`, so the frame is georeferenced at **rotation 0**. + The docs at `:51` state 180 is "correct for most BC film", so the image lands 180° out. +- A factor with ≥ 361 levels would give `n_shifts >= 4` and the `subscript out of bounds` + the round-2 fix was written to prevent — swallowed by the `tryCatch` at `:270`, naming + neither the column nor the value. Exotic, but it is the same door. + +**Why the new test misses it:** `test-fly_georef_digital.R:229` uses `factor("180")`, whose +level *parses*, so `user_rot` is non-`NA` and the `:252` branch handles it — `:258` is never +reached. The test's own comment says *"Converting in two places is how those come apart"*; +they are still converted in two places. + +**Fix:** `:258` should read `user_rot[j]` too, or the branch should be restructured so the +raw column is read exactly once. Note the branch is *also* the auto/bearing path (where +`photos_sf$rotation` was written by `bearing_to_rotation()` and is a clean integer), so the +two cases need separating rather than the read merely being swapped. + +--- + +### 2. **[bug]** `tests/testthat/test-fly_georef_digital.R:185-209` — the regression test for the round-1 fix does not exercise the gate, and its comment states a false number + +The test is named "a film scan that is not exactly square still georeferences" and exists to +pin the `!fly_is_square(fp)` exemption. It uses a **1250 x 1200** image on a 4000 x 4000 m +square ring. That anisotropy is **inside** the guard's tolerance: + +``` +aniso = (4000/1250) / (4000/1200) = 0.96 +abs(log(0.96)) = 0.040822 +log(1.05) = 0.048790 <- threshold +0.040822 > 0.048790 = FALSE -> the guard would NOT have refused it +``` + +So the guard never fires on this fixture, gate or no gate. Confirmed by restoring the bug — +`georef_one()` re-derived with the `if (!fly_is_square(fp))` gate stripped, patch proven to +have taken: + +``` +WITH gate (as staged) -> test would PASS +patch took? body has no fly_is_square: TRUE +WITHOUT gate (bug back) -> test would PASS +``` + +The comment at `:200` asserts the opposite as its premise: + +```r +# 0.96 — outside the guard's 5% tolerance, and identical at every rotation. +``` + +`0.96` is **inside** the tolerance. The assertion two lines down +(`expect_equal(fly_gcp_anisotropy(...), 0.96)`) is correct and passes; it is the +interpretation that is wrong, and it is the interpretation a future editor will trust when +deciding whether the exemption is still needed. + +This is `CLAUDE.md`'s "a fixture that cannot reach the failure mode is not validation", +with the aggravating factor that a stated measurement licenses the claim. + +Dimensions that would actually exercise it (same 4000 m square ring): + +``` + 1250x1250 aniso=1.0000 refuses=FALSE <- cannot exercise + 1250x1200 aniso=0.9600 refuses=FALSE <- cannot exercise (what the test uses) + 1250x1180 aniso=0.9440 refuses=TRUE <- exercises the gate + 9600x9000 aniso=0.9375 refuses=TRUE <- exercises the gate; round 1's own example +``` + +Round 1's finding was motivated by the 9600x9000 full-resolution scan. Using that, or +1250x1180, makes the test discriminating and costs nothing. + +*(The gate is genuinely exercised elsewhere — `:178-181`, the 2000 x 2000 m square ring at +anisotropy 2.0, does fail with the gate removed. So the fix is covered; it is this +specifically-named regression test that is decoration.)* + +--- + +### 3. **[bug]** `R/fly_georef.R:317` — the gate switches the guard off for a digital frame sized through the documented `format_size` escape hatch + +The gate reasons that a square footprint has "no pairing to get wrong". True of the *corner +mapping*. But `!fly_is_square(fp)` is a proxy, and it also stands for **"was sized as +film"** — which for a digital frame is precisely the mismatch the guard's own comment at +`:306-308` says it exists to catch ("a frame sized from an inferred format that does not +match the camera that actually took it"). + +`format_size` is the documented escape hatch "for a camera `fly` does not know" +(`R/fly_footprint.R` roxygen), it takes a single width, and a single width produces a +**square** footprint. So every unknown digital camera routed through the documented escape +hatch now bypasses the guard. Measured on bundled digital frame 19 with a portrait +1063 x 1654 image: + +``` +media: Digital - Colour +fly_footprint(dig, format_size = c("Digital - Colour" = 9)) +footprint square? TRUE +guard runs? NO (gated off) +written? TRUE, 0 warnings +anisotropy the mapping applies: 1.556 (guard refuses at >1.05 or <0.952) +``` + +A 1.556x stretch, written silently, with `success = TRUE` — the exact failure the code +comment describes as "a valid GeoTIFF, in the right CRS, over the right ground — just +squashed … which nothing downstream would report". Before the gate this frame was refused. + +Related documentation gap: the roxygen at `:81-84` still promises, without qualification, +that "A frame whose delivered image aspect does not pair with its footprint edges is +**skipped with a warning**". After the gate that holds only for non-square footprints. The +`@param rotation` text was updated for the square/non-square split; this paragraph was not. + +**Shape of a fix** (round 1 suggested it and it survives this finding): keep the guard +running on every frame, but compare `aniso` against the **best** anisotropy achievable over +the four rotations, refusing only when another pairing would be materially better. That +measures mispairing directly instead of using shape as a stand-in, so it accepts the +off-square film scan *and* still refuses the square-footprint digital frame above. + +--- + +### 4. **[fragile]** `R/fly_georef.R:194-200` — the new validation cannot refuse a value that does not parse + +`bad <- !is.na(user_rot) & !user_rot %in% c(0L, 90L, 180L, 270L)` excludes `NA`, and +`as.integer(as.character(x))` produces `NA` for anything unparseable. So a typo is +indistinguishable from a deliberate `NA` and is silently downgraded to the fallback: + +``` +rotation = c("ninety", NA, "180", NA) -> ACCEPTED, no error +rotations applied: 180, 180, 180, 180 +bearing/auto would have given: 90, 270, 270, 90 +``` + +Row 3's `"180"` was honoured; row 1's `"ninety"` was discarded without a word. The batch +mixes honoured and ignored values with nothing reporting which. The block is titled +"refused by name, not by GDAL" — an unparseable value is exactly the case that most needs +naming, and it is the one case that gets through. + +Cheap close: compute `parsed_na <- is.na(user_rot) & !is.na(photos_sf[["rotation"]])` and +fold it into `bad`. + +Note this also feeds finding 1: it is `user_rot[j]` being `NA` here that routes a factor to +the raw read at `:258`. + +--- + +### 5. **[fragile]** `R/fly_georef.R:66-71` — the documented `case_when` escape hatch does not do what its comment says, and this diff promotes it + +The roxygen example says: + +```r +photos$rotation <- dplyr::case_when( + photos$film_roll == "bc5282" ~ 270, + .default = NA # fall through to auto +) +``` + +`NA` does **not** fall through to auto for a film frame. Presence of the column sets +`has_rotation_col` at `:181`, which suppresses the whole auto/bearing block at `:204`, and +an `NA` row then lands on `if (auto_rotation) 180L`. Measured on the same four film frames: + +``` +rotation column all NA : 180, 180, 180, 180 +no rotation column : 90, 270, 270, 90 <- what "fall through to auto" promises +``` + +So a user who overrides one roll and leaves the rest `NA` silently loses bearing-derived +rotation on **every other roll in the batch** — 180° wrong on two of these four frames. + +The behaviour is pre-existing, but the diff makes it load-bearing and more prominent: the +new `@param rotation` text at `:31-34` instructs users to "set it to `NA` for those rows", +and the new digital test at `:79-84` pins `NA` → fall-through as correct **for digital** +(where it genuinely is, because `non_square[j]` catches it before the film branch). The doc +now describes one behaviour that is true for digital and false for film, in a paragraph that +says "for both". + +Either fix the fall-through (test `all(is.na(user_rot))` rather than column presence when +deciding whether to run the auto block) or correct the comment to say `NA` means 180 for +film. Whichever, the two halves should stop contradicting. + +--- + +## Checked and found sound + +- **Index alignment** — re-derived independently rather than trusting round 1. + `fly_footprint()` preserves row count and order (`nrow` 4 → 4, `airp_id` identical), so + `footprints`, `empty_fp` and `non_square` are all in `photos_sf` row order. + `j <- which(photos_sf[["airp_id"]] == …)[1]` indexes `photos_sf`, and `user_rot` is + built from `photos_sf[["rotation"]]` before any reassignment, so `user_rot[j]`, + `empty_fp[j]` and `non_square[j]` are all keyed on the same thing. The one reassignment + of `photos_sf` (`fly_bearing()` at `:206`) is unreachable when a user column exists, so + it cannot desynchronise `user_rot`. +- **Film behaviour bit-for-bit unchanged when no digital frame is present** — the + anisotropy guard is new in this diff, so a square footprint skipping it is exactly the + pre-diff path. Rotation selection for film re-traced and unchanged. The one behaviour + change is the new `stop()`, below. +- **`fly_is_square()` called per frame inside `georef_one()`** — correct and consistent + with `non_square[j]`: it is per-feature (`vapply` over `seq_along(g)`), uses + `all.equal(max(d), min(d))` with the default relative tolerance, and has no batch-shared + state, so `fly_is_square(footprints)[j]` and `fly_is_square(footprints[j, ])` cannot + disagree. No batch dependence. Cost is one `st_transform()` of a single already-3005 + feature per frame — negligible beside the two GDAL calls that follow. + `st_transform()` would error on a missing CRS, but `footprints` is transformed to 3005 at + `:143` before the loop, so that is unreachable from `fly_georef()`. +- **`fly_is_square()` TRUE-on-EMPTY inside `georef_one()`** — harmless. An empty ring would + die earlier at `sf::st_coordinates(fp)[1:4, ]`, and `fly_georef()` returns at + `if (empty_fp[j]) next` before `fp` is ever subset. +- **Fix 3, the `withCallingHandlers` warning counts** — these can fail, in both directions. + `test-fly_georef_digital.R:114-122` counts warnings matching `"have no footprint"` and + asserts `expect_identical(n, 1L)`; an implementation warning zero times fails, and one + warning per frame fails. The text matches `fly_warn_unsized()`'s + `" frames have no footprint and are excluded from "`. Muffling every warning is correct + here — the fixture also raises `fly_footprint()`'s unknown-format warning, which + `expect_warning()` would re-raise as a test WARNING. + `:131-139` (`expect_warning` "no flight bearing" on one frame, `expect_no_warning` on six) + discriminates a guard that never fires and one that always does. Suite runs `WARN 0`. +- **New `test-fly_georef_gcps.R` tests** — the handedness test is real (reversing the ring + flips the sign, asserted). The anisotropy test is not vacuous despite `wrong` being + hand-computed from the same 3175/2040/1063/1654 numbers `ultracam_ring()` encodes, + because `expect_equal(wrong, 2.4217, tolerance = 1e-4)` pins the value independently. + `aniso(90)` and `aniso(270)` both being 1 is acknowledged in the file header as + geometrically unavoidable. +- **`stop()` reachability** — the new validation is before the loop and before + `dir.create()`'s side effects matter, so it cannot leave a half-written batch. It *is* a + behaviour change: a film batch carrying `-90` or `45` in a `rotation` column previously + ran (silently at rotation 0) and now aborts the entire call rather than skipping the bad + rows. That is the right direction — fail loud — and is worth a NEWS line if there is not + one, but it is not a defect. +- No new hardcoded absolute paths, no secrets, no shell quoting, no file-writing side + effects outside `dest_dir` / `tempfile()`. + +## Not raised + +- The `# unname() matters` comment at `:329-332` still describes the `gcp_args` loop it now + sits two blocks above — readability only, carried over from round 1. +- `DESCRIPTION` `Date:` vs `NEWS.md` 0.7.0 date, carried over from round 1. diff --git a/planning/active/review-round3.md b/planning/active/review-round3.md new file mode 100644 index 0000000..528c34c --- /dev/null +++ b/planning/active/review-round3.md @@ -0,0 +1,294 @@ +# Code check — round 3 — fly#38 staged diff + +Scope: `git diff --cached`. Brief: treat **round 2's fixes** as the prime suspects, verify +they are complete and introduced nothing. + +**Verdict: round 2's fixes 1, 3 and 4 (the `user_rot` write-back, the discriminating film +fixture, the tolerance-pin test) are sound and verified. Fix 2 — replacing the shape gate +with a single tolerance of 1.10 — is off by 0.4% and leaves one of the two bundled cameras +outside the guard, in exactly the case the tolerance was chosen to cover.** + +Verified empirically, not by reading: + +- full suite under `pkgload::load_all()` + `test_dir()`: **FAIL 0 | WARN 0 | SKIP 0 | + PASS 1251**. Every finding below is invisible to it. +- `devtools::document()` re-run: **no drift** — `git status` unchanged afterwards. +- `lintr::lint("R/fly_georef.R")`: **1** lint against **5** at `HEAD`. The one remaining + (`no visible global function definition for 'fly_is_square'`) is the installed-namespace + artifact CLAUDE.md documents — confirmed: + `exists("fly_is_square", asNamespace("fly"))` is `FALSE` on the installed copy. +- Restore-the-bug run on a **new** assertion, both bindings patched, patch proven to have + taken. Result below. +- `data-raw/georef_calibrate-corner_mapping.R` exists; `^planning$` is in `.Rbuildignore`, + so `review-round*.md` do not ship. New test filename is portable. No new examples. + +--- + +## Findings + +### 1. **[bug]** `R/fly_georef.R:483` (`fly_gcp_stretch_max() = 1.10`) — the tolerance is 0.4% too loose, and a Leica DMC II frame on a `format_size` square footprint slips through the guard silently + +This is round 2's fix 2 and the case round 2's finding 3 was written about. The shape gate +was removed in favour of a tolerance specifically so that a digital frame sized through +`format_size` — which takes a single width and therefore produces a **square** footprint — +would still be checked. The new note says so in as many words +(`inst/notes/georeferencing.md`, "The stretch tolerance, and why it is not a shape gate"): + +> `format_size` sizes a frame from a single width, so a digital frame from a camera `fly` +> does not know lands on a **square** footprint — and that is precisely the case the guard +> exists for. A shape gate switches it off there. A tolerance wide enough for the rebate +> does not. + +For the **UltraCam** it does. For the **DMC II** it does not, by 0.4%: + +``` +threshold |log(1.10)| = 0.0953102 +DMC II |log(15552/14144)| = 0.0948987 <- shipped px counts +DMC II |log(972/884)| = 0.0948987 <- delivered thumbnail, aspect test :58 + 0.0948987 > 0.0953102 -> FALSE (NOT refused) +``` + +Confirmed end to end on a bundled DMC II frame, real `fly_footprint()` + real +`georef_one()` + real GDAL: + +``` +media: Digital - Colour +fly_footprint(one, format_size = c("Digital - Colour" = 87)) +footprint square? TRUE +image 1414 x 1555 px (DMC II aspect) + rot 0 aniso 1.099717 |log| 0.095053 refused=FALSE + rot 90 aniso 1.099717 |log| 0.095053 refused=FALSE + rot 180 aniso 1.099717 |log| 0.095053 refused=FALSE + rot 270 aniso 1.099717 |log| 0.095053 refused=FALSE + +georef_one returned: TRUE file written: TRUE warnings: 0 +``` + +A ~10% stretch, written silently, `success = TRUE` — the exact failure the code comment at +`:323-327` describes as "a valid GeoTIFF, in the right CRS, over the right ground — just +squashed … which nothing downstream would report". + +**Every row of `inst/extdata/camera_formats.csv` was checked**, `calib_file` and +`focal_length` alike (19 rows). Two columns matter and only one of them is in the shipped +table/comment: + +| case | quantity | tightest row | refused? | +|---|---|---|---| +| non-square footprint, 90-off pairing | `\|log(r²)\|` | DMC II, 0.18980 | **yes**, all 19 | +| **square footprint (`format_size`)** | `\|log(r)\|` | **DMC II, 0.09490** | **no** — the only row that slips | + +Every other row is refused in both columns (next-tightest square case is the UltraCam Eagle +at 0.42515). So the answer to the brief's question is: **yes, one shipped camera produces a +mispairing inside the tolerance, and it is the square-footprint case the tolerance replaced +the shape gate to cover.** + +**The two tables that argue for 1.10 both omit this row.** `georef_one():335-344` and the +note both list the square-footprint case as `0.442 — a portrait digital frame on a square +footprint`, which is the UltraCam. The DMC II's square case at 0.0949 sits *below* the +`0.095 the tolerance` line and appears in neither. + +**There is an admissible value; 1.10 is just outside it.** The two constraints are +`> 0.06454` (the 9600x9000 film rebate) and `< 0.09490` (DMC II square): + +``` +admissible fly_gcp_stretch_max() in (1.0667, 1.0995) +shipped 1.10 <- outside, by 0.4% +e.g. 1.08 -> |log| 0.07696 1.19x headroom below, 1.23x above +``` + +The band is narrow, and that narrowness is itself worth recording beside the constant — +the current comment claims "roughly 1.5x of headroom below and 2x above", which is true of +the `r²` column only. + +Note also the margin is 0.00026 in log space. A DMC II thumbnail delivered at a slightly +different rounding, or a full-res frame vs a thumbnail, flips this either way — so the +behaviour is not merely wrong, it is unstable at the boundary. + +--- + +### 2. **[fragile]** `tests/testthat/test-fly_georef_digital.R:235` — the tolerance test's third number picks the lenient bundled camera; written for the other one the same assertion fails + +Answering the brief's item 3 directly: the test **can** fail (proven below), and none of its +three numbers is arithmetically wrong. But the third assertion is named for a property that +does not hold: + +```r +expect_gt(abs(log(1654 / 1063)), tol) # portrait frame on a square footprint +``` + +`1654/1063` is the **UltraCam Eagle M3**. Substituting the other bundled camera — the same +claim, the same fixture set, one line — fails: + +``` +UltraCam Eagle M3 portrait frame on a square footprint: 0.4421015 > 0.09531 ? TRUE +Leica DMC II portrait frame on a square footprint: 0.0948987 > 0.09531 ? FALSE + +── Failure: a DMC II frame on a square footprint is also refused ── +Expected `abs(log(972/884))` > `log(fly:::fly_gcp_stretch_max())`. +Actual comparison: 0.09490 <= 0.09531 +``` + +So the test asserting "the tolerance sits between the film rebate and the tightest +mispairing" is asserting it against the second-tightest case. This is the `CLAUDE.md` +"fixture that cannot reach the failure mode" shape, at one remove: the fixture *set* has +both cameras and the assertion reaches only one of them. + +The same omission is in `georef_one()`'s comment table, the `@details` roxygen at +`R/fly_georef.R:90-95`, and the note's table — all four say the square-footprint case is +0.442 without saying it is camera-specific. + +**The test is otherwise real.** Restore-the-bug, tolerance reverted to round 2's 1.05, both +bindings patched: + +``` +proof patch took: namespace value = 1.05 / search-path value = 1.05 +proof it is REACHED: georef_one deparse mentions fly_gcp_stretch_max: TRUE + +unpatched FAIL=0 PASS=45 +patched 1.05 FAIL=10 PASS=35 +failing tests: + - a film scan carrying the negative's rebate still georeferences + - the stretch tolerance sits between the film rebate and the tightest mispairing +``` + +Both go red, so round 2's fix 3 (the 1250x1172 fixture replacing the non-discriminating +1250x1200) genuinely closed the hole round 2 named. That half is verified. + +--- + +### 3. **[fragile]** `R/fly_georef.R:69-70` and `man/fly_georef.Rd:93` — the new `case_when` comment is true only for `rotation = "auto"` + +Answering the brief's item 4. The claim: + +```r +.default = NA # square footprints fall through to 180, non-square to the + # digital mapping — NOT back to the per-photo bearing +``` + +Traced over the whole product of (footprint shape) x (`rotation` argument), with an all-`NA` +column present, `georef_one()` mocked at the boundary: + +``` + film NA column, rotation=auto -> 180, 180, 180, 180 <- claim holds + film NA column, rotation=0 -> 0, 0, 0, 0 <- claim FALSE + film NA column, rotation=90 -> 90, 90, 90, 90 <- claim FALSE + film NA column, rotation=270 -> 270, 270, 270, 270 <- claim FALSE + + digital NA column, rotation=auto -> 270, 270, 270, 270 <- claim holds + digital NA column, rotation=0 -> 270, 270, 270, 270 <- claim holds + digital NA column, rotation=90 -> 270, 270, 270, 270 <- claim holds + digital NA column, rotation=270 -> 270, 270, 270, 270 <- claim holds +``` + +The non-square half is correct in all four. The square half is correct in one of four — +`R/fly_georef.R:279` is `if (auto_rotation) 180L else rotation`, so an explicit `rotation` +argument is what an `NA` row falls through to, not 180. + +The behaviour is defensible (the caller gets what they asked for) and the fix is one +clause: "…fall through to the `rotation` argument, 180 by default". Raised because round 2's +finding 5 was about this exact comment being wrong, and the replacement is wrong in a +narrower way. Note also that `test-fly_georef_digital.R:88-97` pins the *digital* half of +this and nothing pins the square half. + +--- + +## Round 2's fixes — verified complete + +### Fix 1 — one parse into `user_rot`, unparseable refused, written back + +**Sound on every count the brief asked about.** + +- **Write-back safety across sf shapes.** `photos_sf[["rotation"]] <- user_rot` run over + `centroid_shapes()` (plain / tibble / grouped / `bcdc_sf`): value correct, `nrow` + preserved, `sf_column` attribute intact, class **set** preserved in all four. Only the + class *order* moves on `bcdc_sf` (`bcdc_sf, sf, …` → `sf, bcdc_sf, …`), which CLAUDE.md + already records as expected and which cannot matter here because `photos_sf` is a local + copy that is never returned. End to end with a factor column, all four shapes give + identical rotations (`180, 270, 180, 270`). +- **Nothing downstream is disturbed.** `fly_footprint()` runs at `:154`, *before* the + write-back at `:219`, so footprints are built from the caller's original column. + `fly_bearing()` at `:225` is inside `if (auto_rotation && !has_rotation_col)` and + `has_rotation_col` is `TRUE` whenever a user column exists, so it is unreachable on the + write-back path and cannot desynchronise anything. The only later reads are + `names()`, `[["airp_id"]]` and `[["rotation"]][j]`. +- **No remaining raw reads.** `:277` (`as.integer(photos_sf[["rotation"]][j])`) is the site + round 2's finding 1 named; it now reads the written-back integer vector, so `as.integer()` + is the identity there. Round 2's escape route — an unparseable factor level reaching `:277` + and yielding a level code — is closed twice over: the value is refused at `:212` before the + loop, *and* the write-back would have made `:277` correct anyway. Measured: + + ``` + factor c(180, north) STOP: `photos_sf$rotation` must be NA or one of 0, 90, 180, 270. Got: north. + factor('north') STOP: ... Got: north. + chr c(ninety,NA,180,NA) STOP: ... Got: ninety. <- round 2 finding 4, closed + int 360 STOP: ... Got: 360. + int -90 STOP: ... Got: -90. + all NA logical applied: 180, 180, 180, 180 + numeric 180.0 / chr '180' / list col applied: 180, ... + ``` + + The residual fragility is structural rather than live: the conversion still happens in two + places (`:206` and `:277`) and the second is correct only because of the write-back, which + is the hazard `:216-218`'s own comment warns about. Moving `:277` to `user_rot[j]` would + remove the dependence; today nothing is wrong. + +### Fix 3 — the discriminating film fixture + +Verified by restore-the-bug above (goes red at 1.05, green at 1.10). The premise assertions +at `:218-219` are correct: `|log(aniso)|` for 1250x1172 on a 4000 m square ring is +0.064412, matching the pinned `0.0644` and the 9600x9000 scan it stands in for +(`|log(9000/9600)| = 0.064539`). + +### Fix 4 — the `case_when` roxygen + +See finding 3 — the non-square half is right, the square half is right in one of four +combinations. + +--- + +## Checked independently and found sound + +- **Film is bit-for-bit unchanged for a square footprint and a square scan.** The GCP + construction was pulled from `HEAD`'s bytes (`git show HEAD:R/fly_georef.R`, sourced into + its own environment — not rewritten from memory) and compared against the staged one over + all 20 bundled film footprints x 4 rotations x 2 square scan sizes: **160 cases, 0 + differences**. Rotation *selection* re-traced over every combination of (`rotation` arg, + column present, column `NA`, `film_roll`/`frame_number` present) and is unchanged, except + where a previously-silent invalid column value now `stop()`s — a deliberate, NEWS-worthy + behaviour change, and the safe direction. The new anisotropy guard is vacuous at + `aniso = 1`, which `test-fly_georef_aspect.R:110` asserts. +- **The guard cannot be reached with a non-finite value.** `!is.finite(aniso)` + short-circuits before `log()`; `w == 0` gives `log(0) = -Inf`, which exceeds the threshold + and refuses rather than errors; degenerate dimensions return `NA_real_` + (`test-fly_georef_gcps.R:162`). +- **No fixture makes both sides of a comparison identical.** `test-fly_georef_aspect.R:50` + compares the shipped CSV against hand-measured thumbnail dimensions — genuinely + independent numbers. `test-fly_georef_digital.R:67` computes the film expectation from + `bearing_to_rotation(fly_bearing(...))`, the same functions the code calls, but the + bundled bearings give `90, 270, 270, 90` and the two `90`s discriminate against a + "270 everywhere" implementation; the test's comment states the limitation. The digital + premise at `:40-41` is the load-bearing one and it is asserted (bearings ~343° map to 0, + not 270). +- **The `stop()` cannot leave a half-written batch** — it fires before the loop. It does run + after `dir.create()` and after `fly_footprint()`, so an invalid column costs a footprint + computation before aborting; no correctness consequence. +- **`R CMD check` surface.** New file `tests/testthat/test-fly_georef_digital.R` — portable + name, no spaces. `planning/` (both new review files) is excluded by `^planning$` in + `.Rbuildignore`. `data-raw/` likewise, so the `data-raw/georef_calibrate-corner_mapping.R` + pointer in NEWS/roxygen names a file that is not installed — consistent with existing house + style, not raised. No new `@examples`; the existing block is unchanged and its + network dependence is pre-existing. `devtools::document()` produces no drift. + `testthat (>= 3.2.0)` bump is required by `local_mocked_bindings(.package =)` and is + present. +- **`empty_fp` / `non_square` three-way classification** — unchanged from round 1's + verification; `if (empty_fp[j]) next` precedes the only `footprints[j, ]` subset. +- No hardcoded absolute paths, no secrets, no shell quoting, no writes outside `dest_dir` / + `tempfile()`. + +## Not raised + +- `# unname() matters` at `R/fly_georef.R:355-358` still describes the `gcp_args` loop it now + sits above rather than beside — carried from rounds 1 and 2, readability only. +- `DESCRIPTION` `Date: 2026-08-30` now agrees with NEWS — the round-1/2 note is stale. +- `as.integer("180.5")` truncates to `180` and is accepted. Consistent with `as.integer()` + semantics, and 269.9 is refused by name, so no silent wrong rotation is reachable. diff --git a/planning/active/task_plan.md b/planning/active/task_plan.md index 06c4ee2..d6debb7 100644 --- a/planning/active/task_plan.md +++ b/planning/active/task_plan.md @@ -127,18 +127,18 @@ All three measurements below use public data, so the derivation lives in fly. ## Phase 4: Land the constant and remove the exclusion -- [ ] `fly_georef()` applies the constant to non-square footprints and does **not** apply +- [x] `fly_georef()` applies the constant to non-square footprints and does **not** apply `bearing_to_rotation()` to them — the ring already carries the bearing. -- [ ] Branch on `fly_is_square()` (`R/fly_footprint.R:165`), computed for every row before +- [x] Branch on `fly_is_square()` (`R/fly_footprint.R:165`), computed for every row before the per-row loop. Never on `half_cross`/`half_along`, which are NA by construction until a sizing route fills them (CLAUDE.md, fly#32's three-round trap). -- [ ] Keep the `rotation`-column override for non-square frames; document the +- [x] Keep the `rotation`-column override for non-square frames; document the carried-column hazard in `@param rotation`. -- [ ] Delete the exclusion block (`R/fly_georef.R:122-143`) and the `rotated[fp_idx[1]]` +- [x] Delete the exclusion block (`R/fly_georef.R:122-143`) and the `rotated[fp_idx[1]]` skip (`:190`). Keep `fly_is_square()` itself — `test-fly_camera_format.R:266` uses it. -- [ ] Keep the empty-geometry skip (`:195`) and `fly_warn_unsized()` (`:120`) untouched; +- [x] Keep the empty-geometry skip (`:195`) and `fly_warn_unsized()` (`:120`) untouched; they cover a different exclusion. -- [ ] Rewrite the **Rotation** `@details` section (`:38-65`), which currently documents +- [x] Rewrite the **Rotation** `@details` section (`:38-65`), which currently documents only the film scheme, and the `@param rotation` text that says `"auto"` applies everywhere. @@ -148,27 +148,29 @@ bearing rotation and digital the constant, asserted per row; no non-square warni ## Phase 5: Notes, docs, release -- [ ] `inst/notes/georeferencing.md`, companion to `terrain-correction.md` and +- [x] `inst/notes/georeferencing.md`, companion to `terrain-correction.md` and `camera-formats.md`: the ring-order contract, why the bearing must not be applied twice, what the QA measured and against what, and what the aspect invariant can and cannot catch. Ground truth referenced obliquely — no repo, endpoint or database named. -- [ ] Update CLAUDE.md: the Key Decisions entry for #30 and the NEWS line +- [x] Update CLAUDE.md: the Key Decisions entry for #30 and the NEWS line "`fly_georef()` excludes rotated footprints" are both now stale. -- [ ] `devtools::document()` — **read its output**. `Writing '.Rd'` or a +- [x] `devtools::document()` — **read its output**. `Writing '.Rd'` or a falling `grep -c "^export(" NAMESPACE` is the fly#30 rebind. -- [ ] `lintr::lint_package()` against the `HEAD` baseline, not against zero. -- [ ] NEWS entry; version 0.6.0 -> 0.7.0 as the **final** commit of the branch. +- [x] `lintr::lint_package()` against the `HEAD` baseline, not against zero. `R/fly_georef.R` + 1 against a baseline of 5; `data-raw/` script cleaned to 0 +- [x] NEWS entry; version 0.6.0 -> 0.7.0 as the **final** commit of the branch. ## Validation -- [ ] `devtools::test()` passes -- [ ] `/code-check` clean on each commit -- [ ] No fly file names the private ortho repo, its endpoint, or its database: +- [x] `devtools::test()` passes — 1255, 0 failures +- [x] `/code-check` — 3 rounds, 1 + 3 + 3 findings; every round found a defect inside the + previous round's fix. See `review-round{1,2,3}.md` +- [x] No fly file names the private ortho repo, its endpoint, or its database: `git diff main --stat` reviewed, plus a grep of the branch diff for the repo name - [x] Restore-the-bug check on the Phase 2 invariant: it goes red against the rejected rotations, patched in **both** `asNamespace("fly")` and `as.environment("package:fly")`, with a printed value proving the patch took -- [ ] PWF checkboxes match landed work; `/planning-archive` on completion +- [x] PWF checkboxes match landed work; `/planning-archive` on completion ## Open, deliberately not decided yet diff --git a/tests/testthat/test-fly_georef_digital.R b/tests/testthat/test-fly_georef_digital.R new file mode 100644 index 0000000..0ff5d00 --- /dev/null +++ b/tests/testthat/test-fly_georef_digital.R @@ -0,0 +1,289 @@ +# What rotation each frame actually receives, observed at the boundary rather than +# inferred from the output. `georef_one()` is mocked so these need no network and no +# GDAL — the question is which rotation was chosen and why, not what GDAL did with it. + +# Record the rotation `georef_one()` is called with, per source file. +# `.env` is the environment the mock is CLEANED UP with, not the one it is installed +# in — `.package` names the target. Passing `asNamespace("fly")` to `.env` installs the +# mock correctly and then never removes it, because a namespace does not exit: every +# later test in the run keeps the stub. That leaks silently in the direction that reads +# as success, since the stub returns TRUE. +capture_rotations <- function(expr) { + seen <- list() + testthat::local_mocked_bindings( + georef_one = function(src, fp, out_file, srcnodata = "0", rotation = 180) { + seen[[basename(src)]] <<- rotation + TRUE + }, + .package = "fly", + .env = parent.frame() # unwinds with the calling test, not with this helper + ) + force(expr) + seen +} + +fake_fetch <- function(photos) { + files <- file.path(tempdir(), paste0("f", photos$airp_id, ".jpg")) + for (f in files) if (!file.exists(f)) writeLines("x", f) + dplyr::tibble(airp_id = photos$airp_id, dest = files, success = TRUE) +} + +digital_photos <- function(i = 19:24) { + sf::st_read(testdata_path("photo_centroids_digital.gpkg"), quiet = TRUE)[i, ] +} + + +test_that("a non-square footprint gets the measured digital rotation, not the bearing", { + photos <- digital_photos() + # The premise: without this the test would pass for a bearing rule that happened to + # agree. These frames fly ~343 degrees, which `bearing_to_rotation()` maps to 0. + expect_equal(mean(fly_bearing(photos)$bearing, na.rm = TRUE), 343, tolerance = 1) + expect_false(all(bearing_to_rotation(fly_bearing(photos)$bearing) == fly_digital_rotation())) + + seen <- capture_rotations( + suppressWarnings(fly_georef(fake_fetch(photos), photos, dest_dir = tempfile())) + ) + expect_length(seen, nrow(photos)) + expect_equal(unique(unlist(seen)), fly_digital_rotation()) +}) + + +test_that("a mixed batch rotates film by bearing and digital by the constant", { + film <- sf::st_transform(sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE), 3005) + dig <- sf::st_transform(digital_photos(), 3005) + cols <- intersect(names(film), names(dig)) + mix <- rbind(film[1:4, cols], dig[1:3, cols]) + + fp <- suppressWarnings(fly_footprint(mix)) + expect_equal(fly_is_square(fp), c(rep(TRUE, 4), rep(FALSE, 3))) # premise + + seen <- capture_rotations( + suppressWarnings(fly_georef(fake_fetch(mix), mix, dest_dir = tempfile())) + ) + got <- unlist(seen)[paste0("f", mix$airp_id, ".jpg")] + + expect_equal(unname(got[5:7]), rep(fly_digital_rotation(), 3)) + # Film keeps the pre-existing rule exactly, whatever it returns for these bearings. + expect_equal(unname(got[1:4]), bearing_to_rotation(fly_bearing(mix)$bearing[1:4])) +}) + + +test_that("a user rotation column overrides the digital constant", { + photos <- digital_photos() + photos$rotation <- 180L + seen <- capture_rotations( + suppressWarnings(fly_georef(fake_fetch(photos), photos, dest_dir = tempfile())) + ) + expect_equal(unique(unlist(seen)), 180L) + + # NA in the column falls through to the constant rather than to the film default. + photos$rotation <- NA_integer_ + seen <- capture_rotations( + suppressWarnings(fly_georef(fake_fetch(photos), photos, dest_dir = tempfile())) + ) + expect_equal(unique(unlist(seen)), fly_digital_rotation()) +}) + + +test_that("the `rotation` argument does not silently override a digital frame", { + # It is documented as applying to square footprints only. Asserted so that a future + # change making it apply everywhere has to disagree with a test. + photos <- digital_photos() + seen <- capture_rotations( + suppressWarnings(fly_georef(fake_fetch(photos), photos, dest_dir = tempfile(), + rotation = 0)) + ) + expect_equal(unique(unlist(seen)), fly_digital_rotation()) +}) + + +test_that("frames with no footprint are still skipped, and warned about once", { + photos <- mixed_media_fixture() + fp <- suppressWarnings(fly_footprint(photos)) + unsized <- sf::st_is_empty(sf::st_geometry(fp)) + expect_true(any(unsized)) # premise + + # Counted, not merely suppressed. The test is named for the warning and previously + # only checked how many frames reached `georef_one()`, which passes just as happily + # for zero warnings as for one per frame. + # + # `expect_warning()` is the wrong instrument here: this fixture also raises + # `fly_footprint()`'s unknown-format warning, and testthat re-raises every warning it + # did not match, so the unmatched one surfaces as a test WARNING. Handling them all + # and counting the one under test keeps the assertion exact and the run clean. + n <- 0L + seen <- withCallingHandlers( + capture_rotations(fly_georef(fake_fetch(photos), photos, dest_dir = tempfile())), + warning = function(w) { + if (grepl("have no footprint", conditionMessage(w))) n <<- n + 1L + invokeRestart("muffleWarning") + } + ) + expect_identical(n, 1L) + expect_length(seen, sum(!unsized)) +}) + + +test_that("a non-square footprint with no bearing is warned about", { + # `fly_bearing()` needs a neighbour, so one frame on its own is the ordinary way to + # reach this — not an exotic case. + one <- digital_photos(19) + expect_warning( + capture_rotations(fly_georef(fake_fetch(one), one, dest_dir = tempfile())), + "no flight bearing" + ) + # And it does not fire when the bearing is available. + many <- digital_photos() + expect_no_warning( + capture_rotations(fly_georef(fake_fetch(many), many, dest_dir = tempfile())) + ) +}) + + +test_that("a mapping that would stretch the image is refused, not written squashed", { + skip_if_no_terra() + + # A synthetic portrait image and a footprint of matching aspect. Written locally so + # this needs no network — the guard is about geometry, not about pixels. + src <- tempfile(fileext = ".tif") + r <- terra::rast(nrows = 200, ncols = 100, vals = seq_len(20000)) + terra::writeRaster(r, src, overwrite = TRUE) + + ring <- function(hc, ha) { + sf::st_sf(geometry = sf::st_sfc(sf::st_polygon(list(matrix( + c(-hc, -ha, hc, -ha, hc, ha, -hc, ha, -hc, -ha), + ncol = 2, byrow = TRUE + ) + rep(c(1e6, 1e6), each = 5))), crs = 3005)) + } + + # 100 x 200 px onto a 1000 x 2000 m footprint: the height axis is the long one on both + # sides, so rotation 270 pairs them isotropically and the file is written. + ok <- ring(hc = 1000, ha = 500) + out <- tempfile(fileext = ".tif") + expect_true(georef_one(src, ok, out, rotation = fly_digital_rotation())) + expect_true(file.exists(out)) + + # The same image on the same footprint at rotation 0 pairs 100 px with the 2000 m edge + # and 200 px with the 1000 m one — a 4x stretch. Refused. + out2 <- tempfile(fileext = ".tif") + expect_warning(res <- georef_one(src, ok, out2, rotation = 0), "Skipped rather than") + expect_false(res) + expect_false(file.exists(out2)) + + # A square footprint is NOT exempt. It has no pairing to get wrong, but a gross + # mismatch between image and footprint still means they disagree about the frame — + # and exempting it by shape would switch the guard off for a digital frame sized + # through `format_size` into a square footprint, which is the unknown-camera case the + # guard exists for. 100 x 200 px on a 2000 x 2000 m footprint is a 2x stretch. + sq <- ring(hc = 1000, ha = 1000) + out3 <- tempfile(fileext = ".tif") + expect_warning(res3 <- georef_one(src, sq, out3, rotation = 0), "Skipped rather than") + expect_false(res3) + expect_false(file.exists(out3)) + + # A square image on that square footprint is consistent, and is written. + square_image <- tempfile(fileext = ".tif") + terra::writeRaster(terra::rast(nrows = 100, ncols = 100, vals = seq_len(10000)), + square_image, overwrite = TRUE) + out4 <- tempfile(fileext = ".tif") + expect_no_warning(res4 <- georef_one(square_image, sq, out4, rotation = 0)) + expect_true(res4) +}) + + +test_that("a film scan carrying the negative's rebate still georeferences", { + skip_if_no_terra() + # The regression the tolerance is sized for. Every bundled film thumbnail is exactly + # 1250 x 1250, so the fixture set cannot reach this — a full-resolution 9-inch scan + # including the rebate is the ordinary case and lands a few percent off square. 1250 x + # 1172 has the same aspect as the 9600 x 9000 scan that motivated it, at 1/60th the + # pixels. + # + # This is the closest legitimate case to the threshold, which is what makes it worth + # pinning: 0.064 against a tolerance of 0.095. An earlier draft used 1250 x 1200, + # which sits at 0.041 and would have passed a 5% tolerance too — a fixture that could + # not reach the failure it was written for. + src <- tempfile(fileext = ".tif") + terra::writeRaster(terra::rast(nrows = 1172, ncols = 1250, vals = seq_len(1465000)), + src, overwrite = TRUE) + + sq <- sf::st_sf(geometry = sf::st_sfc(sf::st_polygon(list(matrix( + c(-2000, -2000, 2000, -2000, 2000, 2000, -2000, 2000, -2000, -2000), + ncol = 2, byrow = TRUE + ) + rep(c(1e6, 1e6), each = 5))), crs = 3005)) + expect_true(fly_is_square(sq)) # premise + + gcp <- fly_georef_gcps(1250L, 1172L, sf::st_coordinates(sq)[1:4, , drop = FALSE], 0) + aniso <- fly_gcp_anisotropy(gcp, 1250L, 1172L) + expect_equal(abs(log(aniso)), 0.0644, tolerance = 1e-3) # premise + expect_lt(abs(log(aniso)), log(fly_gcp_stretch_max())) # premise + + for (rot in c(0, 90, 180, 270)) { + out <- tempfile(fileext = ".tif") + expect_no_warning(res <- georef_one(src, sq, out, rotation = rot)) + expect_true(res) + } +}) + + +test_that("the stretch tolerance clears every shipped camera, not just the bundled two", { + # The whole argument for the tolerance, computed from the shipped table rather than + # from remembered numbers. An earlier version asserted the square-footprint case using + # the UltraCam at 0.442 — the most eccentric camera, which any tolerance clears. The + # binding case is the LEAST eccentric one, and picking the lenient example is how the + # tolerance came to be set 0.4% too loose and let a DMC II frame through. + tol <- log(fly_gcp_stretch_max()) + cf <- utils::read.csv(system.file("extdata/camera_formats.csv", package = "fly")) + aspect <- ifelse(!is.na(cf$px_cross) & !is.na(cf$px_along), + cf$px_cross / cf$px_along, cf$width_mm / cf$height_mm) + expect_gt(length(aspect), 10) # premise + + # Below: the worst disagreement a legitimate frame produces — a full-resolution + # 9-inch scan carrying the negative's rebate. + expect_lt(abs(log(9000 / 9600)), tol) + + # Above: EVERY shipped camera, on a square footprint (the `format_size` route) and + # mispaired on its own footprint. No row may slip through. + expect_true(all(abs(log(aspect)) > tol)) + expect_true(all(abs(log(aspect^2)) > tol)) + + # And the margin is real on both sides rather than incidental. + expect_equal(min(abs(log(aspect))), 0.0949, tolerance = 1e-3) # DMC II + expect_gt(min(abs(log(aspect))) / tol, 1.15) + expect_gt(tol / abs(log(9000 / 9600)), 1.15) +}) + + +test_that("an invalid `rotation` column is refused by name, not by GDAL", { + photos <- digital_photos() + + # 360 shifts by four and indexes past the ring; without this check it surfaces from + # inside `tryCatch` as "subscript out of bounds", naming neither the column nor the + # value. 45 and -90 are worse — they shift by zero and georeference silently wrong. + for (bad in list(360L, 45L, -90L)) { + photos$rotation <- bad + expect_error(fly_georef(fake_fetch(photos), photos, dest_dir = tempfile()), + "must be NA or one of", info = paste("rotation =", bad)) + } + + # A factor is the realistic way to get an unexpected type here — a `rotation` column + # read from a CSV. `as.integer()` on it returns the LEVEL CODE, so a naive read both + # validates 180 as 1 and applies it as 1. The value is normalised once and read from + # that one place, so the factor is honoured as 180 rather than either erroring or + # silently becoming rotation 0. + photos$rotation <- factor("180") + seen <- capture_rotations( + suppressWarnings(fly_georef(fake_fetch(photos), photos, dest_dir = tempfile())) + ) + expect_equal(unique(unlist(seen)), 180L) + + # The level code is 1, so a version that converted in two places would land here. + expect_false(identical(unique(unlist(seen)), 1L)) + + # And the valid values still get through. + photos$rotation <- 90L + seen <- capture_rotations( + suppressWarnings(fly_georef(fake_fetch(photos), photos, dest_dir = tempfile())) + ) + expect_equal(unique(unlist(seen)), 90L) +}) diff --git a/tests/testthat/test-fly_georef_gcps.R b/tests/testthat/test-fly_georef_gcps.R index 372aeca..8552548 100644 --- a/tests/testthat/test-fly_georef_gcps.R +++ b/tests/testthat/test-fly_georef_gcps.R @@ -112,3 +112,60 @@ test_that("rotation 0 and 180 put image width on the cross-track edge, 90 and 27 expect_equal(width_edge_m(90), 2040) expect_equal(width_edge_m(270), 2040) }) + + +test_that("the ground quad keeps one handedness at every rotation", { + # The mirror the aspect invariant cannot see. Mapping the pixel corners onto the ring + # traversed the other way leaves every edge length identical, the aspect ratio + # identical, and the output reflected. Signed area is what separates them, and it is a + # property of the whole quad rather than of any one corner. + signed_area <- function(g) { + x <- g[, "ground_x"]; y <- g[, "ground_y"] + k <- c(2, 3, 4, 1) + sum(x * y[k] - x[k] * y) / 2 + } + + for (ring in list(film_ring(), ultracam_ring())) { + s <- vapply(c(0, 90, 180, 270), function(r) { + sign(signed_area(fly_georef_gcps(1063L, 1654L, ring, r))) + }, numeric(1)) + expect_equal(s, rep(-1, 4)) + } + + # And it can fail: reverse the ring and the sign flips. + reversed <- ultracam_ring()[c(4, 3, 2, 1), , drop = FALSE] + expect_equal(sign(signed_area(fly_georef_gcps(1063L, 1654L, reversed, 0))), 1) +}) + + +test_that("fly_gcp_anisotropy measures the stretch a mapping would apply", { + # 3175 x 2040 m footprint, 1063 x 1654 px portrait image. The measured mapping is + # isotropic; the pair it rejects stretches by the footprint's aspect squared. + aniso <- function(rot) { + fly_gcp_anisotropy(fly_georef_gcps(1063L, 1654L, ultracam_ring(), rot), 1063L, 1654L) + } + # 1e-3, not 1e-6: a thumbnail is a rounded downscale of the sensor, so even the + # correct pairing lands at 0.99974 rather than exactly 1. That 0.03% is the headroom + # the guard's 5% tolerance has to clear, against a 21% minimum for a wrong pairing + # (the DMC II, the least eccentric camera in the shipped table). + expect_equal(aniso(fly_digital_rotation()), 1, tolerance = 1e-3) + expect_equal(aniso(90), 1, tolerance = 1e-3) + # (m/px across) / (m/px along) when the axes are paired the wrong way round. Not + # exactly the footprint's aspect squared, because the thumbnail is a rounded + # downscale of the sensor rather than an exact one. + wrong <- (3175 / 1063) / (2040 / 1654) + expect_equal(wrong, 2.4217, tolerance = 1e-4) + expect_equal(aniso(0), wrong, tolerance = 1e-6) + expect_equal(aniso(180), wrong, tolerance = 1e-6) + + # Degenerate dimensions give NA rather than a number the guard would act on. + expect_true(is.na(fly_gcp_anisotropy(fly_georef_gcps(10L, 10L, ultracam_ring(), 0), 0, 10))) +}) + + +test_that("the digital rotation constant is 270 and is not derived at run time", { + # Pinned deliberately. It was measured (inst/notes/georeferencing.md); a future reader + # who recomputes it from the ring geometry will get a plausible wrong answer, because + # the geometry alone cannot separate 270 from 90. + expect_identical(fly_digital_rotation(), 270L) +}) diff --git a/vignettes/airphoto-selection.Rmd b/vignettes/airphoto-selection.Rmd index 4d4e96f..9f274d0 100644 --- a/vignettes/airphoto-selection.Rmd +++ b/vignettes/airphoto-selection.Rmd @@ -302,6 +302,15 @@ argument. They are approximate either way, useful for visual context rather than survey-grade positioning. Metadata from the original centroid data (date, scale, focal length) links back via `airp_id`. +Digital frames georeference too. Their footprints are not square and are +already rotated onto the flight line, so they use a fixed corner mapping rather +than the bearing rule film uses — see the **Rotation** section of +`?fly_georef`. One consequence is worth knowing when experimenting: passing a +single frame gives `fly_bearing()` no neighbour to work from, so its footprint +is drawn axis-aligned and warps as though the flight line ran due north. +`fly_georef()` warns when that happens. Pass consecutive frames from the same +roll to get the real azimuth. + # Terrain-adjusted footprints Everything above sizes each frame from its reported scale. That scale is From 5a77a0a160bd48fcbbb496194789f61e9ed1e4d0 Mon Sep 17 00:00:00 2001 From: almac2022 Date: Sun, 30 Aug 2026 16:12:56 -0700 Subject: [PATCH 6/7] Release v0.7.0 Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01GBKqedyBysV7hB4DuL98ZR --- DESCRIPTION | 6 +++--- NEWS.md | 9 +++++++++ 2 files changed, 12 insertions(+), 3 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index ba7934b..a76dc1b 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: fly Title: Historic Airphoto Footprints, Selection and Georeferencing for British Columbia -Version: 0.6.0 -Date: 2026-08-29 +Version: 0.7.0 +Date: 2026-08-30 Authors@R: c( person("Allan", "Irvine", , "al@newgraphenvironment.com", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-3495-2128")), @@ -32,7 +32,7 @@ Suggests: furrr, future, terra, - testthat (>= 3.0.0), + testthat (>= 3.2.0), knitr, rmarkdown Config/testthat/edition: 3 diff --git a/NEWS.md b/NEWS.md index 5107200..387590b 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,14 @@ # fly (development version) +## 0.7.0 (2026-08-30) + +- `fly_georef()` georeferences digital frames ([#38](https://github.com/NewGraphEnvironment/fly/issues/38)). v0.6.0 gave them footprints and then excluded them from georeferencing with a warning; that exclusion is gone, so the whole post-2010 catalogue is now georeferenceable rather than only sizeable +- **The corner mapping for a non-square footprint is rotation 270 — the top-left pixel maps to the ring's rear-left corner, equivalently image columns run in the flight direction and image rows run flight-right.** It was measured, not reasoned: a digital footprint is already rotated onto its flight line, so `bearing_to_rotation()` is not applied to it, and the geometry that remains cannot distinguish 270 from 90 +- The measurement needed no licence-restricted imagery, which the issue had assumed it would. Three independent public routes agree: the per-frame exterior orientation the catalogue publishes through `patb_georef_url` (the UltraCam Eagle's mount is rigid to 0.18 degrees over 6839 frames spanning the compass); adjacent-frame overlap correlation, which needs no reference imagery at all because consecutive frames check each other (+0.616 and +0.659 against at most +0.43 for the alternatives); and FWA lake darkness. `data-raw/georef_calibrate-corner_mapping.R` reproduces all three, and `inst/notes/georeferencing.md` records them — including the one that disagreed and was wrong +- A frame whose delivered image aspect does not pair with its footprint edges is now **skipped with a warning** rather than written stretched. A wrong pairing produces a valid GeoTIFF in the right CRS over the right ground, squashed by the aspect ratio squared, which nothing downstream would report. This also catches a frame sized from an inferred camera format that does not match the camera that took it +- A non-square footprint built without a flight bearing is drawn axis-aligned and so georeferences as though the flight line ran due north. `fly_bearing()` needs a neighbouring frame, so this is the ordinary result of georeferencing one frame on its own, and it is now warned about rather than left to be noticed in the output +- The `rotation` argument applies to square footprints only, and is documented as such; a `rotation` column in `photos_sf` still overrides per-photo for both. Carrying a film-era `rotation` column into a digital batch therefore overrides the correct mapping — drop the column, or set it to `NA` for those rows + ## 0.6.0 (2026-08-30) - `fly_footprint()` now sizes digital frames, closing the gap #30 made honest but left open ([#32](https://github.com/NewGraphEnvironment/fly/issues/32)). Province-wide that is 223,667 of 1,670,471 frames — the package was quietly film-only for anything after ~2010 From 8fd5caac484b6296281ecfaf74bb772172120c4d Mon Sep 17 00:00:00 2001 From: almac2022 Date: Sun, 30 Aug 2026 16:13:22 -0700 Subject: [PATCH 7/7] Archive PWF for #38 Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01GBKqedyBysV7hB4DuL98ZR --- .../README.md | 26 +++++++++++++++++++ .../findings.md | 0 .../progress.md | 0 .../review-round1.md | 0 .../review-round2.md | 0 .../review-round3.md | 0 .../task_plan.md | 0 7 files changed, 26 insertions(+) create mode 100644 planning/archive/2026-08-issue-38-georeference-digital/README.md rename planning/{active => archive/2026-08-issue-38-georeference-digital}/findings.md (100%) rename planning/{active => archive/2026-08-issue-38-georeference-digital}/progress.md (100%) rename planning/{active => archive/2026-08-issue-38-georeference-digital}/review-round1.md (100%) rename planning/{active => archive/2026-08-issue-38-georeference-digital}/review-round2.md (100%) rename planning/{active => archive/2026-08-issue-38-georeference-digital}/review-round3.md (100%) rename planning/{active => archive/2026-08-issue-38-georeference-digital}/task_plan.md (100%) diff --git a/planning/archive/2026-08-issue-38-georeference-digital/README.md b/planning/archive/2026-08-issue-38-georeference-digital/README.md new file mode 100644 index 0000000..0f232a7 --- /dev/null +++ b/planning/archive/2026-08-issue-38-georeference-digital/README.md @@ -0,0 +1,26 @@ +# fly#38 — Georeference digital frames + +Closed by PR against `main`, v0.7.0. `fly_georef()` had excluded every non-square +footprint, so the digital frames #32 had just sized could not be georeferenced. + +The corner mapping is **rotation 270** — the top-left pixel maps to the footprint ring's +rear-left corner. The issue assumed establishing it needed licence-restricted orthophoto +imagery; it did not. The catalogue publishes per-frame exterior orientation through +`patb_georef_url`, a column `fly_fetch()` already supported, and consecutive frames +overlap enough to check each other. Three public routes agree, and the derivation ships +in `data-raw/georef_calibrate-corner_mapping.R` with the record in +`inst/notes/georeferencing.md`. + +Two things worth carrying forward, both in `progress.md` in full: + +- **The most convincing measurement was the wrong one.** Read naively the exterior + orientation puts the Leica DMC II at rotation 90, on 97.6% agreement across four pooled + projects. That camera's project flies east and west only, so its own data separates the + rigid-mount and reflected-frame hypotheses at 97.6% against 98.4% — not at all. +- **Every one of three code-check rounds found a defect inside the previous round's + fix**, and the suite was green throughout. Round 1 caught a guard that refused + full-resolution film scans; round 2 caught the shape gate that fixed it switching the + guard off for the exact case it existed for, plus an incomplete factor normalisation; + round 3 caught the replacement tolerance being 0.4% too loose, guarded by a test that + asserted against the most eccentric camera instead of the least. Convergence never + arrived within three rounds — the reviews are kept beside this README. diff --git a/planning/active/findings.md b/planning/archive/2026-08-issue-38-georeference-digital/findings.md similarity index 100% rename from planning/active/findings.md rename to planning/archive/2026-08-issue-38-georeference-digital/findings.md diff --git a/planning/active/progress.md b/planning/archive/2026-08-issue-38-georeference-digital/progress.md similarity index 100% rename from planning/active/progress.md rename to planning/archive/2026-08-issue-38-georeference-digital/progress.md diff --git a/planning/active/review-round1.md b/planning/archive/2026-08-issue-38-georeference-digital/review-round1.md similarity index 100% rename from planning/active/review-round1.md rename to planning/archive/2026-08-issue-38-georeference-digital/review-round1.md diff --git a/planning/active/review-round2.md b/planning/archive/2026-08-issue-38-georeference-digital/review-round2.md similarity index 100% rename from planning/active/review-round2.md rename to planning/archive/2026-08-issue-38-georeference-digital/review-round2.md diff --git a/planning/active/review-round3.md b/planning/archive/2026-08-issue-38-georeference-digital/review-round3.md similarity index 100% rename from planning/active/review-round3.md rename to planning/archive/2026-08-issue-38-georeference-digital/review-round3.md diff --git a/planning/active/task_plan.md b/planning/archive/2026-08-issue-38-georeference-digital/task_plan.md similarity index 100% rename from planning/active/task_plan.md rename to planning/archive/2026-08-issue-38-georeference-digital/task_plan.md