navcore_route_find/

land_cache.rs

//! A precomputed, ahead-of-time cache of a chart's land geometry, and a
//! raster derived from it, so a route query does not have to fetch,
//! reconcile and simplify `LNDARE` polygons from the chart before
//! testing anything against them.
//!
//! Land is the one hazard class worth this: [`crate::visibility::dijkstra`]
//! tests a candidate edge against it up to [`crate::visibility::MAX_VERTICES`]
//! squared times, and those [`crate::hazards::Hazards::intersects_band`]
//! calls -- each a [`geo::Intersects`] test against a simplified but
//! still real coastline -- are the dominant cost of a hard case's
//! wall-clock time. [`crate::hazards::Hazards::with_simplified_land`]'s
//! per-lattice-node test pays the same real-geometry cost on every route
//! this crate finds, not only the visibility-graph fallback's hard cases.
//!
//! A raster answers the same question -- does this query footprint touch
//! land -- in time proportional to the footprint's size, not to how many
//! points the charted coastline is sampled at. Restricted areas and
//! depth areas are not cached here: both are comparatively few and
//! cheap, and depth is already dropped from the hot screening path (see
//! [`crate::hazards::Hazards::coarsened_for_screening`]'s own doc).
//!
//! Written once, ahead of time, by [`bake`] (driven by `enc-check bake`);
//! read lazily and cached for the life of the process by [`load`]. A
//! missing or stale cache is never a correctness risk, only a speed one:
//! every caller in [`crate::hazards`] falls back to the live path this
//! module exists to avoid paying for on every query.
//!
//! A cache built from over-simplified geometry would itself be a
//! correctness risk: a small islet simplified down to nothing by
//! [`crate::hazards::SCREENING_SIMPLIFY_DEG`]'s tolerance would never
//! reach the raster, and a query small enough for [`LandCache::covers`]
//! to trust the cache would read open water where `check_leg`, reading
//! the chart directly, still finds land. [`bake`] rasterizes from the
//! unsimplified parts instead, kept apart from the
//! [`simplify_land`]-simplified ones [`LandCache::land_geometry`] still
//! stores for the coarse lattice/visibility use that tolerance is for --
//! see that function's own doc for the full reasoning.

use std::fs;
use std::path::{Path, PathBuf};

use geo::{BoundingRect, Coord, Geometry, MultiPolygon, Polygon, Rect};
use rusqlite::Connection;
use wkt::{ToWkt, TryFromWkt};

use enc_store::{ChartStore, Feature, StoreError};

use crate::hazards::{fetch_preferring_detail, simplify_land};

/// The S-57 class this module caches. Mirrors `hazards::ALWAYS_UNSAFE`'s
/// one entry -- not imported from there, the same "answers a different
/// question, stayed simple enough to duplicate" reasoning this crate
/// already applies to its own small `polygon_parts` helpers.
const LAND_CLASS: &str = "LNDARE";

/// The S-57 class marking a chart's own real, surveyed coverage -- the
/// authoritative source [`detect_coverage_area`] uses for "what area is
/// this chart actually for", so nobody baking a cache has to name a
/// region by hand and keep it in sync with whatever a deployment's own
/// chart-ingest pipeline adds over time.
const COVERAGE_CLASS: &str = "M_COVR";

/// How far apart two [`COVERAGE_CLASS`] features' centres may be, in
/// degrees, and still count as the same region in
/// [`detect_coverage_area`]'s clustering: generous enough to bridge an
/// ordinary gap between adjacent chart cells or between separately
/// surveyed stretches of the same coastline, tight enough that a
/// genuinely different part of the world stays its own cluster. Some
/// charts carry `M_COVR` rows far outside their real coverage area,
/// background rows a broader chart product carries alongside the region
/// it was actually produced for.
const CLUSTER_DISTANCE_DEG: f64 = 3.0;

/// The raster's own cell size, in degrees -- finer than
/// [`crate::hazards::SCREENING_SIMPLIFY_DEG`]'s 0.005 degrees, so no
/// coastline detail the simplified geometry still carries is lost to a
/// coarser grid on top of it. A starting point, not a settled constant --
/// tune against a real chart's own bake time and cache size the same way
/// [`crate::visibility::MAX_VERTICES`] was tuned, by measurement rather
/// than guesswork, if either turns out to matter in practice.
const CELL_DEG: f64 = 0.0025;

/// How far beyond the fetched land geometry's own bounding box the raster
/// extends, in degrees -- headroom for a query footprint that grazes the
/// chart's own edge without every such query falling outside the raster
/// and silently reading as clear water for the wrong reason.
const RASTER_PADDING_DEG: f64 = 0.1;

