diff --git a/.Rbuildignore b/.Rbuildignore index a7c9991..bd653dd 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -12,3 +12,5 @@ ^Meta$ ^CLAUDE\.md$ ^\.claude$ +^Rplots\.pdf$ +^planning$ diff --git a/.gitignore b/.gitignore index cc0b0d7..8798d74 100644 --- a/.gitignore +++ b/.gitignore @@ -6,3 +6,4 @@ inst/doc docs /doc/ /Meta/ +Rplots.pdf diff --git a/DESCRIPTION b/DESCRIPTION index e03b26e..e69b278 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: fly Title: Airphoto Footprint Estimation and Coverage Selection -Version: 0.4.0 -Date: 2026-08-28 +Version: 0.5.0 +Date: 2026-08-29 Authors@R: c( person("Allan", "Irvine", , "al@newgraphenvironment.com", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-3495-2128")), @@ -29,6 +29,7 @@ Suggests: bookdown, furrr, future, + terra, testthat (>= 3.0.0), knitr, rmarkdown diff --git a/NEWS.md b/NEWS.md index bc5411b..049536f 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,15 @@ # fly (development version) +## 0.5.0 (2026-08-29) + +- `fly_footprint()` gains a `dem` argument, sizing each frame from its height above ground instead of the reported scale ([#9](https://github.com/NewGraphEnvironment/fly/issues/9)). On the bundled Upper Bulkley AOI the reported scale understates footprint **area by a median 14%, ranging to 26%** — and always in the same direction, because the scale is referenced to an elevation above the valley floor the photos cover. This is a datum offset, not the slope effect the issue described +- `FLYING_HEIGHT` is metres above sea level, so subtracting terrain elevation is what turns it into the height ground coverage scales with. Elevation is the mean under the whole footprint, not a reading at the centroid — the two differ by up to 140 m on a 7.2 km frame. It is measured in two passes, because the footprint being averaged over is itself what the correction changes +- New `footprint_terrain`, `height_agl` and `dem_coverage` columns record which terrain treatment each frame received, the height it was sized from, and how much of its footprint the DEM actually covered — measured against the cells the footprint should have covered, so that a footprint running past the edge of a cropped DEM is reported rather than counted as complete. `footprint_basis` is unchanged: it is already matched by value downstream, so encoding terrain into it would break caller filters +- `dem` is accepted by `fly_coverage()`, `fly_overlap()`, `fly_filter()`, `fly_select()` and `fly_georef()`, so the correction is reachable from every function that builds a footprint +- Frames the DEM cannot correct — outside its coverage, or with unusable `flying_height` / `focal_length` — fall back to nominal scale with a warning rather than being dropped. A frame the DEM covers only partly is still corrected, from the mean of the covered part, and warns below 95% coverage +- New `inst/testdata/dem.tif`, a clip of NRCan's MRDEM-30. `terra` added to Suggests; it is only needed when a `dem` is supplied +- Footprints remain axis-aligned rectangles under a nadir assumption. Per-corner ray-casting measures ~2% against the 14% the DEM addresses, and is deferred + ## 0.4.0 (2026-08-28) - `fly_footprint()` sizes each frame from its `media` value instead of applying a fixed 9-inch negative to everything. A digital frame has no negative, and the catalogue mixes film and digital in one layer ([#30](https://github.com/NewGraphEnvironment/fly/issues/30)) diff --git a/R/fly_coverage.R b/R/fly_coverage.R index 9408618..a2a5ce7 100644 --- a/R/fly_coverage.R +++ b/R/fly_coverage.R @@ -6,6 +6,9 @@ #' @param photos_sf An sf point object with a `scale` column. #' @param aoi_sf An sf polygon to check coverage against. #' @param by Column name to group by (default `"photo_year"`). +#' @param dem Optional elevation raster passed to [fly_footprint()], sizing each +#' frame from its height above ground instead of the reported scale. See the +#' **Terrain** section of [fly_footprint()]. #' @return A tibble with the grouping column, `n_photos`, `covered_km2`, #' and `coverage_pct`. #' @@ -15,7 +18,7 @@ #' fly_coverage(centroids, aoi, by = "scale") #' #' @export -fly_coverage <- function(photos_sf, aoi_sf, by = "photo_year") { +fly_coverage <- function(photos_sf, aoi_sf, by = "photo_year", dem = NULL) { sf::sf_use_s2(FALSE) on.exit(sf::sf_use_s2(TRUE)) @@ -24,7 +27,7 @@ fly_coverage <- function(photos_sf, aoi_sf, by = "photo_year") { sf::st_make_valid() aoi_area <- as.numeric(sf::st_area(aoi_albers)) - footprints <- fly_footprint(photos_sf) |> sf::st_transform(3005) + footprints <- fly_footprint(photos_sf, dem = dem) |> sf::st_transform(3005) fly_warn_unsized(footprints, "coverage") photos_with_fp <- photos_sf diff --git a/R/fly_filter.R b/R/fly_filter.R index 2579e60..5e6e242 100644 --- a/R/fly_filter.R +++ b/R/fly_filter.R @@ -9,6 +9,9 @@ #' @param method One of `"footprint"` (default) or `"centroid"`. #' @param buffer Buffer distance in metres added to the AOI before testing #' intersection (default 0). Applied in BC Albers (EPSG:3005). +#' @param dem Optional elevation raster passed to [fly_footprint()], sizing each +#' frame from its height above ground instead of the reported scale. See the +#' **Terrain** section of [fly_footprint()]. #' @return A subset of `photos_sf` that intersects the AOI. #' #' @examples @@ -20,7 +23,8 @@ #' nrow(fp_result) >= nrow(ct_result) #' #' @export -fly_filter <- function(photos_sf, aoi_sf, method = c("footprint", "centroid"), buffer = 0) { +fly_filter <- function(photos_sf, aoi_sf, method = c("footprint", "centroid"), + buffer = 0, dem = NULL) { method <- match.arg(method) aoi_3005 <- sf::st_transform(aoi_sf, 3005) |> @@ -36,7 +40,7 @@ fly_filter <- function(photos_sf, aoi_sf, method = c("footprint", "centroid"), b if (method == "centroid") { hits <- sf::st_intersects(photos_sf, aoi_test, sparse = FALSE)[, 1] } else { - footprints <- fly_footprint(photos_sf) + footprints <- fly_footprint(photos_sf, dem = dem) fly_warn_unsized(footprints, "this filter") hits <- sf::st_intersects(footprints, aoi_test, sparse = FALSE)[, 1] } diff --git a/R/fly_footprint.R b/R/fly_footprint.R index c072e7e..32cdf93 100644 --- a/R/fly_footprint.R +++ b/R/fly_footprint.R @@ -25,6 +25,105 @@ fly_warn_unsized <- function(footprints, operation) { invisible(footprints) } +# Coverage below which a footprint's mean elevation stops being trustworthy. +# +# Reprojecting a DEM leaves NA slivers along its edges, so a frame near the +# margin is routinely a percent or two short through no fault of the caller — +# warning on any missing cell at all would fire on a bundled frame that is +# 99.96% covered, and a guard that noisy stops being read. A frame missing more +# than a twentieth of its footprint is a different thing: that is enough for the +# covered part to sit systematically higher or lower than the whole. +fly_dem_coverage_min <- function() 0.95 + +# The DEM-aligned grid a single footprint is counted against. +# +# Named and separate so the "one frame at a time" invariant can be asserted +# rather than merely intended: the size of this grid is the whole difference +# between a bounded allocation and one scaled to the distance between photos. +fly_dem_grid <- function(dem, geom) { + terra::rast( + terra::align(terra::ext(terra::vect(geom)), dem), + resolution = terra::res(dem), + crs = terra::crs(dem) + ) +} + +# Mean ground elevation under each rectangle, with the fraction of the +# rectangle the DEM actually described. +# +# `covered` matters because averaging with na.rm = TRUE cannot tell a +# fully-sampled frame from one hanging half off the edge of the data — both +# yield a number, and only one of them means what it appears to. +fly_dem_sample <- function(dem, rects) { + elev <- rep(NA_real_, length(rects)) + covered <- rep(NA_real_, length(rects)) + ok <- !sf::st_is_empty(rects) + if (!any(ok)) { + return(list(elev = elev, covered = covered)) + } + in_dem <- sf::st_transform(sf::st_sf(geometry = rects[ok]), + sf::st_crs(terra::crs(dem))) + v <- terra::vect(in_dem) + + cells <- terra::extract(dem, v) + # split() keys on the ID column, whose values are 1..n in ascending order, so + # the results come back aligned with `rects[ok]`. + per_frame <- split(cells[, 2], cells[, 1]) + elev[ok] <- vapply(per_frame, function(x) mean(x, na.rm = TRUE), numeric(1)) + got <- vapply(per_frame, function(x) sum(!is.na(x)), numeric(1)) + + # The denominator has to be counted the way the numerator is. extract() takes + # a cell when its *centre* falls inside the polygon, so dividing by the + # footprint's area in cell units compares two different measurements and runs + # about 2/k low for a footprint k cells wide — on a DEM with nothing missing + # that reported 91% coverage and warned. + # + # So count cells on a grid aligned to the DEM's own, per frame. align() snaps + # to the DEM's cell boundaries, which is what puts both counts on the same + # centres. + # + # Per frame, not once over their union: the union's bounding box spans the + # whole photo set, and one frame away from the rest sizes the template to the + # gap between them. On a fixture in this package's own suite that is 243 + # million cells against 16 thousand for the same two frames counted + # separately. Each frame's own template is bounded by one footprint. + expected <- vapply(seq_len(nrow(in_dem)), function(i) { + vi <- terra::vect(in_dem[i, ]) + tmpl <- fly_dem_grid(dem, in_dem[i, ]) + terra::values(tmpl) <- 1L + sum(!is.na(terra::extract(tmpl, vi)[, 2])) + }, numeric(1)) + + covered[ok] <- ifelse(expected > 0, pmin(1, got / expected), 0) + + elev[is.nan(elev)] <- NA_real_ + list(elev = elev, covered = covered) +} + +# Build axis-aligned squares of `half_side` metres about each coordinate pair. +# A half-side that is NA or non-finite yields an empty polygon — the #30 +# contract for a frame whose recording format could not be resolved, and the +# only safe answer for a frame whose metadata divides by zero. An infinite +# rectangle is not merely meaningless: it reaches the sampler, which then tries +# to size a raster to it. +fly_rectangles <- function(coords, half_side) { + sf::st_sfc(lapply(seq_len(nrow(coords)), function(i) { + w <- half_side[i] + if (!is.finite(w)) { + return(sf::st_polygon()) + } + cx <- coords[i, 1] + cy <- coords[i, 2] + sf::st_polygon(list(matrix(c( + cx - w, cy - w, + cx + w, cy - w, + cx + w, cy + w, + cx - w, cy + w, + cx - w, cy - w + ), ncol = 2, byrow = TRUE))) + }), crs = 3005) +} + #' Estimate photo footprint polygons from centroids and scale #' #' Creates rectangular polygons representing the estimated ground coverage @@ -39,8 +138,17 @@ fly_warn_unsized <- function(footprints, operation) { #' @param format_size Named numeric vector of recording-format widths in inches, #' keyed by `media` value, merged over the shipped film defaults. Supply this #' to size frames whose format `fly` does not know — see Details. +#' @param dem Optional elevation raster used to size each frame from its true +#' height above ground rather than the reported scale. A `terra::SpatRaster`, +#' a file path, or a `/vsicurl/` URL. Requires `flying_height` and +#' `focal_length` columns, and the `terra` package. `NULL` (default) keeps the +#' flat-terrain behaviour — see **Terrain** below. #' @return An sf polygon object in the same CRS as input, with footprint -#' rectangles and a `footprint_basis` column recording how each was sized. +#' rectangles, a `footprint_basis` column recording how each was sized, a +#' `footprint_terrain` column recording which terrain treatment was applied, +#' `height_agl` giving the metres above ground each footprint was sized from, +#' and `dem_coverage` giving the fraction of each footprint the DEM actually +#' covered (`0` where it covered none, `NA` only where there is no footprint). #' Frames whose format could not be resolved get an empty geometry. #' #' @details @@ -88,13 +196,84 @@ fly_warn_unsized <- function(footprints, operation) { #' latter returns all `NULL`, which reads as missing data rather than a wrong #' field name. #' -#' **Flat-terrain assumption:** footprints are estimated assuming flat ground -#' beneath the aircraft. In reality terrain slope changes the actual ground -#' coverage — downhill slopes increase the true footprint (ground falls away -#' from the camera), while uphill slopes reduce it. In steep terrain typical -#' of BC valleys, true footprints may differ meaningfully from these estimates. -#' Coverage and overlap calculations downstream (e.g. [fly_coverage()], -#' [fly_overlap()]) inherit this limitation. +#' @section Terrain: +#' +#' Without `dem`, footprints are sized from the reported scale, which assumes +#' flat ground at whatever elevation the scale was computed for. That assumption +#' costs more than it looks: on the bundled Upper Bulkley AOI the reported scale +#' **understates footprint area by a median 14%, ranging to 26%** — and always in +#' the same direction, because the scale is referenced to an elevation above the +#' valley floor the photos actually cover. +#' +#' Supplying `dem` removes that bias. `FLYING_HEIGHT` is metres above sea level, +#' not height above ground, so subtracting terrain elevation is what turns it +#' into the height ground coverage actually scales with: +#' +#' ``` +#' height above ground = flying_height - terrain elevation +#' ground width = format width * (height above ground / focal length) +#' ``` +#' +#' Elevation is the **mean under the whole footprint**, not a reading at the +#' centroid — on a 7.2 km wide 1:31680 frame the two differ by up to 140 m. +#' That is measured in two passes, because the footprint being averaged over is +#' itself what the correction changes: the first pass averages over the +#' nominal-scale rectangle, the second over the rectangle the first produced. +#' The second pass is a refinement rather than the substance — it moves area by +#' at most 0.5% against the correction's own 14% — and a third moves it by +#' 0.03%, so two is where this settles. +#' +#' `footprint_terrain` records what happened to each frame: +#' +#' \describe{ +#' \item{`"nominal_scale"`}{sized from the reported scale (no `dem`, or a +#' fallback — see below)} +#' \item{`"dem_agl"`}{sized from height above ground} +#' \item{`"no_dem_coverage"`}{`dem` supplied but does not cover the frame} +#' \item{`NA`}{no footprint to place — see `footprint_basis`} +#' } +#' +#' A frame the DEM cannot correct falls back to nominal scale with a warning, +#' rather than being dropped. The same applies where the DEM puts terrain at or +#' above the aircraft, which means `flying_height` is not in metres ASL. +#' +#' **Still assumed, with or without a DEM:** the camera points straight down. +#' The BC catalogue carries no tilt, roll or crab, so footprints stay +#' axis-aligned rectangles and corner rays are not projected individually. On +#' this AOI that per-corner refinement is worth roughly 2%, against the 14% the +#' DEM addresses. +#' +#' **DEM sources.** Any raster `terra` can open works. Three that suit BC: +#' +#' \itemize{ +#' \item **MRDEM-30** — NRCan's 30 m bare-earth DTM, all of Canada, public +#' and unauthenticated. A good default, and what the bundled `dem.tif` is +#' cut from: +#' `/vsicurl/https://canelevation-dem.s3.ca-central-1.amazonaws.com/mrdem-30/mrdem-30-dtm.tif` +#' \item **LidarBC** — sub-10 m where coverage exists; query the +#' `stac-dem-bc` STAC catalogue and pass an item's COG URL. +#' \item **BC TRIM** — 25 m provincial DEM via the `bcdata` CLI +#' (`bcdata get-dem`). +#' } +#' +#' Resolution matters less here than extent. A 30 m DEM resolves a 2.7 km +#' footprint's mean elevation perfectly well; a DEM that stops short of the +#' frame edges does not, and this is the ordinary failure rather than an exotic +#' one — a DEM cropped to an AOI simply stops. `no_dem_coverage` is reached only +#' when a footprint finds no elevation at all. A footprint that is merely +#' truncated is still corrected, from the mean of the part the DEM described, +#' and warns once that falls below 95%. `dem_coverage` reports the fraction per +#' frame — measured against the cells the footprint should have covered, not the +#' cells that came back — so a truncated footprint can be filtered rather than +#' merely noticed. +#' +#' Buffer past the **corner** of the widest footprint, not its half-side: the +#' far point of a square is `half_side * sqrt(2)`, which at 1:31680 is 5.1 km +#' rather than 3.6 km. Allow more again for the correction itself, which +#' enlarges footprints before the second pass samples them. +#' +#' Coverage and overlap downstream (e.g. [fly_coverage()], [fly_overlap()]) +#' accept the same `dem` argument and inherit whichever basis you give them. #' #' @examples #' centroids <- sf::st_read(system.file("testdata/photo_centroids.gpkg", package = "fly")) @@ -108,8 +287,22 @@ fly_warn_unsized <- function(footprints, operation) { #' sized <- footprints[footprints$footprint_basis != "unknown_format", ] #' nrow(sized) #' +#' # Terrain-adjusted: size each frame from its height above ground instead of +#' # the reported scale. On this AOI every footprint grows, by a median 14%. +#' # terra is Suggests-only, so the DEM path is guarded here. +#' if (requireNamespace("terra", quietly = TRUE)) { +#' terrain <- fly_footprint( +#' centroids, +#' dem = system.file("testdata/dem.tif", package = "fly") +#' ) +#' print(round(100 * (as.numeric(sf::st_area(sf::st_transform(terrain, 3005))) / +#' as.numeric(sf::st_area(sf::st_transform(footprints, 3005))) - 1), 1)) +#' print(table(terrain$footprint_terrain)) +#' } +#' #' @export -fly_footprint <- function(centroids_sf, negative_size = 9, format_size = NULL) { +fly_footprint <- function(centroids_sf, negative_size = 9, format_size = NULL, + dem = NULL) { if (!inherits(centroids_sf, "sf")) { stop("`centroids_sf` must be an sf object.", call. = FALSE) } @@ -124,6 +317,31 @@ fly_footprint <- function(centroids_sf, negative_size = 9, format_size = NULL) { } } + if (!is.null(dem)) { + rlang::check_installed("terra", "for terrain-adjusted footprints.") + missing_cols <- setdiff(c("flying_height", "focal_length"), names(centroids_sf)) + if (length(missing_cols)) { + stop( + "`dem` needs ", paste0("`", missing_cols, "`", collapse = " and "), + " on `centroids_sf`. Ground coverage scales with height above ground, ", + "which is `flying_height` (metres above sea level) minus terrain ", + "elevation \u2014 without it there is nothing for the DEM to correct.", + call. = FALSE + ) + } + if (!inherits(dem, "SpatRaster")) { + dem <- terra::rast(dem) + } + if (!nzchar(terra::crs(dem))) { + stop( + "`dem` has no CRS, so its cells cannot be located against the photo ", + "centroids. Set one with `terra::crs(dem) <- \"EPSG:3005\"` (or ", + "whichever it is) before passing it.", + call. = FALSE + ) + } + } + input_crs <- sf::st_crs(centroids_sf) pts_3005 <- sf::st_transform(centroids_sf, 3005) coords <- sf::st_coordinates(pts_3005) @@ -158,27 +376,116 @@ fly_footprint <- function(centroids_sf, negative_size = 9, format_size = NULL) { half_side <- width_in * scale_num * 0.0254 / 2 - polys <- lapply(seq_len(nrow(coords)), function(i) { - w <- half_side[i] - if (is.na(w)) { - return(sf::st_polygon()) + # Keyed on half_side, not width_in: an unparseable `scale` also leaves a frame + # with no footprint, and a frame with no footprint has had no terrain + # treatment to report. + terrain <- ifelse(is.na(half_side), NA_character_, "nominal_scale") + height_agl <- rep(NA_real_, n) + dem_coverage <- rep(NA_real_, n) + + if (!is.null(dem)) { + focal_m <- centroids_sf$focal_length / 1000 + resize <- function(e) { + width_in * ((centroids_sf$flying_height - e) / focal_m) * 0.0254 / 2 } - cx <- coords[i, 1] - cy <- coords[i, 2] - corners <- matrix(c( - cx - w, cy - w, - cx + w, cy - w, - cx + w, cy + w, - cx - w, cy + w, - cx - w, cy - w - ), ncol = 2, byrow = TRUE) - sf::st_polygon(list(corners)) - }) + + # Two passes. The first averages the DEM over the nominal-scale rectangle, + # which yields a height above ground and so a better rectangle; the second + # averages over that one, since the window being averaged is itself what + # the correction changes. + # + # Keep the size of this in proportion. The correction as a whole moves area + # by a median 14%; the second pass moves it by at most a further 0.53%, and + # a third by 0.03%. It is worth one more extract, not the emphasis a bare + # "iterates until it converges" would imply. + first <- fly_dem_sample(dem, fly_rectangles(coords, half_side)) + second <- fly_dem_sample(dem, fly_rectangles(coords, resize(first$elev))) + + # Keep the first pass wherever the second could not improve on it, so a + # frame is never lost to the resize alone. + elev <- ifelse(is.na(second$elev), first$elev, second$elev) + + sized <- !is.na(half_side) + agl <- centroids_sf$flying_height - elev + candidate <- resize(elev) + + # Classify on the half-side we would actually use, not on the inputs that + # feed it. An NA or zero `focal_length`, or an NA `flying_height`, yields a + # non-finite half-side that would otherwise become an empty geometry — and + # an empty geometry here is indistinguishable from one whose recording + # format was never resolved, so the frame would vanish under a warning + # pointing at `format_size` rather than at its own metadata. + corrected <- sized & is.finite(candidate) & candidate > 0 + uncovered <- sized & !corrected & is.na(elev) + unusable <- sized & !corrected & !uncovered + + # Coverage has to describe the footprint that is actually returned. A + # corrected frame ships the second pass's rectangle, so it takes the second + # pass's coverage; a frame that fell back ships the nominal one, so it takes + # the first pass's. Reading the second pass for a fallback frame is not a + # rounding difference — terrain above the aircraft gives a negative + # half-side, which draws a mirrored rectangle somewhere else entirely, and + # that measured 100% coverage for a footprint only 30% covered. + covered <- ifelse(corrected, second$covered, first$covered) + + # Every fallback keeps the frame at nominal scale rather than dropping it: + # a frame we cannot correct is still a frame. + if (any(uncovered)) { + warning( + sum(uncovered), " of ", sum(sized), " frames fall outside the DEM's ", + "coverage and were sized from nominal scale instead. See ", + "`footprint_terrain`.", + call. = FALSE + ) + } + if (any(unusable)) { + warning( + sum(unusable), " of ", sum(sized), " frames have `flying_height` or ", + "`focal_length` values that give no usable height above ground \u2014 ", + "missing, zero, or terrain at or above the aircraft. Check that ", + "`flying_height` is metres above sea level. Sized from nominal scale ", + "instead. See `footprint_terrain`.", + call. = FALSE + ) + } + + # A footprint hanging off the edge of the DEM still yields a mean, taken + # from whichever part had data. That is the best estimate available and is + # kept — but it is not the full-frame mean it would otherwise be taken for, + # so it is reported rather than passed off silently. + partial <- corrected & !is.na(covered) & covered < fly_dem_coverage_min() + if (any(partial)) { + warning( + sum(partial), " of ", sum(corrected), " corrected frames are less than ", + round(100 * fly_dem_coverage_min()), "% covered by the DEM (as little ", + "as ", round(100 * min(covered[partial])), "% of one footprint). Their ", + "ground elevation is the mean of the covered part, which need not ", + "represent the whole. Buffer the DEM past the corner of the widest ", + "footprint \u2014 half its width times sqrt(2), not half its width. ", + "See `dem_coverage`.", + call. = FALSE + ) + } + + half_side[corrected] <- candidate[corrected] + height_agl[corrected] <- agl[corrected] + # Reported for every frame that had a footprint to sample, not only the + # corrected ones: `no_dem_coverage` is a measured zero, and leaving it NA + # makes the documented "filter on dem_coverage" workflow impossible. + dem_coverage[sized] <- covered[sized] + terrain[sized] <- "nominal_scale" + terrain[corrected] <- "dem_agl" + terrain[uncovered] <- "no_dem_coverage" + terrain[unusable] <- "nominal_scale" + } result <- sf::st_sf( sf::st_drop_geometry(pts_3005), footprint_basis = basis, - geometry = sf::st_sfc(polys, crs = 3005) + footprint_terrain = terrain, + height_agl = height_agl, + dem_coverage = dem_coverage, + geometry = fly_rectangles(coords, half_side) ) sf::st_transform(result, input_crs) diff --git a/R/fly_georef.R b/R/fly_georef.R index 7f0f0bb..dd757bc 100644 --- a/R/fly_georef.R +++ b/R/fly_georef.R @@ -18,6 +18,9 @@ #' 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 dem Optional elevation raster passed to [fly_footprint()], sizing each +#' 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 — @@ -77,10 +80,12 @@ #' 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. +#' **Accuracy:** footprints assume a nadir camera angle, and without `dem` +#' they also assume flat terrain. Passing `dem` sizes each frame from its +#' height above ground, which on steep ground is the larger of the two error +#' terms — but the images stay approximate either way, useful for visual +#' context rather than survey-grade positioning. See the **Terrain** section +#' of [fly_footprint()]. #' #' @examples #' centroids <- sf::st_read(system.file("testdata/photo_centroids.gpkg", package = "fly")) @@ -95,7 +100,7 @@ #' @export fly_georef <- function(fetch_result, photos_sf, dest_dir = "georef", overwrite = FALSE, - srcnodata = "0", rotation = "auto") { + srcnodata = "0", rotation = "auto", dem = NULL) { if (!all(c("airp_id", "dest", "success") %in% names(fetch_result))) { stop("`fetch_result` must be output from `fly_fetch()`.", call. = FALSE) } @@ -111,7 +116,7 @@ fly_georef <- function(fetch_result, photos_sf, dir.create(dest_dir, recursive = TRUE, showWarnings = FALSE) # Build footprints in BC Albers - footprints <- fly_footprint(photos_sf) |> sf::st_transform(3005) + footprints <- fly_footprint(photos_sf, dem = dem) |> sf::st_transform(3005) fly_warn_unsized(footprints, "georeferencing") # Match fetch results to photos by airp_id diff --git a/R/fly_overlap.R b/R/fly_overlap.R index 9833f10..9969501 100644 --- a/R/fly_overlap.R +++ b/R/fly_overlap.R @@ -4,10 +4,17 @@ #' area and the percentage of each photo's footprint that overlaps. #' Most useful on same-scale photos from the same flight. #' -#' Overlap percentages are estimates based on flat-terrain footprints from -#' [fly_footprint()]. See that function for details on terrain limitations. +#' Overlap percentages are estimates from [fly_footprint()], and inherit +#' whatever basis it was given. Without `dem` that is the reported scale, which +#' assumes flat ground and understates footprint area wherever the terrain sits +#' below the elevation the scale was computed for. Pass `dem` to size each frame +#' from its height above ground instead — see the **Terrain** section of +#' [fly_footprint()]. #' #' @param photos_sf An sf point object with a `scale` column. +#' @param dem Optional elevation raster passed to [fly_footprint()], sizing each +#' frame from its height above ground instead of the reported scale. See the +#' **Terrain** section of [fly_footprint()]. #' @return A tibble with columns `photo_a`, `photo_b`, `overlap_km2`, #' `pct_of_a`, and `pct_of_b`. Only pairs with non-zero overlap are returned. #' @@ -19,11 +26,11 @@ #' fly_overlap(selected) #' #' @export -fly_overlap <- function(photos_sf) { +fly_overlap <- function(photos_sf, dem = NULL) { sf::sf_use_s2(FALSE) on.exit(sf::sf_use_s2(TRUE)) - footprints <- fly_footprint(photos_sf) |> sf::st_transform(3005) + footprints <- fly_footprint(photos_sf, dem = dem) |> sf::st_transform(3005) fly_warn_unsized(footprints, "overlap") n <- nrow(footprints) diff --git a/R/fly_select.R b/R/fly_select.R index c8f513b..16796c8 100644 --- a/R/fly_select.R +++ b/R/fly_select.R @@ -11,6 +11,9 @@ #' (every photo touching the AOI). #' @param target_coverage Stop when this fraction is reached (default 0.95). #' Only used when `mode = "minimal"`. +#' @param dem Optional elevation raster passed to [fly_footprint()], sizing each +#' frame from its height above ground instead of the reported scale. See the +#' **Terrain** section of [fly_footprint()]. #' @param component_ensure If `TRUE` (default `FALSE`), guarantee that every #' polygon component of `aoi_sf` is covered by at least one photo before #' running the greedy selection. Useful for multi-polygon AOIs (e.g. patchy @@ -36,22 +39,22 @@ #' @export fly_select <- function(photos_sf, aoi_sf, mode = "minimal", target_coverage = 0.95, - component_ensure = FALSE) { + component_ensure = FALSE, dem = NULL) { mode <- match.arg(mode, c("minimal", "all")) if (mode == "all") { - return(fly_select_all(photos_sf, aoi_sf)) + return(fly_select_all(photos_sf, aoi_sf, dem)) } - fly_select_minimal(photos_sf, aoi_sf, target_coverage, component_ensure) + fly_select_minimal(photos_sf, aoi_sf, target_coverage, component_ensure, dem) } #' @noRd -fly_select_all <- function(photos_sf, aoi_sf) { +fly_select_all <- function(photos_sf, aoi_sf, dem = NULL) { sf::sf_use_s2(FALSE) on.exit(sf::sf_use_s2(TRUE)) - footprints <- fly_footprint(photos_sf) + footprints <- fly_footprint(photos_sf, dem = dem) fly_warn_unsized(footprints, "this selection") aoi_union <- sf::st_transform(aoi_sf, sf::st_crs(footprints)) |> sf::st_union() |> @@ -100,7 +103,7 @@ ensure_component_coverage <- function(footprints, aoi_albers) { #' @noRd fly_select_minimal <- function(photos_sf, aoi_sf, target_coverage, - component_ensure) { + component_ensure, dem = NULL) { sf::sf_use_s2(FALSE) on.exit(sf::sf_use_s2(TRUE)) @@ -109,7 +112,7 @@ fly_select_minimal <- function(photos_sf, aoi_sf, target_coverage, sf::st_make_valid() aoi_area <- as.numeric(sf::st_area(aoi_albers)) - footprints <- fly_footprint(photos_sf) |> sf::st_transform(3005) + footprints <- fly_footprint(photos_sf, dem = dem) |> sf::st_transform(3005) fly_warn_unsized(footprints, "this selection") footprints$photo_idx <- seq_len(nrow(footprints)) diff --git a/data-raw/make_testdata.R b/data-raw/make_testdata.R index 864b64b..0bb575c 100644 --- a/data-raw/make_testdata.R +++ b/data-raw/make_testdata.R @@ -7,6 +7,11 @@ # Dual-scale coverage: 1:12000 and 1:31680 (1968). # # Source: diggs cached data (BC Data Catalogue + flooded VCA output) +# ... except dem.tif, which is fetched over the network from NRCan's MRDEM-30 +# (see the DEM section at the foot of this script). That step needs outbound +# HTTPS to canelevation-dem.s3.ca-central-1.amazonaws.com; everything else is +# local. Takes a few seconds. +# # Run from fly repo root: Rscript data-raw/make_testdata.R library(sf) @@ -116,4 +121,59 @@ if (nrow(streams_clip) > 0) { message("lakes.gpkg: ", nrow(lakes), " lake(s)") } +# --- DEM: MRDEM-30 clip for terrain-adjusted footprints (#9) -------------- +# +# MRDEM-30 is NRCan's Medium-Resolution Digital Elevation Model: a 30 m +# bare-earth DTM covering all of Canada as a single ~84 GB Cloud-Optimized +# GeoTIFF in EPSG:3979, public and unauthenticated. `/vsicurl/` range-reads +# only the bytes intersecting our AOI, so nothing near 84 GB is transferred. +# +# Same product `flooded::fl_dem_aoi()` defaults to, which is why fly reaches +# for it rather than a second elevation source. Compared head to head against +# elevatr z=10 over this AOI, the two agree to within 0.42 percentage points on +# the resulting footprint-area correction (fly#9) — so this choice is about +# provenance and dependencies, not accuracy. +# +# Buffer is 5.4 km. The widest footprint is 1:31680, 7.24 km across, so its +# half-side is 3.62 km — but a square's *corner* is the far point, at +# half_side * sqrt(2) = 5.12 km. Buffering by the half-side leaves the four +# corners of every edge frame hanging over no-data, which is a partial mean +# rather than a clean fallback. 5.4 km clears that with a little margin. +# +# The terrain correction also enlarges footprints by up to ~25% here, and the +# second sampling pass averages over the *enlarged* rectangle, so the margin +# has to cover that too. +# +# Crop in the source CRS, project after. Never reproject the whole COG. + +message("Fetching MRDEM-30 clip (network step) ...") +terra::setGDALconfig("GDAL_HTTP_MAX_RETRY", "3") +terra::setGDALconfig("GDAL_HTTP_RETRY_DELAY", "2") + +mrdem_url <- paste0( + "/vsicurl/https://canelevation-dem.s3.ca-central-1.amazonaws.com/", + "mrdem-30/mrdem-30-dtm.tif" +) + +dem_aoi <- st_sf(geometry = st_union(st_buffer(st_transform(test_photos, 3005), 5400))) +dem_src <- terra::rast(mrdem_url) +dem_clip <- terra::crop( + dem_src, + terra::vect(st_transform(dem_aoi, st_crs(terra::crs(dem_src)))), + snap = "out" +) +dem_clip <- terra::project(dem_clip, "EPSG:3005") + +# INT2S: elevations are metres, and a sub-metre DEM read is not meaningful at +# 30 m posting. Halves the file against FLT4S for no loss that matters here. +dem_path <- file.path(outdir, "dem.tif") +terra::writeRaster( + dem_clip, dem_path, overwrite = TRUE, datatype = "INT2S", + gdal = c("COMPRESS=DEFLATE", "PREDICTOR=2", "TILED=YES") +) +message("dem.tif: ", paste(dim(dem_clip)[1:2], collapse = "x"), " cells at ", + round(terra::res(dem_clip)[1], 1), " m, ", + paste(round(as.vector(terra::minmax(dem_clip, compute = TRUE))), collapse = "-"), + " m, ", round(file.size(dem_path) / 1024), " KB") + message("\nDone. Test data in: ", outdir) diff --git a/inst/testdata/dem.tif b/inst/testdata/dem.tif new file mode 100644 index 0000000..4ddd5eb Binary files /dev/null and b/inst/testdata/dem.tif differ diff --git a/man/fly_coverage.Rd b/man/fly_coverage.Rd index 14fb5c1..14de78d 100644 --- a/man/fly_coverage.Rd +++ b/man/fly_coverage.Rd @@ -4,7 +4,7 @@ \alias{fly_coverage} \title{Check photo coverage of an AOI by group} \usage{ -fly_coverage(photos_sf, aoi_sf, by = "photo_year") +fly_coverage(photos_sf, aoi_sf, by = "photo_year", dem = NULL) } \arguments{ \item{photos_sf}{An sf point object with a \code{scale} column.} @@ -12,6 +12,10 @@ fly_coverage(photos_sf, aoi_sf, by = "photo_year") \item{aoi_sf}{An sf polygon to check coverage against.} \item{by}{Column name to group by (default \code{"photo_year"}).} + +\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 +\strong{Terrain} section of \code{\link[=fly_footprint]{fly_footprint()}}.} } \value{ A tibble with the grouping column, \code{n_photos}, \code{covered_km2}, diff --git a/man/fly_filter.Rd b/man/fly_filter.Rd index ecc7aec..a25a332 100644 --- a/man/fly_filter.Rd +++ b/man/fly_filter.Rd @@ -4,7 +4,13 @@ \alias{fly_filter} \title{Filter photos by spatial relationship with an AOI} \usage{ -fly_filter(photos_sf, aoi_sf, method = c("footprint", "centroid"), buffer = 0) +fly_filter( + photos_sf, + aoi_sf, + method = c("footprint", "centroid"), + buffer = 0, + dem = NULL +) } \arguments{ \item{photos_sf}{An sf point object with a \code{scale} column.} @@ -15,6 +21,10 @@ fly_filter(photos_sf, aoi_sf, method = c("footprint", "centroid"), buffer = 0) \item{buffer}{Buffer distance in metres added to the AOI before testing intersection (default 0). Applied in BC Albers (EPSG:3005).} + +\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 +\strong{Terrain} section of \code{\link[=fly_footprint]{fly_footprint()}}.} } \value{ A subset of \code{photos_sf} that intersects the AOI. diff --git a/man/fly_footprint.Rd b/man/fly_footprint.Rd index 1b96b2b..1cc4574 100644 --- a/man/fly_footprint.Rd +++ b/man/fly_footprint.Rd @@ -4,7 +4,7 @@ \alias{fly_footprint} \title{Estimate photo footprint polygons from centroids and scale} \usage{ -fly_footprint(centroids_sf, negative_size = 9, format_size = NULL) +fly_footprint(centroids_sf, negative_size = 9, format_size = NULL, dem = NULL) } \arguments{ \item{centroids_sf}{An sf point object with a \code{scale} column (e.g. "1:31680"). @@ -18,10 +18,20 @@ recording format per frame when present.} \item{format_size}{Named numeric vector of recording-format widths in inches, keyed by \code{media} value, merged over the shipped film defaults. Supply this to size frames whose format \code{fly} does not know — see Details.} + +\item{dem}{Optional elevation raster used to size each frame from its true +height above ground rather than the reported scale. A \code{terra::SpatRaster}, +a file path, or a \verb{/vsicurl/} URL. Requires \code{flying_height} and +\code{focal_length} columns, and the \code{terra} package. \code{NULL} (default) keeps the +flat-terrain behaviour — see \strong{Terrain} below.} } \value{ An sf polygon object in the same CRS as input, with footprint -rectangles and a \code{footprint_basis} column recording how each was sized. +rectangles, a \code{footprint_basis} column recording how each was sized, a +\code{footprint_terrain} column recording which terrain treatment was applied, +\code{height_agl} giving the metres above ground each footprint was sized from, +and \code{dem_coverage} giving the fraction of each footprint the DEM actually +covered (\code{0} where it covered none, \code{NA} only where there is no footprint). Frames whose format could not be resolved get an empty geometry. } \description{ @@ -71,15 +81,87 @@ Filter on \code{footprint_basis} to keep only frames sized from a known format. the bundled test data. Note the field is \code{SCALE}, not \code{PHOTO_SCALE} — the latter returns all \code{NULL}, which reads as missing data rather than a wrong field name. +} +\section{Terrain}{ + + +Without \code{dem}, footprints are sized from the reported scale, which assumes +flat ground at whatever elevation the scale was computed for. That assumption +costs more than it looks: on the bundled Upper Bulkley AOI the reported scale +\strong{understates footprint area by a median 14\%, ranging to 26\%} — and always in +the same direction, because the scale is referenced to an elevation above the +valley floor the photos actually cover. + +Supplying \code{dem} removes that bias. \code{FLYING_HEIGHT} is metres above sea level, +not height above ground, so subtracting terrain elevation is what turns it +into the height ground coverage actually scales with: + +\if{html}{\out{
}}\preformatted{height above ground = flying_height - terrain elevation +ground width = format width * (height above ground / focal length) +}\if{html}{\out{
}} + +Elevation is the \strong{mean under the whole footprint}, not a reading at the +centroid — on a 7.2 km wide 1:31680 frame the two differ by up to 140 m. +That is measured in two passes, because the footprint being averaged over is +itself what the correction changes: the first pass averages over the +nominal-scale rectangle, the second over the rectangle the first produced. +The second pass is a refinement rather than the substance — it moves area by +at most 0.5\% against the correction's own 14\% — and a third moves it by +0.03\%, so two is where this settles. + +\code{footprint_terrain} records what happened to each frame: + +\describe{ +\item{\code{"nominal_scale"}}{sized from the reported scale (no \code{dem}, or a +fallback — see below)} +\item{\code{"dem_agl"}}{sized from height above ground} +\item{\code{"no_dem_coverage"}}{\code{dem} supplied but does not cover the frame} +\item{\code{NA}}{no footprint to place — see \code{footprint_basis}} +} + +A frame the DEM cannot correct falls back to nominal scale with a warning, +rather than being dropped. The same applies where the DEM puts terrain at or +above the aircraft, which means \code{flying_height} is not in metres ASL. + +\strong{Still assumed, with or without a DEM:} the camera points straight down. +The BC catalogue carries no tilt, roll or crab, so footprints stay +axis-aligned rectangles and corner rays are not projected individually. On +this AOI that per-corner refinement is worth roughly 2\%, against the 14\% the +DEM addresses. -\strong{Flat-terrain assumption:} footprints are estimated assuming flat ground -beneath the aircraft. In reality terrain slope changes the actual ground -coverage — downhill slopes increase the true footprint (ground falls away -from the camera), while uphill slopes reduce it. In steep terrain typical -of BC valleys, true footprints may differ meaningfully from these estimates. -Coverage and overlap calculations downstream (e.g. \code{\link[=fly_coverage]{fly_coverage()}}, -\code{\link[=fly_overlap]{fly_overlap()}}) inherit this limitation. +\strong{DEM sources.} Any raster \code{terra} can open works. Three that suit BC: + +\itemize{ +\item \strong{MRDEM-30} — NRCan's 30 m bare-earth DTM, all of Canada, public +and unauthenticated. A good default, and what the bundled \code{dem.tif} is +cut from: +\verb{/vsicurl/https://canelevation-dem.s3.ca-central-1.amazonaws.com/mrdem-30/mrdem-30-dtm.tif} +\item \strong{LidarBC} — sub-10 m where coverage exists; query the +\code{stac-dem-bc} STAC catalogue and pass an item's COG URL. +\item \strong{BC TRIM} — 25 m provincial DEM via the \code{bcdata} CLI +(\verb{bcdata get-dem}). } + +Resolution matters less here than extent. A 30 m DEM resolves a 2.7 km +footprint's mean elevation perfectly well; a DEM that stops short of the +frame edges does not, and this is the ordinary failure rather than an exotic +one — a DEM cropped to an AOI simply stops. \code{no_dem_coverage} is reached only +when a footprint finds no elevation at all. A footprint that is merely +truncated is still corrected, from the mean of the part the DEM described, +and warns once that falls below 95\%. \code{dem_coverage} reports the fraction per +frame — measured against the cells the footprint should have covered, not the +cells that came back — so a truncated footprint can be filtered rather than +merely noticed. + +Buffer past the \strong{corner} of the widest footprint, not its half-side: the +far point of a square is \code{half_side * sqrt(2)}, which at 1:31680 is 5.1 km +rather than 3.6 km. Allow more again for the correction itself, which +enlarges footprints before the second pass samples them. + +Coverage and overlap downstream (e.g. \code{\link[=fly_coverage]{fly_coverage()}}, \code{\link[=fly_overlap]{fly_overlap()}}) +accept the same \code{dem} argument and inherit whichever basis you give them. +} + \examples{ centroids <- sf::st_read(system.file("testdata/photo_centroids.gpkg", package = "fly")) footprints <- fly_footprint(centroids) @@ -92,4 +174,17 @@ table(footprints$footprint_basis) sized <- footprints[footprints$footprint_basis != "unknown_format", ] nrow(sized) +# Terrain-adjusted: size each frame from its height above ground instead of +# the reported scale. On this AOI every footprint grows, by a median 14\%. +# terra is Suggests-only, so the DEM path is guarded here. +if (requireNamespace("terra", quietly = TRUE)) { + terrain <- fly_footprint( + centroids, + dem = system.file("testdata/dem.tif", package = "fly") + ) + print(round(100 * (as.numeric(sf::st_area(sf::st_transform(terrain, 3005))) / + as.numeric(sf::st_area(sf::st_transform(footprints, 3005))) - 1), 1)) + print(table(terrain$footprint_terrain)) +} + } diff --git a/man/fly_georef.Rd b/man/fly_georef.Rd index 255fbf4..a12010e 100644 --- a/man/fly_georef.Rd +++ b/man/fly_georef.Rd @@ -10,7 +10,8 @@ fly_georef( dest_dir = "georef", overwrite = FALSE, srcnodata = "0", - rotation = "auto" + rotation = "auto", + dem = NULL ) } \arguments{ @@ -39,6 +40,10 @@ 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.} + +\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 +\strong{Terrain} section of \code{\link[=fly_footprint]{fly_footprint()}}.} } \value{ A tibble with columns \code{airp_id}, \code{source}, \code{dest}, and \code{success}. @@ -100,10 +105,12 @@ shadow detail is minimal. For full-resolution scans where shadow 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 images are approximate — useful for visual context, -not survey-grade positioning. See \code{\link[=fly_footprint]{fly_footprint()}} for details on -limitations. +\strong{Accuracy:} footprints assume a nadir camera angle, and without \code{dem} +they also assume flat terrain. Passing \code{dem} sizes each frame from its +height above ground, which on steep ground is the larger of the two error +terms — but the images stay approximate either way, useful for visual +context rather than survey-grade positioning. See the \strong{Terrain} section +of \code{\link[=fly_footprint]{fly_footprint()}}. } \examples{ centroids <- sf::st_read(system.file("testdata/photo_centroids.gpkg", package = "fly")) diff --git a/man/fly_overlap.Rd b/man/fly_overlap.Rd index e9f8177..39718f7 100644 --- a/man/fly_overlap.Rd +++ b/man/fly_overlap.Rd @@ -4,10 +4,14 @@ \alias{fly_overlap} \title{Compute pairwise overlap between photo footprints} \usage{ -fly_overlap(photos_sf) +fly_overlap(photos_sf, dem = NULL) } \arguments{ \item{photos_sf}{An sf point object with a \code{scale} column.} + +\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 +\strong{Terrain} section of \code{\link[=fly_footprint]{fly_footprint()}}.} } \value{ A tibble with columns \code{photo_a}, \code{photo_b}, \code{overlap_km2}, @@ -19,8 +23,12 @@ area and the percentage of each photo's footprint that overlaps. Most useful on same-scale photos from the same flight. } \details{ -Overlap percentages are estimates based on flat-terrain footprints from -\code{\link[=fly_footprint]{fly_footprint()}}. See that function for details on terrain limitations. +Overlap percentages are estimates from \code{\link[=fly_footprint]{fly_footprint()}}, and inherit +whatever basis it was given. Without \code{dem} that is the reported scale, which +assumes flat ground and understates footprint area wherever the terrain sits +below the elevation the scale was computed for. Pass \code{dem} to size each frame +from its height above ground instead — see the \strong{Terrain} section of +\code{\link[=fly_footprint]{fly_footprint()}}. } \examples{ centroids <- sf::st_read(system.file("testdata/photo_centroids.gpkg", package = "fly")) diff --git a/man/fly_select.Rd b/man/fly_select.Rd index 55a546a..a09a532 100644 --- a/man/fly_select.Rd +++ b/man/fly_select.Rd @@ -9,7 +9,8 @@ fly_select( aoi_sf, mode = "minimal", target_coverage = 0.95, - component_ensure = FALSE + component_ensure = FALSE, + dem = NULL ) } \arguments{ @@ -29,6 +30,10 @@ polygon component of \code{aoi_sf} is covered by at least one photo before running the greedy selection. Useful for multi-polygon AOIs (e.g. patchy floodplain fragments) where small components might otherwise get zero coverage. Only used when \code{mode = "minimal"}.} + +\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 +\strong{Terrain} section of \code{\link[=fly_footprint]{fly_footprint()}}.} } \value{ An sf object (subset of \code{photos_sf}). For \code{mode = "minimal"}, diff --git a/planning/archive/2026-08-issue-9-dem-terrain/README.md b/planning/archive/2026-08-issue-9-dem-terrain/README.md new file mode 100644 index 0000000..0667c26 --- /dev/null +++ b/planning/archive/2026-08-issue-9-dem-terrain/README.md @@ -0,0 +1,36 @@ +# Issue #9 — DEM-based terrain-adjusted footprints + +`fly_footprint()` sized every frame from the reported nominal scale. Measurement on the +bundled Upper Bulkley AOI reframed the issue twice over: `FLYING_HEIGHT` turned out to be +metres **above sea level** rather than height above ground, and the resulting error is a +**datum offset** rather than the slope effect the issue described — reported scale +understates footprint area by a median 14%, ranging to 27%, and always in the same +direction, because the scale is referenced to an elevation above the valley floor. + +Built the true-scale rectangle: size each frame from `flying_height - terrain elevation`, +sampled in two passes (centroid, then mean under the resulting rectangle — they differ by +up to 130 m on a 7.2 km frame). Geometry stays rectangular, so downstream consumers were +unaffected. Per-corner ray-casting measured ~2% against the 14% and was deferred. + +DEM source is MRDEM-30, NRCan's 30 m bare-earth DTM, chosen over `elevatr` after a +head-to-head: the two agree to within 0.42 percentage points, so the choice rested on +MRDEM needing no dependency beyond `terra` and being the product `flooded` already uses. + +Two problems surfaced that the issue did not name. `dem` had to be threaded through all six +internal `fly_footprint()` call sites or the correction was unreachable from every function +a user actually calls. And an NA or zero `focal_length` made the corrected half-side +non-finite, producing an **empty geometry while `footprint_terrain` claimed `"dem_agl"`** — +downstream that reads as an unresolved recording format and sends the user to `format_size` +for a metadata problem. Fixed by classifying on the computed half-side rather than its +inputs. + +Deviated from the issue on one point: terrain went into a new `footprint_terrain` column +rather than into `footprint_basis` as suggested, because that column is already matched by +value downstream. Issue body edited to record it. + +`/code-check`'s subagent rounds were not run — the session barred the Agent tool. Reviewed +against the checklist directly, which is what caught the empty-geometry bug. + +Closed by PR (`Fixes #9`). Suite 176 pass / 0 fail; `R CMD check` 0 errors, 0 warnings, +2 pre-existing NOTEs; vignette rebuilds. Released 0.5.0. +Follow-up: #10 (tilt/roll), ray-cast footprints still open. diff --git a/planning/archive/2026-08-issue-9-dem-terrain/findings.md b/planning/archive/2026-08-issue-9-dem-terrain/findings.md new file mode 100644 index 0000000..3273d4b --- /dev/null +++ b/planning/archive/2026-08-issue-9-dem-terrain/findings.md @@ -0,0 +1,135 @@ +# Findings — DEM-based terrain-adjusted footprints (#9) + +## `flying_height` is metres ASL, not height above ground + +Measured on `inst/testdata/photo_centroids.gpkg` (n = 20): + +| scale | focal (mm) | flying_height (m) | = feet | H_agl from scale (m) | implied ground elev (m) | +|---|---|---|---|---|---| +| 1:12000 | 153 | 2438 | 7999 | 1836 | 602 | +| 1:12000 | 153 | 2591 | 8501 | 1836 | 755 | +| 1:31680 | 153 | 5944 | 19501 | 4847 | 1097 | +| 1:31680 | 153 | 6096 | 20000 | 4847 | 1249 | + +The round-feet values are the tell. Implied ground elevations (602-1249 m) bracket real +Houston-area terrain, confirming the field is ASL. This is what makes a DEM necessary +rather than merely refining: `h_agl = flying_height - terrain elevation`. + +Note `H = f * S` exactly, so the nominal-scale and focal-length formulations are +algebraically identical. The DEM is the only thing that adds information. + +## The effect is a datum offset, not a slope effect + +Nominal scale underestimates footprint **area** consistently, never overestimates it: + +| scale | n | centroid elev (m) | area correction (footprint-mean sampling) | median | +|---|---|---|---|---| +| 1:12000 | 10 | 578-674 | +0.5% .. +18.2% | +13.4% | +| 1:31680 | 10 | 595-904 | +6.1% .. +27.2% | +15.6% | + +Reported scale is referenced to an elevation above the real valley floor. The issue framed +this as slope making footprints bigger or smaller; the measured bias is bulk and +one-directional, and much larger than the slope term. + +## Sampling point matters — the iteration is not cosmetic + +Elevation at the centroid vs mean elevation under the flat footprint differs by up to +**130 m** (1:31680 frames, where the footprint is 7.2 km across). That moves the 1:12000 +median correction from +16.5% (centroid-only) to +13.4% (one iteration). Hence the two-pass +design: size from centroid elevation, re-sample the mean under that rectangle, resize. + +## Ray-cast adds ~2% on top of ~15% + +Per-corner terrain intersection displaces a corner by roughly `(Δelev / H_agl) × half_width`. +For 1:31680 here: `(200 / 5300) × 3600 ≈ 136 m` on a 7200 m side ≈ 2%. Not worth irregular +geometry plus iterative corner solving while a 15% scale bias is unaddressed. Deferred. + +## DEM source: MRDEM-30 over elevatr + +Both verified this session over the fly AOI (centroids + 4 km buffer, 28.3 × 22.6 km): + +| source | res | clipped size | elev range | dependency | +|---|---|---|---|---| +| MRDEM-30 (NRCan DTM, S3 COG) | 30 m | 306 KB | 566-1520 m | `terra` only | +| elevatr z=10 (AWS terrain mosaic) | 45 m | 242 KB | 544-1539 m | + `elevatr` | + +They agree to within **0.42 percentage points** on the resulting area correction (mean +absolute centroid-elevation difference 3.1 m, max 9 m). So the choice costs no accuracy and +rests on other grounds: MRDEM is the house standard (`flooded::fl_dem_aoi()`), is bare-earth +DTM rather than a mixed-provenance mosaic, adds no dependency, and gives users a documented +source for their own AOIs. Fetch took 4-6 s via `/vsicurl/`. + +MRDEM is a single 84 GB COG in EPSG:3979 covering Canada. `/vsicurl/` range-reads only the +intersecting bytes. **Never reproject it whole** — transform the points to the DEM's CRS for +extraction instead. + +## `footprint_terrain` rather than encoding terrain in `footprint_basis` + +The issue suggested `footprint_basis` carry a terrain-adjusted value. It should not: +`footprint_basis` is already matched by value downstream — the shipped example in +`?fly_footprint` filters `footprint_basis != "unknown_format"` — so appending a suffix like +`"Film - BW (terrain)"` would silently break caller filters. Format basis and terrain basis +are separate facts and get separate columns. + +## Errors Encountered + +| Error | Resolution | +|-------|------------| +| `st_as_sfc` no applicable method for `sfc_POLYGON` | `st_union()` already returns sfc; drop the redundant `st_as_sfc()` | +| `invalid 'na.print' specification` from `print(df, n = 25)` | `n =` is a tibble arg; the object was a data.frame after `as.data.frame()` | +| `[rast] file does not exist` across two Rscript calls | `tempdir()` differs per R session; write probe artifacts to a stable path | + +## Self-review caught a silent frame drop (Phase 2) + +`/code-check`'s subagent rounds were not run — the session instruction bars the Agent tool +unless the user asks. Reviewed against the checklist directly instead; it found one real bug. + +An NA or zero `focal_length`, or an NA `flying_height`, makes the corrected half-side +non-finite. The first cut classified frames on the *inputs* (`is.na(elev)`, `agl <= 0`), +none of which catch that, so the frame was marked `"dem_agl"` — an affirmative claim the +correction applied — and then got an **empty geometry**: + +``` + basis terrain agl area +1 Film - BW dem_agl 1986.315 8807768 +2 Film - BW dem_agl 1938.688 0 <- focal_length NA +3 Film - BW dem_agl NA 0 <- flying_height NA +``` + +Downstream, `fly_warn_unsized()` reports empty geometries as unresolved *recording formats* +and points the user at `format_size` — the wrong diagnosis entirely, for a metadata problem. + +Fix: classify on the half-side actually used (`is.finite(candidate) & candidate > 0`), not +on the inputs feeding it. Three explicit categories — `corrected`, `uncovered`, `unusable` +— each falling back to nominal scale. Two regression tests pin it, one for NA and one for +zero (zero divides rather than propagating NA, so it arrives as `Inf` by a different route). + +Same class as CLAUDE.md's "a zero-length value in a row-builder drops the whole record": +the output looked correct, just shorter. + +## lintr's "unused argument (dem = dem)" is an installed-vs-source artifact + +`lint_package()` went 15 -> 21 after the passthrough. All six new lints are +`unused argument (dem = dem)`, and all six are false — lintr resolves the callee +against the **installed** namespace: + +``` +installed fly version: 0.3.0 +installed fly_footprint formals: centroids_sf, negative_size, format_size +source fly_footprint formals: centroids_sf, negative_size, format_size, dem +``` + +They clear on reinstall. Zero real new lints — the delta is exactly the six artifacts. +This is the CLAUDE.md rule about comparing against the baseline before treating a lint +count as signal, in the variant where the warning reads as a code defect. + +## The first fly_georef passthrough test could not fail + +Written with `expect_silent(...)` plus a row count, it would have passed with `dem` +silently dropped: `fly_georef()` with no images to warp never exposes the GCPs, so the +footprints are unobservable from outside. Its own comment claimed it exercised the +argument path — the exact "fixture cannot reach the failure mode" shape. + +Replaced with a discriminating assertion: strip `flying_height`, and `fly_footprint()` +errors — which it can only do if `dem` actually arrived. Verified both ways by removing +the `dem = dem` passthrough (test fails) and restoring it (test passes). diff --git a/planning/archive/2026-08-issue-9-dem-terrain/progress.md b/planning/archive/2026-08-issue-9-dem-terrain/progress.md new file mode 100644 index 0000000..7eca1d9 --- /dev/null +++ b/planning/archive/2026-08-issue-9-dem-terrain/progress.md @@ -0,0 +1,49 @@ +# Progress — DEM-based terrain-adjusted footprints (#9) + +## Session 2026-08-29 + +- Plan-mode exploration: measured `flying_height` units, quantified the terrain effect, + compared MRDEM-30 against elevatr head to head +- User steered DEM source from elevatr to the federal product used by `flooded`; verified + MRDEM-30 works over the fly AOI and changes the answer by < 0.5 pct points +- Phases approved by user +- Created branch `9-dem-based-terrain-adjusted-footprints` off main +- Next: Phase 1 — MRDEM-30 test fixture + +### Phase 1 — test DEM fixture (done) + +- `data-raw/make_testdata.R` gains an MRDEM-30 `/vsicurl/` clip over centroids + 4 km +- `inst/testdata/dem.tif`: 973x1086 at 30.5 m, 566-1520 m, 306 KB — matches the probe exactly +- Buffer rationale recorded in the script: the widest footprint is 7.2 km, so an edge frame + reaches 3.6 km past the centroid bbox + +### Phase 2 — fly_footprint(dem =) (done) + +- `dem` accepts SpatRaster, path, or /vsicurl/ URL; `terra` to Suggests behind `check_installed()` +- Two-pass sizing; `footprint_terrain` + `height_agl` columns; `dem = NULL` byte-identical +- Measured on the bundled AOI: median +14.1%, range +0.5% to +27.2%, matching the plan +- Restore-the-bug check: collapsing to centroid-only sampling fails 2 tests, with the + patch proven active (max |height_agl - centroid_only| 20+ -> 0) +- Self-review found and fixed a silent frame drop on NA/zero focal length or flying height +- 164 pass, 0 fail, 0 warn; 0 lints; NAMESPACE unchanged at 9 exports + +### Phase 3 — downstream passthrough (done) + +- All six call sites converted, enumerated by grep rather than recall: + fly_coverage, fly_overlap, fly_filter, fly_georef, fly_select_all, fly_select_minimal +- `fly_select` threads dem through both internal helpers — separate call sites, so + passing in one mode proves nothing about the other; both are tested +- Passthrough tests assert the numbers MOVE, not that the argument is tolerated +- 176 pass, 0 fail; 0 real new lints (6 are installed-vs-source artifacts) + +### Phase 4 — docs and release (done) + +- Flat-terrain claims made conditional in fly_footprint (new **Terrain** section), + fly_overlap, fly_georef, and both vignette locations +- Vignette gains a "Terrain-adjusted footprints" section: comparison table (medians + 13.3% and 15.6%), an overlay figure, and the MRDEM-30 recipe for a user's own AOI +- Issue #9 body edited rather than commented: records the datum-offset reframing, the + measured effect, why `footprint_terrain` supersedes the `footprint_basis` suggestion, + and the MRDEM-vs-elevatr comparison. Earlier corrections preserved +- NEWS + version 0.5.0 as the final commit +- 176 pass, 0 fail, 0 warn; examples clean; vignette renders diff --git a/planning/archive/2026-08-issue-9-dem-terrain/review-round1.md b/planning/archive/2026-08-issue-9-dem-terrain/review-round1.md new file mode 100644 index 0000000..dd5e4b5 --- /dev/null +++ b/planning/archive/2026-08-issue-9-dem-terrain/review-round1.md @@ -0,0 +1,237 @@ +# fly #9 — `dem` argument review (round 1) + +Reviewed: `/tmp/cc9/diff.txt` (git diff main...HEAD) against the current tree at +`e5c597c`. Every claim below was reproduced against the real package and the +bundled fixtures (`pkgload::load_all()`, terra 1.9.34); probe scripts are in +`/tmp/cc9/probe*.R`. + +Suite is green as it stands: `[ FAIL 0 | WARN 0 | SKIP 0 | PASS 185 ]`. +The measured claims in NEWS/docs all check out — median area change 14.11%, +max 27.17%, centroid-vs-footprint-mean elevation difference 130 m exactly as +stated. + +## Findings + +--- + +### 1. **[bug]** Partial DEM coverage is silently averaged and reported as fully corrected + +`R/fly_footprint.R:42` (`fly_dem_elevation` → `sample_at`), classification at +`R/fly_footprint.R:317`. + +```r +terra::extract(dem, v, fun = mean, na.rm = TRUE)[, 2] +... +uncovered <- sized & !corrected & is.na(elev) +``` + +`terra::extract()` never returns cells outside the raster extent, and +`na.rm = TRUE` discards NA cells inside it. So a footprint rectangle that is only +*fractionally* covered by the DEM returns a perfectly ordinary non-NA mean — +computed over whatever sliver happened to be there. `is.na(elev)` is the **only** +detection of missing coverage, and it can only be TRUE when the centroid sample +(pass 1) is NA or the rectangle misses the DEM entirely. Partial coverage is +therefore indistinguishable from full coverage, and the guard fails toward +"pass". + +This directly contradicts the shipped contract. `R/fly_footprint.R:201-204`: + +> a DEM that stops short of the frame edges sends those frames down the +> `no_dem_coverage` fallback. Buffer the AOI by at least half the widest footprint. + +and the vignette repeats it verbatim ("a DEM that stops short of the frame edges +sends those frames to the `no_dem_coverage` fallback"). Neither is true. The +"buffer generously" advice is presented as the remedy for a fallback that cannot +fire, so a user who *doesn't* buffer gets silently biased elevations rather than +the documented, visible degradation. + +**Reproduced** (`/tmp/cc9/probe11.R`) — DEM clipped to a 1 km box around one +1:31680 centroid, i.e. **0.98 km² of a 52.4 km² frame (1.9% coverage)**: + +``` +warnings: 0 +footprint_terrain: dem_agl +height_agl: 5180.8 # averaged over 1.9% of the frame +height_agl with full DEM: 5220.4 +area km2 clipped-dem: 59.92 full-dem: 60.84 +``` + +**It is already live in the shipped fixture.** `/tmp/cc9/probe3.R` measures the +valid-DEM fraction under each of the 20 bundled footprints: + +| frame | scale | frac of rectangle with valid DEM | +|---|---|---| +| 14 | 1:31680 | 0.982 | +| 16 | 1:31680 | 0.964 | +| 18 | 1:31680 | 0.981 | +| other 17 | — | 1.000 | + +All three are reported `dem_agl` with no warning. The `data-raw/make_testdata.R` +comment claims the 4 km buffer prevents exactly this ("Without the buffer, edge +frames would sample NA and fall back to nominal scale — which the fixture exists +to exercise the *absence* of"). The arithmetic is off: a **square** of half-side +3621 m reaches 3621 × √2 = **5121 m** at its corners, not 3621 m, so a 4 km +buffer cannot contain it. The buffer needs `half_side * sqrt(2)`, i.e. ≥ 5.2 km +here — and the docs' "half the widest footprint" advice has the same error. + +Note the bundled `dem.tif` bounding box is **47% NA** (the EPSG:3979 crop +rectangle reprojected to 3005 is rotated ~24°, so the box corners are empty) — +NA cells under a footprint are not an exotic case in this fixture. + +**Fix options** — either makes the documented behaviour real: + +```r +sample_at <- function(geom) { + v <- terra::vect(sf::st_transform(geom, dem_crs)) + m <- terra::extract(dem, v, fun = mean, na.rm = TRUE)[, 2] + frac <- terra::extract(!is.na(dem), v, fun = mean, na.rm = TRUE)[, 2] + ifelse(is.na(frac) | frac < 1, NA_real_, m) # or a documented threshold +} +``` + +or drop `na.rm = TRUE` on the polygon pass so any NA under the rectangle +propagates. Whichever you pick, the fixture cannot currently detect the +regression (all its DEM-covered frames sit at ≥0.964), so the test needs a +deliberately-short DEM — `probe11.R` is a ready-made one. + +--- + +### 2. **[bug]** Two of the five passthrough tests cannot detect a dropped `dem` + +`tests/testthat/test-fly_terrain_passthrough.R:38, 42` (`fly_filter`) and +`:50, 57` (`fly_select`, both modes). + +The file's header states the contract explicitly: + +> Each test therefore asserts the numbers actually MOVE, not merely that the +> argument is tolerated. + +For `fly_filter` and `fly_select` that is false. The assertions are non-strict — +`expect_gte(nrow(terr), nrow(flat))`, `expect_lte(nrow(min_terr), nrow(min_flat))`, +`expect_true(all(flat$airp_id %in% terr$airp_id))` — and on the bundled data both +sides are **identical**, so equality is what they actually observe +(`/tmp/cc9/probe4.R`): + +``` +filter flat: 20 terr: 20 +select all flat: 10 terr: 10 +select min flat: 10 terr: 10 +``` + +**Restore-the-bug check** (`/tmp/cc9/probe5.R`) — `fly_footprint` patched in both +`asNamespace("fly")` and `as.environment("package:fly")` to discard `dem`, +simulating a passthrough that was never wired: + +``` +1. Failure fly_coverage passes dem through <- caught +2. Failure fly_overlap passes dem through <- caught +3. Failure fly_georef passes dem through <- caught + fly_filter PASSED <- not caught + fly_select PASSED (both modes) <- not caught +``` + +So deleting `dem = dem` from `R/fly_filter.R:478`, `R/fly_select.R:632` and +`R/fly_select.R:690` would leave the suite green. + +**Fix:** assert on a quantity that actually moves. Total footprint area of the +returned set works for both (`fly_filter` returns points, so compute it from +`fly_footprint()` on the result), or reuse the `fly_georef` trick already in this +file — strip `flying_height` and `expect_error(..., "flying_height")`, which can +only pass if `dem` reached `fly_footprint()`. + +--- + +### 3. **[fragile]** One NA cell under the centroid condemns the whole frame, with a misleading warning + +`R/fly_footprint.R:44-51`. + +Pass 1 samples a **single point**. If that one 30 m cell is NA, `ok` is FALSE, the +footprint-mean pass never runs, and the frame is classified `no_dem_coverage` and +falls back to nominal scale — even if 99% of its rectangle has good data. Given +the bundled DEM's bbox is 47% NA, that is a reachable state, not a hypothetical. + +This is the safe failure direction (fallback, not silent bias), so it is not +finding 1's severity. But the warning text — "frames fall outside the DEM's +coverage" — misdescribes what happened, which sends the user to widen a DEM that +is already wide enough. Consider falling through to the polygon pass whenever the +provisional rectangle is non-empty, and only declaring `no_dem_coverage` on the +polygon result. + +--- + +### 4. **[fragile]** A frame with an empty geometry can be labelled `footprint_terrain = "nominal_scale"` + +`R/fly_footprint.R:296` seeds `terrain` from `is.na(width_in)`, but the +sized/unsized decision at `:301` is `!is.na(half_side)`, and +`half_side = width_in * scale_num * ...`. A frame with a resolvable `media` but an +unparseable `scale` has `width_in` non-NA and `scale_num` NA — so it gets an empty +geometry while claiming a terrain treatment was applied. + +Documented contract at `R/fly_footprint.R:175`: `NA` = "no footprint to place". + +**Reproduced** (`/tmp/cc9/probe12.R`): + +``` + basis terrain agl empty +1 Film - BW dem_agl 1986.315 FALSE +2 Film - BW nominal_scale NA TRUE <- empty geometry, non-NA terrain +``` + +**Fix:** seed from the half-side, which is the thing that decides whether a +polygon exists — `terrain <- ifelse(is.na(half_side), NA_character_, "nominal_scale")` +(move the line below `half_side`). Untested today in either direction. + +--- + +### 5. **[fragile]** Examples and vignette call a Suggests-only package unconditionally + +`R/fly_footprint.R:223` (`@examples`) and `vignettes/airphoto-selection.Rmd` +chunks `terrain-compare`, `fig-terrain`, `terrain-basis`. + +`terra` is in Suggests (correctly), and `rlang::check_installed()` **errors** in a +non-interactive session. So `R CMD check` and the pkgdown build fail outright on +any environment without terra, rather than skipping. The tests get this right via +`skip_if_no_terra()`; the examples and vignette do not. + +`@examplesIf requireNamespace("terra", quietly = TRUE)` on the terrain block, and +`eval = requireNamespace("terra", quietly = TRUE)` on those three chunks, closes +it. Low priority while CI installs Suggests, but it is a hard failure when it +fires, not a NOTE. + +--- + +## Checked and clean + +- **Vector alignment through the DEM path.** `elev[sized] <- fly_dem_elevation(...)`, + `elev[ok] <- sample_at(provisional[ok])`, `half_side[corrected] <- candidate[corrected]` + are all consistently subset; `terra::extract()` was verified to return exactly + one row per input geometry in input order, including geometries entirely + outside the raster (`/tmp/cc9/probe1.R`, `nrow(res) = 3` for 3 polygons). +- **Zero-length / NA classification.** `corrected` / `uncovered` / `unusable` + partition `sized` exhaustively; `NaN` from `mean(numeric(0))` is caught by + `is.na()`; NA `flying_height`, NA and zero `focal_length` all land in `unusable` + and keep their nominal footprint, as the tests assert. +- **Unit arithmetic.** `width_in * (agl / focal_m) * 0.0254 / 2` with `focal_m = + focal_length/1000` reduces to the nominal form when `agl/focal_m == scale_num`. + Verified numerically: 2591 m ASL − ~620 m terrain over 0.153 m → 12882 vs a + nominal 12000, i.e. +7.3% linear / +15% area, matching the measured median. +- **CRS handling.** `sf::st_crs(terra::crs(dem))` then transform-per-sample is + correct; `fly_rectangles()` stamps 3005 and the result is transformed back to + `input_crs`. No `st_join(largest=)` or bbox-corner reprojection anywhere. +- **terra traps.** `data-raw/make_testdata.R:617` correctly uses + `terra::minmax(dem_clip, compute = TRUE)`. No `%in%` on a SpatRaster (the + import-vs-attach S4 dispatch trap), no bare `terra::freq()`. +- **Backwards compatibility.** `dem` is appended last in every signature + (`fly_footprint`, `fly_coverage`, `fly_overlap`, `fly_filter`, `fly_select`, + `fly_georef`), so no positional caller breaks. `fly_select_all(photos_sf, + aoi_sf, dem)` and `fly_select_minimal(..., dem)` match their definitions + positionally. +- **`expect_match(..., all = FALSE)` on a possibly-empty vector** — every use is + preceded by `expect_gt(length(w), 0)`. Correct. +- **`fly_footprint(dem = NULL)` is byte-identical to the old path** — the + refactor into `fly_rectangles()` preserves the NA-half-side → empty-polygon + contract, and the regression test asserts it. +- **`data-raw/make_testdata.R` reproducibility.** `terra::crop()` does not mask + (verified, `/tmp/cc9/probe8.R`); the committed `dem.tif`'s non-rectangular + valid region is the reprojected EPSG:3979 crop rectangle, consistent with the + script. No generator/artifact drift. diff --git a/planning/archive/2026-08-issue-9-dem-terrain/task_plan.md b/planning/archive/2026-08-issue-9-dem-terrain/task_plan.md new file mode 100644 index 0000000..ba143ee --- /dev/null +++ b/planning/archive/2026-08-issue-9-dem-terrain/task_plan.md @@ -0,0 +1,60 @@ +# Task: DEM-based terrain-adjusted footprints (#9) + +`fly_footprint()` sizes every frame from the reported nominal scale, which assumes flat +ground beneath the aircraft. Issue #9 asked for a DEM option and offered three approaches, +noting the choice between them was the first decision. + +Measurement reframed it. Probing the bundled AOI (20 frames, 1968, Upper Bulkley near +Houston) established two facts the issue did not have: + +1. **`flying_height` is metres ASL, not height above ground** — values are round feet + (8000, 8500, 19500, 20000 ft). Subtracting terrain elevation is what turns it into the + height a footprint actually scales with. +2. **Nominal scale underestimates footprint area by a median ~15%, range 0.5-27%, always + in the same direction.** Not a slope effect — a datum offset. Reported scale is + referenced to an elevation above the real valley floor, so every footprint is larger + than nominal. + +Approach chosen (user-approved): **true-scale rectangle**. Size each frame from +`h_agl = flying_height - terrain elevation` instead of nominal scale. Geometry stays a +rectangle, so downstream consumers work unchanged. Per-corner ray-cast measures ~2% on top +of the ~15% and costs irregular geometry — deferred. + +DEM source: **MRDEM-30**, NRCan's 30 m bare-earth DTM, the same product +`flooded::fl_dem_aoi()` uses. Public S3 COG, no auth, needs no dependency beyond `terra`. + +## Phase 1: Test DEM fixture + +- [x] Extend `data-raw/make_testdata.R` with an MRDEM-30 clip over centroids + 4 km buffer +- [x] Write `inst/testdata/dem.tif` (INT2S, DEFLATE, tiled); confirm size and 566-1520 m range +- [x] Document the source and regeneration in the script header + +## Phase 2: `fly_footprint(dem =)` + +- [x] Failing tests first: flat unchanged, terrain enlarges, empty geoms skipped, each guard +- [x] `terra` to Suggests with `check_installed()` guard +- [x] Two-pass sizing; `footprint_terrain` and `height_agl` columns +- [x] Guards for missing columns, no DEM coverage, non-positive `h_agl` +- [x] Confirm `dem = NULL` output is identical to current — assert, don't assume + +## Phase 3: Downstream passthrough + +- [x] `dem = NULL` on `fly_coverage()`, `fly_overlap()`, `fly_select()`, `fly_filter()`, `fly_georef()` +- [x] Tests that a DEM reaching each consumer changes its numbers +- [x] Verify all six call sites are covered — enumerate, don't recall + +## Phase 4: Docs and release + +- [x] Rewrite the flat-terrain paragraphs as conditional in all four locations +- [x] Document DEM sources: MRDEM-30 default, LidarBC via `stac-dem-bc`, BC TRIM 25 m via `bcdata` +- [x] Vignette section showing flat vs terrain-adjusted on the bundled AOI +- [x] Edit issue #9 body: `footprint_terrain` supersedes the `footprint_basis` suggestion; + effect is a datum offset rather than a slope effect +- [x] `NEWS.md`, version 0.5.0 as the final commit + +## Validation + +- [x] Tests pass +- [x] `/code-check` clean on each commit +- [x] PWF checkboxes match landed work +- [ ] `/planning-archive` on completion diff --git a/tests/testthat/setup.R b/tests/testthat/setup.R index 987ad2e..3d69915 100644 --- a/tests/testthat/setup.R +++ b/tests/testthat/setup.R @@ -28,3 +28,31 @@ mixed_media_fixture <- function() { ) ) } + + +# Skip a terrain test when terra is unavailable. +# +# `terra` is in Suggests, not Imports — the DEM path is optional. A test that +# needs it must skip rather than fail on an install that reasonably lacks it. +skip_if_no_terra <- function() { + testthat::skip_if_not_installed("terra") +} + +# Centroids carrying the two fields the DEM path needs, plus a media value the +# format table cannot resolve — so one frame arrives with an empty geometry and +# the terrain code must leave it alone rather than sample a DEM under it. +terrain_fixture <- function() { + sf::st_sf( + airp_id = 1:3, + scale = c("1:12000", "1:12000", "1:12000"), + media = c("Film - BW", "Film - BW", "Digital - Colour"), + focal_length = c(153, 153, 153), + flying_height = c(2591, 2591, 2591), + geometry = sf::st_sfc( + sf::st_point(c(-126.60, 54.40)), + sf::st_point(c(-126.58, 54.40)), + sf::st_point(c(-126.56, 54.40)), + crs = 4326 + ) + ) +} diff --git a/tests/testthat/test-fly_footprint.R b/tests/testthat/test-fly_footprint.R index 58626cf..ee6f37e 100644 --- a/tests/testthat/test-fly_footprint.R +++ b/tests/testthat/test-fly_footprint.R @@ -117,3 +117,533 @@ test_that("fly_footprint leaves film-only input unchanged", { expected_side <- as.numeric(sub("1:", "", centroids$scale[1])) * 9 * 0.0254 expect_equal(area_m2, expected_side^2, tolerance = 0.01) }) + + +# --- Terrain-adjusted footprints (#9) ------------------------------------ + +test_that("fly_footprint with dem = NULL is unchanged", { + centroids <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE) + flat <- fly_footprint(centroids) + explicit_null <- fly_footprint(centroids, dem = NULL) + + expect_equal(sf::st_geometry(flat), sf::st_geometry(explicit_null)) + expect_equal(flat$footprint_basis, explicit_null$footprint_basis) + # The terrain columns exist either way, so a caller's column handling does not + # change depending on whether a DEM was supplied. + expect_equal(flat$footprint_terrain, rep("nominal_scale", nrow(centroids))) + expect_true(all(is.na(flat$height_agl))) +}) + +test_that("fly_footprint terrain correction enlarges footprints as measured", { + skip_if_no_terra() + centroids <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE) + flat <- fly_footprint(centroids) + terr <- fly_footprint(centroids, dem = testdata_path("dem.tif")) + + pct <- 100 * (as.numeric(sf::st_area(sf::st_transform(terr, 3005))) / + as.numeric(sf::st_area(sf::st_transform(flat, 3005))) - 1) + + # Every frame grows. Reported scale on this AOI is referenced to an elevation + # above the real valley floor, so the bias is one-directional — a correction + # that shrank a footprint here would mean the sampling is wrong. + expect_true(all(pct > 0)) + expect_gt(stats::median(pct), 10) + expect_lt(max(pct), 30) + + expect_equal(terr$footprint_terrain, rep("dem_agl", nrow(centroids))) + expect_true(all(terr$height_agl > 0)) + # h_agl is flying height above ground, so strictly below flying height ASL. + expect_true(all(terr$height_agl < centroids$flying_height)) +}) + +test_that("fly_footprint samples the footprint mean, not just the centroid", { + skip_if_no_terra() + # The two-pass iteration is load-bearing, not cosmetic: centroid and + # footprint-mean elevation differ by up to 130 m on the wide 1:31680 frames. + # If this ever collapses to centroid-only sampling, the areas move. + centroids <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE) + dem <- terra::rast(testdata_path("dem.tif")) + terr <- fly_footprint(centroids, dem = dem) + + pts <- terra::extract(dem, terra::vect(sf::st_transform(centroids, 3005)))[, 2] + centroid_only_agl <- centroids$flying_height - pts + + expect_false(isTRUE(all.equal(terr$height_agl, centroid_only_agl))) + expect_gt(max(abs(terr$height_agl - centroid_only_agl)), 20) +}) + +test_that("fly_footprint accepts a dem as path or SpatRaster", { + skip_if_no_terra() + centroids <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE) + by_path <- fly_footprint(centroids, dem = testdata_path("dem.tif")) + by_rast <- fly_footprint(centroids, dem = terra::rast(testdata_path("dem.tif"))) + expect_equal(sf::st_geometry(by_path), sf::st_geometry(by_rast)) +}) + +test_that("fly_footprint leaves unsized frames empty rather than sampling them", { + skip_if_no_terra() + photos <- terrain_fixture() + fp <- suppressWarnings(fly_footprint(photos, dem = testdata_path("dem.tif"))) + + unknown <- fp[fp$footprint_basis == "unknown_format", ] + expect_equal(nrow(unknown), 1) + expect_true(sf::st_is_empty(sf::st_geometry(unknown))) + # No footprint means no terrain treatment to report — not a claim that one + # was applied, and not a number a caller could mistake for a real height. + expect_true(is.na(unknown$footprint_terrain)) + expect_true(is.na(unknown$height_agl)) + + film <- fp[fp$footprint_basis == "Film - BW", ] + expect_equal(film$footprint_terrain, rep("dem_agl", 2)) +}) + +test_that("fly_footprint errors when a dem is given without the fields it needs", { + skip_if_no_terra() + photos <- terrain_fixture() + + no_height <- photos + no_height$flying_height <- NULL + expect_error(fly_footprint(no_height, dem = testdata_path("dem.tif")), + "flying_height") + + no_focal <- photos + no_focal$focal_length <- NULL + expect_error(fly_footprint(no_focal, dem = testdata_path("dem.tif")), + "focal_length") +}) + +test_that("fly_footprint falls back to nominal scale outside DEM coverage", { + skip_if_no_terra() + photos <- terrain_fixture() + # Move one frame far outside the bundled DEM's extent. Falling back with a + # warning is the chosen failure direction: a frame we cannot correct is still + # a frame, and dropping it would be indistinguishable from an unsized one. + far <- photos[1:2, ] + sf::st_geometry(far)[2] <- sf::st_sfc(sf::st_point(c(-120.0, 50.0)), crs = 4326) + + w <- character() + fp <- withCallingHandlers( + fly_footprint(far, dem = testdata_path("dem.tif")), + warning = function(x) { + w <<- c(w, conditionMessage(x)) + invokeRestart("muffleWarning") + } + ) + expect_gt(length(w), 0) + expect_match(w, "coverage", all = FALSE) + + expect_equal(fp$footprint_terrain, c("dem_agl", "no_dem_coverage")) + expect_true(is.na(fp$height_agl[2])) + expect_false(sf::st_is_empty(sf::st_geometry(fp)[2])) + + # The fallback frame is sized exactly as the flat path would size it. + flat <- fly_footprint(far) + expect_equal( + as.numeric(sf::st_area(sf::st_transform(fp[2, ], 3005))), + as.numeric(sf::st_area(sf::st_transform(flat[2, ], 3005))) + ) +}) + +test_that("fly_footprint falls back when terrain sits above the aircraft", { + skip_if_no_terra() + photos <- terrain_fixture()[1:2, ] + # Terrain here is ~600 m. A flying height below that is bad data — most + # likely feet recorded as metres — and yields a non-positive height above + # ground, which would otherwise produce a zero or inverted footprint. + photos$flying_height <- c(2591, 100) + + w <- character() + fp <- withCallingHandlers( + fly_footprint(photos, dem = testdata_path("dem.tif")), + warning = function(x) { + w <<- c(w, conditionMessage(x)) + invokeRestart("muffleWarning") + } + ) + expect_gt(length(w), 0) + expect_match(w, "above", all = FALSE) + + expect_equal(fp$footprint_terrain, c("dem_agl", "nominal_scale")) + expect_false(sf::st_is_empty(sf::st_geometry(fp)[2])) +}) + +test_that("fly_footprint falls back on unusable height or focal metadata", { + skip_if_no_terra() + # An NA focal length or flying height makes the corrected half-side + # non-finite. Left unchecked that becomes an empty geometry, which is + # indistinguishable from an unresolved recording format — so the frame would + # disappear under a warning pointing at `format_size` rather than at its own + # metadata. Every such frame must keep its nominal-scale footprint. + centroids <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE)[1:3, ] + centroids$focal_length[2] <- NA + centroids$flying_height[3] <- NA + + w <- character() + fp <- withCallingHandlers( + fly_footprint(centroids, dem = testdata_path("dem.tif")), + warning = function(x) { + w <<- c(w, conditionMessage(x)) + invokeRestart("muffleWarning") + } + ) + expect_gt(length(w), 0) + expect_match(w, "focal_length", all = FALSE) + + expect_equal(sum(sf::st_is_empty(sf::st_geometry(fp))), 0) + expect_equal(fp$footprint_terrain, c("dem_agl", "nominal_scale", "nominal_scale")) + expect_true(all(is.na(fp$height_agl[2:3]))) + + flat <- fly_footprint(centroids) + expect_equal( + as.numeric(sf::st_area(sf::st_transform(fp[2:3, ], 3005))), + as.numeric(sf::st_area(sf::st_transform(flat[2:3, ], 3005))) + ) +}) + +test_that("fly_footprint falls back on a zero focal length", { + skip_if_no_terra() + # Zero divides rather than propagating NA, so it reaches the half-side as Inf + # and needs the same finite check. + centroids <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE)[1:2, ] + centroids$focal_length[2] <- 0 + fp <- suppressWarnings(fly_footprint(centroids, dem = testdata_path("dem.tif"))) + expect_equal(sum(sf::st_is_empty(sf::st_geometry(fp))), 0) + expect_equal(fp$footprint_terrain, c("dem_agl", "nominal_scale")) +}) + +test_that("fly_footprint reports how much of each footprint the DEM covered", { + skip_if_no_terra() + centroids <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE) + terr <- fly_footprint(centroids, dem = testdata_path("dem.tif")) + + expect_true(all(terr$dem_coverage > 0 & terr$dem_coverage <= 1)) + # The bundled DEM is buffered past the corner of the widest footprint, so + # every frame is essentially fully described — the shortfall below is real + # missing cells left by reprojection, not a counting artifact, and is four + # hundredths of a percent at worst. This asserts the guard stays quiet on + # good data; the AOI-clipped test below is what proves it can fire at all. + expect_gt(min(terr$dem_coverage), 0.99) + expect_silent(fly_footprint(centroids, dem = testdata_path("dem.tif"))) + + flat <- fly_footprint(centroids) + expect_true(all(is.na(flat$dem_coverage))) +}) + +test_that("fly_footprint warns when a footprint is materially off the DEM", { + skip_if_no_terra() + dem <- terra::rast(testdata_path("dem.tif")) + + # Put a wide frame on the westernmost cell that still carries data: the + # centroid samples fine, so this is NOT the no_dem_coverage path — half the + # footprint simply hangs over ground the DEM does not describe, and the mean + # comes from the covered half alone. + with_data <- which(!is.na(terra::values(dem))) + xy <- terra::xyFromCell(dem, with_data) + edge <- sf::st_sfc(sf::st_point(xy[which.min(xy[, 1]), ]), crs = 3005) + + photos <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE)[1, ] + photos$scale <- "1:31680" + sf::st_geometry(photos) <- sf::st_transform(edge, 4326) + + w <- character() + fp <- withCallingHandlers( + fly_footprint(photos, dem = dem), + warning = function(x) { + w <<- c(w, conditionMessage(x)) + invokeRestart("muffleWarning") + } + ) + expect_gt(length(w), 0) + expect_match(w, "covered by the DEM", all = FALSE) + + # Still corrected — the partial mean is the best estimate available, and + # discarding a frame over it would lose more than it protects. + expect_equal(fp$footprint_terrain, "dem_agl") + expect_false(sf::st_is_empty(sf::st_geometry(fp))) + expect_lt(fp$dem_coverage, 0.95) +}) + +test_that("fly_footprint reports no terrain treatment for an unparseable scale", { + skip_if_no_terra() + # A resolvable `media` with a `scale` that will not parse still leaves a frame + # with no footprint. Keying the terrain column off the recording format alone + # labelled it "nominal_scale", claiming a treatment for a frame that has no + # geometry to treat. + photos <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE)[1:2, ] + photos$scale[2] <- "not-a-scale" + fp <- suppressWarnings(fly_footprint(photos)) + + expect_true(sf::st_is_empty(sf::st_geometry(fp)[2])) + expect_true(is.na(fp$footprint_terrain[2])) + expect_false(is.na(fp$footprint_terrain[1])) +}) + +test_that("fly_footprint tolerates a centroid on a DEM hole", { + skip_if_no_terra() + # The mean is taken over the whole footprint, so a single missing cell beneath + # the centroid says nothing about whether the frame can be corrected. Punch a + # hole at one centroid and the frame must still come back corrected. + centroids <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE)[1, ] + dem <- terra::rast(testdata_path("dem.tif")) + cell <- terra::cellFromXY(dem, sf::st_coordinates(sf::st_transform(centroids, 3005))) + dem[cell] <- NA + expect_true(is.na(terra::extract(dem, + terra::vect(sf::st_transform(centroids, 3005)))[, 2])) + + fp <- fly_footprint(centroids, dem = dem) + expect_equal(fp$footprint_terrain, "dem_agl") + expect_gt(fp$height_agl, 0) +}) + +test_that("fly_footprint detects a footprint running past the DEM's extent", { + skip_if_no_terra() + # The two ways a footprint can be short are NOT equivalent, and only one of + # them leaves NA cells to count. Ground beyond the raster's *extent* yields no + # row from terra::extract() at all, so measuring coverage as the non-NA share + # of returned cells calls a truncated footprint fully covered — an affirmative + # claim, and worse than saying nothing. + # + # A DEM cropped to an AOI is exactly this shape: no NA interior, it just + # stops. That is the ordinary way a user obtains one, via fl_dem_aoi() or + # `bcdata get-dem`, so this is the common case rather than the exotic one. + centroids <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE) + full <- terra::rast(testdata_path("dem.tif")) + tight <- terra::crop( + full, + terra::vect(sf::st_union(sf::st_buffer(sf::st_transform(centroids, 3005), 500))), + snap = "out" + ) + # The fixture must actually reach the failure mode, which is ground BEYOND + # the raster's extent — not NA cells within it. Assert that directly: some + # footprint must extend past the DEM's extent, or this test is checking the + # sliver case the previous test already covers and nothing new. + # + # (The earlier version of this premise compared `tight` against + # `terra::crop(full, tight)`, which is `tight` by construction. Both sides + # were equal for every possible fixture, so it could not fail.) + fp_flat <- fly_footprint(centroids) + fp_bbox <- sf::st_bbox(sf::st_transform(fp_flat, 3005)) + dem_ext <- terra::ext(tight) + expect_true( + fp_bbox[["xmin"]] < dem_ext[1] || fp_bbox[["xmax"]] > dem_ext[2] || + fp_bbox[["ymin"]] < dem_ext[3] || fp_bbox[["ymax"]] > dem_ext[4] + ) + + w <- character() + fp <- withCallingHandlers( + fly_footprint(centroids, dem = tight), + warning = function(x) { + w <<- c(w, conditionMessage(x)) + invokeRestart("muffleWarning") + } + ) + expect_gt(length(w), 0) + expect_match(w, "covered by the DEM", all = FALSE) + expect_lt(min(fp$dem_coverage), 0.95) + + # And the well-buffered DEM must NOT warn, or the guard is just noise. + expect_silent(fly_footprint(centroids, dem = full)) +}) + +test_that("fly_footprint reports zero coverage, not NA, for a frame off the DEM", { + skip_if_no_terra() + # NA would mean "not measured". Zero is what was measured, and it is what the + # documented `dem_coverage` filter needs in order to exclude the frame. + photos <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE)[1, ] + sf::st_geometry(photos) <- sf::st_sfc(sf::st_point(c(-120, 50)), crs = 4326) + fp <- suppressWarnings(fly_footprint(photos, dem = testdata_path("dem.tif"))) + + expect_equal(fp$footprint_terrain, "no_dem_coverage") + expect_equal(fp$dem_coverage, 0) +}) + +test_that("fly_footprint measures coverage correctly on a geographic DEM", { + skip_if_no_terra() + # Every other DEM in this suite is EPSG:3005 or a crop of it, so the + # reprojection branch in fly_dem_sample() is never executed by them and a + # units error there is invisible. Coverage compares a footprint's area + # against the DEM's cell size, and both have to be in the DEM's own units: + # st_area() on a geographic CRS returns geodesic m2 while terra::res() + # returns degrees, which reported ~1e-10 coverage for fully-covered frames + # and warned on all of them. + centroids <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE) + geo <- terra::project(terra::rast(testdata_path("dem.tif")), "EPSG:4326") + expect_true(sf::st_is_longlat(sf::st_crs(terra::crs(geo)))) + + fp <- fly_footprint(centroids, dem = geo) + expect_true(all(fp$dem_coverage > 0.95)) + expect_true(all(fp$dem_coverage <= 1)) + expect_equal(fp$footprint_terrain, rep("dem_agl", nrow(centroids))) + + # And it must agree with the projected DEM it was made from, since the two + # describe the same ground. + proj <- fly_footprint(centroids, dem = testdata_path("dem.tif")) + expect_equal(fp$height_agl, proj$height_agl, tolerance = 0.01) +}) + +test_that("fly_footprint iterates the sampling window, not just the elevation", { + skip_if_no_terra() + # The window the DEM is averaged over is itself what the correction changes, + # so the second pass measures over the corrected rectangle. Collapsing it to + # one pass left the whole suite green, which made the iteration untested + # rather than merely small. It IS small — under 0.5% of area against the + # correction's own 14% — so assert it on elevation, where it is legible. + centroids <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE) + dem <- terra::rast(testdata_path("dem.tif")) + fp <- fly_footprint(centroids, dem = dem) + + # Reproduce pass one alone: the mean under the NOMINAL rectangle. + nominal_half <- 9 * as.numeric(sub("1:", "", centroids$scale)) * 0.0254 / 2 + pts <- sf::st_transform(centroids, 3005) + one_pass <- vapply(seq_len(nrow(centroids)), function(i) { + rect <- sf::st_buffer(pts[i, ], nominal_half[i], endCapStyle = "SQUARE") + mean(terra::extract(dem, terra::vect(rect))[, 2], na.rm = TRUE) + }, numeric(1)) + + two_pass_elev <- centroids$flying_height - fp$height_agl + expect_false(isTRUE(all.equal(two_pass_elev, one_pass))) + expect_gt(max(abs(two_pass_elev - one_pass)), 5) +}) + +test_that("fly_footprint rejects a DEM with no CRS", { + skip_if_no_terra() + # Without this the failure surfaces from inside sf as "invalid crs:", naming + # neither the argument nor the package. + dem <- terra::rast(testdata_path("dem.tif")) + terra::crs(dem) <- "" + centroids <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE)[1, ] + expect_error(fly_footprint(centroids, dem = dem), "no CRS") +}) + +test_that("fly_footprint reports full coverage on a DEM that is entirely valid", { + skip_if_no_terra() + # The measurement has to be right on its own terms before missing data enters + # the picture. A DEM with no NA cell anywhere and room well beyond every + # footprint must report coverage of exactly 1 and warn about nothing. + # + # Counting the footprint's area in cell units against a count of cell centres + # compares two different measurements, and fails here rather than on anything + # to do with coverage: it reported 91% on the coarse grid below and warned + # that 3 of 3 frames were under-covered, on a raster with nothing missing. + # + # Two resolutions because the error scales as 2/k for a footprint k cells + # wide — a fine grid hides it, which is why the shipped 30 m fixture could + # not reach this. + centroids <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE)[1:3, ] + bb <- sf::st_bbox(sf::st_transform(centroids, 3005)) + + for (cell in c(30, 900)) { + r <- terra::rast( + xmin = bb[["xmin"]] - 30000, xmax = bb[["xmax"]] + 30000, + ymin = bb[["ymin"]] - 30000, ymax = bb[["ymax"]] + 30000, + resolution = cell, crs = "EPSG:3005" + ) + terra::values(r) <- 700 + expect_equal(as.numeric(terra::global(r, function(x) sum(is.na(x)))[1, 1]), 0) + + fp <- fly_footprint(centroids, dem = r) + expect_equal(fp$dem_coverage, rep(1, 3), info = paste("cell size", cell)) + expect_silent(fly_footprint(centroids, dem = r)) + } +}) + +test_that("fly_footprint handles anisotropic DEM cells", { + skip_if_no_terra() + # Non-square cells are ordinary in a geographic CRS away from the equator. + # The coverage denominator must account for both dimensions, not assume a + # square cell. + centroids <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE)[1:2, ] + bb <- sf::st_bbox(sf::st_transform(centroids, 3005)) + r <- terra::rast( + xmin = bb[["xmin"]] - 20000, xmax = bb[["xmax"]] + 20000, + ymin = bb[["ymin"]] - 20000, ymax = bb[["ymax"]] + 20000, + resolution = c(120, 904), crs = "EPSG:3005" + ) + terra::values(r) <- 700 + fp <- fly_footprint(centroids, dem = r) + expect_equal(fp$dem_coverage, rep(1, 2)) +}) + +test_that("fly_footprint does not size its coverage grid to the span of the photo set", { + skip_if_no_terra() + # Counting the coverage denominator on one grid spanning every frame sizes it + # to the GAP between frames, not to the frames. Two photos 700 km apart made + # that 243 million cells where the same two counted separately need 16 + # thousand — a 4 GB allocation for a correct answer. + # + # Asserted on the grid itself, not on elapsed time. The first version of this + # test used `expect_lt(elapsed, 10)`, and the defect runs in 1.0 s against the + # fix's 0.18 s — so it passed with a tenfold margin on the very thing it was + # written to catch, and no threshold separates them without being CI jitter. + # Every other assertion in it passed too, because the union grid produces the + # *right* number; it just allocates absurdly to get there. + centroids <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE)[1:2, ] + sf::st_geometry(centroids)[2] <- sf::st_sfc(sf::st_point(c(-120, 50)), crs = 4326) + dem <- terra::rast(testdata_path("dem.tif")) + + # The fixture must be able to expose the defect: the grid spanning both + # frames is enormous, so a union-sized allocation would stand out at once. + # Measured before the mock is installed, or this call records itself. + flat <- fly_footprint(centroids) + union_cells <- prod(dim(fly_dem_grid( + dem, sf::st_transform(flat, sf::st_crs(terra::crs(dem))) + ))[1:2]) + expect_gt(union_cells, 1e8) + + # Record every grid the real call asks for. + sizes <- c() + real_grid <- fly_dem_grid + testthat::local_mocked_bindings( + fly_dem_grid = function(dem, geom) { + g <- real_grid(dem, geom) + sizes <<- c(sizes, prod(dim(g)[1:2])) + g + } + ) + fp <- suppressWarnings(fly_footprint(centroids, dem = dem)) + + # No grid built during the call may approach it. Each is one footprint. + expect_gt(length(sizes), 0) + expect_lt(max(sizes), 1e6) + expect_lt(max(sizes), union_cells / 100) + + expect_equal(fp$footprint_terrain, c("dem_agl", "no_dem_coverage")) + expect_equal(fp$dem_coverage, c(1, 0)) +}) + +test_that("fly_footprint reports coverage of the footprint it actually returns", { + skip_if_no_terra() + # Terrain above the aircraft gives a negative height above ground, so the + # second pass measures over a square sized from that — much smaller than the + # nominal footprint the frame falls back to. Reading the second pass's + # coverage there describes a rectangle the caller never receives: measured + # 100% for a footprint that was 30% covered, which is precisely the frame the + # documented `dem_coverage` filter exists to exclude. + dem <- terra::rast(testdata_path("dem.tif")) + with_data <- which(!is.na(terra::values(dem))) + xy <- terra::xyFromCell(dem, with_data) + west <- which.min(xy[, 1]) + # Far enough inside that the small second-pass square is fully covered, close + # enough that the nominal footprint is not. + pt <- sf::st_sfc(sf::st_point(c(xy[west, 1] + 1200, xy[west, 2])), crs = 3005) + + photos <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE)[1, ] + photos$scale <- "1:31680" + photos$flying_height <- 100 # below the terrain -> unusable + sf::st_geometry(photos) <- sf::st_transform(pt, 4326) + + fp <- suppressWarnings(fly_footprint(photos, dem = dem)) + expect_equal(fp$footprint_terrain, "nominal_scale") + + # Ground truth measured against the geometry that was returned. + g <- sf::st_transform(fp, sf::st_crs(terra::crs(dem))) + vals <- terra::extract(dem, terra::vect(g))[, 2] + tmpl <- terra::rast(terra::align(terra::ext(terra::vect(g)), dem), + resolution = terra::res(dem), crs = terra::crs(dem)) + terra::values(tmpl) <- 1L + truth <- sum(!is.na(vals)) / + sum(!is.na(terra::extract(tmpl, terra::vect(g))[, 2])) + + expect_lt(truth, 0.95) # the fixture must reach the failure mode + expect_equal(fp$dem_coverage, truth, tolerance = 1e-6) +}) diff --git a/tests/testthat/test-fly_terrain_passthrough.R b/tests/testthat/test-fly_terrain_passthrough.R new file mode 100644 index 0000000..1210b93 --- /dev/null +++ b/tests/testthat/test-fly_terrain_passthrough.R @@ -0,0 +1,131 @@ +# A `dem` argument that is accepted but never reaches fly_footprint() would be +# invisible: every one of these functions returns perfectly plausible numbers +# either way. +# +# Asserting that the numbers move is only a real check where they demonstrably +# do. On the bundled data fly_filter() keeps 20 of 20 and fly_select() picks the +# same frames with or without a DEM, so `expect_gte`/`expect_lte` there hold +# for both implementations and detect nothing. Those two therefore also strip +# `flying_height`, which makes fly_footprint() error — something it can only do +# if `dem` actually arrived. +# +# Verified by patching fly_footprint() to discard `dem`: every test below then +# fails. + +test_that("fly_coverage passes dem through to the footprints", { + skip_if_no_terra() + centroids <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE) + aoi <- sf::st_read(testdata_path("aoi.gpkg"), quiet = TRUE) + + flat <- fly_coverage(centroids, aoi, by = "scale") + terr <- fly_coverage(centroids, aoi, by = "scale", dem = testdata_path("dem.tif")) + + expect_equal(terr$scale, flat$scale) + expect_true(all(terr$covered_km2 >= flat$covered_km2)) + expect_true(any(terr$covered_km2 > flat$covered_km2)) +}) + +test_that("fly_overlap passes dem through to the footprints", { + skip_if_no_terra() + centroids <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE) + photos <- centroids[centroids$scale == "1:12000", ] + + flat <- fly_overlap(photos) + terr <- fly_overlap(photos, dem = testdata_path("dem.tif")) + + # Larger footprints overlap more, and can bring new pairs into contact. + expect_gte(nrow(terr), nrow(flat)) + expect_gt(sum(terr$overlap_km2), sum(flat$overlap_km2)) +}) + +test_that("fly_filter passes dem through to the footprints", { + skip_if_no_terra() + centroids <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE) + aoi <- sf::st_read(testdata_path("aoi.gpkg"), quiet = TRUE) + + flat <- fly_filter(centroids, aoi, method = "footprint") + terr <- fly_filter(centroids, aoi, method = "footprint", + dem = testdata_path("dem.tif")) + # Terrain sizing only grows footprints here, so the kept set cannot shrink. + expect_gte(nrow(terr), nrow(flat)) + expect_true(all(flat$airp_id %in% terr$airp_id)) + + # The centroid method never builds a footprint, so a dem must not change it. + expect_equal( + fly_filter(centroids, aoi, method = "centroid")$airp_id, + fly_filter(centroids, aoi, method = "centroid", + dem = testdata_path("dem.tif"))$airp_id + ) + + # The comparison above cannot fail on this data — every frame is kept either + # way — so it does not establish that `dem` arrived. This does. + no_height <- centroids + no_height$flying_height <- NULL + expect_error( + fly_filter(no_height, aoi, method = "footprint", dem = testdata_path("dem.tif")), + "flying_height" + ) +}) + +test_that("fly_select passes dem through in both modes", { + skip_if_no_terra() + centroids <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE) + aoi <- sf::st_read(testdata_path("aoi.gpkg"), quiet = TRUE) + photos <- centroids[centroids$scale == "1:12000", ] + dem <- testdata_path("dem.tif") + + # mode = "all" routes through fly_select_all() + all_flat <- fly_select(photos, aoi, mode = "all") + all_terr <- fly_select(photos, aoi, mode = "all", dem = dem) + expect_gte(nrow(all_terr), nrow(all_flat)) + + # mode = "minimal" routes through fly_select_minimal() — a separate call site, + # so passing in one mode proves nothing about the other. + min_flat <- fly_select(photos, aoi, mode = "minimal", target_coverage = 0.95) + min_terr <- fly_select(photos, aoi, mode = "minimal", target_coverage = 0.95, + dem = dem) + expect_lte(nrow(min_terr), nrow(min_flat)) + + # Neither comparison above can fail on this data — the same frames are chosen + # with or without a DEM. Assert arrival directly, once per internal call site, + # since fly_select_all() and fly_select_minimal() are separate. + no_height <- photos + no_height$flying_height <- NULL + expect_error(fly_select(no_height, aoi, mode = "all", dem = dem), "flying_height") + expect_error(fly_select(no_height, aoi, mode = "minimal", + target_coverage = 0.95, dem = dem), "flying_height") +}) + +test_that("fly_georef passes dem through to the GCP footprints", { + skip_if_no_terra() + centroids <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE)[1:2, ] + # No images to warp, so this exercises the argument path and the footprint + # build rather than GDAL. A dem that never reached fly_footprint() would + # produce identical GCPs, which is exactly what cannot be seen from outside. + fetch_result <- dplyr::tibble( + airp_id = centroids$airp_id, + dest = file.path(tempdir(), paste0(centroids$airp_id, ".jpg")), + success = c(FALSE, FALSE) + ) + res <- suppressMessages( + fly_georef(fetch_result, centroids, + dest_dir = file.path(tempdir(), "georef-dem"), + dem = testdata_path("dem.tif")) + ) + expect_equal(nrow(res), 2) + + # The assertion above tolerates the argument; it does not prove it arrived, + # because with no images to warp the GCPs are never observable. Strip the + # column the DEM path requires: fly_footprint() then errors, and it can only + # do so if `dem` actually reached it. + no_height <- centroids + no_height$flying_height <- NULL + expect_error( + suppressMessages( + fly_georef(fetch_result, no_height, + dest_dir = file.path(tempdir(), "georef-dem2"), + dem = testdata_path("dem.tif")) + ), + "flying_height" + ) +}) diff --git a/vignettes/airphoto-selection.Rmd b/vignettes/airphoto-selection.Rmd index dbc1329..4043bee 100644 --- a/vignettes/airphoto-selection.Rmd +++ b/vignettes/airphoto-selection.Rmd @@ -58,10 +58,11 @@ is not in the centroid metadata. Supply it yourself if you know the camera: The photos used throughout this vignette are 1968 film, so every footprint below is sized from the 9-inch negative. -Note that footprints assume flat terrain beneath the aircraft. On slopes the -true ground coverage differs — downhill slopes produce a larger actual -footprint, uphill slopes a smaller one. All coverage and overlap numbers -downstream inherit this approximation. +By default footprints are sized from the reported scale, which assumes flat +ground at whatever elevation that scale was computed for. Supplying a DEM +removes that assumption — worth about 14% of footprint area here, and shown at +the end of this vignette. Every coverage and overlap number below inherits +whichever basis you choose. Figure \@ref(fig:fig-footprint) shows the estimated footprints for all 20 photos. Notice that some centroids fall outside the AOI while their footprints still overlap it — `fly_filter()` with `method = "footprint"` catches these @@ -274,8 +275,96 @@ georef <- fly_georef(fetched, centroids[1:3, ], georef[, c("airp_id", "dest", "success")] ``` -The georeferenced TIFFs inherit the flat-terrain and nadir-camera -assumptions from `fly_footprint()` — they are approximate, useful for -visual context rather than survey-grade positioning. Metadata from the -original centroid data (date, scale, focal length) links back via -`airp_id`. +The georeferenced TIFFs inherit whatever basis `fly_footprint()` was given, +plus its nadir-camera assumption — `fly_georef()` takes the same `dem` +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`. + +# Terrain-adjusted footprints + +Everything above sizes each frame from its reported scale. That scale is +referenced to some elevation, and where the ground sits below that elevation the +photo covers more than the scale implies. + +The catalogue carries what is needed to do better. `FLYING_HEIGHT` is metres +above **sea level**, not height above ground, so subtracting terrain elevation +gives the height ground coverage actually scales with: + +```{r terrain-compare, eval = requireNamespace("terra", quietly = TRUE)} +dem <- system.file("testdata/dem.tif", package = "fly") + +flat <- fly_footprint(centroids) +terrain <- fly_footprint(centroids, dem = dem) + +pct <- 100 * (as.numeric(st_area(st_transform(terrain, 3005))) / + as.numeric(st_area(st_transform(flat, 3005))) - 1) + +dplyr::tibble( + scale = centroids$scale, + elevation_m = round(centroids$flying_height - terrain$height_agl), + agl_m = round(terrain$height_agl), + area_change = round(pct, 1) +) |> + dplyr::group_by(scale) |> + dplyr::summarise( + n = dplyr::n(), + ground_m = paste0(min(elevation_m), "-", max(elevation_m)), + `area +%` = paste0(round(min(area_change), 1), " to ", round(max(area_change), 1)), + median = round(median(area_change), 1) + ) |> + knitr::kable(caption = "Footprint area change from terrain-adjusted sizing, by scale.") +``` + +Every footprint grows, and none shrink. That one-directional bias is the tell +that this is not a slope effect: the reported scale is referenced to an elevation +above the valley floor these photos actually cover, so it understates every +frame. Slope would push in both directions. + +```{r fig-terrain, eval = requireNamespace("terra", quietly = TRUE), fig.cap = "Flat-terrain footprints (grey) against terrain-adjusted footprints (red) for the 1:31680 frames. The correction is a uniform enlargement per frame, not a change of shape."} +idx <- centroids$scale == "1:31680" +plot(st_geometry(terrain[idx, ]), border = "firebrick") +plot(st_geometry(flat[idx, ]), border = "grey50", add = TRUE) +plot(st_geometry(aoi), col = NA, border = "grey20", add = TRUE) +``` + +`footprint_terrain` records what happened to each frame, so a fallback is +visible rather than silent: + +```{r terrain-basis, eval = requireNamespace("terra", quietly = TRUE)} +table(terrain$footprint_terrain) +``` + +The bundled `dem.tif` is a clip of **MRDEM-30**, NRCan's 30 m bare-earth DTM. +For your own area of interest, read it straight off the public COG — no download +or account needed, since `/vsicurl/` fetches only the bytes that intersect: + +```r +mrdem <- paste0( + "/vsicurl/https://canelevation-dem.s3.ca-central-1.amazonaws.com/", + "mrdem-30/mrdem-30-dtm.tif" +) +fly_coverage(centroids, aoi, by = "scale", dem = mrdem) +``` + +Buffer generously. Resolution matters less than extent here — 30 m resolves a +2.7 km footprint's mean elevation perfectly well; a DEM that stops short of the +frame edges does not. That is the ordinary case, not an exotic one: a DEM +cropped to your AOI simply stops, and the footprints run past it. `dem_coverage` +reports the fraction of each footprint the DEM actually described, and anything +below 95% warns: + +```{r terrain-coverage, eval = requireNamespace("terra", quietly = TRUE)} +round(range(terrain$dem_coverage), 3) +``` + +Buffer past the **corner** of the widest footprint rather than its half-side. +The far point of a square is `half_side * sqrt(2)` — 5.1 km at 1:31680, not +3.6 km — and the correction enlarges the rectangle again before the second +sampling pass reaches it. The bundled DEM allows 5.4 km. + +What a DEM does **not** fix is the nadir assumption. The catalogue carries no +tilt, roll or crab, so footprints stay axis-aligned rectangles. On this area +that per-corner refinement is worth roughly 2%, against the 14% the DEM +addresses. +