|
| 1 | +#' Georeference airphoto images to footprint polygons |
| 2 | +#' |
| 3 | +#' Warps images to their estimated ground footprint using GCPs (ground control |
| 4 | +#' points) derived from [fly_footprint()]. Produces georeferenced GeoTIFFs in |
| 5 | +#' BC Albers (EPSG:3005). Works with thumbnails and full-resolution scans. |
| 6 | +#' |
| 7 | +#' @param fetch_result A tibble returned by [fly_fetch()], with columns |
| 8 | +#' `airp_id`, `dest`, and `success`. |
| 9 | +#' @param photos_sf The same sf object passed to `fly_fetch()`, with a |
| 10 | +#' `scale` column for footprint estimation. If a `rotation` column is |
| 11 | +#' present, per-photo rotation values are used (see **Rotation** below). |
| 12 | +#' @param dest_dir Directory for output GeoTIFFs. Created if it does not |
| 13 | +#' exist. |
| 14 | +#' @param overwrite If `FALSE` (default), skip files that already exist. |
| 15 | +#' @param srcnodata Source nodata value passed to GDAL warp. Black pixels |
| 16 | +#' matching this value are treated as transparent (alpha=0 for RGB, |
| 17 | +#' nodata for grayscale). Default `"0"` masks camera frame borders and |
| 18 | +#' film holder edges at the cost of losing real black pixels — acceptable |
| 19 | +#' for thumbnails but may need adjustment for full-resolution scans. |
| 20 | +#' Set to `NULL` to disable source nodata detection entirely. |
| 21 | +#' @param rotation Image rotation in degrees clockwise. One of `"auto"`, |
| 22 | +#' `0`, `90`, `180`, or `270`. `"auto"` (default) computes flight line |
| 23 | +#' bearing from consecutive centroids and derives rotation per-photo — |
| 24 | +#' requires `film_roll` and `frame_number` columns. Fixed values apply |
| 25 | +#' the same rotation to all photos. Overridden per-photo if `photos_sf` |
| 26 | +#' contains a `rotation` column. |
| 27 | +#' @return A tibble with columns `airp_id`, `source`, `dest`, and `success`. |
| 28 | +#' |
| 29 | +#' @details |
| 30 | +#' Each image's four corners are mapped to the corresponding footprint |
| 31 | +#' polygon corners computed by [fly_footprint()] in BC Albers. GDAL |
| 32 | +#' translates the image with GCPs then warps to the target CRS using |
| 33 | +#' bilinear resampling. |
| 34 | +#' |
| 35 | +#' **Rotation:** Aerial photos may appear rotated in their footprints |
| 36 | +#' because the camera orientation relative to north varies by flight |
| 37 | +#' direction, camera mounting, and scanner orientation. The `rotation` |
| 38 | +#' parameter rotates the GCP corner mapping: |
| 39 | +#' \itemize{ |
| 40 | +#' \item `0` — top of image maps to north edge of footprint (original behavior) |
| 41 | +#' \item `90` — top of image maps to east edge (90° clockwise) |
| 42 | +#' \item `180` — top of image maps to south edge (default, correct for most BC photos) |
| 43 | +#' \item `270` — top of image maps to west edge |
| 44 | +#' } |
| 45 | +#' |
| 46 | +#' When `rotation = "auto"`, the bearing-to-rotation formula is: |
| 47 | +#' `floor((bearing + 91) / 90) * 90 %% 360`. This was calibrated on |
| 48 | +#' BC aerial photos spanning 1968–2019 across multiple camera systems |
| 49 | +#' and scanners. Photos on diagonal flight lines (~45° off cardinal) |
| 50 | +#' may be imperfect — check visually and override with a `rotation` |
| 51 | +#' column if needed. |
| 52 | +#' |
| 53 | +#' Within a film roll, consecutive flight legs alternate direction |
| 54 | +#' (back-and-forth pattern), so different frames on the same roll may |
| 55 | +#' need different rotations. This is why `"auto"` computes per-photo, |
| 56 | +#' not per-roll. To override, add a `rotation` column to `photos_sf`: |
| 57 | +#' ``` |
| 58 | +#' photos$rotation <- dplyr::case_when( |
| 59 | +#' photos$film_roll == "bc5282" ~ 270, |
| 60 | +#' .default = NA # fall through to auto |
| 61 | +#' ) |
| 62 | +#' ``` |
| 63 | +#' |
| 64 | +#' **Nodata handling:** Two sources of unwanted black pixels are masked: |
| 65 | +#' |
| 66 | +#' 1. **Warp fill** — GDAL creates black pixels outside the rotated source |
| 67 | +#' frame. RGB images get an alpha band (`-dstalpha`); grayscale use |
| 68 | +#' `dstnodata=0`. |
| 69 | +#' 2. **Camera frame borders** — film holder edges, fiducial marks, and |
| 70 | +#' scanning artifacts produce black (value 0) pixels within the source |
| 71 | +#' image. The `srcnodata` parameter (default `"0"`) tells GDAL to treat |
| 72 | +#' these as transparent before warping. |
| 73 | +#' |
| 74 | +#' **Tradeoff:** `srcnodata = "0"` also masks real black pixels (deep |
| 75 | +#' shadows). At thumbnail resolution (~1250x1250) this is acceptable — |
| 76 | +#' shadow detail is minimal. For full-resolution scans where shadow |
| 77 | +#' detail matters, set `srcnodata = NULL` and handle frame masking |
| 78 | +#' downstream (e.g., circle detection). |
| 79 | +#' |
| 80 | +#' **Accuracy:** footprints assume flat terrain and nadir camera angle. |
| 81 | +#' The georeferenced images are approximate — useful for visual context, |
| 82 | +#' not survey-grade positioning. See [fly_footprint()] for details on |
| 83 | +#' limitations. |
| 84 | +#' |
| 85 | +#' @examples |
| 86 | +#' centroids <- sf::st_read(system.file("testdata/photo_centroids.gpkg", package = "fly")) |
| 87 | +#' |
| 88 | +#' # Fetch and georeference with auto rotation (uses bearing from centroids) |
| 89 | +#' fetched <- fly_fetch(centroids[1:2, ], type = "thumbnail", |
| 90 | +#' dest_dir = tempdir()) |
| 91 | +#' georef <- fly_georef(fetched, centroids[1:2, ], |
| 92 | +#' dest_dir = tempdir()) |
| 93 | +#' georef |
| 94 | +#' |
| 95 | +#' @export |
| 96 | +fly_georef <- function(fetch_result, photos_sf, |
| 97 | + dest_dir = "georef", overwrite = FALSE, |
| 98 | + srcnodata = "0", rotation = "auto") { |
| 99 | + if (!all(c("airp_id", "dest", "success") %in% names(fetch_result))) { |
| 100 | + stop("`fetch_result` must be output from `fly_fetch()`.", call. = FALSE) |
| 101 | + } |
| 102 | + |
| 103 | + auto_rotation <- identical(rotation, "auto") |
| 104 | + if (!auto_rotation) { |
| 105 | + rotation <- as.integer(rotation) |
| 106 | + if (!rotation %in% c(0L, 90L, 180L, 270L)) { |
| 107 | + stop("`rotation` must be one of \"auto\", 0, 90, 180, 270.", call. = FALSE) |
| 108 | + } |
| 109 | + } |
| 110 | + |
| 111 | + dir.create(dest_dir, recursive = TRUE, showWarnings = FALSE) |
| 112 | + |
| 113 | + # Build footprints in BC Albers |
| 114 | + footprints <- fly_footprint(photos_sf) |> sf::st_transform(3005) |
| 115 | + |
| 116 | + # Match fetch results to photos by airp_id |
| 117 | + ids <- fetch_result$airp_id |
| 118 | + |
| 119 | + # Per-photo rotation: column overrides auto/default |
| 120 | + |
| 121 | + has_rotation_col <- "rotation" %in% names(photos_sf) |
| 122 | + |
| 123 | + # Auto-compute bearing → rotation when needed |
| 124 | + if (auto_rotation && !has_rotation_col) { |
| 125 | + if (all(c("film_roll", "frame_number") %in% names(photos_sf))) { |
| 126 | + photos_sf <- fly_bearing(photos_sf) |
| 127 | + photos_sf$rotation <- bearing_to_rotation(photos_sf$bearing) |
| 128 | + has_rotation_col <- TRUE |
| 129 | + } else { |
| 130 | + message("No film_roll/frame_number columns for auto rotation, using 180") |
| 131 | + rotation <- 180L |
| 132 | + auto_rotation <- FALSE |
| 133 | + } |
| 134 | + } |
| 135 | + |
| 136 | + results <- dplyr::tibble( |
| 137 | + airp_id = ids, |
| 138 | + source = fetch_result$dest, |
| 139 | + dest = NA_character_, |
| 140 | + success = FALSE |
| 141 | + ) |
| 142 | + |
| 143 | + for (i in seq_len(nrow(results))) { |
| 144 | + if (!fetch_result$success[i]) next |
| 145 | + src <- results$source[i] |
| 146 | + if (is.na(src) || !file.exists(src)) next |
| 147 | + |
| 148 | + out_file <- file.path(dest_dir, |
| 149 | + sub("\\.[^.]+$", ".tif", basename(src))) |
| 150 | + results$dest[i] <- out_file |
| 151 | + |
| 152 | + if (!overwrite && file.exists(out_file)) { |
| 153 | + results$success[i] <- TRUE |
| 154 | + next |
| 155 | + } |
| 156 | + |
| 157 | + # Find matching footprint |
| 158 | + fp_idx <- which(photos_sf[["airp_id"]] == results$airp_id[i]) |
| 159 | + if (length(fp_idx) == 0) next |
| 160 | + fp <- footprints[fp_idx[1], ] |
| 161 | + |
| 162 | + # Per-photo rotation from column, or default |
| 163 | + rot <- if (has_rotation_col) { |
| 164 | + val <- as.integer(photos_sf[["rotation"]][fp_idx[1]]) |
| 165 | + if (is.na(val)) { |
| 166 | + if (auto_rotation) 180L else rotation |
| 167 | + } else val |
| 168 | + } else { |
| 169 | + rotation |
| 170 | + } |
| 171 | + |
| 172 | + results$success[i] <- tryCatch( |
| 173 | + georef_one(src, fp, out_file, srcnodata = srcnodata, rotation = rot), |
| 174 | + error = function(e) { |
| 175 | + message("Failed to georef ", basename(src), ": ", e$message) |
| 176 | + FALSE |
| 177 | + } |
| 178 | + ) |
| 179 | + } |
| 180 | + |
| 181 | + n_ok <- sum(results$success) |
| 182 | + message("Georeferenced ", n_ok, " of ", nrow(results), " images") |
| 183 | + results |
| 184 | +} |
| 185 | + |
| 186 | +#' Georeference a single image to a footprint polygon |
| 187 | +#' @noRd |
| 188 | +georef_one <- function(src, fp, out_file, srcnodata = "0", rotation = 180) { |
| 189 | + # Get footprint corner coordinates |
| 190 | + # fly_footprint builds: BL, BR, TR, TL, BL (closing) |
| 191 | + coords <- sf::st_coordinates(fp)[1:4, , drop = FALSE] |
| 192 | + # coords: [1]=BL, [2]=BR, [3]=TR, [4]=TL |
| 193 | + |
| 194 | + # Read image dimensions and band count via GDAL |
| 195 | + info <- sf::gdal_utils("info", source = src, quiet = TRUE) |
| 196 | + dims <- regmatches(info, regexpr("Size is \\d+, \\d+", info)) |
| 197 | + if (length(dims) == 0) return(FALSE) |
| 198 | + px <- as.integer(strsplit(sub("Size is ", "", dims), ", ")[[1]]) |
| 199 | + ncol_px <- px[1] |
| 200 | + nrow_px <- px[2] |
| 201 | + |
| 202 | + # Count bands from "Band N" lines |
| 203 | + n_bands <- length(gregexpr("Band \\d+", info)[[1]]) |
| 204 | + is_rgb <- n_bands >= 3 |
| 205 | + |
| 206 | + # Pixel corners: TL, TR, BR, BL |
| 207 | + pixel_corners <- list( |
| 208 | + c(0, 0), # TL |
| 209 | + c(ncol_px, 0), # TR |
| 210 | + c(ncol_px, nrow_px), # BR |
| 211 | + c(0, nrow_px) # BL |
| 212 | + ) |
| 213 | + |
| 214 | + # Footprint corners in same order: TL, TR, BR, BL |
| 215 | + fp_corners <- list( |
| 216 | + coords[4, 1:2], # TL |
| 217 | + coords[3, 1:2], # TR |
| 218 | + coords[2, 1:2], # BR |
| 219 | + coords[1, 1:2] # BL |
| 220 | + ) |
| 221 | + |
| 222 | + # Rotation: shift the footprint corner mapping |
| 223 | + |
| 224 | + # rotation=0: pixel TL → footprint TL (north-up, original behavior) |
| 225 | + # rotation=90: pixel TL → footprint TR (top of image = east) |
| 226 | + # rotation=180: pixel TL → footprint BR (top of image = south) |
| 227 | + # rotation=270: pixel TL → footprint BL (top of image = west) |
| 228 | + n_shifts <- rotation %/% 90 |
| 229 | + if (n_shifts > 0) { |
| 230 | + fp_corners <- c( |
| 231 | + fp_corners[(n_shifts + 1):4], |
| 232 | + fp_corners[1:n_shifts] |
| 233 | + ) |
| 234 | + } |
| 235 | + |
| 236 | + # Build GCP args mapping pixel corners to (rotated) footprint corners |
| 237 | + gcp_args <- character(0) |
| 238 | + for (j in seq_along(pixel_corners)) { |
| 239 | + gcp_args <- c(gcp_args, |
| 240 | + "-gcp", pixel_corners[[j]][1], pixel_corners[[j]][2], |
| 241 | + fp_corners[[j]][1], fp_corners[[j]][2] |
| 242 | + ) |
| 243 | + } |
| 244 | + |
| 245 | + # Step 1: translate with GCPs |
| 246 | + tmp_file <- tempfile(fileext = ".tif") |
| 247 | + on.exit(unlink(tmp_file), add = TRUE) |
| 248 | + |
| 249 | + sf::gdal_utils("translate", |
| 250 | + source = src, |
| 251 | + destination = tmp_file, |
| 252 | + options = c("-a_srs", "EPSG:3005", gcp_args) |
| 253 | + ) |
| 254 | + |
| 255 | + # Step 2: warp to target CRS with nodata handling |
| 256 | + # srcnodata: masks black source pixels (camera frame borders) |
| 257 | + # RGB: alpha band (-dstalpha) for transparent fill in mosaics |
| 258 | + # Grayscale: dstnodata=0 for nodata metadata |
| 259 | + warp_opts <- c("-t_srs", "EPSG:3005", "-r", "bilinear") |
| 260 | + if (!is.null(srcnodata)) { |
| 261 | + src_val <- if (is_rgb) { |
| 262 | + paste(rep(srcnodata, n_bands), collapse = " ") |
| 263 | + } else { |
| 264 | + srcnodata |
| 265 | + } |
| 266 | + warp_opts <- c(warp_opts, "-srcnodata", src_val) |
| 267 | + } |
| 268 | + if (is_rgb) { |
| 269 | + warp_opts <- c(warp_opts, "-dstalpha") |
| 270 | + } else { |
| 271 | + warp_opts <- c(warp_opts, "-dstnodata", "0") |
| 272 | + } |
| 273 | + |
| 274 | + sf::gdal_utils("warp", |
| 275 | + source = tmp_file, |
| 276 | + destination = out_file, |
| 277 | + options = warp_opts |
| 278 | + ) |
| 279 | + |
| 280 | + file.exists(out_file) && file.size(out_file) > 0 |
| 281 | +} |
| 282 | + |
| 283 | +#' Convert flight bearing to GCP rotation |
| 284 | +#' |
| 285 | +#' Formula calibrated on BC aerial photos (1968–2019). |
| 286 | +#' @param bearing Numeric vector of bearings (degrees, 0–360). |
| 287 | +#' @return Integer vector of rotations (0, 90, 180, or 270). NA bearings |
| 288 | +#' return 180 (most common default). |
| 289 | +#' @noRd |
| 290 | +bearing_to_rotation <- function(bearing) { |
| 291 | + rot <- (floor((bearing + 91) / 90) * 90L) %% 360L |
| 292 | + rot[is.na(rot)] <- 180L |
| 293 | + as.integer(rot) |
| 294 | +} |
0 commit comments