/// The cache format's own version, bumped whenever [`bake`]'s output or
/// [`load`]'s expectations of it change shape -- a version mismatch is
/// read the same as a missing file, not an error.
const FORMAT_VERSION: i64 = 1;

/// Everything that can go wrong baking a chart's land cache.
#[derive(Debug)]
pub enum LandCacheError {
    /// Reading the source chart failed.
    Store(StoreError),
    /// Writing or reading the cache file itself failed.
    Database(rusqlite::Error),
    /// The source chart's own file metadata could not be read.
    Io(std::io::Error),
    /// The chart carries no `COVERAGE_CLASS` features at all, so
    /// `detect_coverage_area` has nothing to detect a region from.
    NoCoverage,
}

impl std::fmt::Display for LandCacheError {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        match self {
            Self::Store(error) => write!(f, "chart: {error}"),
            Self::Database(error) => write!(f, "land cache: {error}"),
            Self::Io(error) => write!(f, "chart file: {error}"),
            Self::NoCoverage => write!(f, "chart carries no M_COVR coverage features to detect a bake region from"),
        }
    }
}

impl std::error::Error for LandCacheError {}

impl From<StoreError> for LandCacheError {
    fn from(error: StoreError) -> Self {
        Self::Store(error)
    }
}

impl From<rusqlite::Error> for LandCacheError {
    fn from(error: rusqlite::Error) -> Self {
        Self::Database(error)
    }
}

impl From<std::io::Error> for LandCacheError {
    fn from(error: std::io::Error) -> Self {
        Self::Io(error)
    }
}

/// The sidecar path [`bake`] writes to and [`load`] reads from for a chart
/// opened from `chart_path` -- appended to the chart's own name rather
/// than replacing its extension, so it is unambiguous which chart file a
/// cache belongs to even if a directory holds charts that differ only in
/// extension.
fn cache_path(chart_path: &Path) -> PathBuf {
    let mut name = chart_path.as_os_str().to_owned();
    name.push(".hazcache");
    PathBuf::from(name)
}

/// A fingerprint of a chart file cheap enough to check on every
/// [`load`] -- not a content hash, just enough to catch "this chart was
/// replaced since the cache was baked" without reading the whole file
/// again to confirm it.
fn fingerprint(chart_path: &Path) -> Result<(u64, i64), std::io::Error> {
    let metadata = fs::metadata(chart_path)?;
    let modified = metadata
        .modified()?
        .duration_since(std::time::UNIX_EPOCH)
        .map(|duration| duration.as_secs() as i64)
        .unwrap_or(0);
    Ok((metadata.len(), modified))
}

/// A precomputed land raster and the simplified geometry it was rasterized
/// from, loaded once and kept for the life of the process -- see
/// [`crate::hazards`]'s process-wide cache of this type.
pub(crate) struct LandCache {
    origin: Coord<f64>,
    cell_deg: f64,
    cols: usize,
    rows: usize,
    raster: Vec<bool>,
    land_geometry: Geometry<f64>,
}

impl LandCache {
    /// The simplified land geometry this cache was built from -- what
    /// [`crate::hazards::Hazards::land`] hands back when a cache is
    /// loaded, identical in shape to what
    /// [`crate::hazards::Hazards::with_simplified_land`] computes live,
    /// just precomputed. Always a [`Geometry::MultiPolygon`], to match
    /// the type every other land entry in that iterator already is.
    pub(crate) fn land_geometry(&self) -> &Geometry<f64> {
        &self.land_geometry
    }

    /// Whether `area` falls entirely within this cache's own baked
    /// extent.
    ///
    /// A chart is one file, but nothing says every query against it
    /// stays inside whatever region it was baked for: a chart can carry
    /// real, detailed coastline far from the area actually baked.
    /// Trusting the cache there regardless would not just be slow:
    /// [`Self::cell_of`] reads "outside the raster" as "no land here",
    /// the same conclusion an empty live fetch would reach, so real
    /// charted coastline outside the baked extent would go completely
    /// unreported rather than merely uncached. [`Hazards::fetch`] checks
    /// this before trusting a loaded cache for a given query; a query
    /// that fails it falls back to the live fetch, the same fallback a
    /// missing or stale cache takes.
    pub(crate) fn covers(&self, area: Rect<f64>) -> bool {
        #[allow(clippy::cast_precision_loss)]
        let (max_x, max_y) = (self.origin.x + self.cell_deg * self.cols as f64, self.origin.y + self.cell_deg * self.rows as f64);
        area.min().x >= self.origin.x
            && area.min().y >= self.origin.y
            && area.max().x <= max_x
            && area.max().y <= max_y
    }

