From e16960eeecc10623ee01c1a4a35fc1df5e2240b9 Mon Sep 17 00:00:00 2001 From: almac2022 Date: Thu, 5 Mar 2026 11:39:36 -0800 Subject: [PATCH] Add fly_overlap, fly_select modes, and document footprint assumptions (#8) - fly_select() gains mode param: "minimal" (fewest photos, default) and "all" (every photo touching AOI) - New fly_overlap() reports pairwise overlap % between photo footprints - Document 9x9 negative assumption and flat-terrain limitation in fly_footprint(), fly_overlap(), and vignette - Guard st_area() against empty geometry from full coverage - Simplify vignette with both modes, overlap analysis, and multi-scale workflow that sorts scales finest-first - 64 tests all passing, 0 errors 0 warnings on R CMD check Closes #8 Co-Authored-By: Claude Opus 4.6 --- NAMESPACE | 1 + R/fly_footprint.R | 18 ++++ R/fly_overlap.R | 81 +++++++++++++++ R/fly_select.R | 55 +++++++++-- man/fly-package.Rd | 2 + man/fly_footprint.Rd | 18 ++++ man/fly_overlap.Rd | 32 ++++++ man/fly_select.Rd | 27 +++-- tests/testthat/test-fly_overlap.R | 55 +++++++++++ tests/testthat/test-fly_select.R | 77 +++++++++++++-- vignettes/airphoto-selection.Rmd | 157 +++++++++++++++++++----------- 11 files changed, 443 insertions(+), 80 deletions(-) create mode 100644 R/fly_overlap.R create mode 100644 man/fly_overlap.Rd create mode 100644 tests/testthat/test-fly_overlap.R diff --git a/NAMESPACE b/NAMESPACE index 1cc995e..162e4e1 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -3,6 +3,7 @@ export(fly_coverage) export(fly_filter) export(fly_footprint) +export(fly_overlap) export(fly_query_habitat) export(fly_query_lakes) export(fly_select) diff --git a/R/fly_footprint.R b/R/fly_footprint.R index c4c3296..b0c7acf 100644 --- a/R/fly_footprint.R +++ b/R/fly_footprint.R @@ -12,6 +12,24 @@ #' per side. Rectangles are constructed in BC Albers (EPSG:3005) for accurate #' metric distances, then transformed back to the input CRS. #' +#' The scale denominator is parsed from the `scale` column string (e.g. +#' `"1:12000"` becomes `12000`). +#' +#' **9x9 assumption:** the default `negative_size = 9` (inches) reflects the +#' standard 228 mm format used by BC aerial survey cameras (e.g. Wild RC-10, +#' Zeiss RMK). The BC Air Photo Database records camera focal length per roll +#' (Type 02 field 3.2.2) but this is not available in the simplified centroid +#' data from the catalogue. If working with non-standard format photography, +#' override `negative_size` accordingly. +#' +#' **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. +#' #' @examples #' centroids <- sf::st_read(system.file("testdata/photo_centroids.gpkg", package = "fly")) #' footprints <- fly_footprint(centroids) diff --git a/R/fly_overlap.R b/R/fly_overlap.R new file mode 100644 index 0000000..879c9f2 --- /dev/null +++ b/R/fly_overlap.R @@ -0,0 +1,81 @@ +#' Compute pairwise overlap between photo footprints +#' +#' For each pair of photos whose footprints intersect, computes the overlap +#' 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. +#' +#' @param photos_sf An sf point object with a `scale` column. +#' @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. +#' +#' @examples +#' centroids <- sf::st_read(system.file("testdata/photo_centroids.gpkg", package = "fly")) +#' aoi <- sf::st_read(system.file("testdata/aoi.gpkg", package = "fly")) +#' photos_12k <- centroids[centroids$scale == "1:12000", ] +#' selected <- fly_select(photos_12k, aoi, mode = "all") +#' fly_overlap(selected) +#' +#' @export +fly_overlap <- function(photos_sf) { + sf::sf_use_s2(FALSE) + on.exit(sf::sf_use_s2(TRUE)) + + footprints <- fly_footprint(photos_sf) |> sf::st_transform(3005) + n <- nrow(footprints) + + if (n < 2) { + return(dplyr::tibble( + photo_a = integer(0), photo_b = integer(0), + overlap_km2 = numeric(0), pct_of_a = numeric(0), pct_of_b = numeric(0) + )) + } + + fp_areas <- as.numeric(sf::st_area(footprints)) + pairs <- sf::st_intersects(footprints) + + ids <- if ("airp_id" %in% names(footprints)) { + footprints$airp_id + } else { + seq_len(n) + } + + results <- list() + for (i in seq_len(n)) { + neighbors <- pairs[[i]] + neighbors <- neighbors[neighbors > i] + if (length(neighbors) == 0) next + + for (j in neighbors) { + overlap_geom <- tryCatch( + sf::st_intersection(sf::st_geometry(footprints[i, ]), + sf::st_geometry(footprints[j, ])) |> + sf::st_make_valid(), + error = function(e) NULL + ) + if (is.null(overlap_geom) || length(overlap_geom) == 0) next + + overlap_area <- as.numeric(sf::st_area(overlap_geom)) + if (overlap_area <= 0) next + + results <- c(results, list(dplyr::tibble( + photo_a = ids[i], + photo_b = ids[j], + overlap_km2 = round(overlap_area / 1e6, 3), + pct_of_a = round(overlap_area / fp_areas[i] * 100, 1), + pct_of_b = round(overlap_area / fp_areas[j] * 100, 1) + ))) + } + } + + if (length(results) == 0) { + return(dplyr::tibble( + photo_a = integer(0), photo_b = integer(0), + overlap_km2 = numeric(0), pct_of_a = numeric(0), pct_of_b = numeric(0) + )) + } + + dplyr::bind_rows(results) +} diff --git a/R/fly_select.R b/R/fly_select.R index 11d7fab..f32e47c 100644 --- a/R/fly_select.R +++ b/R/fly_select.R @@ -1,23 +1,60 @@ -#' Select minimum photo set to cover an AOI (greedy set cover) +#' Select photos covering an AOI #' -#' Iteratively picks the photo whose footprint covers the most uncovered -#' area until the target coverage is reached. +#' Two modes: `"minimal"` picks the fewest photos to reach target coverage +#' (greedy set-cover); `"all"` returns every photo whose footprint intersects +#' the AOI. #' #' @param photos_sf An sf point object with a `scale` column #' (pre-filtered to target year/scale). #' @param aoi_sf An sf polygon to cover. +#' @param mode Either `"minimal"` (fewest photos to reach target) or `"all"` +#' (every photo touching the AOI). #' @param target_coverage Stop when this fraction is reached (default 0.95). -#' @return An sf object (subset of `photos_sf`) with added columns -#' `selection_order` and `cumulative_coverage_pct`. +#' Only used when `mode = "minimal"`. +#' @return An sf object (subset of `photos_sf`). For `mode = "minimal"`, +#' includes `selection_order` and `cumulative_coverage_pct` columns. #' #' @examples #' centroids <- sf::st_read(system.file("testdata/photo_centroids.gpkg", package = "fly")) #' aoi <- sf::st_read(system.file("testdata/aoi.gpkg", package = "fly")) -#' selected <- fly_select(centroids, aoi, target_coverage = 0.80) -#' selected[, c("airp_id", "scale", "selection_order", "cumulative_coverage_pct")] +#' +#' # Fewest photos to reach 80% coverage +#' fly_select(centroids, aoi, mode = "minimal", target_coverage = 0.80) +#' +#' # All photos touching the AOI +#' fly_select(centroids, aoi, mode = "all") #' #' @export -fly_select <- function(photos_sf, aoi_sf, target_coverage = 0.95) { +fly_select <- function(photos_sf, aoi_sf, mode = "minimal", + target_coverage = 0.95) { + mode <- match.arg(mode, c("minimal", "all")) + + if (mode == "all") { + return(fly_select_all(photos_sf, aoi_sf)) + } + + fly_select_minimal(photos_sf, aoi_sf, target_coverage) +} + +#' @noRd +fly_select_all <- function(photos_sf, aoi_sf) { + sf::sf_use_s2(FALSE) + on.exit(sf::sf_use_s2(TRUE)) + + footprints <- fly_footprint(photos_sf) + aoi_union <- sf::st_transform(aoi_sf, sf::st_crs(footprints)) |> + sf::st_union() |> + sf::st_make_valid() + + touches <- sf::st_intersects(footprints, aoi_union, sparse = FALSE)[, 1] + result <- photos_sf[touches, ] + message("Selected ", nrow(result), " of ", nrow(photos_sf), + " photos intersecting the AOI") + result +} + +#' @noRd +fly_select_minimal <- function(photos_sf, aoi_sf, target_coverage) { sf::sf_use_s2(FALSE) on.exit(sf::sf_use_s2(TRUE)) @@ -66,7 +103,7 @@ fly_select <- function(photos_sf, aoi_sf, target_coverage = 0.95) { error = function(e) aoi_albers ) - pct <- as.numeric(sf::st_area(covered_in_aoi)) / aoi_area + pct <- sum(as.numeric(sf::st_area(covered_in_aoi))) / aoi_area coverage_pcts <- c(coverage_pcts, pct) if (length(selected_idx) %% 10 == 0 || pct >= target_coverage) { diff --git a/man/fly-package.Rd b/man/fly-package.Rd index e10b8e1..b0882e1 100644 --- a/man/fly-package.Rd +++ b/man/fly-package.Rd @@ -6,6 +6,8 @@ \alias{fly-package} \title{fly: Airphoto Footprint Estimation and Coverage Selection} \description{ +\if{html}{\figure{logo.png}{options: style='float: right' alt='logo' width='120'}} + Estimate ground footprints from airphoto centroids and scale, compute coverage of areas of interest, and select minimum photo sets using greedy set-cover. Includes optional helpers for querying fish habitat streams from bcfishpass. } \seealso{ diff --git a/man/fly_footprint.Rd b/man/fly_footprint.Rd index 2522c6d..a9c9265 100644 --- a/man/fly_footprint.Rd +++ b/man/fly_footprint.Rd @@ -22,6 +22,24 @@ of each airphoto, based on film negative dimensions and the reported scale. Ground coverage is computed as \code{negative_size * scale_number * 0.0254} metres per side. Rectangles are constructed in BC Albers (EPSG:3005) for accurate metric distances, then transformed back to the input CRS. + +The scale denominator is parsed from the \code{scale} column string (e.g. +\code{"1:12000"} becomes \code{12000}). + +\strong{9x9 assumption:} the default \code{negative_size = 9} (inches) reflects the +standard 228 mm format used by BC aerial survey cameras (e.g. Wild RC-10, +Zeiss RMK). The BC Air Photo Database records camera focal length per roll +(Type 02 field 3.2.2) but this is not available in the simplified centroid +data from the catalogue. If working with non-standard format photography, +override \code{negative_size} accordingly. + +\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. } \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 new file mode 100644 index 0000000..e9f8177 --- /dev/null +++ b/man/fly_overlap.Rd @@ -0,0 +1,32 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/fly_overlap.R +\name{fly_overlap} +\alias{fly_overlap} +\title{Compute pairwise overlap between photo footprints} +\usage{ +fly_overlap(photos_sf) +} +\arguments{ +\item{photos_sf}{An sf point object with a \code{scale} column.} +} +\value{ +A tibble with columns \code{photo_a}, \code{photo_b}, \code{overlap_km2}, +\code{pct_of_a}, and \code{pct_of_b}. Only pairs with non-zero overlap are returned. +} +\description{ +For each pair of photos whose footprints intersect, computes the overlap +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. +} +\examples{ +centroids <- sf::st_read(system.file("testdata/photo_centroids.gpkg", package = "fly")) +aoi <- sf::st_read(system.file("testdata/aoi.gpkg", package = "fly")) +photos_12k <- centroids[centroids$scale == "1:12000", ] +selected <- fly_select(photos_12k, aoi, mode = "all") +fly_overlap(selected) + +} diff --git a/man/fly_select.Rd b/man/fly_select.Rd index e156733..3aa718a 100644 --- a/man/fly_select.Rd +++ b/man/fly_select.Rd @@ -2,9 +2,9 @@ % Please edit documentation in R/fly_select.R \name{fly_select} \alias{fly_select} -\title{Select minimum photo set to cover an AOI (greedy set cover)} +\title{Select photos covering an AOI} \usage{ -fly_select(photos_sf, aoi_sf, target_coverage = 0.95) +fly_select(photos_sf, aoi_sf, mode = "minimal", target_coverage = 0.95) } \arguments{ \item{photos_sf}{An sf point object with a \code{scale} column @@ -12,20 +12,29 @@ fly_select(photos_sf, aoi_sf, target_coverage = 0.95) \item{aoi_sf}{An sf polygon to cover.} -\item{target_coverage}{Stop when this fraction is reached (default 0.95).} +\item{mode}{Either \code{"minimal"} (fewest photos to reach target) or \code{"all"} +(every photo touching the AOI).} + +\item{target_coverage}{Stop when this fraction is reached (default 0.95). +Only used when \code{mode = "minimal"}.} } \value{ -An sf object (subset of \code{photos_sf}) with added columns -\code{selection_order} and \code{cumulative_coverage_pct}. +An sf object (subset of \code{photos_sf}). For \code{mode = "minimal"}, +includes \code{selection_order} and \code{cumulative_coverage_pct} columns. } \description{ -Iteratively picks the photo whose footprint covers the most uncovered -area until the target coverage is reached. +Two modes: \code{"minimal"} picks the fewest photos to reach target coverage +(greedy set-cover); \code{"all"} returns every photo whose footprint intersects +the AOI. } \examples{ centroids <- sf::st_read(system.file("testdata/photo_centroids.gpkg", package = "fly")) aoi <- sf::st_read(system.file("testdata/aoi.gpkg", package = "fly")) -selected <- fly_select(centroids, aoi, target_coverage = 0.80) -selected[, c("airp_id", "scale", "selection_order", "cumulative_coverage_pct")] + +# Fewest photos to reach 80\% coverage +fly_select(centroids, aoi, mode = "minimal", target_coverage = 0.80) + +# All photos touching the AOI +fly_select(centroids, aoi, mode = "all") } diff --git a/tests/testthat/test-fly_overlap.R b/tests/testthat/test-fly_overlap.R new file mode 100644 index 0000000..7bda253 --- /dev/null +++ b/tests/testthat/test-fly_overlap.R @@ -0,0 +1,55 @@ +test_that("fly_overlap returns tibble with expected columns", { + centroids <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE) + result <- fly_overlap(centroids) + expect_s3_class(result, "tbl_df") + expect_true(all(c("photo_a", "photo_b", "overlap_km2", + "pct_of_a", "pct_of_b") %in% names(result))) +}) + +test_that("fly_overlap finds overlapping same-scale photos", { + centroids <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE) + photos_12k <- centroids[centroids$scale == "1:12000", ] + result <- fly_overlap(photos_12k) + # adjacent flight-line photos should have some overlap + + expect_gt(nrow(result), 0) + expect_true(all(result$overlap_km2 > 0)) + expect_true(all(result$pct_of_a >= 0 & result$pct_of_a <= 100)) + expect_true(all(result$pct_of_b >= 0 & result$pct_of_b <= 100)) +}) + +test_that("fly_overlap uses airp_id when available", { + centroids <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE) + photos_12k <- centroids[centroids$scale == "1:12000", ] + result <- fly_overlap(photos_12k) + if (nrow(result) > 0) { + expect_true(all(result$photo_a %in% photos_12k$airp_id)) + expect_true(all(result$photo_b %in% photos_12k$airp_id)) + } +}) + +test_that("fly_overlap returns empty tibble for single photo", { + centroids <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE) + result <- fly_overlap(centroids[1, ]) + expect_equal(nrow(result), 0) + expect_s3_class(result, "tbl_df") +}) + +test_that("fly_overlap pairs are unique (no duplicates)", { + centroids <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE) + result <- fly_overlap(centroids) + if (nrow(result) > 0) { + pair_keys <- paste(result$photo_a, result$photo_b, sep = "-") + expect_equal(length(pair_keys), length(unique(pair_keys))) + } +}) + +test_that("fly_overlap larger scale has larger overlaps", { + centroids <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE) + overlap_12k <- fly_overlap(centroids[centroids$scale == "1:12000", ]) + overlap_31k <- fly_overlap(centroids[centroids$scale == "1:31680", ]) + # 1:31680 footprints are ~7x larger so overlap area should be larger + if (nrow(overlap_12k) > 0 && nrow(overlap_31k) > 0) { + expect_gt(max(overlap_31k$overlap_km2), max(overlap_12k$overlap_km2)) + } +}) diff --git a/tests/testthat/test-fly_select.R b/tests/testthat/test-fly_select.R index 18101bd..f2a14e3 100644 --- a/tests/testthat/test-fly_select.R +++ b/tests/testthat/test-fly_select.R @@ -1,25 +1,90 @@ -test_that("fly_select returns sf with selection columns", { +test_that("fly_select minimal mode returns sf with selection columns", { centroids <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE) aoi <- sf::st_read(testdata_path("aoi.gpkg"), quiet = TRUE) - result <- fly_select(centroids, aoi, target_coverage = 0.50) + result <- fly_select(centroids, aoi, mode = "minimal", target_coverage = 0.50) expect_s3_class(result, "sf") expect_true("selection_order" %in% names(result)) expect_true("cumulative_coverage_pct" %in% names(result)) }) -test_that("fly_select selects fewer photos than input", { +test_that("fly_select minimal mode selects fewer photos than input", { centroids <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE) aoi <- sf::st_read(testdata_path("aoi.gpkg"), quiet = TRUE) - result <- fly_select(centroids, aoi, target_coverage = 0.50) + result <- fly_select(centroids, aoi, mode = "minimal", target_coverage = 0.50) expect_lt(nrow(result), nrow(centroids)) }) -test_that("fly_select coverage increases monotonically", { +test_that("fly_select minimal coverage increases monotonically", { centroids <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE) aoi <- sf::st_read(testdata_path("aoi.gpkg"), quiet = TRUE) - result <- fly_select(centroids, aoi, target_coverage = 0.50) + result <- fly_select(centroids, aoi, mode = "minimal", target_coverage = 0.50) if (nrow(result) > 1) { diffs <- diff(result$cumulative_coverage_pct) expect_true(all(diffs >= 0)) } }) + +test_that("fly_select defaults to minimal mode", { + centroids <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE) + aoi <- sf::st_read(testdata_path("aoi.gpkg"), quiet = TRUE) + result_default <- fly_select(centroids, aoi, target_coverage = 0.50) + result_explicit <- fly_select(centroids, aoi, mode = "minimal", target_coverage = 0.50) + expect_equal(nrow(result_default), nrow(result_explicit)) + expect_equal(result_default$selection_order, result_explicit$selection_order) +}) + +test_that("fly_select all mode returns every photo touching AOI", { + centroids <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE) + aoi <- sf::st_read(testdata_path("aoi.gpkg"), quiet = TRUE) + result <- fly_select(centroids, aoi, mode = "all") + expect_s3_class(result, "sf") + # all mode should return at least as many as minimal + result_min <- fly_select(centroids, aoi, mode = "minimal", target_coverage = 0.99) + expect_gte(nrow(result), nrow(result_min)) +}) + +test_that("fly_select all mode does not add selection columns", { + centroids <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE) + aoi <- sf::st_read(testdata_path("aoi.gpkg"), quiet = TRUE) + result <- fly_select(centroids, aoi, mode = "all") + expect_false("selection_order" %in% names(result)) + expect_false("cumulative_coverage_pct" %in% names(result)) +}) + +test_that("fly_select all mode only returns photos intersecting AOI", { + centroids <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE) + aoi <- sf::st_read(testdata_path("aoi.gpkg"), quiet = TRUE) + result <- fly_select(centroids, aoi, mode = "all") + # verify every selected photo footprint actually intersects the AOI + fp <- fly_footprint(result) + aoi_t <- sf::st_transform(aoi, sf::st_crs(fp)) + touches <- sf::st_intersects(fp, aoi_t, sparse = FALSE)[, 1] + expect_true(all(touches)) +}) + +test_that("fly_select all on single scale returns subset", { + centroids <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE) + aoi <- sf::st_read(testdata_path("aoi.gpkg"), quiet = TRUE) + photos_12k <- centroids[centroids$scale == "1:12000", ] + result <- fly_select(photos_12k, aoi, mode = "all") + expect_s3_class(result, "sf") + expect_lte(nrow(result), nrow(photos_12k)) + expect_true(all(result$scale == "1:12000")) +}) + +test_that("fly_select rejects invalid mode", { + centroids <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE) + aoi <- sf::st_read(testdata_path("aoi.gpkg"), quiet = TRUE) + expect_error(fly_select(centroids, aoi, mode = "bogus")) +}) + +test_that("fly_select minimal handles full coverage without error", { + centroids <- sf::st_read(testdata_path("photo_centroids.gpkg"), quiet = TRUE) + aoi <- sf::st_read(testdata_path("aoi.gpkg"), quiet = TRUE) + # target 1.0 forces the algorithm to exhaust all photos or hit 100% + result <- fly_select(centroids, aoi, mode = "minimal", target_coverage = 1.0) + expect_s3_class(result, "sf") + expect_gt(nrow(result), 0) + # coverage should be scalar, not length 0 + expect_length(result$cumulative_coverage_pct, nrow(result)) +}) diff --git a/vignettes/airphoto-selection.Rmd b/vignettes/airphoto-selection.Rmd index f167a68..baab93e 100644 --- a/vignettes/airphoto-selection.Rmd +++ b/vignettes/airphoto-selection.Rmd @@ -31,8 +31,18 @@ aoi <- st_read(system.file("testdata/aoi.gpkg", package = "fly"), quiet = TRUE) ## Footprint estimation `fly_footprint()` converts point centroids into rectangular polygons -representing estimated ground coverage. The standard 9" x 9" negative -produces a footprint width of `9 * scale * 0.0254` metres. +representing estimated ground coverage. The standard 9" x 9" (228 mm) negative +used by BC aerial survey cameras produces a footprint width of +`9 * scale * 0.0254` metres. The 9x9 format is the default (`negative_size = 9`) +and is correct for standard BC catalogue photography. The full BC Air Photo +Database records camera focal length per roll, but this field is not available +in the simplified centroid data — override `negative_size` if working with +non-standard format photography. + +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. ```{r footprint} footprints <- fly_footprint(centroids) @@ -71,91 +81,126 @@ footprints, grouped by any column. fly_coverage(centroids, aoi, by = "scale") ``` -## Greedy photo selection +## Photo selection -`fly_select()` picks the minimum set of photos needed to reach a target -coverage, using a greedy set-cover algorithm that iteratively selects the -photo covering the most uncovered area. +`fly_select()` has two modes: -```{r select} -selected <- fly_select(centroids, aoi, target_coverage = 0.80) +- `mode = "minimal"` — fewest photos to reach target coverage +- `mode = "all"` — every photo whose footprint touches the AOI + +### Minimal selection + +```{r select-minimal} +selected <- fly_select(centroids, aoi, mode = "minimal", target_coverage = 0.80) selected[, c("airp_id", "scale", "selection_order", "cumulative_coverage_pct")] ``` -```{r plot-selected} +```{r plot-minimal} sel_fp <- fly_footprint(selected) plot(st_geometry(aoi), col = "lightyellow", border = "grey40", - main = paste(nrow(selected), "photos selected (greedy)")) + main = paste(nrow(selected), "photos (minimal selection)")) plot(st_geometry(sel_fp), border = "steelblue", col = adjustcolor("steelblue", 0.15), add = TRUE) plot(st_geometry(selected), pch = 20, col = "red", add = TRUE) ``` -## Priority selection: best resolution first +### All photos touching AOI + +```{r select-all} +all_in_aoi <- fly_select(centroids, aoi, mode = "all") +cat(nrow(all_in_aoi), "photos intersect the AOI\n") +``` + +## Overlap analysis -In practice you often want **all** photos at the best available resolution, -then backfill uncovered area with coarser scales. This combines `fly_select()` -with a priority loop over scales. +`fly_overlap()` reports pairwise overlap between photo footprints. +Run it on same-scale subsets to understand coverage quality. + +```{r overlap} +photos_12k <- centroids[centroids$scale == "1:12000", ] +overlap_12k <- fly_overlap(photos_12k) +overlap_12k +``` + +```{r overlap-summary} +if (nrow(overlap_12k) > 0) { + cat("1:12000 overlap range:", + round(min(overlap_12k$pct_of_a), 1), "% -", + round(max(overlap_12k$pct_of_a), 1), "%\n") +} -```{r priority} +photos_31k <- centroids[centroids$scale == "1:31680", ] +overlap_31k <- fly_overlap(photos_31k) +if (nrow(overlap_31k) > 0) { + cat("1:31680 overlap range:", + round(min(overlap_31k$pct_of_a), 1), "% -", + round(max(overlap_31k$pct_of_a), 1), "%\n") +} +``` + +## Multi-scale workflow: best resolution first + +In practice you want the finest-scale photos first, then fill gaps with +coarser scales. Sort scales finest-first by parsing the numeric denominator: + +```{r multi-scale} sf_use_s2(FALSE) -scales_priority <- c("1:12000", "1:31680") -target_coverage <- 0.80 +# Sort scales finest-first +scales <- sort(unique(as.numeric(gsub("1:", "", centroids$scale)))) +cat("Scales (finest first):", paste0("1:", scales), "\n") + +target_coverage <- 0.80 aoi_albers <- st_transform(aoi, 3005) |> st_union() |> st_make_valid() aoi_area <- as.numeric(st_area(aoi_albers)) selected_all <- NULL remaining_aoi <- aoi_albers -for (i in seq_along(scales_priority)) { - sc <- scales_priority[i] +for (sc_num in scales) { + sc <- paste0("1:", sc_num) photos_sc <- centroids[centroids$scale == sc, ] if (nrow(photos_sc) == 0) next - if (i == 1) { - # Best resolution: take ALL photos - cat(sc, ": taking all", nrow(photos_sc), "photos\n") - sel <- photos_sc - sel$selection_order <- seq_len(nrow(sel)) - sel$cumulative_coverage_pct <- NA_real_ - } else { - # Coarser scales: greedy select on remaining uncovered area - remaining_sf <- st_sf(geometry = st_geometry(remaining_aoi)) |> - st_transform(4326) |> st_make_valid() - sel <- fly_select(photos_sc, remaining_sf, target_coverage = target_coverage) - } - - if (nrow(sel) > 0) { - fp <- fly_footprint(sel) |> st_transform(3005) - fp_union <- st_union(fp) |> st_make_valid() - remaining_aoi <- tryCatch( - st_difference(remaining_aoi, fp_union) |> st_make_valid(), - error = function(e) remaining_aoi - ) - covered_pct <- 1 - as.numeric(st_area(remaining_aoi)) / aoi_area - cat(" -> ", nrow(sel), " photos at ", sc, - " (cumulative: ", round(covered_pct * 100, 1), "%)\n") - sel$priority_scale <- sc - selected_all <- rbind(selected_all, sel) - } + # Take all photos at this scale that touch the (remaining) AOI + remaining_sf <- st_sf(geometry = st_geometry(remaining_aoi)) |> + st_transform(4326) |> st_make_valid() + sel <- fly_select(photos_sc, remaining_sf, mode = "all") + if (nrow(sel) == 0) next + + # Update remaining uncovered area + fp <- fly_footprint(sel) |> st_transform(3005) + fp_union <- st_union(fp) |> st_make_valid() + remaining_aoi <- tryCatch( + st_difference(remaining_aoi, fp_union) |> st_make_valid(), + error = function(e) remaining_aoi + ) + remaining_area <- as.numeric(st_area(remaining_aoi)) + covered_pct <- 1 - sum(remaining_area) / aoi_area + + cat(sc, ":", nrow(sel), "photos (cumulative coverage:", + round(covered_pct * 100, 1), "%)\n") + + sel$priority_scale <- sc + selected_all <- rbind(selected_all, sel) + + if (covered_pct >= target_coverage) break } cat("\nTotal:", nrow(selected_all), "photos\n") as.data.frame(table(selected_all$priority_scale)) ``` -```{r plot-priority} +```{r plot-multi-scale} sel_fp <- fly_footprint(selected_all) plot(st_geometry(aoi), col = "lightyellow", border = "grey40", - main = paste(nrow(selected_all), "photos (priority selection)")) -# Color by scale -cols <- ifelse(selected_all$priority_scale == "1:12000", "steelblue", "darkorange") -border_cols <- cols -fill_cols <- adjustcolor(cols, 0.15) + main = paste(nrow(selected_all), "photos (best resolution first)")) +scale_labels <- sort(unique(selected_all$priority_scale)) +palette <- c("steelblue", "darkorange", "forestgreen", "firebrick") +cols <- palette[match(selected_all$priority_scale, scale_labels)] for (j in seq_len(nrow(sel_fp))) { - plot(st_geometry(sel_fp[j, ]), border = border_cols[j], - col = fill_cols[j], add = TRUE) + plot(st_geometry(sel_fp[j, ]), border = cols[j], + col = adjustcolor(cols[j], 0.15), add = TRUE) } -legend("topright", legend = scales_priority, - fill = adjustcolor(c("steelblue", "darkorange"), 0.3), - border = c("steelblue", "darkorange"), bty = "n") +legend("topright", legend = scale_labels, + fill = adjustcolor(palette[seq_along(scale_labels)], 0.3), + border = palette[seq_along(scale_labels)], bty = "n") ```