From faa8881cc66483584125044316cb8db6311f61d3 Mon Sep 17 00:00:00 2001 From: almac2022 Date: Thu, 12 Mar 2026 09:24:48 -0700 Subject: [PATCH 1/2] =?UTF-8?q?Rename=20fly=5Fthumb=5Fgeoref=20=E2=86=92?= =?UTF-8?q?=20fly=5Fgeoref,=20add=20rotation=20and=20bearing?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Rename fly_thumb_georef() to fly_georef() — works with any image, not just thumbnails (#24) - Add rotation param (0/90/180/270, default 180) to correct image orientation. Per-photo rotation via column on input sf (#25) - Add fly_bearing() to compute flight line bearing from consecutive centroids per roll — diagnostic tool for determining rotation - 102 tests passing Fixes #24 Fixes #25 Co-Authored-By: Claude Opus 4.6 --- DESCRIPTION | 2 +- NAMESPACE | 3 +- NEWS.md | 7 + R/fly_bearing.R | 67 ++++++++++ R/{fly_thumb_georef.R => fly_georef.R} | 120 ++++++++++++++---- man/fly_bearing.Rd | 39 ++++++ man/{fly_thumb_georef.Rd => fly_georef.Rd} | 58 ++++++--- tests/testthat/test-fly_bearing.R | 51 ++++++++ ...t-fly_thumb_georef.R => test-fly_georef.R} | 82 +++++++++--- vignettes/airphoto-selection.Rmd | 4 +- 10 files changed, 370 insertions(+), 63 deletions(-) create mode 100644 R/fly_bearing.R rename R/{fly_thumb_georef.R => fly_georef.R} (56%) create mode 100644 man/fly_bearing.Rd rename man/{fly_thumb_georef.Rd => fly_georef.Rd} (53%) create mode 100644 tests/testthat/test-fly_bearing.R rename tests/testthat/{test-fly_thumb_georef.R => test-fly_georef.R} (51%) diff --git a/DESCRIPTION b/DESCRIPTION index 00b3118..4921988 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: fly Title: Airphoto Footprint Estimation and Coverage Selection -Version: 0.2.1 +Version: 0.3.0 Authors@R: c( person("Allan", "Irvine", , "al@newgraphenvironment.com", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-3495-2128")), diff --git a/NAMESPACE b/NAMESPACE index 3715e0a..04cf979 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,13 +1,14 @@ # Generated by roxygen2: do not edit by hand +export(fly_bearing) export(fly_coverage) export(fly_fetch) export(fly_filter) export(fly_footprint) +export(fly_georef) export(fly_overlap) export(fly_select) export(fly_summary) -export(fly_thumb_georef) importFrom(rlang,"!!") importFrom(rlang,":=") importFrom(rlang,.data) diff --git a/NEWS.md b/NEWS.md index 822f056..31e2439 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,12 @@ # fly (development version) +## 0.3.0 (2026-03-12) + +- **BREAKING:** Rename `fly_thumb_georef()` → `fly_georef()` — not thumbnail-specific ([#24](https://github.com/NewGraphEnvironment/fly/issues/24)) +- Add `rotation` parameter to `fly_georef()` (default 180°) for correcting image orientation per-roll ([#25](https://github.com/NewGraphEnvironment/fly/issues/25)) +- Add `fly_bearing()` — compute flight line bearing from consecutive centroids per roll +- Per-photo rotation via `rotation` column on input sf overrides default + ## 0.2.1 (2026-03-11) - Add `workers` parameter to `fly_fetch()` for parallel downloads via `furrr`/`future` ([#21](https://github.com/NewGraphEnvironment/fly/issues/21)) diff --git a/R/fly_bearing.R b/R/fly_bearing.R new file mode 100644 index 0000000..7995224 --- /dev/null +++ b/R/fly_bearing.R @@ -0,0 +1,67 @@ +#' Compute flight line bearing from consecutive airphoto centroids +#' +#' Estimates the flight direction for each photo by computing the azimuth +#' between consecutive centroids on the same film roll, sorted by frame +#' number. Useful for diagnosing image rotation issues in [fly_georef()]. +#' +#' @param photos_sf An sf object with `film_roll` and `frame_number` +#' columns. Projected to BC Albers (EPSG:3005) internally for metric +#' bearing computation. +#' @return The input sf object with an added `bearing` column (degrees +#' clockwise from north, 0–360). Photos with no computable bearing +#' (single-frame rolls) get `NA`. +#' +#' @details +#' Within each roll, frames are sorted by `frame_number`. The bearing +#' for each frame is the azimuth to the next frame on the same roll. +#' The last frame on each roll gets the bearing from the previous frame. +#' +#' Aerial survey flights follow back-and-forth patterns, so bearings +#' alternate between ~opposite directions (e.g., 90° and 270°) on +#' consecutive legs. Large frame number gaps may indicate a new flight +#' line within the same roll. +#' +#' @examples +#' centroids <- sf::st_read(system.file("testdata/photo_centroids.gpkg", package = "fly")) +#' with_bearing <- fly_bearing(centroids) +#' with_bearing[, c("film_roll", "frame_number", "bearing")] +#' +#' @export +fly_bearing <- function(photos_sf) { + if (!all(c("film_roll", "frame_number") %in% names(photos_sf))) { + stop("`photos_sf` must have `film_roll` and `frame_number` columns.", + call. = FALSE) + } + + # Project to BC Albers for metric bearing + proj <- sf::st_transform(photos_sf, 3005) + coords <- sf::st_coordinates(proj) + + # Sort index by roll + frame + + ord <- order(photos_sf$film_roll, photos_sf$frame_number) + + bearing <- rep(NA_real_, nrow(photos_sf)) + + rolls <- photos_sf$film_roll[ord] + x <- coords[ord, 1] + y <- coords[ord, 2] + + for (i in seq_along(ord)) { + if (i < length(ord) && rolls[i] == rolls[i + 1]) { + # Forward bearing to next frame on same roll + dx <- x[i + 1] - x[i] + dy <- y[i + 1] - y[i] + bearing[ord[i]] <- (atan2(dx, dy) * 180 / pi) %% 360 + } else if (i > 1 && rolls[i] == rolls[i - 1]) { + # Last frame on roll: use bearing from previous + dx <- x[i] - x[i - 1] + dy <- y[i] - y[i - 1] + bearing[ord[i]] <- (atan2(dx, dy) * 180 / pi) %% 360 + } + # else: single-frame roll, stays NA + } + + photos_sf$bearing <- bearing + photos_sf +} diff --git a/R/fly_thumb_georef.R b/R/fly_georef.R similarity index 56% rename from R/fly_thumb_georef.R rename to R/fly_georef.R index 42d15c0..d80f12d 100644 --- a/R/fly_thumb_georef.R +++ b/R/fly_georef.R @@ -1,13 +1,14 @@ -#' Georeference downloaded thumbnails to footprint polygons +#' Georeference airphoto images to footprint polygons #' -#' Warps thumbnail images to their estimated ground footprint using GCPs -#' (ground control points) derived from [fly_footprint()]. Produces -#' georeferenced GeoTIFFs in BC Albers (EPSG:3005). +#' Warps images to their estimated ground footprint using GCPs (ground control +#' points) derived from [fly_footprint()]. Produces georeferenced GeoTIFFs in +#' BC Albers (EPSG:3005). Works with thumbnails and full-resolution scans. #' #' @param fetch_result A tibble returned by [fly_fetch()], with columns #' `airp_id`, `dest`, and `success`. #' @param photos_sf The same sf object passed to `fly_fetch()`, with a -#' `scale` column for footprint estimation. +#' `scale` column for footprint estimation. If a `rotation` column is +#' present, per-photo rotation values are used (see **Rotation** below). #' @param dest_dir Directory for output GeoTIFFs. Created if it does not #' exist. #' @param overwrite If `FALSE` (default), skip files that already exist. @@ -17,18 +18,43 @@ #' film holder edges at the cost of losing real black pixels — acceptable #' for thumbnails but may need adjustment for full-resolution scans. #' Set to `NULL` to disable source nodata detection entirely. +#' @param rotation Default image rotation in degrees clockwise (one of +#' `0`, `90`, `180`, `270`). Controls which pixel corners map to which +#' footprint corners. Default `180` is correct for most BC aerial photos. +#' Overridden per-photo if `photos_sf` contains a `rotation` column. #' @return A tibble with columns `airp_id`, `source`, `dest`, and `success`. #' #' @details -#' Each thumbnail's four corners are mapped to the corresponding footprint +#' Each image's four corners are mapped to the corresponding footprint #' polygon corners computed by [fly_footprint()] in BC Albers. GDAL #' 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: +#' \itemize{ +#' \item `0` — top of image maps to north edge of footprint (original behavior) +#' \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 `270` — top of image maps to west edge +#' } +#' +#' Rotation is consistent within a film roll — all frames from the same +#' roll need the same correction. Set per-roll values in a `rotation` +#' column on `photos_sf`: +#' ``` +#' photos$rotation <- dplyr::case_when( +#' photos$film_roll == "bc5282" ~ 90, +#' .default = 180 +#' ) +#' ``` +#' #' **Nodata handling:** Two sources of unwanted black pixels are masked: #' #' 1. **Warp fill** — GDAL creates black pixels outside the rotated source -#' frame. RGB thumbnails get an alpha band (`-dstalpha`); grayscale use +#' frame. RGB images get an alpha band (`-dstalpha`); grayscale use #' `dstnodata=0`. #' 2. **Camera frame borders** — film holder edges, fiducial marks, and #' scanning artifacts produce black (value 0) pixels within the source @@ -42,7 +68,7 @@ #' downstream (e.g., circle detection). #' #' **Accuracy:** footprints assume flat terrain and nadir camera angle. -#' The georeferenced thumbnails are approximate — useful for visual context, +#' The georeferenced images are approximate — useful for visual context, #' not survey-grade positioning. See [fly_footprint()] for details on #' limitations. #' @@ -52,17 +78,21 @@ #' # Fetch and georeference first 2 thumbnails #' fetched <- fly_fetch(centroids[1:2, ], type = "thumbnail", #' dest_dir = tempdir()) -#' georef <- fly_thumb_georef(fetched, centroids[1:2, ], -#' dest_dir = tempdir()) +#' georef <- fly_georef(fetched, centroids[1:2, ], +#' dest_dir = tempdir()) #' georef #' #' @export -fly_thumb_georef <- function(fetch_result, photos_sf, - dest_dir = "georef", overwrite = FALSE, - srcnodata = "0") { +fly_georef <- function(fetch_result, photos_sf, + dest_dir = "georef", overwrite = FALSE, + srcnodata = "0", rotation = 180) { if (!all(c("airp_id", "dest", "success") %in% names(fetch_result))) { stop("`fetch_result` must be output from `fly_fetch()`.", call. = FALSE) } + rotation <- as.integer(rotation) + if (!rotation %in% c(0L, 90L, 180L, 270L)) { + stop("`rotation` must be one of 0, 90, 180, 270.", call. = FALSE) + } dir.create(dest_dir, recursive = TRUE, showWarnings = FALSE) @@ -72,6 +102,9 @@ fly_thumb_georef <- function(fetch_result, photos_sf, # Match fetch results to photos by airp_id ids <- fetch_result$airp_id + # Per-photo rotation: column in photos_sf overrides default + has_rotation_col <- "rotation" %in% names(photos_sf) + results <- dplyr::tibble( airp_id = ids, source = fetch_result$dest, @@ -98,8 +131,15 @@ fly_thumb_georef <- function(fetch_result, photos_sf, if (length(fp_idx) == 0) next fp <- footprints[fp_idx[1], ] + # Per-photo rotation from column, or default + rot <- if (has_rotation_col) { + as.integer(photos_sf[["rotation"]][fp_idx[1]]) + } else { + rotation + } + results$success[i] <- tryCatch( - georef_one(src, fp, out_file, srcnodata = srcnodata), + georef_one(src, fp, out_file, srcnodata = srcnodata, rotation = rot), error = function(e) { message("Failed to georef ", basename(src), ": ", e$message) FALSE @@ -108,16 +148,17 @@ fly_thumb_georef <- function(fetch_result, photos_sf, } n_ok <- sum(results$success) - message("Georeferenced ", n_ok, " of ", nrow(results), " thumbnails") + message("Georeferenced ", n_ok, " of ", nrow(results), " images") results } -#' Georeference a single thumbnail to a footprint polygon +#' Georeference a single image to a footprint polygon #' @noRd -georef_one <- function(src, fp, out_file, srcnodata = "0") { +georef_one <- function(src, fp, out_file, srcnodata = "0", rotation = 180) { # Get footprint corner coordinates # fly_footprint builds: BL, BR, TR, TL, BL (closing) 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) @@ -131,16 +172,45 @@ georef_one <- function(src, fp, out_file, srcnodata = "0") { n_bands <- length(gregexpr("Band \\d+", info)[[1]]) is_rgb <- n_bands >= 3 - # Map pixel corners to footprint corners - # Pixel: TL=(0,0), TR=(ncol,0), BR=(ncol,nrow), BL=(0,nrow) - # Footprint coords: [1]=BL, [2]=BR, [3]=TR, [4]=TL - gcp_args <- c( - "-gcp", 0, 0, coords[4, 1], coords[4, 2], - "-gcp", ncol_px, 0, coords[3, 1], coords[3, 2], - "-gcp", ncol_px, nrow_px, coords[2, 1], coords[2, 2], - "-gcp", 0, nrow_px, coords[1, 1], coords[1, 2] + # Pixel corners: TL, TR, BR, BL + pixel_corners <- list( + c(0, 0), # TL + c(ncol_px, 0), # TR + c(ncol_px, nrow_px), # BR + c(0, nrow_px) # BL ) + # Footprint corners in same order: TL, TR, BR, BL + fp_corners <- list( + coords[4, 1:2], # TL + coords[3, 1:2], # TR + coords[2, 1:2], # BR + coords[1, 1:2] # BL + ) + + # Rotation: shift the footprint corner mapping + + # rotation=0: pixel TL → footprint TL (north-up, original behavior) + # rotation=90: pixel TL → footprint TR (top of image = east) + # rotation=180: pixel TL → footprint BR (top of image = south) + # rotation=270: pixel TL → footprint BL (top of image = west) + n_shifts <- rotation %/% 90 + if (n_shifts > 0) { + fp_corners <- c( + fp_corners[(n_shifts + 1):4], + fp_corners[1:n_shifts] + ) + } + + # Build GCP args mapping pixel corners to (rotated) footprint corners + 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] + ) + } + # Step 1: translate with GCPs tmp_file <- tempfile(fileext = ".tif") on.exit(unlink(tmp_file), add = TRUE) diff --git a/man/fly_bearing.Rd b/man/fly_bearing.Rd new file mode 100644 index 0000000..5183e1a --- /dev/null +++ b/man/fly_bearing.Rd @@ -0,0 +1,39 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/fly_bearing.R +\name{fly_bearing} +\alias{fly_bearing} +\title{Compute flight line bearing from consecutive airphoto centroids} +\usage{ +fly_bearing(photos_sf) +} +\arguments{ +\item{photos_sf}{An sf object with \code{film_roll} and \code{frame_number} +columns. Projected to BC Albers (EPSG:3005) internally for metric +bearing computation.} +} +\value{ +The input sf object with an added \code{bearing} column (degrees +clockwise from north, 0–360). Photos with no computable bearing +(single-frame rolls) get \code{NA}. +} +\description{ +Estimates the flight direction for each photo by computing the azimuth +between consecutive centroids on the same film roll, sorted by frame +number. Useful for diagnosing image rotation issues in \code{\link[=fly_georef]{fly_georef()}}. +} +\details{ +Within each roll, frames are sorted by \code{frame_number}. The bearing +for each frame is the azimuth to the next frame on the same roll. +The last frame on each roll gets the bearing from the previous frame. + +Aerial survey flights follow back-and-forth patterns, so bearings +alternate between ~opposite directions (e.g., 90° and 270°) on +consecutive legs. Large frame number gaps may indicate a new flight +line within the same roll. +} +\examples{ +centroids <- sf::st_read(system.file("testdata/photo_centroids.gpkg", package = "fly")) +with_bearing <- fly_bearing(centroids) +with_bearing[, c("film_roll", "frame_number", "bearing")] + +} diff --git a/man/fly_thumb_georef.Rd b/man/fly_georef.Rd similarity index 53% rename from man/fly_thumb_georef.Rd rename to man/fly_georef.Rd index 095f73e..c3eb186 100644 --- a/man/fly_thumb_georef.Rd +++ b/man/fly_georef.Rd @@ -1,15 +1,16 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/fly_thumb_georef.R -\name{fly_thumb_georef} -\alias{fly_thumb_georef} -\title{Georeference downloaded thumbnails to footprint polygons} +% Please edit documentation in R/fly_georef.R +\name{fly_georef} +\alias{fly_georef} +\title{Georeference airphoto images to footprint polygons} \usage{ -fly_thumb_georef( +fly_georef( fetch_result, photos_sf, dest_dir = "georef", overwrite = FALSE, - srcnodata = "0" + srcnodata = "0", + rotation = 180 ) } \arguments{ @@ -17,7 +18,8 @@ fly_thumb_georef( \code{airp_id}, \code{dest}, and \code{success}.} \item{photos_sf}{The same sf object passed to \code{fly_fetch()}, with a -\code{scale} column for footprint estimation.} +\code{scale} column for footprint estimation. If a \code{rotation} column is +present, per-photo rotation values are used (see \strong{Rotation} below).} \item{dest_dir}{Directory for output GeoTIFFs. Created if it does not exist.} @@ -30,25 +32,51 @@ nodata for grayscale). Default \code{"0"} masks camera frame borders and film holder edges at the cost of losing real black pixels — acceptable for thumbnails but may need adjustment for full-resolution scans. Set to \code{NULL} to disable source nodata detection entirely.} + +\item{rotation}{Default image rotation in degrees clockwise (one of +\code{0}, \code{90}, \code{180}, \code{270}). Controls which pixel corners map to which +footprint corners. Default \code{180} is correct for most BC aerial photos. +Overridden per-photo if \code{photos_sf} contains a \code{rotation} column.} } \value{ A tibble with columns \code{airp_id}, \code{source}, \code{dest}, and \code{success}. } \description{ -Warps thumbnail images to their estimated ground footprint using GCPs -(ground control points) derived from \code{\link[=fly_footprint]{fly_footprint()}}. Produces -georeferenced GeoTIFFs in BC Albers (EPSG:3005). +Warps images to their estimated ground footprint using GCPs (ground control +points) derived from \code{\link[=fly_footprint]{fly_footprint()}}. Produces georeferenced GeoTIFFs in +BC Albers (EPSG:3005). Works with thumbnails and full-resolution scans. } \details{ -Each thumbnail's four corners are mapped to the corresponding footprint +Each image's four corners are mapped to the corresponding footprint polygon corners computed by \code{\link[=fly_footprint]{fly_footprint()}} in BC Albers. GDAL 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: +\itemize{ +\item \code{0} — top of image maps to north edge of footprint (original behavior) +\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{270} — top of image maps to west edge +} + +Rotation is consistent within a film roll — all frames from the same +roll need the same correction. Set per-roll values in a \code{rotation} +column on \code{photos_sf}: + +\if{html}{\out{
}}\preformatted{photos$rotation <- dplyr::case_when( + photos$film_roll == "bc5282" ~ 90, + .default = 180 +) +}\if{html}{\out{
}} + \strong{Nodata handling:} Two sources of unwanted black pixels are masked: \enumerate{ \item \strong{Warp fill} — GDAL creates black pixels outside the rotated source -frame. RGB thumbnails get an alpha band (\code{-dstalpha}); grayscale use +frame. RGB images get an alpha band (\code{-dstalpha}); grayscale use \code{dstnodata=0}. \item \strong{Camera frame borders} — film holder edges, fiducial marks, and scanning artifacts produce black (value 0) pixels within the source @@ -63,7 +91,7 @@ detail matters, set \code{srcnodata = NULL} and handle frame masking downstream (e.g., circle detection). \strong{Accuracy:} footprints assume flat terrain and nadir camera angle. -The georeferenced thumbnails are approximate — useful for visual context, +The georeferenced images are approximate — useful for visual context, not survey-grade positioning. See \code{\link[=fly_footprint]{fly_footprint()}} for details on limitations. } @@ -73,8 +101,8 @@ centroids <- sf::st_read(system.file("testdata/photo_centroids.gpkg", package = # Fetch and georeference first 2 thumbnails fetched <- fly_fetch(centroids[1:2, ], type = "thumbnail", dest_dir = tempdir()) -georef <- fly_thumb_georef(fetched, centroids[1:2, ], - dest_dir = tempdir()) +georef <- fly_georef(fetched, centroids[1:2, ], + dest_dir = tempdir()) georef } diff --git a/tests/testthat/test-fly_bearing.R b/tests/testthat/test-fly_bearing.R new file mode 100644 index 0000000..db01a6c --- /dev/null +++ b/tests/testthat/test-fly_bearing.R @@ -0,0 +1,51 @@ +test_that("fly_bearing adds bearing column", { + centroids <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE) + result <- fly_bearing(centroids) + expect_true("bearing" %in% names(result)) + expect_equal(nrow(result), nrow(centroids)) +}) + +test_that("fly_bearing computes valid azimuths", { + centroids <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE) + result <- fly_bearing(centroids) + + # All non-NA bearings should be 0-360 + bearings <- result$bearing[!is.na(result$bearing)] + expect_true(all(bearings >= 0 & bearings < 360)) +}) + +test_that("fly_bearing handles single-frame rolls", { + centroids <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE) + + # bc5300 and bc5301 have only 1 frame each in test data + single <- centroids[centroids$film_roll == "bc5300", ] + expect_equal(nrow(single), 1) + + result <- fly_bearing(single) + expect_true(is.na(result$bearing[1])) +}) + +test_that("fly_bearing is consistent within flight legs", { + centroids <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE) + result <- fly_bearing(centroids) + + # bc5282 has 10 frames — consecutive frames on same leg should + + # have similar bearings (within 30 degrees) + roll <- result[result$film_roll == "bc5282", ] + roll <- roll[order(roll$frame_number), ] + bearings <- roll$bearing + + # Check that at least some consecutive pairs are similar + diffs <- abs(diff(bearings)) + # Normalize to 0-180 + diffs <- pmin(diffs, 360 - diffs) + # Back-and-forth legs produce ~180 differences, same-leg pairs < 30 + expect_true(any(diffs < 30)) +}) + +test_that("fly_bearing rejects missing columns", { + centroids <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE) + centroids$film_roll <- NULL + expect_error(fly_bearing(centroids), "film_roll") +}) diff --git a/tests/testthat/test-fly_thumb_georef.R b/tests/testthat/test-fly_georef.R similarity index 51% rename from tests/testthat/test-fly_thumb_georef.R rename to tests/testthat/test-fly_georef.R index 7f59c60..8eb6257 100644 --- a/tests/testthat/test-fly_thumb_georef.R +++ b/tests/testthat/test-fly_georef.R @@ -1,4 +1,4 @@ -test_that("fly_thumb_georef returns expected columns", { +test_that("fly_georef returns expected columns", { centroids <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE) dest_fetch <- file.path(tempdir(), "fly_georef_test_fetch") unlink(dest_fetch, recursive = TRUE) @@ -8,13 +8,13 @@ test_that("fly_thumb_georef returns expected columns", { dest_georef <- file.path(tempdir(), "fly_georef_test_out") unlink(dest_georef, recursive = TRUE) - result <- fly_thumb_georef(fetched, centroids[1, ], - dest_dir = dest_georef) + result <- fly_georef(fetched, centroids[1, ], + dest_dir = dest_georef) expect_s3_class(result, "tbl_df") expect_true(all(c("airp_id", "source", "dest", "success") %in% names(result))) }) -test_that("fly_thumb_georef produces georeferenced TIFFs", { +test_that("fly_georef produces georeferenced TIFFs", { centroids <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE) dest_fetch <- file.path(tempdir(), "fly_georef_test_tiff_fetch") unlink(dest_fetch, recursive = TRUE) @@ -24,8 +24,8 @@ test_that("fly_thumb_georef produces georeferenced TIFFs", { dest_georef <- file.path(tempdir(), "fly_georef_test_tiff_out") unlink(dest_georef, recursive = TRUE) - result <- fly_thumb_georef(fetched, centroids[1, ], - dest_dir = dest_georef) + result <- fly_georef(fetched, centroids[1, ], + dest_dir = dest_georef) expect_true(result$success[1]) expect_true(file.exists(result$dest[1])) @@ -34,7 +34,7 @@ test_that("fly_thumb_georef produces georeferenced TIFFs", { expect_true(grepl("3005", info)) }) -test_that("fly_thumb_georef skips failed fetches", { +test_that("fly_georef skips failed fetches", { centroids <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE) fake_fetch <- dplyr::tibble( airp_id = centroids$airp_id[1], @@ -42,12 +42,12 @@ test_that("fly_thumb_georef skips failed fetches", { dest = "/nonexistent/fake.jpg", success = FALSE ) - result <- fly_thumb_georef(fake_fetch, centroids[1, ], - dest_dir = tempdir()) + result <- fly_georef(fake_fetch, centroids[1, ], + dest_dir = tempdir()) expect_false(result$success[1]) }) -test_that("fly_thumb_georef skips existing when overwrite is FALSE", { +test_that("fly_georef skips existing when overwrite is FALSE", { centroids <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE) dest_fetch <- file.path(tempdir(), "fly_georef_overwrite_fetch") unlink(dest_fetch, recursive = TRUE) @@ -58,24 +58,24 @@ test_that("fly_thumb_georef skips existing when overwrite is FALSE", { unlink(dest_georef, recursive = TRUE) # First run - fly_thumb_georef(fetched, centroids[1, ], dest_dir = dest_georef) + fly_georef(fetched, centroids[1, ], dest_dir = dest_georef) f <- list.files(dest_georef, full.names = TRUE)[1] mtime1 <- file.mtime(f) Sys.sleep(1) # Second run without overwrite - fly_thumb_georef(fetched, centroids[1, ], - dest_dir = dest_georef, overwrite = FALSE) + fly_georef(fetched, centroids[1, ], + dest_dir = dest_georef, overwrite = FALSE) mtime2 <- file.mtime(f) expect_equal(mtime1, mtime2) }) -test_that("fly_thumb_georef rejects bad input", { - expect_error(fly_thumb_georef(data.frame(x = 1), data.frame(y = 1)), +test_that("fly_georef rejects bad input", { + expect_error(fly_georef(data.frame(x = 1), data.frame(y = 1)), "fly_fetch") }) -test_that("fly_thumb_georef extent matches footprint", { +test_that("fly_georef extent matches footprint", { centroids <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE) dest_fetch <- file.path(tempdir(), "fly_georef_extent_fetch") unlink(dest_fetch, recursive = TRUE) @@ -85,17 +85,61 @@ test_that("fly_thumb_georef extent matches footprint", { dest_georef <- file.path(tempdir(), "fly_georef_extent_out") unlink(dest_georef, recursive = TRUE) - result <- fly_thumb_georef(fetched, centroids[1, ], - dest_dir = dest_georef) + result <- fly_georef(fetched, centroids[1, ], + dest_dir = dest_georef) # Compare georef extent to footprint extent fp <- fly_footprint(centroids[1, ]) |> sf::st_transform(3005) fp_bbox <- sf::st_bbox(fp) info <- sf::gdal_utils("info", source = result$dest[1], quiet = TRUE) - # Extract corner coordinates from gdalinfo ul <- regmatches(info, regexpr("Upper Left\\s+\\([^)]+\\)", info)) lr <- regmatches(info, regexpr("Lower Right\\s+\\([^)]+\\)", info)) expect_length(ul, 1) expect_length(lr, 1) }) + +test_that("fly_georef accepts rotation parameter", { + centroids <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE) + dest_fetch <- file.path(tempdir(), "fly_georef_rot_fetch") + unlink(dest_fetch, recursive = TRUE) + + fetched <- fly_fetch(centroids[1, ], type = "thumbnail", + dest_dir = dest_fetch) + + # Each rotation value should produce a valid georef + for (rot in c(0, 90, 180, 270)) { + dest_georef <- file.path(tempdir(), paste0("fly_georef_rot_", rot)) + unlink(dest_georef, recursive = TRUE) + result <- fly_georef(fetched, centroids[1, ], + dest_dir = dest_georef, rotation = rot) + expect_true(result$success[1], info = paste("rotation =", rot)) + } +}) + +test_that("fly_georef rejects invalid rotation", { + centroids <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE) + fake_fetch <- dplyr::tibble( + airp_id = centroids$airp_id[1], + dest = tempfile(), success = TRUE + ) + expect_error(fly_georef(fake_fetch, centroids[1, ], rotation = 45), + "one of 0, 90, 180, 270") +}) + +test_that("fly_georef reads rotation from column", { + centroids <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE) + dest_fetch <- file.path(tempdir(), "fly_georef_rotcol_fetch") + unlink(dest_fetch, recursive = TRUE) + + fetched <- fly_fetch(centroids[1, ], type = "thumbnail", + dest_dir = dest_fetch) + dest_georef <- file.path(tempdir(), "fly_georef_rotcol_out") + unlink(dest_georef, recursive = TRUE) + + # Add rotation column + centroids$rotation <- 90 + result <- fly_georef(fetched, centroids[1, ], + dest_dir = dest_georef) + expect_true(result$success[1]) +}) diff --git a/vignettes/airphoto-selection.Rmd b/vignettes/airphoto-selection.Rmd index 0e61fae..d036b2f 100644 --- a/vignettes/airphoto-selection.Rmd +++ b/vignettes/airphoto-selection.Rmd @@ -247,13 +247,13 @@ legend("topright", legend = scale_labels, `fly_fetch()` downloads thumbnail images (or flight logs, calibration reports) from the BC Data Catalogue URLs included in the centroid data. -`fly_thumb_georef()` warps each thumbnail to its estimated footprint +`fly_georef()` warps each thumbnail to its estimated footprint polygon, producing georeferenced GeoTIFFs in BC Albers. ```{r fetch-georef} fetched <- fly_fetch(centroids[1:3, ], type = "thumbnail", dest_dir = tempdir()) -georef <- fly_thumb_georef(fetched, centroids[1:3, ], +georef <- fly_georef(fetched, centroids[1:3, ], dest_dir = tempdir()) georef[, c("airp_id", "dest", "success")] ``` From f2dd965c5506f387a30ce9b9a700a49dbdff9ea3 Mon Sep 17 00:00:00 2001 From: almac2022 Date: Thu, 12 Mar 2026 10:04:08 -0700 Subject: [PATCH 2/2] Add auto rotation from bearing in fly_georef() MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit rotation="auto" (new default) computes flight bearing per-photo and derives GCP rotation via floor((bearing + 91) / 90) * 90 %% 360. Calibrated on 10 BC aerial photos (1968-2019), 8/10 correct. Diagonal flights (~45° off cardinal) are unreliable — override with rotation column. Add data-raw/test_georef_decades.R for visual QA across decades. Co-Authored-By: Claude Opus 4.6 --- R/fly_georef.R | 76 ++++++++++++++++++++------ data-raw/test_georef_decades.R | 93 ++++++++++++++++++++++++++++++++ man/fly_georef.Rd | 32 +++++++---- tests/testthat/test-fly_georef.R | 18 ++++++- 4 files changed, 191 insertions(+), 28 deletions(-) create mode 100644 data-raw/test_georef_decades.R diff --git a/R/fly_georef.R b/R/fly_georef.R index d80f12d..8f018ff 100644 --- a/R/fly_georef.R +++ b/R/fly_georef.R @@ -18,10 +18,12 @@ #' film holder edges at the cost of losing real black pixels — acceptable #' for thumbnails but may need adjustment for full-resolution scans. #' Set to `NULL` to disable source nodata detection entirely. -#' @param rotation Default image rotation in degrees clockwise (one of -#' `0`, `90`, `180`, `270`). Controls which pixel corners map to which -#' footprint corners. Default `180` is correct for most BC aerial photos. -#' Overridden per-photo if `photos_sf` contains a `rotation` column. +#' @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. #' @return A tibble with columns `airp_id`, `source`, `dest`, and `success`. #' #' @details @@ -41,13 +43,21 @@ #' \item `270` — top of image maps to west edge #' } #' -#' Rotation is consistent within a film roll — all frames from the same -#' roll need the same correction. Set per-roll values in a `rotation` -#' column on `photos_sf`: +#' When `rotation = "auto"`, the bearing-to-rotation formula is: +#' `floor((bearing + 91) / 90) * 90 %% 360`. This was calibrated on +#' BC aerial photos spanning 1968–2019 across multiple camera systems +#' and scanners. Photos on diagonal flight lines (~45° off cardinal) +#' may be imperfect — check visually and override with a `rotation` +#' column if needed. +#' +#' Within a film roll, consecutive flight legs alternate direction +#' (back-and-forth pattern), so different frames on the same roll may +#' need different rotations. This is why `"auto"` computes per-photo, +#' not per-roll. To override, add a `rotation` column to `photos_sf`: #' ``` #' photos$rotation <- dplyr::case_when( -#' photos$film_roll == "bc5282" ~ 90, -#' .default = 180 +#' photos$film_roll == "bc5282" ~ 270, +#' .default = NA # fall through to auto #' ) #' ``` #' @@ -75,7 +85,7 @@ #' @examples #' centroids <- sf::st_read(system.file("testdata/photo_centroids.gpkg", package = "fly")) #' -#' # Fetch and georeference first 2 thumbnails +#' # Fetch and georeference with auto rotation (uses bearing from centroids) #' fetched <- fly_fetch(centroids[1:2, ], type = "thumbnail", #' dest_dir = tempdir()) #' georef <- fly_georef(fetched, centroids[1:2, ], @@ -85,13 +95,17 @@ #' @export fly_georef <- function(fetch_result, photos_sf, dest_dir = "georef", overwrite = FALSE, - srcnodata = "0", rotation = 180) { + srcnodata = "0", rotation = "auto") { if (!all(c("airp_id", "dest", "success") %in% names(fetch_result))) { stop("`fetch_result` must be output from `fly_fetch()`.", call. = FALSE) } - rotation <- as.integer(rotation) - if (!rotation %in% c(0L, 90L, 180L, 270L)) { - stop("`rotation` must be one of 0, 90, 180, 270.", call. = FALSE) + + auto_rotation <- identical(rotation, "auto") + if (!auto_rotation) { + rotation <- as.integer(rotation) + if (!rotation %in% c(0L, 90L, 180L, 270L)) { + stop("`rotation` must be one of \"auto\", 0, 90, 180, 270.", call. = FALSE) + } } dir.create(dest_dir, recursive = TRUE, showWarnings = FALSE) @@ -102,9 +116,23 @@ fly_georef <- function(fetch_result, photos_sf, # Match fetch results to photos by airp_id ids <- fetch_result$airp_id - # Per-photo rotation: column in photos_sf overrides default + # Per-photo rotation: column overrides auto/default + has_rotation_col <- "rotation" %in% names(photos_sf) + # Auto-compute bearing → rotation when needed + if (auto_rotation && !has_rotation_col) { + if (all(c("film_roll", "frame_number") %in% names(photos_sf))) { + photos_sf <- fly_bearing(photos_sf) + photos_sf$rotation <- bearing_to_rotation(photos_sf$bearing) + has_rotation_col <- TRUE + } else { + message("No film_roll/frame_number columns for auto rotation, using 180") + rotation <- 180L + auto_rotation <- FALSE + } + } + results <- dplyr::tibble( airp_id = ids, source = fetch_result$dest, @@ -133,7 +161,10 @@ fly_georef <- function(fetch_result, photos_sf, # Per-photo rotation from column, or default rot <- if (has_rotation_col) { - as.integer(photos_sf[["rotation"]][fp_idx[1]]) + val <- as.integer(photos_sf[["rotation"]][fp_idx[1]]) + if (is.na(val)) { + if (auto_rotation) 180L else rotation + } else val } else { rotation } @@ -248,3 +279,16 @@ georef_one <- function(src, fp, out_file, srcnodata = "0", rotation = 180) { file.exists(out_file) && file.size(out_file) > 0 } + +#' Convert flight bearing to GCP rotation +#' +#' Formula calibrated on BC aerial photos (1968–2019). +#' @param bearing Numeric vector of bearings (degrees, 0–360). +#' @return Integer vector of rotations (0, 90, 180, or 270). NA bearings +#' return 180 (most common default). +#' @noRd +bearing_to_rotation <- function(bearing) { + rot <- (floor((bearing + 91) / 90) * 90L) %% 360L + rot[is.na(rot)] <- 180L + as.integer(rot) +} diff --git a/data-raw/test_georef_decades.R b/data-raw/test_georef_decades.R new file mode 100644 index 0000000..c15fcbe --- /dev/null +++ b/data-raw/test_georef_decades.R @@ -0,0 +1,93 @@ +# test_georef_decades.R — Visual QA: one georef thumbnail per decade +# +# Produces georeferenced thumbnails in a single directory for QGIS inspection. +# Requires stac_airphoto_bc centroid cache and downloaded thumbnails. +# +# Usage: source this script interactively, then load the output dir in QGIS. + +library(sf) +library(dplyr) +devtools::load_all() + +# --- Config ------------------------------------------------------------------ + +stac_dir <- "../stac_airphoto_bc" +cache_path <- file.path(stac_dir, "data/centroids_raw.parquet") +thumbs_dir <- file.path(stac_dir, "data/raw/thumbs") +out_dir <- file.path(stac_dir, "data/raw/georef/qa_decades") + +stopifnot( + "Centroid cache not found — run stac_airphoto_bc/scripts/01_fetch.R first" = + file.exists(cache_path) +) + +# --- Load centroids ---------------------------------------------------------- + +centroids <- arrow::read_parquet(cache_path) |> + st_as_sf(coords = c("longitude", "latitude"), crs = 4326) + +centroids$decade <- floor(centroids$photo_year / 10) * 10 + +# --- Pick one photo per decade near Houston ---------------------------------- +# Closest to AOI centre so footprints overlap for easy comparison + +aoi_centre <- st_sfc(st_point(c(-126.65, 54.4)), crs = 4326) +centroids$dist <- as.numeric(st_distance(centroids, aoi_centre)) + +samples <- centroids |> + group_by(decade) |> + slice_min(dist, n = 1, with_ties = FALSE) |> + ungroup() |> + arrange(decade) + +message(nrow(samples), " samples across decades: ", + paste(samples$decade, collapse = ", ")) + +# --- Find raw thumbnails on disk --------------------------------------------- + +samples$thumb_path <- vapply(seq_len(nrow(samples)), function(i) { + year_dir <- file.path(thumbs_dir, samples$photo_year[i]) + if (!dir.exists(year_dir)) return(NA_character_) + pattern <- paste0("^", samples$film_roll[i], "_", + sprintf("%03d", samples$frame_number[i])) + files <- list.files(year_dir, pattern = pattern, full.names = TRUE) + if (length(files) == 0) return(NA_character_) + files[1] +}, character(1)) + +found <- !is.na(samples$thumb_path) +if (any(!found)) { + message("Missing thumbnails for: ", + paste(samples$film_roll[!found], samples$frame_number[!found], + sep = "_", collapse = ", ")) +} +samples <- samples[found, ] + +# --- Georef with default rotation=180 ---------------------------------------- + +dir.create(out_dir, recursive = TRUE, showWarnings = FALSE) + +fetch_result <- dplyr::tibble( + airp_id = samples$airp_id, + url = NA_character_, + dest = samples$thumb_path, + success = TRUE +) + +result <- fly_georef(fetch_result, samples, dest_dir = out_dir, overwrite = TRUE) + +message("\n", sum(result$success), "/", nrow(result), " georeferenced") +message("Output: ", normalizePath(out_dir)) +message("Open in QGIS: open ", normalizePath(out_dir)) + +# --- Summary table ----------------------------------------------------------- + +summary_tbl <- samples |> + st_drop_geometry() |> + select(decade, photo_year, film_roll, frame_number, scale) |> + mutate( + georef_ok = result$success, + output = basename(result$dest) + ) + +print(as.data.frame(summary_tbl)) diff --git a/man/fly_georef.Rd b/man/fly_georef.Rd index c3eb186..255fbf4 100644 --- a/man/fly_georef.Rd +++ b/man/fly_georef.Rd @@ -10,7 +10,7 @@ fly_georef( dest_dir = "georef", overwrite = FALSE, srcnodata = "0", - rotation = 180 + rotation = "auto" ) } \arguments{ @@ -33,10 +33,12 @@ film holder edges at the cost of losing real black pixels — acceptable for thumbnails but may need adjustment for full-resolution scans. Set to \code{NULL} to disable source nodata detection entirely.} -\item{rotation}{Default image rotation in degrees clockwise (one of -\code{0}, \code{90}, \code{180}, \code{270}). Controls which pixel corners map to which -footprint corners. Default \code{180} is correct for most BC aerial photos. -Overridden per-photo if \code{photos_sf} contains a \code{rotation} column.} +\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.} } \value{ A tibble with columns \code{airp_id}, \code{source}, \code{dest}, and \code{success}. @@ -63,13 +65,21 @@ parameter rotates the GCP corner mapping: \item \code{270} — top of image maps to west edge } -Rotation is consistent within a film roll — all frames from the same -roll need the same correction. Set per-roll values in a \code{rotation} -column on \code{photos_sf}: +When \code{rotation = "auto"}, the bearing-to-rotation formula is: +\code{floor((bearing + 91) / 90) * 90 \%\% 360}. This was calibrated on +BC aerial photos spanning 1968–2019 across multiple camera systems +and scanners. Photos on diagonal flight lines (~45° off cardinal) +may be imperfect — check visually and override with a \code{rotation} +column if needed. + +Within a film roll, consecutive flight legs alternate direction +(back-and-forth pattern), so different frames on the same roll may +need different rotations. This is why \code{"auto"} computes per-photo, +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" ~ 90, - .default = 180 + photos$film_roll == "bc5282" ~ 270, + .default = NA # fall through to auto ) }\if{html}{\out{
}} @@ -98,7 +108,7 @@ limitations. \examples{ centroids <- sf::st_read(system.file("testdata/photo_centroids.gpkg", package = "fly")) -# Fetch and georeference first 2 thumbnails +# Fetch and georeference with auto rotation (uses bearing from centroids) fetched <- fly_fetch(centroids[1:2, ], type = "thumbnail", dest_dir = tempdir()) georef <- fly_georef(fetched, centroids[1:2, ], diff --git a/tests/testthat/test-fly_georef.R b/tests/testthat/test-fly_georef.R index 8eb6257..74d4bb9 100644 --- a/tests/testthat/test-fly_georef.R +++ b/tests/testthat/test-fly_georef.R @@ -124,7 +124,23 @@ test_that("fly_georef rejects invalid rotation", { dest = tempfile(), success = TRUE ) expect_error(fly_georef(fake_fetch, centroids[1, ], rotation = 45), - "one of 0, 90, 180, 270") + "one of") +}) + +test_that("fly_georef auto rotation uses bearing", { + centroids <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE) + dest_fetch <- file.path(tempdir(), "fly_georef_auto_fetch") + unlink(dest_fetch, recursive = TRUE) + + fetched <- fly_fetch(centroids[1, ], type = "thumbnail", + dest_dir = dest_fetch) + dest_georef <- file.path(tempdir(), "fly_georef_auto_out") + unlink(dest_georef, recursive = TRUE) + + # Default is "auto" — should work with film_roll + frame_number + result <- fly_georef(fetched, centroids[1, ], + dest_dir = dest_georef) + expect_true(result$success[1]) }) test_that("fly_georef reads rotation from column", {