    /// The raster cell `coord` falls in, or `None` outside the raster's
    /// own coverage -- which [`Self::blocked`] reads as "no land the
    /// chart this cache was baked from ever reported here", the same
    /// conclusion an empty fetch would already reach live.
    fn cell_of(&self, coord: Coord<f64>) -> Option<(usize, usize)> {
        let col = ((coord.x - self.origin.x) / self.cell_deg).floor();
        let row = ((coord.y - self.origin.y) / self.cell_deg).floor();
        if col < 0.0 || row < 0.0 {
            return None;
        }
        let (col, row) = (col as usize, row as usize);
        (col < self.cols && row < self.rows).then_some((row, col))
    }

    fn cell_blocked(&self, row: usize, col: usize) -> bool {
        self.raster[row * self.cols + col]
    }

    /// Whether `polygon` -- a query footprint standing in for the water a
    /// candidate step or leg needs clear, the same role
    /// [`crate::hazards::Hazards::intersects_band`]'s own `band` plays --
    /// touches a raster cell marked as land.
    ///
    /// Samples the polygon's exterior ring at [`CELL_DEG`]-spaced points
    /// along every edge rather than rasterizing its full interior: exact
    /// for a shape no wider than a handful of cells (a vessel's corridor
    /// is rarely more than a couple of tenths of a nautical mile either
    /// side of the line it is testing), and for a wider one this can in
    /// principle miss a small obstacle entirely inside the footprint.
    /// That gap is no wider than [`crate::hazards::SCREENING_SIMPLIFY_DEG`]'s
    /// own simplification tolerance for the same coarse-screening
    /// purpose, and this cache is never the last word regardless:
    /// [`crate::smooth::straighten`]'s `check_leg` calls read the
    /// unsimplified chart directly and catch whatever this test misses.
    #[must_use]
    pub(crate) fn blocked(&self, polygon: &Polygon<f64>) -> bool {
        let exterior = polygon.exterior();
        let vertices = exterior.0.as_slice();
        if vertices.is_empty() {
            return false;
        }
        if self.point_blocked(vertices[0]) {
            return true;
        }
        vertices.windows(2).any(|edge| self.segment_blocked(edge[0], edge[1]))
    }

    fn point_blocked(&self, point: Coord<f64>) -> bool {
        self.cell_of(point).is_some_and(|(row, col)| self.cell_blocked(row, col))
    }

    fn segment_blocked(&self, start: Coord<f64>, end: Coord<f64>) -> bool {
        let length = ((end.x - start.x).powi(2) + (end.y - start.y).powi(2)).sqrt();
        let steps = (length / self.cell_deg).ceil().max(1.0) as usize;
        (1..=steps).any(|step| {
            let t = f64::from(u32::try_from(step).unwrap_or(u32::MAX)) / f64::from(u32::try_from(steps).unwrap_or(1));
            let point = Coord { x: start.x + (end.x - start.x) * t, y: start.y + (end.y - start.y) * t };
            self.point_blocked(point)
        })
    }
}

/// The individual polygon parts inside `geometry` -- see
/// `hazards::polygon_parts`'s own doc for why this crate duplicates
/// rather than shares this helper across modules.
fn polygon_parts(geometry: &Geometry<f64>) -> Vec<Polygon<f64>> {
    match geometry {
        Geometry::Polygon(polygon) => vec![polygon.clone()],
        Geometry::MultiPolygon(multi) => multi.iter().cloned().collect(),
        _ => Vec::new(),
    }
}

