ras_commander_fast_map/
terrain.rs

1//! Native same-CRS north-up GeoTIFF terrain through terrio 0.4.2 lazy reads.
2//!
3//! Serial decode guarantee: every terrio `read_window` call requests exactly one TIFF block
4//! (tile or strip) clipped to the image. terrio's `decode_in_threads` spawns
5//! `rayon::current_num_threads().min(chunks)` scoped threads; with one chunk that is one, so the
6//! chunk is decoded on the calling worker thread and the outer pool is never multiplied. This
7//! holds regardless of `RAYON_NUM_THREADS` or the calling pool size. Evidence is the pinned
8//! terrio `=0.4.2` source guard (`decode_in_threads` spawns only when `workers > 1`); re-verify it
9//! on any terrio upgrade. The `debug_assert` in `block_window` only checks this crate's own grid.
10//!
11//! Size arithmetic is checked: block bytes are computed once at open with checked
12//! multiplication and absurd layouts fail with a typed `unsupported` error instead of panicking.
13use crate::MapError;
14use crate::cache::{CacheStats, TileCache};
15use crate::model::{CrsIdentity, NodataPolicy, QueryConvention};
16use serde::Serialize;
17use std::sync::atomic::{AtomicU64, Ordering};
18use std::sync::{Arc, Mutex};
19use std::time::Instant;
20use terrio::{LazyRaster, Window, WindowOptions};
21
22pub struct NativeTerrain {
23    pub path: String,
24    pub gt: [f64; 6],
25    pub width: usize,
26    pub height: usize,
27    pub block: (usize, usize),
28    pub nodata: Option<f64>,
29    pub epsg: u32,
30    pub crs_identity: &'static str,
31    pub crs_wkt_fnv1a64: Option<String>,
32    pub query: QueryConvention,
33    pub policy: NodataPolicy,
34    /// Bytes of one full decoded block (checked at open).
35    block_bytes: u64,
36    readers: Mutex<Vec<LazyRaster>>,
37    readers_opened: AtomicU64,
38    decode_nanos: AtomicU64,
39    pub cache: TileCache,
40}
41
42#[derive(Serialize, Debug, Clone)]
43pub struct TerrainInfo {
44    pub path: String,
45    pub geo_transform: [f64; 6],
46    pub width: usize,
47    pub height: usize,
48    pub block: (usize, usize),
49    pub blocks_x: usize,
50    pub blocks_y: usize,
51    pub nodata: Option<f64>,
52    pub epsg: Option<u32>,
53    pub data_type: String,
54    pub band_count: usize,
55    pub crs: terrio::crs::SpatialRef,
56    pub crs_wkt_fnv1a64: Option<String>,
57}
58
59#[derive(Serialize, Debug, Clone)]
60pub struct TerrainStats {
61    /// Decoded *resident* tile accounting (the enforced cache bound).
62    pub cache: CacheStats,
63    pub readers_opened: u64,
64    pub decode_seconds: f64,
65    /// Estimate (not enforced, not measured) of decoder transients outside the cache budget:
66    /// `peak_concurrent_loads x decoder_transient_bytes_per_load_estimate`.
67    pub decoder_transient_bytes_per_load_estimate: u64,
68    pub decoder_transient_peak_estimate_bytes: u64,
69}
70
71/// Index of a terrain pixel hit by a sample point.
72#[derive(Clone, Copy, Debug, PartialEq)]
73pub enum Hit {
74    Outside,
75    Pixel {
76        col: usize,
77        row: usize,
78        near_boundary: bool,
79    },
80}
81
82pub fn inspect(path: &str) -> Result<TerrainInfo, MapError> {
83    let r = open_lazy(path)?;
84    let block = r.block_size().unwrap_or((0, 0));
85    Ok(TerrainInfo {
86        path: path.to_string(),
87        geo_transform: r.geo_transform(),
88        width: r.width(),
89        height: r.height(),
90        block,
91        blocks_x: if block.0 > 0 {
92            r.width().div_ceil(block.0)
93        } else {
94            0
95        },
96        blocks_y: if block.1 > 0 {
97            r.height().div_ceil(block.1)
98        } else {
99            0
100        },
101        nodata: r.nodata(),
102        epsg: r.crs().epsg,
103        data_type: format!("{:?}", r.data_type()),
104        band_count: r.band_count(),
105        crs: r.crs(),
106        crs_wkt_fnv1a64: wkt_fnv(&r.crs()),
107    })
108}
109
110/// FNV-1a 64 of terrio's WKT for the terrain CRS (identity pin for caller-verified CRS).
111pub fn wkt_fnv(crs: &terrio::crs::SpatialRef) -> Option<String> {
112    crs.to_wkt()
113        .map(|w| format!("{:016x}", crate::sink::fnv1a64(w.as_bytes())))
114}
115
116fn open_lazy(path: &str) -> Result<LazyRaster, MapError> {
117    let r = LazyRaster::open_tiff(path)
118        .map_err(|e| MapError::Terrain(format!("{path}: {e}")))?
119        .ok_or_else(|| {
120            MapError::Unsupported(format!(
121                "{path}: not a TIFF (VRT/HDF/other terrain unsupported)"
122            ))
123        })?;
124    if !r.is_lazy() {
125        return Err(MapError::Unsupported(format!(
126            "{path}: terrio did not open a lazy GeoTIFF"
127        )));
128    }
129    Ok(r)
130}
131
132impl NativeTerrain {
133    pub fn open(
134        path: &str,
135        expected_epsg: u32,
136        identity: &CrsIdentity,
137        query: QueryConvention,
138        policy: NodataPolicy,
139        cache_budget: u64,
140    ) -> Result<Self, MapError> {
141        let r = open_lazy(path)?;
142        if r.band_count() != 1 {
143            return Err(MapError::Unsupported(format!(
144                "{path}: {} bands; single-band required",
145                r.band_count()
146            )));
147        }
148        if r.data_type() != terrio::DataType::F32 {
149            return Err(MapError::Unsupported(format!(
150                "{path}: data type {:?}; float32 required",
151                r.data_type()
152            )));
153        }
154        let gt = r.geo_transform();
155        if !gt.iter().all(|v| v.is_finite())
156            || gt[2] != 0.0
157            || gt[4] != 0.0
158            || gt[1] <= 0.0
159            || gt[5] >= 0.0
160        {
161            return Err(MapError::Unsupported(format!(
162                "{path}: geotransform {gt:?} is not north-up axis-aligned (rotated/south-up unsupported)"
163            )));
164        }
165        let crs = r.crs();
166        let crs_identity = match (crs.epsg, identity) {
167            (Some(e), _) if e != expected_epsg => {
168                return Err(MapError::Unsupported(format!(
169                    "{path}: terrain EPSG:{e} differs from session EPSG:{expected_epsg}; no reprojection is performed"
170                )));
171            }
172            (Some(_), CrsIdentity::Epsg) => "terrio_identified_epsg",
173            // A supplied WKT pin is always verified, even when terrio also identifies the EPSG.
174            (Some(_), CrsIdentity::CallerVerifiedWkt { fnv1a64 }) => {
175                check_wkt_pin(path, &crs, fnv1a64)?;
176                "terrio_identified_epsg_and_caller_verified_wkt_fnv1a64"
177            }
178            (None, CrsIdentity::Epsg) => {
179                return Err(MapError::Unsupported(format!(
180                    "{path}: terrio does not identify an EPSG code for this terrain; use crs_identity caller_verified_wkt after independent verification"
181                )));
182            }
183            (None, CrsIdentity::CallerVerifiedWkt { fnv1a64 }) => {
184                check_wkt_pin(path, &crs, fnv1a64)?;
185                "caller_verified_wkt_fnv1a64"
186            }
187        };
188        let block = r
189            .block_size()
190            .filter(|b| b.0 > 0 && b.1 > 0)
191            .ok_or_else(|| MapError::Unsupported(format!("{path}: no block layout")))?;
192        let (width, height) = (r.width(), r.height());
193        if width == 0 || height == 0 || width > u32::MAX as usize || height > u32::MAX as usize {
194            return Err(MapError::Unsupported(format!(
195                "{path}: unsupported size {width}x{height}"
196            )));
197        }
198        let block_bytes = checked_block_bytes(block)
199            .filter(|b| *b <= MAX_BLOCK_BYTES)
200            .ok_or_else(|| {
201                MapError::Unsupported(format!(
202                    "{path}: block {}x{} float32 exceeds {MAX_BLOCK_BYTES} bytes or overflows",
203                    block.0, block.1
204                ))
205            })?;
206        let nodata = r.nodata();
207        Ok(Self {
208            path: path.to_string(),
209            gt,
210            width,
211            height,
212            block,
213            nodata,
214            epsg: expected_epsg,
215            crs_identity,
216            crs_wkt_fnv1a64: wkt_fnv(&crs),
217            query,
218            policy,
219            block_bytes,
220            readers: Mutex::new(vec![r]),
221            readers_opened: AtomicU64::new(1),
222            decode_nanos: AtomicU64::new(0),
223            cache: TileCache::new(cache_budget),
224        })
225    }
226
227    pub fn blocks(&self) -> (usize, usize) {
228        (
229            self.width.div_ceil(self.block.0),
230            self.height.div_ceil(self.block.1),
231        )
232    }
233
234    /// Pixel window of block (bx, by), clipped to the image. Always within one TIFF chunk.
235    pub fn block_window(&self, bx: usize, by: usize) -> Window {
236        let (bw, bh) = self.block;
237        let x = bx * bw;
238        let y = by * bh;
239        let w = bw.min(self.width - x);
240        let h = bh.min(self.height - y);
241        debug_assert!(x / bw == (x + w - 1) / bw && y / bh == (y + h - 1) / bh);
242        Window::new(x, y, w, h)
243    }
244
245    pub fn max_block_bytes(&self) -> u64 {
246        self.block_bytes
247    }
248
249    /// Estimated transient bytes of one concurrent block decode, outside the cache budget:
250    /// terrio window buffer, tiff decoded chunk, compressed chunk (assumed <= decoded size) and
251    /// the `Arc<[f32]>` copy (4 x block bytes), plus terrio's per-read copies of the full chunk
252    /// offset and byte-count tables (16 bytes per chunk). From source reading; not measured.
253    pub fn decoder_transient_bytes_per_load(&self) -> u64 {
254        let (bx, by) = self.blocks();
255        let tables = (bx as u64).saturating_mul(by as u64).saturating_mul(16);
256        self.block_bytes.saturating_mul(4).saturating_add(tables)
257    }
258
259    /// Shared decoded block through the cache.
260    pub fn get_block(&self, bx: usize, by: usize) -> Result<Arc<[f32]>, MapError> {
261        let w = self.block_window(bx, by);
262        // At most block_bytes, which was checked at open.
263        let bytes = (w.width * w.height * 4) as u64;
264        self.cache
265            .get_or_load((bx as u32, by as u32), bytes, || self.decode(w))
266    }
267
268    fn decode(&self, w: Window) -> Result<Vec<f32>, MapError> {
269        let reader = self.readers.lock().unwrap_or_else(|p| p.into_inner()).pop();
270        let mut reader = match reader {
271            Some(r) => r,
272            None => {
273                self.readers_opened.fetch_add(1, Ordering::Relaxed);
274                open_lazy(&self.path)?
275            }
276        };
277        let start = Instant::now();
278        let options = WindowOptions {
279            bands: vec![1],
280            overview: None,
281        };
282        let ds = reader.read_window(w, &options).map_err(|e| {
283            MapError::Terrain(format!("{}: block read {:?}: {e}", self.path, (w.x, w.y)))
284        })?;
285        self.decode_nanos
286            .fetch_add(start.elapsed().as_nanos() as u64, Ordering::Relaxed);
287        // Return the reader only after a successful read.
288        self.readers
289            .lock()
290            .unwrap_or_else(|p| p.into_inner())
291            .push(reader);
292        let terrio::raster::RasterData::F32(values) = ds.raster.data else {
293            return Err(MapError::Unsupported(
294                "terrio returned non-float32 block".into(),
295            ));
296        };
297        if values.len() != w.width * w.height {
298            return Err(MapError::Terrain(format!(
299                "block has {} values, expected {}",
300                values.len(),
301                w.width * w.height
302            )));
303        }
304        Ok(values)
305    }
306
307    /// Explicit sampling formula: col = floor((xs - gt0)/gt1), row = floor((ys - gt3)/gt5).
308    #[inline]
309    pub fn hit(&self, xs: f64, ys: f64) -> Hit {
310        let fc = (xs - self.gt[0]) / self.gt[1];
311        let fr = (ys - self.gt[3]) / self.gt[5];
312        if !(fc.is_finite() && fr.is_finite()) {
313            return Hit::Outside;
314        }
315        let (c, r) = (fc.floor(), fr.floor());
316        if c < 0.0 || r < 0.0 || c >= self.width as f64 || r >= self.height as f64 {
317            return Hit::Outside;
318        }
319        let tol = 1e-9;
320        let near = (fc - c) < tol || (c + 1.0 - fc) < tol || (fr - r) < tol || (r + 1.0 - fr) < tol;
321        Hit::Pixel {
322            col: c as usize,
323            row: r as usize,
324            near_boundary: near,
325        }
326    }
327
328    /// Floor index along one axis, unclamped (used for footprint planning).
329    pub fn floor_index(&self, xs: f64, ys: f64) -> Option<(f64, f64)> {
330        let fc = ((xs - self.gt[0]) / self.gt[1]).floor();
331        let fr = ((ys - self.gt[3]) / self.gt[5]).floor();
332        (fc.is_finite() && fr.is_finite()).then_some((fc, fr))
333    }
334
335    /// Applies the explicit nodata policy to a raw terrain value. NaN is dry. An infinite value
336    /// that is not a dry sentinel is returned as is; the session rejects it with a typed terrain
337    /// error (matching the caller_buffer rejection of infinite terrain).
338    #[inline]
339    pub fn classify(&self, v: f32) -> Option<f64> {
340        if v.is_nan() {
341            return None;
342        }
343        let sentinel = self.nodata.is_some_and(|nd| v == nd as f32);
344        match (sentinel, self.policy) {
345            (true, NodataPolicy::Dry) => None,
346            _ => Some(v as f64),
347        }
348    }
349
350    pub fn is_sentinel(&self, v: f32) -> bool {
351        self.nodata
352            .is_some_and(|nd| v == nd as f32 || (nd.is_nan() && v.is_nan()))
353    }
354
355    pub fn stats(&self) -> TerrainStats {
356        let cache = self.cache.stats();
357        let per = self.decoder_transient_bytes_per_load();
358        TerrainStats {
359            decoder_transient_peak_estimate_bytes: per.saturating_mul(cache.peak_concurrent_loads),
360            decoder_transient_bytes_per_load_estimate: per,
361            cache,
362            readers_opened: self.readers_opened.load(Ordering::Relaxed),
363            decode_seconds: self.decode_nanos.load(Ordering::Relaxed) as f64 * 1e-9,
364        }
365    }
366}
367
368/// Largest single decoded block accepted (1 GiB of float32).
369pub const MAX_BLOCK_BYTES: u64 = 1 << 30;
370
371/// `bw * bh * 4` with overflow detection.
372pub fn checked_block_bytes(block: (usize, usize)) -> Option<u64> {
373    (block.0 as u64)
374        .checked_mul(block.1 as u64)
375        .and_then(|v| v.checked_mul(4))
376}
377
378fn check_wkt_pin(path: &str, crs: &terrio::crs::SpatialRef, fnv1a64: &str) -> Result<(), MapError> {
379    let got =
380        wkt_fnv(crs).ok_or_else(|| MapError::Unsupported(format!("{path}: terrain has no CRS")))?;
381    if got.eq_ignore_ascii_case(fnv1a64) {
382        Ok(())
383    } else {
384        Err(MapError::Unsupported(format!(
385            "{path}: terrain WKT hash {got} differs from caller-verified {fnv1a64}"
386        )))
387    }
388}
389
390/// Every declared nodata is accepted, including NaN and +/-inf (GDAL-like: a pixel is the
391/// sentinel iff `v == nodata as f32`). Whether that cast is exact is reported, not guessed.
392pub fn nodata_f32_exact(nodata: Option<f64>) -> Option<bool> {
393    nodata.map(|nd| nd.is_nan() || (nd as f32) as f64 == nd)
394}
395
396#[cfg(test)]
397mod tests {
398    use super::*;
399
400    #[test]
401    fn absurd_block_sizes_are_detected_not_panics() {
402        assert_eq!(checked_block_bytes((256, 256)), Some(262_144));
403        assert_eq!(checked_block_bytes((usize::MAX, 2)), None);
404        assert_eq!(
405            checked_block_bytes((u32::MAX as usize, u32::MAX as usize)),
406            None
407        );
408        assert!(checked_block_bytes((65_536, 65_536)).unwrap() > MAX_BLOCK_BYTES);
409    }
410
411    #[test]
412    fn nodata_cast_exactness_is_reported() {
413        assert_eq!(nodata_f32_exact(None), None);
414        assert_eq!(nodata_f32_exact(Some(-9999.0)), Some(true));
415        assert_eq!(nodata_f32_exact(Some(f64::NAN)), Some(true));
416        assert_eq!(nodata_f32_exact(Some(f64::NEG_INFINITY)), Some(true));
417        assert_eq!(nodata_f32_exact(Some(-3.4028234663852886e38)), Some(true));
418        assert_eq!(nodata_f32_exact(Some(-3.40282e38)), Some(false));
419    }
420}
ras_commander_fast_map/terrain.rs

ras_commander_fast_map/
terrain.rs
