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_georef.R b/R/fly_georef.R new file mode 100644 index 0000000..8f018ff --- /dev/null +++ b/R/fly_georef.R @@ -0,0 +1,294 @@ +#' Georeference airphoto images to footprint polygons +#' +#' 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. 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. +#' @param srcnodata Source nodata value passed to GDAL warp. Black pixels +#' matching this value are treated as transparent (alpha=0 for RGB, +#' nodata for grayscale). Default `"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 `NULL` to disable source nodata detection entirely. +#' @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 +#' 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 +#' } +#' +#' 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" ~ 270, +#' .default = NA # fall through to auto +#' ) +#' ``` +#' +#' **Nodata handling:** Two sources of unwanted black pixels are masked: +#' +#' 1. **Warp fill** — GDAL creates black pixels outside the rotated source +#' 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 +#' image. The `srcnodata` parameter (default `"0"`) tells GDAL to treat +#' these as transparent before warping. +#' +#' **Tradeoff:** `srcnodata = "0"` also masks real black pixels (deep +#' shadows). At thumbnail resolution (~1250x1250) this is acceptable — +#' shadow detail is minimal. For full-resolution scans where shadow +#' detail matters, set `srcnodata = NULL` and handle frame masking +#' downstream (e.g., circle detection). +#' +#' **Accuracy:** footprints assume flat terrain and nadir camera angle. +#' The georeferenced images are approximate — useful for visual context, +#' not survey-grade positioning. See [fly_footprint()] for details on +#' limitations. +#' +#' @examples +#' centroids <- sf::st_read(system.file("testdata/photo_centroids.gpkg", package = "fly")) +#' +#' # 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, ], +#' dest_dir = tempdir()) +#' georef +#' +#' @export +fly_georef <- function(fetch_result, photos_sf, + dest_dir = "georef", overwrite = FALSE, + 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) + } + + 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) + + # Build footprints in BC Albers + footprints <- fly_footprint(photos_sf) |> sf::st_transform(3005) + + # Match fetch results to photos by airp_id + ids <- fetch_result$airp_id + + # 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, + dest = NA_character_, + success = FALSE + ) + + for (i in seq_len(nrow(results))) { + if (!fetch_result$success[i]) next + src <- results$source[i] + if (is.na(src) || !file.exists(src)) next + + out_file <- file.path(dest_dir, + sub("\\.[^.]+$", ".tif", basename(src))) + results$dest[i] <- out_file + + if (!overwrite && file.exists(out_file)) { + results$success[i] <- TRUE + next + } + + # Find matching footprint + fp_idx <- which(photos_sf[["airp_id"]] == results$airp_id[i]) + if (length(fp_idx) == 0) next + fp <- footprints[fp_idx[1], ] + + # Per-photo rotation from column, or default + rot <- if (has_rotation_col) { + val <- as.integer(photos_sf[["rotation"]][fp_idx[1]]) + if (is.na(val)) { + if (auto_rotation) 180L else rotation + } else val + } else { + rotation + } + + results$success[i] <- tryCatch( + georef_one(src, fp, out_file, srcnodata = srcnodata, rotation = rot), + error = function(e) { + message("Failed to georef ", basename(src), ": ", e$message) + FALSE + } + ) + } + + n_ok <- sum(results$success) + message("Georeferenced ", n_ok, " of ", nrow(results), " images") + results +} + +#' 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) + 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) + dims <- regmatches(info, regexpr("Size is \\d+, \\d+", info)) + if (length(dims) == 0) return(FALSE) + px <- as.integer(strsplit(sub("Size is ", "", dims), ", ")[[1]]) + ncol_px <- px[1] + nrow_px <- px[2] + + # Count bands from "Band N" lines + n_bands <- length(gregexpr("Band \\d+", info)[[1]]) + is_rgb <- n_bands >= 3 + + # Pixel corners: TL, TR, BR, BL + pixel_corners <- list( + c(0, 0), # TL + c(ncol_px, 0), # TR + c(ncol_px, nrow_px), # BR + c(0, nrow_px) # BL + ) + + # Footprint corners in same order: TL, TR, BR, BL + fp_corners <- list( + coords[4, 1:2], # TL + coords[3, 1:2], # TR + coords[2, 1:2], # BR + coords[1, 1:2] # BL + ) + + # Rotation: shift the footprint corner mapping + + # rotation=0: pixel TL → footprint TL (north-up, original behavior) + # rotation=90: pixel TL → footprint TR (top of image = east) + # rotation=180: pixel TL → footprint BR (top of image = south) + # rotation=270: pixel TL → footprint BL (top of image = west) + n_shifts <- rotation %/% 90 + if (n_shifts > 0) { + fp_corners <- c( + fp_corners[(n_shifts + 1):4], + fp_corners[1:n_shifts] + ) + } + + # Build GCP args mapping pixel corners to (rotated) footprint corners + 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) + + sf::gdal_utils("translate", + source = src, + destination = tmp_file, + options = c("-a_srs", "EPSG:3005", gcp_args) + ) + + # Step 2: warp to target CRS with nodata handling + # srcnodata: masks black source pixels (camera frame borders) + # RGB: alpha band (-dstalpha) for transparent fill in mosaics + # Grayscale: dstnodata=0 for nodata metadata + warp_opts <- c("-t_srs", "EPSG:3005", "-r", "bilinear") + if (!is.null(srcnodata)) { + src_val <- if (is_rgb) { + paste(rep(srcnodata, n_bands), collapse = " ") + } else { + srcnodata + } + warp_opts <- c(warp_opts, "-srcnodata", src_val) + } + if (is_rgb) { + warp_opts <- c(warp_opts, "-dstalpha") + } else { + warp_opts <- c(warp_opts, "-dstnodata", "0") + } + + sf::gdal_utils("warp", + source = tmp_file, + destination = out_file, + options = warp_opts + ) + + 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/R/fly_thumb_georef.R b/R/fly_thumb_georef.R deleted file mode 100644 index 42d15c0..0000000 --- a/R/fly_thumb_georef.R +++ /dev/null @@ -1,180 +0,0 @@ -#' Georeference downloaded thumbnails 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). -#' -#' @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. -#' @param dest_dir Directory for output GeoTIFFs. Created if it does not -#' exist. -#' @param overwrite If `FALSE` (default), skip files that already exist. -#' @param srcnodata Source nodata value passed to GDAL warp. Black pixels -#' matching this value are treated as transparent (alpha=0 for RGB, -#' nodata for grayscale). Default `"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 `NULL` to disable source nodata detection entirely. -#' @return A tibble with columns `airp_id`, `source`, `dest`, and `success`. -#' -#' @details -#' Each thumbnail'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. -#' -#' **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 -#' `dstnodata=0`. -#' 2. **Camera frame borders** — film holder edges, fiducial marks, and -#' scanning artifacts produce black (value 0) pixels within the source -#' image. The `srcnodata` parameter (default `"0"`) tells GDAL to treat -#' these as transparent before warping. -#' -#' **Tradeoff:** `srcnodata = "0"` also masks real black pixels (deep -#' shadows). At thumbnail resolution (~1250x1250) this is acceptable — -#' shadow detail is minimal. For full-resolution scans where shadow -#' detail matters, set `srcnodata = NULL` and handle frame masking -#' downstream (e.g., circle detection). -#' -#' **Accuracy:** footprints assume flat terrain and nadir camera angle. -#' The georeferenced thumbnails are approximate — useful for visual context, -#' not survey-grade positioning. See [fly_footprint()] for details on -#' limitations. -#' -#' @examples -#' centroids <- sf::st_read(system.file("testdata/photo_centroids.gpkg", package = "fly")) -#' -#' # 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 -#' -#' @export -fly_thumb_georef <- function(fetch_result, photos_sf, - dest_dir = "georef", overwrite = FALSE, - srcnodata = "0") { - if (!all(c("airp_id", "dest", "success") %in% names(fetch_result))) { - stop("`fetch_result` must be output from `fly_fetch()`.", call. = FALSE) - } - - dir.create(dest_dir, recursive = TRUE, showWarnings = FALSE) - - # Build footprints in BC Albers - footprints <- fly_footprint(photos_sf) |> sf::st_transform(3005) - - # Match fetch results to photos by airp_id - ids <- fetch_result$airp_id - - results <- dplyr::tibble( - airp_id = ids, - source = fetch_result$dest, - dest = NA_character_, - success = FALSE - ) - - for (i in seq_len(nrow(results))) { - if (!fetch_result$success[i]) next - src <- results$source[i] - if (is.na(src) || !file.exists(src)) next - - out_file <- file.path(dest_dir, - sub("\\.[^.]+$", ".tif", basename(src))) - results$dest[i] <- out_file - - if (!overwrite && file.exists(out_file)) { - results$success[i] <- TRUE - next - } - - # Find matching footprint - fp_idx <- which(photos_sf[["airp_id"]] == results$airp_id[i]) - if (length(fp_idx) == 0) next - fp <- footprints[fp_idx[1], ] - - results$success[i] <- tryCatch( - georef_one(src, fp, out_file, srcnodata = srcnodata), - error = function(e) { - message("Failed to georef ", basename(src), ": ", e$message) - FALSE - } - ) - } - - n_ok <- sum(results$success) - message("Georeferenced ", n_ok, " of ", nrow(results), " thumbnails") - results -} - -#' Georeference a single thumbnail to a footprint polygon -#' @noRd -georef_one <- function(src, fp, out_file, srcnodata = "0") { - # Get footprint corner coordinates - # fly_footprint builds: BL, BR, TR, TL, BL (closing) - coords <- sf::st_coordinates(fp)[1:4, , drop = FALSE] - - # Read image dimensions and band count via GDAL - info <- sf::gdal_utils("info", source = src, quiet = TRUE) - dims <- regmatches(info, regexpr("Size is \\d+, \\d+", info)) - if (length(dims) == 0) return(FALSE) - px <- as.integer(strsplit(sub("Size is ", "", dims), ", ")[[1]]) - ncol_px <- px[1] - nrow_px <- px[2] - - # Count bands from "Band N" lines - 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] - ) - - # Step 1: translate with GCPs - tmp_file <- tempfile(fileext = ".tif") - on.exit(unlink(tmp_file), add = TRUE) - - sf::gdal_utils("translate", - source = src, - destination = tmp_file, - options = c("-a_srs", "EPSG:3005", gcp_args) - ) - - # Step 2: warp to target CRS with nodata handling - # srcnodata: masks black source pixels (camera frame borders) - # RGB: alpha band (-dstalpha) for transparent fill in mosaics - # Grayscale: dstnodata=0 for nodata metadata - warp_opts <- c("-t_srs", "EPSG:3005", "-r", "bilinear") - if (!is.null(srcnodata)) { - src_val <- if (is_rgb) { - paste(rep(srcnodata, n_bands), collapse = " ") - } else { - srcnodata - } - warp_opts <- c(warp_opts, "-srcnodata", src_val) - } - if (is_rgb) { - warp_opts <- c(warp_opts, "-dstalpha") - } else { - warp_opts <- c(warp_opts, "-dstnodata", "0") - } - - sf::gdal_utils("warp", - source = tmp_file, - destination = out_file, - options = warp_opts - ) - - file.exists(out_file) && file.size(out_file) > 0 -} 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_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_georef.Rd b/man/fly_georef.Rd new file mode 100644 index 0000000..255fbf4 --- /dev/null +++ b/man/fly_georef.Rd @@ -0,0 +1,118 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/fly_georef.R +\name{fly_georef} +\alias{fly_georef} +\title{Georeference airphoto images to footprint polygons} +\usage{ +fly_georef( + fetch_result, + photos_sf, + dest_dir = "georef", + overwrite = FALSE, + srcnodata = "0", + rotation = "auto" +) +} +\arguments{ +\item{fetch_result}{A tibble returned by \code{\link[=fly_fetch]{fly_fetch()}}, with columns +\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. 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.} + +\item{overwrite}{If \code{FALSE} (default), skip files that already exist.} + +\item{srcnodata}{Source nodata value passed to GDAL warp. Black pixels +matching this value are treated as transparent (alpha=0 for RGB, +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}{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}. +} +\description{ +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 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 +} + +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{