/// Marks every raster cell `part` covers as land, in `bits` (row-major,
/// `cols` wide).
///
/// A proper scanline fill, not "test every cell in the bounding box
/// against the polygon": a single large, many-vertex coastline's own
/// bounding box alone can hold hundreds of thousands of cells, and
/// testing each one against the polygon directly would pay a full
/// polygon-intersects test whose own cost scales with the ring's vertex
/// count -- the same "cost scales with coastline complexity" problem
/// this whole cache exists to get away from, just moved from query time
/// to bake time. Walking one horizontal line per raster row
/// instead -- intersecting it against every ring edge once, even-odd
/// filling the spans between crossings -- costs rows times edges, not
/// cells times edges; `cols` drops out of the cost entirely.
fn rasterize_part(part: &Polygon<f64>, origin: Coord<f64>, cell_deg: f64, cols: usize, rows: usize, bits: &mut [bool]) {
    let Some(bbox) = part.bounding_rect() else { return };
    #[allow(clippy::cast_precision_loss)]
    let (raster_max_x, raster_max_y) = (origin.x + cell_deg * cols as f64, origin.y + cell_deg * rows as f64);
    // A part whose own bounding box does not touch the raster at all --
    // possible even with an area-scoped fetch, since a feature is
    // returned whenever its box merely overlaps the query area, not only
    // when it sits fully inside it. Left unguarded, the clamped indices
    // below would fold such a part onto the raster's own edge row/column
    // instead of skipping it, marking cells as land that are nowhere
    // near it.
    if bbox.max().x < origin.x || bbox.min().x > raster_max_x || bbox.max().y < origin.y || bbox.min().y > raster_max_y {
        return;
    }
    let min_row = (((bbox.min().y - origin.y) / cell_deg).floor().max(0.0) as usize).min(rows.saturating_sub(1));
    let max_row = (((bbox.max().y - origin.y) / cell_deg).ceil() as usize).min(rows.saturating_sub(1));

    let rings: Vec<&geo::LineString<f64>> = std::iter::once(part.exterior()).chain(part.interiors()).collect();

    for row in min_row..=max_row {
        #[allow(clippy::cast_precision_loss)]
        let y = origin.y + (row as f64 + 0.5) * cell_deg;
        let mut crossings: Vec<f64> = rings
            .iter()
            .flat_map(|ring| ring.lines())
            .filter_map(|edge| {
                let (y0, y1) = (edge.start.y, edge.end.y);
                // The standard half-open scanline test: an edge crosses
                // the line at `y` once, not zero or two times, even when
                // `y` passes exactly through a shared vertex between two
                // consecutive edges.
                let crosses = (y0 <= y) != (y1 <= y);
                crosses.then(|| edge.start.x + (y - y0) / (y1 - y0) * (edge.end.x - edge.start.x))
            })
            .collect();
        crossings.sort_by(f64::total_cmp);

        for span in crossings.chunks_exact(2) {
            let min_col = (((span[0] - origin.x) / cell_deg).floor().max(0.0) as usize).min(cols.saturating_sub(1));
            let max_col = (((span[1] - origin.x) / cell_deg).ceil() as usize).min(cols.saturating_sub(1));
            for col in min_col..=max_col {
                bits[row * cols + col] = true;
            }
        }
    }
}

/// The real geographic extent `chart` provides coverage for, detected
/// from its own [`COVERAGE_CLASS`] features rather than named by a
/// caller.
///
/// Not simply the bounding box of every [`COVERAGE_CLASS`] feature in
/// the file: as [`CLUSTER_DISTANCE_DEG`]'s own doc explains, a chart can
/// carry a handful of such rows unrelated to the rest of the file, and a
/// plain bounding box is not robust to even one such outlier. Clustering
/// by proximity first, then keeping only the largest cluster, is robust:
/// a chart's real coverage is one geographically contiguous region, and
/// a background row elsewhere on the planet cannot join that cluster as
/// long as it stays farther than [`CLUSTER_DISTANCE_DEG`] from it. This
/// also makes baking automatic rather than something a caller names by
/// hand: a deployment that adds another country's charts to its ingest
/// pipeline grows this cluster along with it, nothing here to update.
///
/// # Errors
///
/// If the chart cannot be read, or carries no [`COVERAGE_CLASS`]
/// features at all.
fn detect_coverage_area(chart: &ChartStore) -> Result<Rect<f64>, LandCacheError> {
    let world = Rect::new(Coord { x: -180.0, y: -90.0 }, Coord { x: 180.0, y: 90.0 });
    let features = chart.features_in(COVERAGE_CLASS, world)?;
    let boxes: Vec<Rect<f64>> = features.iter().filter_map(|feature| feature.geometry.bounding_rect()).collect();
    largest_cluster_bbox(&boxes).ok_or(LandCacheError::NoCoverage)
}

/// The combined bounding box of the largest cluster among `boxes`,
/// where two boxes join the same cluster whenever their own centres sit
/// within [`CLUSTER_DISTANCE_DEG`] of each other -- transitively, so an
/// elongated coastline's own boxes all end up in one cluster even
/// though its own two ends may be much farther apart than the threshold
/// itself. `None` only for an empty `boxes`.
fn largest_cluster_bbox(boxes: &[Rect<f64>]) -> Option<Rect<f64>> {
    if boxes.is_empty() {
        return None;
    }
    let centres: Vec<Coord<f64>> = boxes.iter().map(|bbox| bbox.center()).collect();

    let mut parent: Vec<usize> = (0..centres.len()).collect();
    for i in 0..centres.len() {
        for j in (i + 1)..centres.len() {
            let (a, b) = (centres[i], centres[j]);
            let close = (a.x - b.x).abs() <= CLUSTER_DISTANCE_DEG && (a.y - b.y).abs() <= CLUSTER_DISTANCE_DEG;
            if close {
                union(&mut parent, i, j);
            }
        }
    }

    let mut cluster_size: std::collections::HashMap<usize, usize> = std::collections::HashMap::new();
    for i in 0..centres.len() {
        *cluster_size.entry(find(&mut parent, i)).or_insert(0) += 1;
    }
    let largest_root = cluster_size
        .into_iter()
        .max_by_key(|&(_, count)| count)
        .map(|(root, _)| root)
        .expect("boxes is non-empty, so at least one cluster exists");

    let mut area: Option<Rect<f64>> = None;
    for (index, &bbox) in boxes.iter().enumerate() {
        if find(&mut parent, index) != largest_root {
            continue;
        }
        area = Some(area.map_or(bbox, |existing| union_rect(existing, bbox)));
    }
    area
}

