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/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 diff --git a/R/fly_georef.R b/R/fly_georef.R index 843217f..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) @@ -238,43 +320,51 @@ 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 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) } - # Build GCP args mapping pixel corners to (rotated) footprint corners + # `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 +417,98 @@ 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 +} + +#' 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 new file mode 100644 index 0000000..f1a4f78 --- /dev/null +++ b/data-raw/georef_calibrate-corner_mapping.R @@ -0,0 +1,199 @@ +# 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", 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) { + 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_query <- bcdata::bcdc_query_geodata("WHSE_BASEMAPPING.FWA_LAKES_POLY") +lakes <- st_transform( + bcdata::collect(bcdata::filter(lakes_query, 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..d632235 --- /dev/null +++ b/inst/notes/georeferencing.md @@ -0,0 +1,150 @@ +# 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 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 +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/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/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/archive/2026-08-issue-38-georeference-digital/findings.md b/planning/archive/2026-08-issue-38-georeference-digital/findings.md new file mode 100644 index 0000000..41d4253 --- /dev/null +++ b/planning/archive/2026-08-issue-38-georeference-digital/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/archive/2026-08-issue-38-georeference-digital/progress.md b/planning/archive/2026-08-issue-38-georeference-digital/progress.md new file mode 100644 index 0000000..78f7197 --- /dev/null +++ b/planning/archive/2026-08-issue-38-georeference-digital/progress.md @@ -0,0 +1,231 @@ +# 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 + +### 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 + +### 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 + +### 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. + +### 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/archive/2026-08-issue-38-georeference-digital/review-round1.md b/planning/archive/2026-08-issue-38-georeference-digital/review-round1.md new file mode 100644 index 0000000..c037db2 --- /dev/null +++ b/planning/archive/2026-08-issue-38-georeference-digital/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/archive/2026-08-issue-38-georeference-digital/review-round2.md b/planning/archive/2026-08-issue-38-georeference-digital/review-round2.md new file mode 100644 index 0000000..d39ffeb --- /dev/null +++ b/planning/archive/2026-08-issue-38-georeference-digital/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/archive/2026-08-issue-38-georeference-digital/review-round3.md b/planning/archive/2026-08-issue-38-georeference-digital/review-round3.md new file mode 100644 index 0000000..528c34c --- /dev/null +++ b/planning/archive/2026-08-issue-38-georeference-digital/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/archive/2026-08-issue-38-georeference-digital/task_plan.md b/planning/archive/2026-08-issue-38-georeference-digital/task_plan.md new file mode 100644 index 0000000..d6debb7 --- /dev/null +++ b/planning/archive/2026-08-issue-38-georeference-digital/task_plan.md @@ -0,0 +1,182 @@ +# 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 + +- [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()`. +- [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). +- [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 +pulled from git, and fails if the extraction changed any coordinate. + +## Phase 2: Pin the aspect invariant (the half needing no imagery) + +- [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`. +- [x] Run it over the bundled digital frames, **both** cameras, from + `inst/testdata/photo_centroids_digital.gpkg` via `digital_fixture()` / + `tests/testthat/setup.R`. +- [x] Assert the invariant **fails** at the rotations it must reject. A test that only + ever passes is decoration. +- [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 +the DMC II / UltraCam distortion factors are recorded so a future fixture change that +drops the UltraCam is visible. + +## 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 + +- [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. +- [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). +- [x] Keep the `rotation`-column override for non-square frames; document the + carried-column hazard in `@param rotation`. +- [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. +- [x] Keep the empty-geometry skip (`:195`) and `fly_warn_unsized()` (`:120`) untouched; + they cover a different exclusion. +- [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. + +**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 + +- [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. +- [x] Update CLAUDE.md: the Key Decisions entry for #30 and the NEWS line + "`fly_georef()` excludes rotated footprints" are both now stale. +- [x] `devtools::document()` — **read its output**. `Writing '.Rd'` or a + falling `grep -c "^export(" NAMESPACE` is the fly#30 rebind. +- [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 + +- [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 +- [x] 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. 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)) + } +}) 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 new file mode 100644 index 0000000..8552548 --- /dev/null +++ b/tests/testthat/test-fly_georef_gcps.R @@ -0,0 +1,171 @@ +# 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) +}) + + +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