ras_commander_fast_map/
terrain_scope.rs

1//! Mapping on the top-priority source of a multi-source terrain, without a
2//! RAS Mapper terrain export (`top-priority-source-footprint/1`).
3//!
4//! RAS Mapper writes stored maps per terrain source, on that source's own grid.
5//! Inside the highest-priority source's footprint (bit 0 of its level-0 tile
6//! `Mask` in the terrain HDF), the effective terrain is that source's original
7//! pixels: on BaldEagleCrkMulti2D p07 (two sources), the RAS Mapper 1x export
8//! equals the priority-1 TIFF at every footprint pixel, and mapping directly on
9//! that TIFF matches native HEC-RAS 7.0.1 stored maps in every valid pixel and
10//! float32 bit (Max and one output time) without HEC-RAS.
11//!
12//! What is not reproduced, and therefore refused rather than approximated:
13//! pixels outside the footprint. RAS Mapper fills them from lower-priority
14//! sources (its export resamples them with a rule close to, but not,
15//! pixel-centre bilinear) and stitches; the footprint also drops valid source
16//! pixels (23,329 on p07, including 229 single-pixel slivers that the export
17//! fills from the coarser source).
18//!
19//! Foreign pixels are therefore allowed only where the result cannot depend on
20//! their terrain: in sloping rendering, pixels of cells that fail the native
21//! cell wet gate are never drawn. Any foreign pixel of a wet cell and any
22//! foreign pixel in horizontal rendering (no cell minimum is supplied) is
23//! refused.
24use crate::{MapError, model, terrain::NativeTerrain};
25use serde::Deserialize;
26use serde_json::{Value, json};
27use std::{
28    fs::File,
29    io::{Read, Seek, SeekFrom},
30    path::Path,
31};
32
33pub const METHOD_ID: &str = "top-priority-source-footprint/1";
34const TILE: usize = 256;
35
36/// Level-0 tile mask of the source (terrain HDF `Terrain/<source>/0/Mask`).
37#[derive(Clone, Debug, Deserialize)]
38#[serde(deny_unknown_fields)]
39pub struct MaskArray {
40    pub path: String,
41    pub dtype: String,
42    /// `[tiles_y * tiles_x, 65536]`, tiles row-major over the source raster.
43    pub shape: [usize; 2],
44}
45
46/// Declares how the supplied `terrain` relates to the registered terrain.
47#[derive(Clone, Debug, Deserialize)]
48#[serde(deny_unknown_fields)]
49pub struct TerrainScope {
50    /// `top_priority_source`: `terrain` is the original TIFF of the
51    /// highest-priority source of a registered terrain with `source_count` sources.
52    pub kind: String,
53    pub source_count: usize,
54    pub mask: MaskArray,
55}
56
57impl TerrainScope {
58    pub fn validate(&self) -> Result<(), MapError> {
59        if self.kind != "top_priority_source"
60            || self.source_count == 0
61            || self.mask.dtype != "|u1"
62            || self.mask.shape[1] != TILE * TILE
63            || self.mask.shape[0] == 0
64        {
65            return Err(MapError::Invalid(
66                "terrain_scope needs kind top_priority_source, source_count >= 1 and a |u1 [tiles, 65536] mask".into(),
67            ));
68        }
69        Ok(())
70    }
71}
72
73pub fn method_record() -> Value {
74    json!({
75        "method_id": METHOD_ID,
76        "terrain": "original float32 TIFF of the highest-priority registered source, 256 x 256 tiles; no export, resampling or mosaic",
77        "footprint": "bit 0 of the source's level-0 terrain-HDF tile mask, with valid source data",
78        "foreign_pixels": "owned pixels outside the footprint; allowed only in sloping cells that fail the native wet gate",
79        "refused": "foreign pixels of wet cells, foreign pixels in horizontal rendering, cell-face sloping",
80        "qualified": "BaldEagleCrkMulti2D p07 (2 sources), native HEC-RAS 7.0.1 plain sloping Max and one output time"
81    })
82}
83
84/// Footprint of the window: mask bit 0 and valid source data, per window pixel.
85fn footprint(
86    scope: &TerrainScope,
87    base: &Path,
88    terrain: &NativeTerrain,
89    window: [usize; 4],
90) -> Result<Vec<bool>, MapError> {
91    let [col, row, width, height] = window;
92    let (tiles_x, tiles_y) = (terrain.width.div_ceil(TILE), terrain.height.div_ceil(TILE));
93    if terrain.block != (TILE, TILE) || scope.mask.shape[0] != tiles_x * tiles_y {
94        return Err(MapError::Unsupported(
95            "terrain_scope needs the original 256 x 256 tiled source and its full level-0 mask"
96                .into(),
97        ));
98    }
99    let path = model::resolve(base, &scope.mask.path);
100    let expected = (scope.mask.shape[0] * TILE * TILE) as u64;
101    if std::fs::metadata(&path)?.len() != expected {
102        return Err(MapError::Invalid(
103            "terrain_scope mask length mismatch".into(),
104        ));
105    }
106    let mut file = File::open(&path)?;
107    let mut inside = vec![false; width * height];
108    let mut tile = vec![0_u8; TILE * TILE];
109    for ty in row / TILE..=(row + height - 1) / TILE {
110        for tx in col / TILE..=(col + width - 1) / TILE {
111            file.seek(SeekFrom::Start(((ty * tiles_x + tx) * TILE * TILE) as u64))?;
112            file.read_exact(&mut tile)?;
113            let block = terrain.get_block(tx, ty)?;
114            let win = terrain.block_window(tx, ty);
115            for r in row.max(win.y)..(row + height).min(win.y + win.height) {
116                for c in col.max(win.x)..(col + width).min(win.x + win.width) {
117                    let (lr, lc) = (r - win.y, c - win.x);
118                    let valid = terrain
119                        .classify(block[lr * win.width + lc])
120                        .is_some_and(f64::is_finite);
121                    inside[(r - row) * width + c - col] = valid && tile[lr * TILE + lc] & 1 != 0;
122                }
123            }
124        }
125    }
126    Ok(inside)
127}
128
129/// Result of the terrain-scope guard on one window.
130pub struct Scoped {
131    pub record: Value,
132    /// 2D-owned window pixels outside the footprint (window order). They may
133    /// be mapped only while their cell fails the wet gate; empty for one source.
134    pub foreign: Vec<usize>,
135}
136
137/// Classifies owned pixels against the top-priority source's footprint.
138pub fn check(
139    scope: &TerrainScope,
140    base: &Path,
141    terrain: &NativeTerrain,
142    window: [usize; 4],
143    owner: &[u32],
144) -> Result<Scoped, MapError> {
145    scope.validate()?;
146    let [_, _, width, height] = window;
147    if owner.len() != width * height {
148        return Err(MapError::Internal("ownership/window size mismatch".into()));
149    }
150    let inside = footprint(scope, base, terrain, window)?;
151    let (owned, foreign) = classify(&inside, owner);
152    let record = json!({"method": method_record(), "source_count": scope.source_count,
153        "owned_pixels": owned, "foreign_2d_pixels": foreign.len()});
154    if scope.source_count == 1 {
155        return Ok(Scoped {
156            record,
157            foreign: Vec::new(),
158        });
159    }
160    Ok(Scoped { record, foreign })
161}
162
163/// Fails when a 2D cell that owns a foreign pixel is wet: the pixel's terrain
164/// would come from a lower-priority source or stitch, which is unqualified.
165pub fn deny_wet_foreign(
166    foreign: &[usize],
167    owner: &[u32],
168    owner_cells: &[[usize; 2]],
169    wet: impl Fn(usize, usize) -> bool,
170) -> Result<(), MapError> {
171    let wet_pixels = foreign
172        .iter()
173        .filter(|&&p| {
174            let [mesh, cell] = owner_cells[owner[p] as usize - 1];
175            wet(mesh, cell)
176        })
177        .count();
178    if wet_pixels > 0 {
179        return Err(MapError::Unsupported(format!(
180            "{wet_pixels} pixels of wet cells lie outside the top-priority source's footprint; lower-priority sources and stitches are not qualified"
181        )));
182    }
183    Ok(())
184}
185
186/// Owned pixel count and 2D-owned pixels outside the footprint (window order).
187fn classify(inside: &[bool], owner: &[u32]) -> (u64, Vec<usize>) {
188    let mut owned = 0_u64;
189    let mut foreign = Vec::new();
190    for (p, _) in owner.iter().enumerate().filter(|(_, label)| **label != 0) {
191        owned += 1;
192        if !inside[p] {
193            foreign.push(p);
194        }
195    }
196    (owned, foreign)
197}
198
199#[cfg(test)]
200mod tests {
201    use super::*;
202
203    fn scope(kind: &str, count: usize, dtype: &str, tiles: usize) -> TerrainScope {
204        TerrainScope {
205            kind: kind.into(),
206            source_count: count,
207            mask: MaskArray {
208                path: "mask.u1".into(),
209                dtype: dtype.into(),
210                shape: [tiles, 65536],
211            },
212        }
213    }
214
215    #[test]
216    fn footprint_classification_of_2d_pixels() {
217        let inside = [true, false, true, false, true, true];
218        let owner = [1, 1, 0, 2, 0, 0];
219        assert_eq!(classify(&inside, &owner), (3, vec![1, 3]));
220        assert_eq!(classify(&inside, &[0; 6]), (0, vec![]));
221    }
222
223    #[test]
224    fn foreign_pixels_are_allowed_only_in_dry_cells() {
225        let owner = [0, 1, 2, 2];
226        let cells = [[0, 7], [1, 3]];
227        // Pixel 1 belongs to mesh 0 cell 7; pixels 2-3 to mesh 1 cell 3.
228        assert!(deny_wet_foreign(&[1, 3], &owner, &cells, |_, _| false).is_ok());
229        let err = deny_wet_foreign(&[1, 3], &owner, &cells, |m, c| (m, c) == (1, 3)).unwrap_err();
230        assert_eq!(err.kind(), "unsupported");
231        assert!(err.to_string().starts_with("unsupported: 1 pixels"));
232        assert!(deny_wet_foreign(&[], &owner, &cells, |_, _| true).is_ok());
233    }
234
235    #[test]
236    fn scope_fields_are_validated() {
237        assert!(scope("top_priority_source", 2, "|u1", 4).validate().is_ok());
238        for bad in [
239            scope("mosaic", 2, "|u1", 4),
240            scope("top_priority_source", 0, "|u1", 4),
241            scope("top_priority_source", 2, "<f4", 4),
242            scope("top_priority_source", 2, "|u1", 0),
243        ] {
244            assert_eq!(bad.validate().unwrap_err().kind(), "invalid");
245        }
246    }
247}
ras_commander_fast_map/terrain_scope.rs

ras_commander_fast_map/
terrain_scope.rs