/// The union-find root of `x`, with path compression.
fn find(parent: &mut [usize], x: usize) -> usize {
    if parent[x] != x {
        parent[x] = find(parent, parent[x]);
    }
    parent[x]
}

/// Joins the clusters `a` and `b` belong to into one.
fn union(parent: &mut [usize], a: usize, b: usize) {
    let (root_a, root_b) = (find(parent, a), find(parent, b));
    if root_a != root_b {
        parent[root_a] = root_b;
    }
}

/// The smallest rectangle containing both `a` and `b`.
fn union_rect(a: Rect<f64>, b: Rect<f64>) -> Rect<f64> {
    Rect::new(
        Coord { x: a.min().x.min(b.min().x), y: a.min().y.min(b.min().y) },
        Coord { x: a.max().x.max(b.max().x), y: a.max().y.max(b.max().y) },
    )
}

/// How large each tile [`fetch_land_tiled`] reconciles `LNDARE`
/// overview coverage against is, in degrees.
///
/// [`detect_coverage_area`] can find a region far larger than any real
/// route's own local query, up to a whole sea basin.
/// `fetch_preferring_detail` only counts an overview feature as covered
/// by detail within the area it was asked about -- correct for a live
/// query at a route's local scale (see that function's own doc for why
/// dropping overview coverage outright is wrong), but not once that area
/// spans a whole sea: detail on one side then does nothing to trim an
/// overview feature's uncovered remainder on the other, and
/// [`crate::visibility::MAX_OBSTACLE_SPAN_DEG`]'s own doc names what a
/// hull that size does to a route search that finds it. Reconciling one
/// small tile at a time, rather than the whole detected region in one
/// call, keeps every reconciliation as local as a real route's query
/// area already is.
const BAKE_TILE_DEG: f64 = 2.0;

/// Every `LNDARE` feature [`fetch_preferring_detail`] returns across
/// `area`, tiled at [`BAKE_TILE_DEG`] rather than fetched in one call --
/// see that constant's own doc for why. A feature straddling a tile
/// boundary is fetched, reconciled and simplified once per tile it
/// overlaps rather than once overall; harmless here, since every
/// caller only ever folds the result into a simplified land shape or a
/// raster, and rasterizing (or re-simplifying) the same small overlap
/// twice changes nothing about either.
///
/// # Errors
///
/// If the chart cannot be read.
fn fetch_land_tiled(chart: &ChartStore, area: Rect<f64>) -> Result<Vec<Feature>, LandCacheError> {
    let mut features = Vec::new();
    let mut min_y = area.min().y;
    while min_y < area.max().y {
        let max_y = (min_y + BAKE_TILE_DEG).min(area.max().y);
        let mut min_x = area.min().x;
        while min_x < area.max().x {
            let max_x = (min_x + BAKE_TILE_DEG).min(area.max().x);
            let tile = Rect::new(Coord { x: min_x, y: min_y }, Coord { x: max_x, y: max_y });
            features.extend(fetch_preferring_detail(chart, tile, LAND_CLASS)?);
            min_x += BAKE_TILE_DEG;
        }
        min_y += BAKE_TILE_DEG;
    }
    Ok(features)
}

