DEM-based terrain-adjusted footprints - #34
Merged
NewGraphEnvironment merged 15 commits intoAug 29, 2026
Conversation
MRDEM-30 is NRCan's 30 m bare-earth DTM, the same product flooded::fl_dem_aoi() defaults to. Fetched via /vsicurl/ so only the intersecting bytes transfer. 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, so the choice is about provenance and dependencies rather than accuracy: MRDEM needs no dependency beyond terra, is bare earth rather than a mixed-provenance mosaic, and is a source users can reach for their own AOIs. 973x1086 cells at 30.5 m, 566-1520 m, 306 KB. Relates to #9 Co-Authored-By: Claude Opus 5 <[email protected]> Claude-Session: https://claude.ai/code/session_01GBKqedyBysV7hB4DuL98ZR
Reported scale assumes flat ground at whatever elevation the scale was computed for. On the bundled Upper Bulkley AOI that understates footprint area by a median 14%, ranging to 27%, and always in the same direction — the scale is referenced to an elevation above the valley floor the photos actually cover. This is a datum offset, not the slope effect #9 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 sampled twice — at the centroid, then as the mean under the resulting rectangle — because on a 7.2 km 1:31680 frame the two differ by up to 130 m. Geometry stays rectangular, so downstream consumers are unaffected. Per-corner ray-casting measures ~2% against the 14% this addresses, and costs irregular geometry; deferred. footprint_terrain records the treatment per frame and height_agl makes it auditable. footprint_basis is left alone deliberately: it is already matched by value downstream, so encoding terrain into it as #9 suggested would break caller filters. Every failure falls back to nominal scale with a warning rather than dropping the frame — outside DEM coverage, and unusable flying_height/focal_length. That second case classifies on the computed half-side rather than its inputs: an NA or zero focal length yields a non-finite half-side that would otherwise become an empty geometry, indistinguishable from an unresolved recording format, sending the user to format_size for a metadata problem. Relates to #9 Co-Authored-By: Claude Opus 5 <[email protected]> Claude-Session: https://claude.ai/code/session_01GBKqedyBysV7hB4DuL98ZR
fly_footprint() gained a dem argument, but all six internal call sites called it bare — so the terrain correction was unreachable from fly_coverage(), fly_overlap(), fly_filter(), fly_select() and fly_georef(), which is every way a user actually reaches a footprint. fly_select() threads it through both fly_select_all() and fly_select_minimal(); they are separate call sites, so passing in one mode would prove nothing about the other. The tests assert the numbers move rather than that the argument is tolerated. A dem accepted and silently dropped is invisible otherwise — each of these returns entirely plausible output either way. fly_georef() has no observable GCPs without images to warp, so it is checked by stripping flying_height and asserting the error, which can only arise if dem reached fly_footprint(). Relates to #9 Co-Authored-By: Claude Opus 5 <[email protected]> Claude-Session: https://claude.ai/code/session_01GBKqedyBysV7hB4DuL98ZR
The flat-terrain caveats in fly_footprint(), fly_overlap(), fly_georef() and the vignette asserted an assumption the package no longer has to make. Each is now conditional on whether a dem was supplied, and fly_footprint() gains a Terrain section carrying the measured effect, the fallback semantics, and three BC DEM sources. The vignette gains a worked comparison on the bundled AOI — medians of 13.3% and 15.6% by scale — plus the MRDEM-30 recipe for a reader's own area. The one-directional result is called out, since that is what distinguishes a datum offset from the slope effect the issue described. Issue #9's body was edited rather than commented on, per convention: it had argued for the slope framing and suggested encoding terrain into footprint_basis, both superseded by measurement. Earlier corrections and the original body stay preserved in the collapsed block. Fixes #9 Co-Authored-By: Claude Opus 5 <[email protected]> Claude-Session: https://claude.ai/code/session_01GBKqedyBysV7hB4DuL98ZR
Co-Authored-By: Claude Opus 5 <[email protected]> Claude-Session: https://claude.ai/code/session_01GBKqedyBysV7hB4DuL98ZR
Committed accidentally by `git add -A` in two earlier commits on this branch. It is a side effect of running examples that plot, not source, and as a top-level file absent from .Rbuildignore it would ship in the built package. Co-Authored-By: Claude Opus 5 <[email protected]> Claude-Session: https://claude.ai/code/session_01GBKqedyBysV7hB4DuL98ZR
R CMD build ships every top-level directory absent from .Rbuildignore, so planning/ was landing in the tarball and therefore in the library of anyone installing from GitHub. Pre-existing — main shipped 11 such entries — but this branch archived a PWF and took it to 16, so it is fixed here rather than left to grow. Verified against the tarball rather than the config: 0 planning entries, with inst/testdata/dem.tif still present. Co-Authored-By: Claude Opus 5 <[email protected]> Claude-Session: https://claude.ai/code/session_01GBKqedyBysV7hB4DuL98ZR
A footprint hanging over the edge of the DEM still produced a mean elevation, taken with na.rm from whichever part had data, and was reported as "dem_agl" — an affirmative claim of a full-frame correction. Measured on a constructed case: a frame 50% off the DEM was corrected from a mean of 677 m where its centroid read 861 m, silently. Same class as the empty-geometry bug earlier on this branch: an operation that partly failed presenting as one that wholly succeeded. The partial mean is kept — it is the best estimate available, and discarding the frame would lose more than it protects — but dem_coverage now reports the fraction per frame so it can be filtered, and anything under 95% warns. The threshold is not decoration. Warning on any missing cell at all fires on 3 of the 20 bundled frames at 96% coverage, because reprojecting a DEM leaves NA slivers along its edges; a guard that noisy stops being read. Missing more than a twentieth of a footprint is a different claim. Also corrects the documentation, which asserted that a DEM stopping short of the frame edges sends those frames to no_dem_coverage. It does not — that path catches only a frame whose centroid has no elevation. Relates to #9 Co-Authored-By: Claude Opus 5 <[email protected]> Claude-Session: https://claude.ai/code/session_01GBKqedyBysV7hB4DuL98ZR
Partial DEM coverage, blind passthrough tests, a wrong buffer geometry, a mislabelled terrain column, and an unguarded optional dependency. - The DEM is now sampled as the mean over the whole footprint in two genuine passes, gated on the rectangle rather than on the centroid. The old centroid sample sized nothing — it was only a coverage gate — so a single NA cell under a centroid condemned a frame that was 99% covered, and the documented "two passes" described an iteration the code never performed. A real second pass moves elevation by up to 10.9 m; a third moves area by under 0.02%. - fly_filter and fly_select passthrough tests could not detect a dropped dem. Their expect_gte/expect_lte comparisons hold for both implementations on the bundled data — 20 of 20 kept either way — so deleting the passthrough left the suite green. Both now also assert arrival directly, fly_select once per internal call site. Verified by patching fly_footprint to discard dem: all five tests fail, where two previously passed. - The test DEM was buffered by the footprint half-side, 3.6 km. A square's far point is its corner, half_side * sqrt(2) = 5.1 km, so every edge frame had its corners over no-data. Regenerated at 5.4 km, and the same error is corrected in the buffering advice given to users. - footprint_terrain was keyed on the recording format, so a frame with a resolvable media and an unparseable scale got "nominal_scale" despite having no geometry. Keyed on the footprint now. - terra is Suggests-only, but the example and four vignette chunks called the DEM path unguarded, so R CMD check and pkgdown would hard-fail without it. Regenerating the fixture moved the measured numbers: median area correction 13.8% (was 14.1%), range 0.5-26.4%, centroid-vs-mean 140 m. Every documented figure updated to match. Relates to #9 Co-Authored-By: Claude Opus 5 <[email protected]> Claude-Session: https://claude.ai/code/session_01GBKqedyBysV7hB4DuL98ZR
Round 2 found the round-1 partial-coverage fix incomplete, and wrong in the more dangerous direction. Coverage was the non-NA share of the cells terra::extract() returned. Ground beyond the raster's *extent* yields no row at all rather than an NA row, so a footprint running off the edge of the data was reported as fully covered. The bundled fixture could only ever reach the other case — NA slivers left inside the extent by reprojection — so both tests written for the fix passed while the real failure stayed invisible. That is the ordinary case, not an exotic one. A DEM cropped to an AOI, which is what fl_dem_aoi() and produce, has no NA interior: it simply stops. Measured on one: all 20 frames reported coverage 1 with no warning, while height_agl was wrong by up to 83 m and footprint area by 3.3%. A reported 1 is an affirmative claim of completeness, so the fix had made this worse than the bug it replaced. Coverage is now the non-NA count over the cells the footprint's area implies, capped at 1 since extract() takes cells by centre. Verified in both directions: the well-buffered bundled DEM reports 0.9945-1 and warns on nothing, while an AOI-clipped one reports down to 0.56 and flags 6 of 20. dem_coverage is also reported for every frame that had a footprint, not only corrected ones. no_dem_coverage means a measured zero, and returning NA there made the documented "filter on dem_coverage" workflow impossible. Three documentation claims corrected: the roxygen Terrain section and the vignette still described the centroid gate removed in the previous commit, and NEWS still described the pre-fix sampling scheme. The new test uses an AOI-clipped DEM and asserts the fixture has no NA interior, so it reaches the extent case rather than the sliver case. Restoring the round-1 version fails 4 tests, with the patch proven active. Relates to #9 Co-Authored-By: Claude Opus 5 <[email protected]> Claude-Session: https://claude.ai/code/session_01GBKqedyBysV7hB4DuL98ZR
Round 3 found the round-2 coverage fix wrong for a geographic DEM, and its fixture unable to notice. Coverage compares a footprint's area against the DEM's cell size. st_area() on a geographic CRS returns geodesic metres squared while terra::res() returns degrees, so the ratio was meaningless: an EPSG:4326 DEM reported coverage of 1.4e-10 for twenty fully-covered frames and warned on all of them. The footprints themselves stayed correct; what broke was the column and the guard, so the documented filter would have dropped everything and the guard became the noise its own comment forbids. Every DEM in the suite was EPSG:3005 or a crop of it, which makes the reprojection inside fly_dem_sample() an identity in 100% of tests, so the whole branch was unexecuted. This is the third round in which a fix's own tests could not reach the failure mode they were written for. The new test uses a reprojected geographic DEM and asserts its height_agl agrees with the projected one, since both describe the same ground. Dropping the CRS forces planar area, which is degrees squared for a geographic DEM and metres squared for a projected one, matching res() in both. Also fixed: the partial-coverage warning told the user to buffer by half the widest footprint, which is the under-buffering corrected two commits ago and would re-trigger the warning it appears in. Every other statement of that advice already says to clear the corner. And a premise assertion that could not fail: it compared the cropped DEM against a crop of itself, so both sides were equal for any fixture, and the premise it stated was false besides. It now asserts what the test actually needs, that some footprint extends past the DEM's extent. Restoring the units defect reads 1.386e-10 where the fix reads ~1, and fails the new test. Relates to #9 Co-Authored-By: Claude Opus 5 <[email protected]> Claude-Session: https://claude.ai/code/session_01GBKqedyBysV7hB4DuL98ZR
Round 4 found the round-3 coverage measure wrong on its own terms, with no missing data involved at all. terra::extract() takes a cell when its centre falls inside the polygon, so the numerator was a count of centres while the denominator was the footprint's area in cell units. Those are different measurements, and the gap runs about 2/k for a footprint k cells wide. On a synthetic DEM with no NA cell anywhere and 30 km of room beyond every frame, coverage reported 0.913 at 900 m cells and warned that 3 of 3 frames were under-covered. Nothing was missing. The denominator is now a centre count too, taken on a grid aligned to the DEM's own via terra::align() and covering only the footprints. Extending the DEM instead would give the same answer and was the obvious move — extend, then count NAs. Measured before taking it: for one frame outside a local DEM that allocates 260 million cells, a 196x blowup, because extend sizes to the union of raster and footprints. The template is only ever as large as the frames. The zero-focal-length test caught a regression from that change: an infinite half-side reached the sampler, which tried to size a raster to an infinite extent and failed with "vector". fly_rectangles() now treats any non-finite half-side as empty, not only NA. Restoring the round-3 measure now fails 4 tests. It failed none when the fix first landed — the shipped 30 m fixture is too fine to expose a 2/k error, so the new tests use 30 m and 900 m grids, plus anisotropic 120x904 m cells. Also from round 4: - The roxygen justified the second sampling pass with the correction's own 25%, which belongs to the first pass. The iteration is worth under 0.5% of area and now says so. Collapsing it to one pass had left the suite green, so it is now asserted against an independently computed one-pass elevation. - A test comment claimed reprojection slivers made coverage less than 1 on the bundled DEM. The measured NA fraction is zero for nineteen of twenty frames; the assertion had been passing on the counting artifact above. - A DEM with no CRS died inside sf as "invalid crs:", naming neither the argument nor the package. Relates to #9 Co-Authored-By: Claude Opus 5 <[email protected]> Claude-Session: https://claude.ai/code/session_01GBKqedyBysV7hB4DuL98ZR
Round 5, two defects. The coverage denominator was counted on one grid spanning every frame, so it was sized by the GAP between frames rather than by the frames. That is the same allocation the previous commit set out to remove, and this package's own test fixture triggers it: two photos 700 km apart built a 243-million-cell template where the same two counted separately need 16 thousand, a factor of 14,950. The suite was peaking around 4.6 GB. The comment claiming the template was "only ever as large as the frames" was simply wrong. Counted per frame now; the distant-frame call goes from seconds and gigabytes to 0.18 s. Second, dem_coverage could describe a rectangle the caller never receives. Terrain at or above the aircraft gives a negative height above ground, so the second pass measures over a square sized from that — far smaller than the nominal footprint the frame then falls back to. Coverage came from the second pass regardless, which reported a fully-covered small square for a footprint that was 30% covered. That frame is exactly what the documented dem_coverage filter exists to exclude, and it passed the filter. Coverage now comes from the pass whose rectangle is actually returned. Both have tests that fail when the defect is restored: the first asserts elapsed time, since the failure is an allocation rather than a wrong number; the second measures ground truth against the returned geometry and asserts the fixture reaches the failure mode before asserting the value. Relates to #9 Co-Authored-By: Claude Opus 5 <[email protected]> Claude-Session: https://claude.ai/code/session_01GBKqedyBysV7hB4DuL98ZR
Round 6, one finding, and it is the guard from the previous commit. The test written for round 5's allocation defect could not fail. It asserted `expect_lt(elapsed, 10)`; 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 existed to catch. Its other two assertions passed under the defect as well, because a union-sized grid produces the *right* denominator — it just allocates 243 million cells and about 4 GB to get there. The previous commit message claimed both round-5 fixes had tests that fail when the defect is restored. For this one that was false. No threshold repairs it: 0.18 s against 1.0 s is CI jitter. So the invariant is now asserted directly. The grid construction is named, `fly_dem_grid()`, and the test mocks it to record every grid the call asks for, then asserts that none approaches the grid spanning both frames. That is the quantity that differs by 14,950x, rather than a timing proxy for it. Verified against the exact prior implementation taken from git rather than rewritten from memory — which matters, because a hand-reconstruction of it failed 4 tests and looked like a working guard. The difference was `length()` versus `sum(!is.na())` on the template extract; the reconstruction was a different bug. Restored properly, the old test suite passes 222/222 while the new test fails. Round 6 confirmed the rest under execution rather than reading: every documented figure recomputes (median 13.77%, range 0.45-26.39%, 140.06 m, second pass 0.529%, third 0.034%), inst/testdata/dem.tif is byte-equivalent to a live MRDEM-30 regeneration, first/second indexing is aligned with `sized` including with holes in the middle, and dem threading is clean through all six downstream call sites with empty geometries present. Relates to #9 Co-Authored-By: Claude Opus 5 <[email protected]> Claude-Session: https://claude.ai/code/session_01GBKqedyBysV7hB4DuL98ZR
NewGraphEnvironment
deleted the
9-dem-based-terrain-adjusted-footprints
branch
August 29, 2026 23:29
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
What
fly_footprint()gains ademargument that sizes each frame from its height above ground instead of the reported map scale, threaded through every function that builds a footprint. Released as 0.5.0.Fixes #9
Why the issue was wrong about its own premise
FLYING_HEIGHTis metres above sea level, not height above ground. The values are round feet — 8000, 8500, 19500, 20000 ft — and the implied ground elevations bracket real Houston-area terrain. Subtracting terrain elevation is what turns the field into the height ground coverage scales with.The error is a datum offset, not a slope effect. Reported scale understates footprint area by a median 13.8%, ranging to 26.4%, and always in the same direction — the scale is referenced to an elevation above the valley floor the photos cover. Slope would push both ways. The issue was filed describing slope; the bulk bias is larger and is what the DEM removes.
Approach
Issue option 3 without the irregular geometry:
h_agl = flying_height − terrain elevation, keeping a rectangle so downstream consumers are unaffected. Per-corner ray-casting measured ~2% against the 14% the scale correction addresses — deferred, not rejected.Elevation is the mean under the whole footprint, not a centroid reading (they differ by up to 140 m on a 7.2 km frame), measured in two passes because the window being averaged is itself what the correction changes.
DEM source: MRDEM-30, NRCan's 30 m bare-earth DTM — public S3 COG, no auth, the product
flooded::fl_dem_aoi()already defaults to. Chosen overelevatrafter a head-to-head: they agree within 0.42 percentage points, so it costs no accuracy and needs no dependency beyondterra(Suggests).New columns
footprint_terrain(nominal_scale/dem_agl/no_dem_coverage/NA),height_agl, anddem_coverage. Terrain deliberately did not go intofootprint_basisas the issue suggested — that column is already matched by value downstream, so a suffix would silently break caller filters.Every failure falls back to nominal scale with a warning rather than dropping the frame.
Review
Five rounds of adversarial review ran. Every round found a real defect in the previous round's fix — fourteen findings total, all fixed, each with a test verified to fail when the defect is restored.
dem; a centroid NA condemned a 99%-covered frame;footprint_terrainkeyed on format not footprint;terraunguarded in example/vignette; fixture buffered by half-side when a square's corner ishalf_side × √2height_aglwas wrong by 83 m;dem_coverageNA where 0 was the truth; 3 stale doc claimsdem_coveragecould describe a rectangle the caller never receives — 100% reported for a 30%-covered footprintThe recurring cause was the same each time: the fixture could not reach the failure mode. The bundled 30 m EPSG:3005 DEM cannot exercise a coarse grid, a geographic CRS, a truncating extent, or a wide photo spread. Tests now cover all four.
Two fixes were rejected after measuring rather than shipped:
terra::extend()(260 M cells for one distant frame) and a content-based partial-coverage threshold at "any missing cell" (fires on reprojection slivers at 99.96%).Verification
R CMD check: 0 errors, 0 warnings, 0 notes —maincarries 2 notes, cleared here by removing a strayRplots.pdfand keepingplanning/out of the buildNote
Commit
13ae908's message reads "fl_dem_aoi() and produce" — a backtick code span in a-mstring was command-substituted. Left as-is; commit messages are immutable history by convention.🤖 Generated with Claude Code
https://claude.ai/code/session_01GBKqedyBysV7hB4DuL98ZR