Estimate photo footprint polygons from centroids and scale

Creates rectangular polygons representing the estimated ground coverage of each airphoto, based on film negative dimensions and the reported scale.

Usage

fly_footprint(centroids_sf, negative_size = 9, format_size = NULL, dem = NULL)

Arguments

centroids_sf: An sf point object with a scale column (e.g. "1:31680"). A media column (e.g. "Film - BW", "Digital - Colour") selects the recording format per frame when present. Geometry must be POINT. sf::st_coordinates() returns one row per vertex for anything else, so a POLYGON input silently multiplies the rows by the vertex count; this is refused rather than coerced, because re-estimating a footprint from the centroid of an estimated footprint is not meaningful.
negative_size: Negative dimension in inches (default 9 for standard 9" x 9"). Applies to film frames — media "Film - BW", "Film - Colour", "Film - BW IR" and "Film - Colour IR" — and to every frame when there is no media column. It never sizes a digital frame — see format_size.
format_size: Named numeric vector of recording-format widths in inches, keyed by media value, merged over the shipped film defaults. Frames it names are sized from the reported scale, as film is, and it takes precedence over the shipped camera table — it is the escape hatch for a camera fly does not know. See Details.
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.

Value

An sf polygon object in the same CRS as input, with footprint rectangles, a footprint_basis column recording how each was sized, a footprint_terrain column recording which terrain treatment was applied, footprint_bearing giving the flight azimuth each rectangle was rotated onto (NA where it was drawn axis-aligned because no bearing could be computed), height_agl giving the metres above ground each was sized from, height_source recording where that height came from or why there is none (see Terrain), dem_coverage giving the fraction of each footprint the DEM actually covered (0 where it covered none, NA only where there is no footprint), dem_elev_sd giving the standard deviation in metres of the DEM values under that footprint — NA wherever dem_coverage is, and also wherever the DEM described fewer than two cells, since a spread needs two. The common case of the second is a frame the DEM does not reach at all, which has dem_coverage of 0 rather than NA and so is the one place the two columns differ. And dem_shortfall_m giving the metres the DEM would have to extend to contain that footprint, where 0 is an answer: the DEM already spans the frame, so any missing cells are nodata inside its extent and no re-crop will recover them. Frames whose format could not be resolved get an empty geometry. Every class the input carries is carried through, so a tibble-backed sf — which is what bcdata::collect() returns — comes back tibble-backed. The order is not preserved: sf::st_transform() moves sf to the front, so a bcdc_sf input returns sf, bcdc_sf, ..., as it always has.

Details

Ground coverage is computed as 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 scale column string (e.g. "1:12000" becomes 12000).

Film and digital are not the same measurement. The 9-inch default reflects the standard 228 mm negative used by BC aerial survey cameras (e.g. Wild RC-10, Zeiss RMK). A digital frame has no negative, and the catalogue mixes the two in one layer — roughly a fifth of frames in a sampled area are Digital - Colour. Ground width scales with the recording format, so applying a negative dimension to a sensor produces a rectangle that is wrong by an unknown factor while still drawing, still overlapping neighbours, and still yielding a coverage percentage.

Each row is therefore sized from its media value, and footprint_basis records the outcome:

the media value: format resolved from the format table: one of the four film values, sized from negative_size (infrared film since fly#89, measured to be the same 9-inch negative); a value named in format_size; or a digital frame whose camera the shipped table resolved by calibration or PAT-B identity
"inferred_format": digital frame with no calibration and no camera named in its PAT-B file, sized from a format inferred from its focal_length
"assumed_default": no media column; negative_size applied
"unknown_format": media present but unknown; empty geometry

Digital frames

A digital frame has no negative, and the catalogue mixes film and digital in one layer — 223,667 of 1,670,471 frames province-wide are Digital - Colour. Sensor dimensions are not in the centroid metadata, but they are recoverable from the calibration report each frame links to through camera_calibration_url, and fly ships them (inst/extdata/camera_formats.csv, built by data-raw/make_camera_formats.R).

Digital frames are sized as pixel count x ground_sample_distance, which needs neither scale nor a DEM. scale is never used for a digital frame fly sized itself. That field is not the true image scale for digital: measured against terrain on 40 UltraCam Eagle frames it gives 34% of true width, because it is a derived nominal figure — the pixel pitch it implies is about 12.5 um for every camera regardless of model, against real pitches of 3.9 to 12 um. ground_sample_distance is in centimetres.

Where a frame carries no camera_calibration_url — about a fifth of digital frames — there are two remaining routes, in this order.

The catalogue publishes the camera that flew the frame in the per-frame georeferencing file it links through patb_georef_url. fly_camera_patb() fetches and parses those and attaches camera_serial, camera_name and patb_gsd; where those columns are present, fly_footprint() sizes the frame from the named sensor on the ordinary px x GSD route, needing no DEM, and width_source records patb_serial= or patb_camera=. patb_gsd is used only where the catalogue's own ground_sample_distance is absent or zero, which it is for whole years of digital imagery. This function never reaches the network itself.

Failing that, the format is inferred from focal_length and footprint_basis records "inferred_format". Sensor width spreads only 1-3% at a given focal length, but pixel count spreads 32-83%, so an inferred frame can only be sized through a DEM (width x height above ground / focal length) and never from its GSD.

width_source names the calibration file, camera or fallback rule per row, so every footprint traces back to a source — and so does every refusal. Calibrations that could not be corroborated are listed in inst/extdata/camera_formats_excluded.csv with the reason, and frames naming one are refused rather than inferred (withheld:). A camera identity the shipped table does not recognise is likewise refused rather than guessed from the model string beside it (unknown_serial:, unknown_camera:, ambiguous_serial:); unlike a withheld calibration, that refusal does not block the focal-length inference, because it says nothing about the sensor's size.

Every footprint is rotated onto the flight line using fly_bearing(), and footprint_bearing records the azimuth each was rotated onto. Where no bearing can be computed — a single-frame roll, or absent film_roll / frame_number columns — the rectangle stays axis-aligned, footprint_bearing is NA and width_source says so.

Digital sensors are not square, running from 1.10:1 (Leica DMC II) to 1.80:1 (Intergraph DMC), so rotating them is visibly load-bearing. Film is square, which is why rotating it was deferred until fly#26 — but a square rotated 45 degrees overlaps its axis-aligned self by only about 83%, so leaving film alone was not neutral. It was reporting coverage of ground the frame does not photograph. Film footprints on off-cardinal headings therefore move as of v0.9.0.

Supply format_size to size a frame fly cannot, or to override it:

fly_footprint(photos, format_size = c("Digital - Colour" = 3.54))

Filter on footprint_basis to keep only frames sized from a known format.

Focal length and flying height are available. FOCAL_LENGTH, FLYING_HEIGHT and SCALE are fully populated in WHSE_IMAGERY_AND_BASE_MAPS.AIMG_PHOTO_CENTROIDS_SP and carried through to the bundled test data. Note the field is SCALE, not PHOTO_SCALE — the latter returns all NULL, which reads as missing data rather than a wrong field name.

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:

"nominal_scale": sized from the reported scale (no dem, or a fallback — see below)
"gsd_scaled": digital frame sized from its pixel count and ground sample distance; used neither the reported scale nor a DEM, so height_agl and dem_coverage are NA
"dem_agl": sized from height above ground
"no_dem_coverage": dem supplied but does not cover the frame
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.

Or that the frame is not where the catalogue puts it. On five roll-heights (bc77026 2042 m, bc77070 1158 m, bc77072 1829 m and 1981 m, bc77087 1158 m) the logbook pages write the catalogued height above sea level for at least 90% of the frames they read, and the photos show the catalogue's step between centroids is longer than the air base (fly#97), so a second warning names any of their frames that land under the terrain, keyed on film_roll, flying_height, focal_length and scale. The label is per roll-height: not every frame on them is read by the logbook, bc77087's page rests on a digit settled by a read that was not blind (fly#101), and bc77070's step margin is at the instrument's resolution. The footprint is unchanged, still drawn at the centroid, and nothing is added to the output. The list is the roll-heights labelled so far, read from inst/extdata/flying_height_image_overlap_keys.csv, not a census: an unlisted roll-height is not thereby clean. Most were never measured, and four that fly#97 measured carry another label. Which frames are named depends on your DEM, and only frames under the terrain are. Of the 311 frames on these roll-heights that fly#93's census holds (those in band above sea level but below it over the ground, and those under the terrain), 112 are under MRDEM's terrain; the other 199 are below the band and refused by the height check with its own warning. Frames in band are sized from the DEM as usual. Without dem, nothing is named.

flying_height is checked before it is believed. Sizing from it means inheriting whatever is wrong with it, and the catalogue's is about 10.76 times (3.28084 squared) too large on 1,589 film frames from 13 rolls flown between 1974 and 2005 — a feet-to-metres conversion applied the wrong way round, which draws a 1:35000 frame 110 km across. On six roll-heights flown 1978-2000 the flight logbooks show the factor is exactly 10 instead (fly#71), and those frames are carried in the roll table below, as are two more that the same roll's adjacent frames settle (fly#74). A film frame states its height above ground twice, as flying_height minus terrain and as scale times focal length, so the two are compared, and height_source records the outcome:

"reported": the two agree to within a factor of 1.6, as 99% of frames do; sized from flying_height as supplied
"corrected_unit_slip": they disagree, and dividing flying_height by 10.76 brings them back into agreement; sized from the corrected height. height_agl is the height used. Your flying_height column is not overwritten, so on these rows flying_height - height_agl is not the ground elevation
"corrected_roll_table": they disagree, and the frame sits on one of the roll-heights in inst/extdata/flying_height_rolls.csv, which the province's flight logbooks, or failing them an adjacent frame on the same roll in an exact named relation, together with the spacing between adjacent frames, settled (fly#60, fly#71, fly#72, fly#74, fly#91): a height with one or two digits dropped (a logbook "20.0", thousands of feet, catalogued as 2,000 ft), a height recorded ten times too large, a leading digit added, or a correct height beside a wrong scale — on either side of the band, including frames around twice their nominal scale whose scale is recorded at half its denominator. Checked before the 10.76 slip above. Sized from the measured height; matched on film_roll, flying_height, focal_length and scale together, so it needs a film_roll column. As with the slip, flying_height is not overwritten
"implausible": they disagree some other way, or flying_height is above 16,000 m or below the terrain. Nothing can say whether the height, the scale or the focal length is the wrong one — a 305 mm lens catalogued as 153 and a scale recorded at half its denominator look the same from here, and the roll table separates them only where it was measured — so the height is not used: a film frame falls back to nominal scale, and a digital frame with no other route has no footprint
NA: no height was judged: no dem, a frame the DEM route does not size (one sized from its ground sample distance, or with no footprint at all for another reason), one the DEM does not cover, or flying_height or focal_length missing — except that a height above 16,000 m is "implausible" even where focal_length is missing

Measured over the whole catalogue, the two corrections for a height that is too large apply to those 1,589 frames and to nothing else: 308 by the roll table, the rest by 10.76. (The table's factor-1 rows correct a scale, not a height, so they are not among them.) A digital frame's scale is a nominal figure, not its image scale, so it is not compared: only the 16,000 m ceiling, or terrain at or above the aircraft, can refuse one. To list the frames worth a second look: dplyr::filter(fp, height_source != "reported").

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:

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 Reading it over /vsicurl/ needs no download: the DEM is sampled through one window per frame, so two frames take seconds rather than the minutes a whole-vector read cost before 0.13.0.
LidarBC — sub-10 m where coverage exists; query the stac-elevation-bc STAC catalogue and pass an item's COG URL.
BC TRIM — 25 m provincial DEM via the bcdata CLI (bcdata get-dem).

Bare earth or surface model. Over forest the camera images the canopy, so a bare-earth DEM draws a frame slightly too wide. Measured against MRDEM's own DSM (mrdem-30-dsm.tif, on the DTM's grid), that is a weighted median 0.17% of width and 0.46% at the 95th percentile, with 2.1% of fine-scale frames over 1% — under what the package acts on, so the DTM stays the recommendation. Either surface works; note that a canopy also feeds the flying_height check, so a frame near the edge of the band can change height_source between them. See inst/notes/terrain-correction.md (fly#80).

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.

Two remedies, and the column that tells them apart. Partial coverage has two causes that dem_coverage reports identically. A DEM whose extent stops short of the frame can be re-cropped and the frame recovered exactly, and dem_shortfall_m says by how much — max(fp$dem_shortfall_m) is the buffer that would contain every footprint in the batch. A DEM with nodata inside its extent cannot be re-cropped at all, and reports 0. Against a remote MRDEM read the shortfall is always 0, because the window is fetched per frame; the first case is what a DEM cropped to an area of interest produces.

What truncation costs, and why two columns report it. Measured over 11,520 truncations of 120 real frames against each frame's own full-coverage answer: the error is exactly the elevation bias divided by the height above ground, so it is a median 0.18% of footprint width at 80-90% coverage, 0.34% at 60-80%, and 1.5% below 20% — but it reaches 24% on broken ground. The median is not what one frame is exposed to: at a fixed dem_coverage the error spans 13 to 52 times between frames, which is why dem_elev_sd ships beside it. Frames above the median spread for their coverage sit 2.3 to 4.1 times further from the full-coverage answer than those below it. Read the pair together — coverage says how much is missing, dem_elev_sd says how much that could be worth — and see inst/notes/terrain-correction.md for the distribution and its limits. Falling back to nominal scale is not the better answer at any coverage: the DEM route beats it in the median even below 20% covered.

A coastal frame is sized from land and sea together. MRDEM carries near-shore sea at about 0.14 m, not nodata, so the mean takes in the sea surface the photo images and dem_coverage stays near 1. Against a ray-cast of the true footprint on bare earth, that is closer on area than the land-only mean, and the land-only mean places the land edge better. The first holds, to first order, under the canopy BC actually has (fly#80), and matched for relief the sea does make the land edge worse than it is inland. LidarBC was mostly nodata over sea where probed, so there the same frame is sized nearer the land-only mean. See inst/notes/terrain-correction.md (fly#65).

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.

A remote DEM is read a window at a time, and each window is one footprint — so cost scales with the number of frames and never with the distance between them. Cropping a /vsicurl/ raster to your AOI first still saves a little where the frames are close together (measured: 1.00 s against 1.23 s for eight contiguous frames), but it allocates for the whole AOI rather than one footprint, and it stops paying as the frames spread out — a crop of the bundled extent costs 4.4 s on its own. Worth doing to work offline; not worth doing for speed.

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"))
#> Reading layer `photo_centroids' from data source 
#>   `/home/runner/work/_temp/Library/fly/testdata/photo_centroids.gpkg' 
#>   using driver `GPKG'
#> Simple feature collection with 20 features and 16 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: -126.7631 ymin: 54.34512 xmax: -126.449 ymax: 54.47635
#> Geodetic CRS:  WGS 84
footprints <- fly_footprint(centroids)
plot(sf::st_geometry(footprints))


# How each footprint was sized
table(footprints$footprint_basis)
#> 
#> Film - BW 
#>        20 

# Keep only frames sized from a known recording format
sized <- footprints[footprints$footprint_basis != "unknown_format", ]
nrow(sized)
#> [1] 20

# 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))
}
#>  [1] 16.9 11.1 15.3  4.9  0.5 13.4 14.0  2.1 18.1 12.8 16.1 16.0 15.0 11.0  9.0
#> [16] 26.4  6.2  9.1 16.3 15.5
#> 
#> dem_agl 
#>      20