/// Bakes a land cache for `chart`, covering whatever area
/// `detect_coverage_area` finds plus `RASTER_PADDING_DEG`, and
/// writes it to `<chart.path()>.hazcache` (or `out`, if given) -- run
/// once, ahead of time, by `enc-check bake`, never at route-query time.
///
/// The region is detected from the chart itself, not named by a caller:
/// see `detect_coverage_area`'s own doc for why a caller-supplied area
/// would either have to be kept in sync by hand with whatever a
/// deployment's own chart-ingest pipeline adds to the file over time, or
/// risk baking a region that quietly stops matching what the chart
/// actually covers.
///
/// # Errors
///
/// If the chart cannot be read, carries no coverage to detect a region
/// from, its own file metadata cannot be read, or the cache file cannot
/// be written.
pub fn bake(chart: &ChartStore, out: Option<&Path>) -> Result<PathBuf, LandCacheError> {
    let target = out.map_or_else(|| cache_path(chart.path()), Path::to_path_buf);
    let (source_size, source_mtime) = fingerprint(chart.path())?;

    let area = detect_coverage_area(chart)?;
    let features = fetch_land_tiled(chart, area)?;

    // The raster is rasterized from the *un*simplified parts, deliberately
    // -- it stands in for `Hazards::blocked`'s own exact-position test
    // (see that method's own doc: "the one place in the search that has
    // to answer for the exact position rather than the lattice's
    // approximation of it"), and `simplify_land`'s tolerance
    // (`SCREENING_SIMPLIFY_DEG`, roughly 0.3 NM) is coarse enough to
    // collapse a small real islet to nothing. Found, not assumed: a real
    // one in the Marano-Grado lagoon vanishes from a simplified raster
    // this way, and `find_gateway` accepts a point sitting on it as
    // clear water -- confirmed against the chart with `check_leg`,
    // which reads geometry this method never touches and catches the
    // same islet without trouble. `land_geometry`
    // stays simplified below -- it feeds `with_simplified_land`'s own
    // coarse lattice-node test and `visibility`'s screening graph, both
    // of which want exactly this tolerance, not the raster's.
    let raw_parts: Vec<Polygon<f64>> = features.iter().flat_map(|feature| polygon_parts(&feature.geometry)).collect();

    let simplified: Vec<Geometry<f64>> = features.iter().map(|feature| simplify_land(&feature.geometry)).collect();
    let simplified_parts: Vec<Polygon<f64>> = simplified.iter().flat_map(polygon_parts).collect();
    let land_geometry = MultiPolygon::new(simplified_parts);

    let (origin, cols, rows) = raster_extent(area);
    let mut bits = vec![false; cols * rows];
    for part in &raw_parts {
        rasterize_part(part, origin, CELL_DEG, cols, rows, &mut bits);
    }

    if let Some(parent) = target.parent() {
        if !parent.as_os_str().is_empty() {
            fs::create_dir_all(parent)?;
        }
    }
    let _ = fs::remove_file(&target);
    let connection = Connection::open(&target)?;
    connection.execute_batch(
        "CREATE TABLE meta (
            format_version INTEGER NOT NULL,
            source_size INTEGER NOT NULL,
            source_mtime INTEGER NOT NULL,
            origin_lon REAL NOT NULL,
            origin_lat REAL NOT NULL,
            cell_deg REAL NOT NULL,
            cols INTEGER NOT NULL,
            rows INTEGER NOT NULL
        );
        CREATE TABLE raster (bits BLOB NOT NULL);
        CREATE TABLE land_geometry (wkt TEXT NOT NULL);",
    )?;
    #[allow(clippy::cast_possible_wrap)]
    connection.execute(
        "INSERT INTO meta (format_version, source_size, source_mtime, origin_lon, origin_lat, cell_deg, cols, rows) \
         VALUES (?1, ?2, ?3, ?4, ?5, ?6, ?7, ?8)",
        rusqlite::params![
            FORMAT_VERSION,
            source_size as i64,
            source_mtime,
            origin.x,
            origin.y,
            CELL_DEG,
            cols as i64,
            rows as i64
        ],
    )?;
    connection.execute("INSERT INTO raster (bits) VALUES (?1)", rusqlite::params![pack_bits(&bits)])?;
    connection.execute(
        "INSERT INTO land_geometry (wkt) VALUES (?1)",
        rusqlite::params![land_geometry.wkt_string()],
    )?;

    Ok(target)
}

/// The raster's own origin (south-west corner) and dimensions: `area`
/// padded by [`RASTER_PADDING_DEG`] on every side, gridded at
/// [`CELL_DEG`].
fn raster_extent(area: Rect<f64>) -> (Coord<f64>, usize, usize) {
    let origin = Coord { x: area.min().x - RASTER_PADDING_DEG, y: area.min().y - RASTER_PADDING_DEG };
    let width = (area.max().x + RASTER_PADDING_DEG) - origin.x;
    let height = (area.max().y + RASTER_PADDING_DEG) - origin.y;
    let cols = ((width / CELL_DEG).ceil() as usize).max(1);
    let rows = ((height / CELL_DEG).ceil() as usize).max(1);
    (origin, cols, rows)
}

fn pack_bits(bits: &[bool]) -> Vec<u8> {
    let mut bytes = vec![0u8; bits.len().div_ceil(8)];
    for (index, &bit) in bits.iter().enumerate() {
        if bit {
            bytes[index / 8] |= 1 << (index % 8);
        }
    }
    bytes
}

fn unpack_bits(bytes: &[u8], count: usize) -> Vec<bool> {
    (0..count).map(|index| bytes[index / 8] & (1 << (index % 8)) != 0).collect()
}

/// Loads the land cache for a chart opened from `chart_path`, if one
/// exists, is readable, and still matches that chart's own file
/// metadata -- `None` for absent, stale, or corrupt alike, all read the
/// same way by every caller in [`crate::hazards`]: fall back to the live
/// path this cache exists to avoid.
pub(crate) fn load(chart_path: &Path) -> Option<LandCache> {
    try_load(chart_path).ok().flatten()
}

fn try_load(chart_path: &Path) -> Result<Option<LandCache>, LandCacheError> {
    let path = cache_path(chart_path);
    if !path.exists() {
        return Ok(None);
    }
    let (source_size, source_mtime) = fingerprint(chart_path)?;

    let connection = Connection::open_with_flags(&path, rusqlite::OpenFlags::SQLITE_OPEN_READ_ONLY)?;
    let meta = connection.query_row(
        "SELECT format_version, source_size, source_mtime, origin_lon, origin_lat, cell_deg, cols, rows FROM meta",
        [],
        |row| {
            Ok((
                row.get::<_, i64>(0)?,
                row.get::<_, i64>(1)?,
                row.get::<_, i64>(2)?,
                row.get::<_, f64>(3)?,
                row.get::<_, f64>(4)?,
                row.get::<_, f64>(5)?,
                row.get::<_, i64>(6)?,
                row.get::<_, i64>(7)?,
            ))
        },
    );
    let Ok((format_version, cached_size, cached_mtime, origin_lon, origin_lat, cell_deg, cols, rows)) = meta else {
        return Ok(None);
    };
    #[allow(clippy::cast_possible_wrap)]
    let size_matches = cached_size == source_size as i64;
    if format_version != FORMAT_VERSION || !size_matches || cached_mtime != source_mtime {
        return Ok(None);
    }

    let Ok(bits_blob) = connection.query_row("SELECT bits FROM raster", [], |row| row.get::<_, Vec<u8>>(0)) else {
        return Ok(None);
    };
    let Ok(wkt_text) = connection.query_row("SELECT wkt FROM land_geometry", [], |row| row.get::<_, String>(0)) else {
        return Ok(None);
    };
    let Ok(land_geometry) = MultiPolygon::<f64>::try_from_wkt_str(&wkt_text) else {
        return Ok(None);
    };

    let (cols, rows) = (cols as usize, rows as usize);
    let raster = unpack_bits(&bits_blob, cols * rows);

    Ok(Some(LandCache {
        origin: Coord { x: origin_lon, y: origin_lat },
        cell_deg,
        cols,
        rows,
        raster,
        land_geometry: Geometry::MultiPolygon(land_geometry),
    }))
}

#[cfg(test)]
mod tests {
    use super::*;
    use geo::LineString;

    fn square(min_lon: f64, min_lat: f64, max_lon: f64, max_lat: f64) -> Polygon<f64> {
        Polygon::new(
            LineString::new(vec![
                Coord { x: min_lon, y: min_lat },
                Coord { x: max_lon, y: min_lat },
                Coord { x: max_lon, y: max_lat },
                Coord { x: min_lon, y: max_lat },
                Coord { x: min_lon, y: min_lat },
            ]),
            Vec::new(),
        )
    }

    fn small_cache() -> LandCache {
        let land = vec![square(13.0, 45.0, 13.1, 45.1)];
        let area = Rect::new(Coord { x: 12.5, y: 44.5 }, Coord { x: 13.5, y: 45.5 });
        let (origin, cols, rows) = raster_extent(area);
        let mut bits = vec![false; cols * rows];
        for part in &land {
            rasterize_part(part, origin, CELL_DEG, cols, rows, &mut bits);
        }
        LandCache {
            origin,
            cell_deg: CELL_DEG,
            cols,
            rows,
            raster: bits,
            land_geometry: Geometry::MultiPolygon(MultiPolygon::new(land)),
        }
    }

    fn rect(min_lon: f64, min_lat: f64, max_lon: f64, max_lat: f64) -> Rect<f64> {
        Rect::new(Coord { x: min_lon, y: min_lat }, Coord { x: max_lon, y: max_lat })
    }

    #[test]
    fn largest_cluster_bbox_ignores_a_far_away_outlier() {
        // A dozen boxes scattered across a real coastline's own scale
        // (a few degrees apart, chained so the whole run is one
        // cluster), plus one far away on its own -- the Antarctica case
        // this function exists for, in miniature.
        let mainland: Vec<Rect<f64>> = (0..12).map(|i| rect(f64::from(i) * 2.0, 45.0, f64::from(i) * 2.0 + 0.5, 45.5)).collect();
        let mut boxes = mainland.clone();
        boxes.push(rect(160.0, -75.0, 160.5, -74.5));

        let area = largest_cluster_bbox(&boxes).expect("a cluster exists");
        assert!(area.max().x < 30.0, "the outlier must not have widened the detected area: {area:?}");
        assert!(area.min().y > 0.0, "the outlier must not have widened the detected area: {area:?}");

        let expected = mainland.into_iter().reduce(union_rect).unwrap();
        assert_eq!(area, expected);
    }

    #[test]
    fn largest_cluster_bbox_of_no_boxes_is_none() {
        assert_eq!(largest_cluster_bbox(&[]), None);
    }

    #[test]
    fn covers_accepts_only_an_area_fully_inside_the_baked_extent() {
        let cache = small_cache();
        assert!(cache.covers(Rect::new(Coord { x: 12.6, y: 44.6 }, Coord { x: 13.4, y: 45.4 })));
        assert!(
            !cache.covers(Rect::new(Coord { x: 13.0, y: 37.0 }, Coord { x: 15.0, y: 38.0 })),
            "an area far outside the baked extent -- Sicily, say, against a cache baked for the \
             Adriatic -- must not be reported as covered"
        );
        assert!(
            !cache.covers(Rect::new(Coord { x: 13.0, y: 45.0 }, Coord { x: 20.0, y: 46.0 })),
            "an area only partly inside the baked extent is not covered either"
        );
    }

    #[test]
    fn a_footprint_touching_land_reads_blocked() {
        let cache = small_cache();
        assert!(cache.blocked(&square(13.04, 45.04, 13.06, 45.06)));
    }

    #[test]
    fn a_footprint_in_clear_water_reads_clear() {
        let cache = small_cache();
        assert!(!cache.blocked(&square(20.0, 50.0, 20.1, 50.1)));
    }

    #[test]
    fn a_long_thin_corridor_crossing_land_is_caught_along_its_edge() {
        // A corridor running the length of the raster, only grazing the
        // land block's own southern edge -- the case a naive bounding-box
        // test over-approximates massively, and a boundary-only test
        // could in principle miss if it only sampled interior points.
        let cache = small_cache();
        let corridor = square(12.9, 45.0, 13.2, 45.005);
        assert!(cache.blocked(&corridor));
    }

    #[test]
    fn a_long_thin_corridor_missing_land_entirely_is_clear() {
        let cache = small_cache();
        let corridor = square(12.9, 45.5, 13.2, 45.505);
        assert!(!cache.blocked(&corridor));
    }

    #[test]
    fn bake_and_load_round_trip_agree_with_the_live_geometry() {
        let dir = std::env::temp_dir().join(format!("land-cache-test-{}", std::process::id()));
        std::fs::create_dir_all(&dir).unwrap();
        let chart_path = dir.join("does-not-need-to-exist-for-this-fingerprint.gpkg");
        std::fs::write(&chart_path, b"not a real chart, just a fingerprint source").unwrap();

        let land = vec![square(13.0, 45.0, 13.1, 45.1)];
        let area = Rect::new(Coord { x: 12.5, y: 44.5 }, Coord { x: 13.5, y: 45.5 });
        let (origin, cols, rows) = raster_extent(area);
        let mut bits = vec![false; cols * rows];
        for part in &land {
            rasterize_part(part, origin, CELL_DEG, cols, rows, &mut bits);
        }
        let (source_size, source_mtime) = fingerprint(&chart_path).unwrap();
        let target = cache_path(&chart_path);
        let connection = Connection::open(&target).unwrap();
        connection
            .execute_batch(
                "CREATE TABLE meta (format_version INTEGER, source_size INTEGER, source_mtime INTEGER, \
                 origin_lon REAL, origin_lat REAL, cell_deg REAL, cols INTEGER, rows INTEGER); \
                 CREATE TABLE raster (bits BLOB); CREATE TABLE land_geometry (wkt TEXT);",
            )
            .unwrap();
        connection
            .execute(
                "INSERT INTO meta VALUES (?1, ?2, ?3, ?4, ?5, ?6, ?7, ?8)",
                rusqlite::params![
                    FORMAT_VERSION,
                    source_size as i64,
                    source_mtime,
                    origin.x,
                    origin.y,
                    CELL_DEG,
                    cols as i64,
                    rows as i64
                ],
            )
            .unwrap();
        connection.execute("INSERT INTO raster VALUES (?1)", rusqlite::params![pack_bits(&bits)]).unwrap();
        connection
            .execute(
                "INSERT INTO land_geometry VALUES (?1)",
                rusqlite::params![MultiPolygon::new(land).wkt_string()],
            )
            .unwrap();
        drop(connection);

        let loaded = load(&chart_path).expect("a freshly baked cache should load");
        assert!(loaded.blocked(&square(13.04, 45.04, 13.06, 45.06)));
        assert!(!loaded.blocked(&square(20.0, 50.0, 20.1, 50.1)));

        std::fs::write(&chart_path, b"a different chart entirely, same path").unwrap();
        assert!(load(&chart_path).is_none(), "a changed chart must invalidate its own stale cache");

        let _ = std::fs::remove_dir_all(&dir);
    }
}