use std::borrow::Cow;
use std::collections::HashMap;
use std::fs;
use std::path::{Path, PathBuf};
use crate::geoid::{egm96_undulation, GeoidError, GeoidGrid};
use crate::terrain::{
self, terrain_grid_candidates, validate_lookup_coordinates, DtedInterpolation,
DtedLookupOptions, DtedTile,
};
use crate::{Error, Result};
const STORE_MAGIC: &[u8; 8] = b"TMMAP001";
const STORE_VERSION: u16 = 1;
const STORE_ALIGNMENT: usize = 4096;
const STORE_HEADER_LEN: usize = 64;
const STORE_INDEX_RECORD_LEN: usize = 80;
const HEADER_VERSION_OFFSET: usize = 8;
const HEADER_DATUM_OFFSET: usize = 10;
const HEADER_TILE_COUNT_OFFSET: usize = 12;
const HEADER_INDEX_OFFSET_OFFSET: usize = 16;
const HEADER_DATA_OFFSET_OFFSET: usize = 24;
const HEADER_TOTAL_LEN_OFFSET: usize = 32;
const INDEX_LAT_OFFSET: usize = 0;
const INDEX_LON_OFFSET: usize = 4;
const INDEX_LON_COUNT_OFFSET: usize = 8;
const INDEX_LAT_COUNT_OFFSET: usize = 12;
const INDEX_DATA_OFFSET_OFFSET: usize = 16;
const INDEX_DATA_LEN_OFFSET: usize = 24;
const INDEX_CHECKSUM_OFFSET: usize = 32;
const INDEX_MIN_LAT_OFFSET: usize = 40;
const INDEX_MIN_LON_OFFSET: usize = 48;
const INDEX_MAX_LAT_OFFSET: usize = 56;
const INDEX_MAX_LON_OFFSET: usize = 64;
const INDEX_DATUM_OFFSET: usize = 72;
const FNV_OFFSET_BASIS: u64 = 0xcbf2_9ce4_8422_2325;
const FNV_PRIME: u64 = 0x0000_0100_0000_01b3;
const EGM96_DAC_REMEDIATION: &str =
"obtain the public NGA EGM96 15-arcminute WW15MGH.DAC file and load it with Egm96FifteenMinuteGeoid::from_ww15mgh_dac_path or Egm96FifteenMinuteGeoid::from_ww15mgh_dac_bytes";
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum VerticalDatum {
Egm96MslOrthometric,
}
impl VerticalDatum {
fn tag(self) -> u8 {
match self {
Self::Egm96MslOrthometric => 1,
}
}
fn from_tag(tag: u8) -> core::result::Result<Self, TerrainStoreError> {
match tag {
1 => Ok(Self::Egm96MslOrthometric),
other => Err(TerrainStoreError::UnsupportedDatum { tag: other }),
}
}
}
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct OrthometricHeightM {
pub value_m: f64,
}
impl OrthometricHeightM {
#[must_use]
pub const fn new(value_m: f64) -> Self {
Self { value_m }
}
#[must_use]
pub const fn metres(self) -> f64 {
self.value_m
}
pub fn to_ellipsoidal_height_deg(
self,
latitude_deg: f64,
longitude_deg: f64,
geoid: TerrainGeoidModel<'_>,
) -> core::result::Result<EllipsoidalHeightM, TerrainDatumError> {
Ok(EllipsoidalHeightM::new(
self.value_m + geoid.undulation_deg(latitude_deg, longitude_deg),
))
}
pub fn to_ellipsoidal_height_rad(
self,
latitude_rad: f64,
longitude_rad: f64,
geoid: TerrainGeoidModel<'_>,
) -> core::result::Result<EllipsoidalHeightM, TerrainDatumError> {
Ok(EllipsoidalHeightM::new(
self.value_m + geoid.undulation_rad(latitude_rad, longitude_rad),
))
}
}
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct EllipsoidalHeightM {
pub value_m: f64,
}
impl EllipsoidalHeightM {
#[must_use]
pub const fn new(value_m: f64) -> Self {
Self { value_m }
}
#[must_use]
pub const fn metres(self) -> f64 {
self.value_m
}
}
#[derive(Clone, Debug, PartialEq)]
pub struct Egm96FifteenMinuteGeoid {
grid: GeoidGrid,
}
impl Egm96FifteenMinuteGeoid {
pub fn from_ww15mgh_dac_bytes(bytes: &[u8]) -> core::result::Result<Self, TerrainDatumError> {
let grid = GeoidGrid::from_egm96_dac(bytes).map_err(TerrainDatumError::Geoid)?;
Ok(Self { grid })
}
pub fn from_ww15mgh_dac_path(
path: impl AsRef<Path>,
) -> core::result::Result<Self, TerrainDatumError> {
let path = path.as_ref();
let bytes = match fs::read(path) {
Ok(bytes) => bytes,
Err(err) if err.kind() == std::io::ErrorKind::NotFound => {
return Err(TerrainDatumError::MissingEgm96Dac {
path: path.to_path_buf(),
remediation: EGM96_DAC_REMEDIATION,
});
}
Err(err) => {
return Err(TerrainDatumError::Io {
path: path.to_path_buf(),
message: err.to_string(),
});
}
};
Self::from_ww15mgh_dac_bytes(&bytes)
}
#[must_use]
pub const fn grid(&self) -> &GeoidGrid {
&self.grid
}
}
#[derive(Clone, Copy, Debug)]
pub enum TerrainGeoidModel<'a> {
Egm96OneDegree,
Egm96FifteenMinute(&'a Egm96FifteenMinuteGeoid),
}
impl TerrainGeoidModel<'_> {
fn undulation_deg(self, latitude_deg: f64, longitude_deg: f64) -> f64 {
match self {
Self::Egm96OneDegree => {
egm96_undulation(latitude_deg.to_radians(), longitude_deg.to_radians())
}
Self::Egm96FifteenMinute(grid) => grid.grid.undulation_deg(latitude_deg, longitude_deg),
}
}
fn undulation_rad(self, latitude_rad: f64, longitude_rad: f64) -> f64 {
match self {
Self::Egm96OneDegree => egm96_undulation(latitude_rad, longitude_rad),
Self::Egm96FifteenMinute(grid) => grid.grid.undulation_rad(latitude_rad, longitude_rad),
}
}
}
#[derive(Debug, Clone, PartialEq)]
pub enum TerrainDatumError {
Terrain(Error),
Geoid(GeoidError),
Io {
path: PathBuf,
message: String,
},
MissingEgm96Dac {
path: PathBuf,
remediation: &'static str,
},
}
impl core::fmt::Display for TerrainDatumError {
fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
match self {
Self::Terrain(err) => write!(f, "terrain lookup failed: {err}"),
Self::Geoid(err) => write!(f, "geoid grid failed: {err}"),
Self::Io { path, message } => {
write!(f, "{} could not be read: {message}", path.display())
}
Self::MissingEgm96Dac { path, remediation } => {
write!(f, "{} is missing; {remediation}", path.display())
}
}
}
}
impl std::error::Error for TerrainDatumError {}
impl From<Error> for TerrainDatumError {
fn from(value: Error) -> Self {
Self::Terrain(value)
}
}
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct TerrainStoreTileIndex {
pub lat_index: i32,
pub lon_index: i32,
pub min_longitude_deg: f64,
pub min_latitude_deg: f64,
pub max_longitude_deg: f64,
pub max_latitude_deg: f64,
pub lon_count: u32,
pub lat_count: u32,
pub data_offset: u64,
pub data_len: u64,
pub checksum64: u64,
pub vertical_datum: VerticalDatum,
}
#[derive(Debug, Clone, PartialEq, Eq)]
pub enum TerrainStoreError {
Io {
path: PathBuf,
message: String,
},
Parse {
reason: String,
},
UnsupportedVersion {
version: u16,
},
UnsupportedDatum {
tag: u8,
},
DuplicateTile {
lat_index: i32,
lon_index: i32,
},
Checksum {
lat_index: i32,
lon_index: i32,
expected: u64,
found: u64,
},
}
impl core::fmt::Display for TerrainStoreError {
fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
match self {
Self::Io { path, message } => write!(f, "{} failed: {message}", path.display()),
Self::Parse { reason } => write!(f, "terrain store parse error: {reason}"),
Self::UnsupportedVersion { version } => {
write!(f, "terrain store version {version} is not supported")
}
Self::UnsupportedDatum { tag } => {
write!(f, "terrain store vertical datum tag {tag} is not supported")
}
Self::DuplicateTile {
lat_index,
lon_index,
} => write!(f, "duplicate terrain tile ({lat_index},{lon_index})"),
Self::Checksum {
lat_index,
lon_index,
expected,
found,
} => write!(
f,
"terrain tile ({lat_index},{lon_index}) checksum expected {expected:#x} but found {found:#x}"
),
}
}
}
impl std::error::Error for TerrainStoreError {}
#[derive(Clone, Debug)]
struct MmapTile {
index: TerrainStoreTileIndex,
}
impl MmapTile {
fn contains(&self, longitude_deg: f64, latitude_deg: f64) -> bool {
latitude_deg >= self.index.min_latitude_deg
&& latitude_deg <= self.index.max_latitude_deg
&& longitude_deg >= self.index.min_longitude_deg
&& longitude_deg <= self.index.max_longitude_deg
}
fn get_elevation(&self, bytes: &[u8], longitude_deg: f64, latitude_deg: f64) -> Result<i16> {
if !self.contains(longitude_deg, latitude_deg) {
return Err(Error::Parse(format!(
"point ({longitude_deg},{latitude_deg}) is outside terrain store tile ({},{})",
self.index.min_longitude_deg, self.index.min_latitude_deg
)));
}
let lat_count = self.index.lat_count as usize;
let lon_count = self.index.lon_count as usize;
let latitude_index = terrain::py_round_to_usize(
(latitude_deg - self.index.min_latitude_deg) * (lat_count - 1) as f64,
)
.map_err(Error::Parse)?;
let longitude_index = terrain::py_round_to_usize(
(longitude_deg - self.index.min_longitude_deg) * (lon_count - 1) as f64,
)
.map_err(Error::Parse)?;
if latitude_index >= lat_count || longitude_index >= lon_count {
return Err(Error::Parse(format!(
"posting index out of bounds lon={longitude_index} lat={latitude_index}"
)));
}
let sample_start =
self.index.data_offset as usize + 2 * (longitude_index * lat_count + latitude_index);
Ok(i16::from_le_bytes([
bytes[sample_start],
bytes[sample_start + 1],
]))
}
}
#[derive(Clone, Debug)]
pub struct MmapTerrain<'a> {
bytes: Cow<'a, [u8]>,
tiles: Vec<MmapTile>,
by_grid: HashMap<(i32, i32), usize>,
tile_index: Vec<TerrainStoreTileIndex>,
vertical_datum: VerticalDatum,
}
impl MmapTerrain<'static> {
pub fn from_vec(bytes: Vec<u8>) -> core::result::Result<Self, TerrainStoreError> {
Self::from_cow(Cow::Owned(bytes))
}
pub fn from_path(path: impl AsRef<Path>) -> core::result::Result<Self, TerrainStoreError> {
let path = path.as_ref();
let bytes = fs::read(path).map_err(|err| TerrainStoreError::Io {
path: path.to_path_buf(),
message: err.to_string(),
})?;
Self::from_vec(bytes)
}
}
impl<'a> MmapTerrain<'a> {
pub fn from_bytes(bytes: &'a [u8]) -> core::result::Result<Self, TerrainStoreError> {
Self::from_cow(Cow::Borrowed(bytes))
}
fn from_cow(bytes: Cow<'a, [u8]>) -> core::result::Result<Self, TerrainStoreError> {
let parsed = parse_store(bytes.as_ref())?;
Ok(Self {
bytes,
tiles: parsed.tiles,
by_grid: parsed.by_grid,
tile_index: parsed.tile_index,
vertical_datum: parsed.vertical_datum,
})
}
#[must_use]
pub fn as_bytes(&self) -> &[u8] {
self.bytes.as_ref()
}
#[must_use]
pub const fn vertical_datum(&self) -> VerticalDatum {
self.vertical_datum
}
#[must_use]
pub fn tile_index(&self) -> &[TerrainStoreTileIndex] {
&self.tile_index
}
#[must_use]
pub fn checksum64(&self) -> u64 {
terrain_store_checksum64(self.bytes.as_ref())
}
#[must_use]
pub fn to_bytes(&self) -> Vec<u8> {
let pending = self
.tiles
.iter()
.map(|tile| PendingTile {
lat_index: tile.index.lat_index,
lon_index: tile.index.lon_index,
min_latitude_deg: tile.index.min_latitude_deg,
min_longitude_deg: tile.index.min_longitude_deg,
max_latitude_deg: tile.index.max_latitude_deg,
max_longitude_deg: tile.index.max_longitude_deg,
lon_count: tile.index.lon_count,
lat_count: tile.index.lat_count,
data: self.tile_payload(tile).to_vec(),
vertical_datum: tile.index.vertical_datum,
})
.collect();
build_store(pending).expect("parsed terrain store can be reserialized")
}
pub fn height_m(&mut self, longitude_deg: f64, latitude_deg: f64) -> Result<f64> {
self.height_m_with_options(longitude_deg, latitude_deg, DtedLookupOptions::default())
}
pub fn height_m_with_options(
&mut self,
longitude_deg: f64,
latitude_deg: f64,
options: DtedLookupOptions,
) -> Result<f64> {
self.orthometric_height_m_with_options(longitude_deg, latitude_deg, options)
.map(OrthometricHeightM::metres)
}
pub fn orthometric_height_m(
&self,
longitude_deg: f64,
latitude_deg: f64,
) -> Result<OrthometricHeightM> {
self.orthometric_height_m_with_options(
longitude_deg,
latitude_deg,
DtedLookupOptions::default(),
)
}
pub fn orthometric_height_m_with_options(
&self,
longitude_deg: f64,
latitude_deg: f64,
options: DtedLookupOptions,
) -> Result<OrthometricHeightM> {
validate_lookup_coordinates(longitude_deg, latitude_deg)?;
let Some(tile_idx) = self.resolve_grid(longitude_deg, latitude_deg) else {
return Ok(OrthometricHeightM::new(0.0));
};
height_from_tile(
self.bytes.as_ref(),
&self.tiles[tile_idx],
longitude_deg,
latitude_deg,
options,
)
.map(OrthometricHeightM::new)
}
pub fn height_batch(
&mut self,
points: &[(f64, f64)],
options: DtedLookupOptions,
) -> Vec<Result<f64>> {
self.orthometric_height_batch(points, options)
.into_iter()
.map(|result| result.map(OrthometricHeightM::metres))
.collect()
}
pub fn orthometric_height_batch(
&self,
points: &[(f64, f64)],
options: DtedLookupOptions,
) -> Vec<Result<OrthometricHeightM>> {
let mut out = Vec::with_capacity(points.len());
let mut current = None;
for &(longitude_deg, latitude_deg) in points {
if let Err(err) = validate_lookup_coordinates(longitude_deg, latitude_deg) {
out.push(Err(err));
continue;
}
let primary_grid = terrain::terrain_grid(longitude_deg, latitude_deg);
if current == Some(primary_grid) {
if let Some(&tile_idx) = self.by_grid.get(&primary_grid) {
let tile = &self.tiles[tile_idx];
if tile.contains(longitude_deg, latitude_deg) {
out.push(
height_from_tile(
self.bytes.as_ref(),
tile,
longitude_deg,
latitude_deg,
options,
)
.map(OrthometricHeightM::new),
);
continue;
}
}
}
match self.resolve_grid(longitude_deg, latitude_deg) {
Some(tile_idx) => {
let tile = &self.tiles[tile_idx];
current = Some((tile.index.lat_index, tile.index.lon_index));
out.push(
height_from_tile(
self.bytes.as_ref(),
tile,
longitude_deg,
latitude_deg,
options,
)
.map(OrthometricHeightM::new),
);
}
None => {
current = None;
out.push(Ok(OrthometricHeightM::new(0.0)));
}
}
}
out
}
pub fn ellipsoidal_height_m(
&self,
longitude_deg: f64,
latitude_deg: f64,
) -> core::result::Result<EllipsoidalHeightM, TerrainDatumError> {
self.ellipsoidal_height_m_with_options(
longitude_deg,
latitude_deg,
DtedLookupOptions::default(),
)
}
pub fn ellipsoidal_height_m_with_options(
&self,
longitude_deg: f64,
latitude_deg: f64,
options: DtedLookupOptions,
) -> core::result::Result<EllipsoidalHeightM, TerrainDatumError> {
self.ellipsoidal_height_m_with_model(
longitude_deg,
latitude_deg,
options,
TerrainGeoidModel::Egm96OneDegree,
)
}
pub fn ellipsoidal_height_m_with_model(
&self,
longitude_deg: f64,
latitude_deg: f64,
options: DtedLookupOptions,
geoid: TerrainGeoidModel<'_>,
) -> core::result::Result<EllipsoidalHeightM, TerrainDatumError> {
let orthometric = self
.orthometric_height_m_with_options(longitude_deg, latitude_deg, options)
.map_err(TerrainDatumError::Terrain)?;
orthometric.to_ellipsoidal_height_deg(latitude_deg, longitude_deg, geoid)
}
fn resolve_grid(&self, longitude_deg: f64, latitude_deg: f64) -> Option<usize> {
for grid_idx in terrain_grid_candidates(longitude_deg, latitude_deg) {
if let Some(&tile_idx) = self.by_grid.get(&grid_idx) {
if self.tiles[tile_idx].contains(longitude_deg, latitude_deg) {
return Some(tile_idx);
}
}
}
None
}
fn tile_payload(&self, tile: &MmapTile) -> &[u8] {
let start = tile.index.data_offset as usize;
let end = start + tile.index.data_len as usize;
&self.bytes[start..end]
}
}
pub fn dted_tree_to_mmap_store(
root: impl AsRef<Path>,
) -> core::result::Result<Vec<u8>, TerrainStoreError> {
let root = root.as_ref();
let mut paths = Vec::new();
collect_dted_tile_paths(root, &mut paths)?;
paths.sort();
let mut pending = Vec::with_capacity(paths.len());
for path in paths {
let tile = DtedTile::from_path(&path).map_err(|reason| TerrainStoreError::Parse {
reason: format!("{}: {reason}", path.display()),
})?;
let decoded =
tile.decoded_postings_lon_major()
.map_err(|reason| TerrainStoreError::Parse {
reason: format!("{}: {reason}", path.display()),
})?;
let mut data = Vec::with_capacity(decoded.len() * 2);
for posting in decoded {
data.extend_from_slice(&posting.to_le_bytes());
}
let lat_index = tile.origin_latitude().floor() as i32;
let lon_index = tile.origin_longitude().floor() as i32;
pending.push(PendingTile {
lat_index,
lon_index,
min_latitude_deg: tile.origin_latitude(),
min_longitude_deg: tile.origin_longitude(),
max_latitude_deg: tile.origin_latitude() + 1.0,
max_longitude_deg: tile.origin_longitude() + 1.0,
lon_count: u32::try_from(tile.lon_count()).map_err(|_| TerrainStoreError::Parse {
reason: format!("{} longitude count exceeds u32", path.display()),
})?,
lat_count: u32::try_from(tile.lat_count()).map_err(|_| TerrainStoreError::Parse {
reason: format!("{} latitude count exceeds u32", path.display()),
})?,
data,
vertical_datum: VerticalDatum::Egm96MslOrthometric,
});
}
build_store(pending)
}
pub fn write_dted_tree_to_mmap_store(
root: impl AsRef<Path>,
output_path: impl AsRef<Path>,
) -> core::result::Result<(), TerrainStoreError> {
let bytes = dted_tree_to_mmap_store(root)?;
let output_path = output_path.as_ref();
fs::write(output_path, &bytes).map_err(|err| TerrainStoreError::Io {
path: output_path.to_path_buf(),
message: err.to_string(),
})
}
#[must_use]
pub fn terrain_store_checksum64(bytes: &[u8]) -> u64 {
fnv1a64(bytes)
}
#[derive(Debug)]
struct PendingTile {
lat_index: i32,
lon_index: i32,
min_latitude_deg: f64,
min_longitude_deg: f64,
max_latitude_deg: f64,
max_longitude_deg: f64,
lon_count: u32,
lat_count: u32,
data: Vec<u8>,
vertical_datum: VerticalDatum,
}
#[derive(Debug)]
struct ParsedStore {
vertical_datum: VerticalDatum,
tiles: Vec<MmapTile>,
by_grid: HashMap<(i32, i32), usize>,
tile_index: Vec<TerrainStoreTileIndex>,
}
fn height_from_tile(
bytes: &[u8],
tile: &MmapTile,
longitude_deg: f64,
latitude_deg: f64,
options: DtedLookupOptions,
) -> Result<f64> {
if options.interpolation == DtedInterpolation::NearestPosting {
return tile
.get_elevation(bytes, longitude_deg, latitude_deg)
.map(|v| v as f64);
}
let postings_per_deg_lon = tile.index.lon_count as usize - 1;
let postings_per_deg_lat = tile.index.lat_count as usize - 1;
let lon_idx = (longitude_deg - tile.index.min_longitude_deg) * postings_per_deg_lon as f64;
let lat_idx = (latitude_deg - tile.index.min_latitude_deg) * postings_per_deg_lat as f64;
let lon_lo = lon_idx.floor() as i64;
let lat_lo = lat_idx.floor() as i64;
let fx = lon_idx - lon_lo as f64;
let fy = lat_idx - lat_lo as f64;
let mut z = 0.0;
for (di, wx) in [(0i64, 1.0 - fx), (1i64, fx)] {
for (dj, wy) in [(0i64, 1.0 - fy), (1i64, fy)] {
let w = wx * wy;
if w == 0.0 {
continue;
}
let posting_lon =
tile.index.min_longitude_deg + (lon_lo + di) as f64 / postings_per_deg_lon as f64;
let posting_lat =
tile.index.min_latitude_deg + (lat_lo + dj) as f64 / postings_per_deg_lat as f64;
z += w * f64::from(tile.get_elevation(bytes, posting_lon, posting_lat)?);
}
}
Ok(z)
}
fn collect_dted_tile_paths(
root: &Path,
out: &mut Vec<PathBuf>,
) -> core::result::Result<(), TerrainStoreError> {
let entries = fs::read_dir(root).map_err(|err| TerrainStoreError::Io {
path: root.to_path_buf(),
message: err.to_string(),
})?;
for entry in entries {
let entry = entry.map_err(|err| TerrainStoreError::Io {
path: root.to_path_buf(),
message: err.to_string(),
})?;
let path = entry.path();
let file_type = entry.file_type().map_err(|err| TerrainStoreError::Io {
path: path.clone(),
message: err.to_string(),
})?;
if file_type.is_dir() {
collect_dted_tile_paths(&path, out)?;
} else if file_type.is_file() && is_dted_tile_path(&path) {
out.push(path);
}
}
Ok(())
}
fn is_dted_tile_path(path: &Path) -> bool {
path.file_name()
.and_then(|name| name.to_str())
.is_some_and(|name| name.ends_with(".dt2"))
}
fn parse_store(bytes: &[u8]) -> core::result::Result<ParsedStore, TerrainStoreError> {
if bytes.len() < STORE_HEADER_LEN {
return Err(TerrainStoreError::Parse {
reason: format!(
"store has {} bytes but needs at least {STORE_HEADER_LEN}",
bytes.len()
),
});
}
if &bytes[..STORE_MAGIC.len()] != STORE_MAGIC {
return Err(TerrainStoreError::Parse {
reason: "missing terrain store magic".to_string(),
});
}
let version = read_u16(bytes, HEADER_VERSION_OFFSET)?;
if version != STORE_VERSION {
return Err(TerrainStoreError::UnsupportedVersion { version });
}
ensure_zero(bytes, 11, 12, "header reserved byte")?;
ensure_zero(bytes, 40, STORE_HEADER_LEN, "header reserved bytes")?;
let vertical_datum = VerticalDatum::from_tag(bytes[HEADER_DATUM_OFFSET])?;
let tile_count = read_u32(bytes, HEADER_TILE_COUNT_OFFSET)? as usize;
let index_offset = read_u64(bytes, HEADER_INDEX_OFFSET_OFFSET)? as usize;
let data_offset = read_u64(bytes, HEADER_DATA_OFFSET_OFFSET)? as usize;
let total_len = read_u64(bytes, HEADER_TOTAL_LEN_OFFSET)? as usize;
if total_len != bytes.len() {
return Err(TerrainStoreError::Parse {
reason: format!(
"header total length {total_len} does not match {}",
bytes.len()
),
});
}
if index_offset != STORE_HEADER_LEN {
return Err(TerrainStoreError::Parse {
reason: format!("index offset must be {STORE_HEADER_LEN}, got {index_offset}"),
});
}
let index_len = tile_count
.checked_mul(STORE_INDEX_RECORD_LEN)
.ok_or_else(|| TerrainStoreError::Parse {
reason: "tile index length overflows usize".to_string(),
})?;
let index_end =
index_offset
.checked_add(index_len)
.ok_or_else(|| TerrainStoreError::Parse {
reason: "tile index end overflows usize".to_string(),
})?;
if index_end > bytes.len() {
return Err(TerrainStoreError::Parse {
reason: "tile index extends past store length".to_string(),
});
}
let expected_data_offset = align_up(index_end, STORE_ALIGNMENT)?;
if data_offset != expected_data_offset {
return Err(TerrainStoreError::Parse {
reason: format!("data offset must be {expected_data_offset}, got {data_offset}"),
});
}
ensure_zero(bytes, index_end, data_offset, "index padding")?;
let mut tiles = Vec::with_capacity(tile_count);
let mut tile_index = Vec::with_capacity(tile_count);
let mut by_grid = HashMap::with_capacity(tile_count);
let mut previous_id = None;
let mut expected_next = data_offset;
for idx in 0..tile_count {
let record_offset = index_offset + idx * STORE_INDEX_RECORD_LEN;
let record = &bytes[record_offset..record_offset + STORE_INDEX_RECORD_LEN];
let lat_index = read_i32(record, INDEX_LAT_OFFSET)?;
let lon_index = read_i32(record, INDEX_LON_OFFSET)?;
let tile_id = (lat_index, lon_index);
if previous_id.is_some_and(|previous| tile_id <= previous) {
return Err(TerrainStoreError::Parse {
reason: "tile index records are not strictly sorted".to_string(),
});
}
previous_id = Some(tile_id);
let lon_count = read_u32(record, INDEX_LON_COUNT_OFFSET)?;
let lat_count = read_u32(record, INDEX_LAT_COUNT_OFFSET)?;
if lon_count < 2 || lat_count < 2 {
return Err(TerrainStoreError::Parse {
reason: format!(
"tile ({lat_index},{lon_index}) has invalid dimensions lon_count={lon_count} lat_count={lat_count}"
),
});
}
let offset = read_u64(record, INDEX_DATA_OFFSET_OFFSET)? as usize;
let data_len = read_u64(record, INDEX_DATA_LEN_OFFSET)? as usize;
let expected_len = (lon_count as usize)
.checked_mul(lat_count as usize)
.and_then(|count| count.checked_mul(2))
.ok_or_else(|| TerrainStoreError::Parse {
reason: format!("tile ({lat_index},{lon_index}) data length overflows usize"),
})?;
if data_len != expected_len {
return Err(TerrainStoreError::Parse {
reason: format!(
"tile ({lat_index},{lon_index}) data length must be {expected_len}, got {data_len}"
),
});
}
let expected_offset = align_up(expected_next, STORE_ALIGNMENT)?;
ensure_zero(bytes, expected_next, expected_offset, "tile padding")?;
if offset != expected_offset {
return Err(TerrainStoreError::Parse {
reason: format!(
"tile ({lat_index},{lon_index}) data offset must be {expected_offset}, got {offset}"
),
});
}
let end = offset
.checked_add(data_len)
.ok_or_else(|| TerrainStoreError::Parse {
reason: format!("tile ({lat_index},{lon_index}) data end overflows usize"),
})?;
if end > bytes.len() {
return Err(TerrainStoreError::Parse {
reason: format!("tile ({lat_index},{lon_index}) data extends past store length"),
});
}
let checksum64 = read_u64(record, INDEX_CHECKSUM_OFFSET)?;
let found = fnv1a64(&bytes[offset..end]);
if found != checksum64 {
return Err(TerrainStoreError::Checksum {
lat_index,
lon_index,
expected: checksum64,
found,
});
}
let min_latitude_deg = read_f64(record, INDEX_MIN_LAT_OFFSET)?;
let min_longitude_deg = read_f64(record, INDEX_MIN_LON_OFFSET)?;
let max_latitude_deg = read_f64(record, INDEX_MAX_LAT_OFFSET)?;
let max_longitude_deg = read_f64(record, INDEX_MAX_LON_OFFSET)?;
for (field, value) in [
("min_latitude_deg", min_latitude_deg),
("min_longitude_deg", min_longitude_deg),
("max_latitude_deg", max_latitude_deg),
("max_longitude_deg", max_longitude_deg),
] {
if !value.is_finite() {
return Err(TerrainStoreError::Parse {
reason: format!("tile ({lat_index},{lon_index}) {field} is not finite"),
});
}
}
let tile_datum = VerticalDatum::from_tag(record[INDEX_DATUM_OFFSET])?;
if tile_datum != vertical_datum {
return Err(TerrainStoreError::Parse {
reason: format!("tile ({lat_index},{lon_index}) datum differs from header"),
});
}
ensure_zero(
record,
INDEX_DATUM_OFFSET + 1,
STORE_INDEX_RECORD_LEN,
"tile index reserved bytes",
)?;
let index = TerrainStoreTileIndex {
lat_index,
lon_index,
min_longitude_deg,
min_latitude_deg,
max_longitude_deg,
max_latitude_deg,
lon_count,
lat_count,
data_offset: offset as u64,
data_len: data_len as u64,
checksum64,
vertical_datum: tile_datum,
};
by_grid.insert(tile_id, tiles.len());
tiles.push(MmapTile { index });
tile_index.push(index);
expected_next = end;
}
if expected_next != bytes.len() {
return Err(TerrainStoreError::Parse {
reason: format!(
"store has trailing bytes: expected length {expected_next}, got {}",
bytes.len()
),
});
}
Ok(ParsedStore {
vertical_datum,
tiles,
by_grid,
tile_index,
})
}
fn build_store(mut tiles: Vec<PendingTile>) -> core::result::Result<Vec<u8>, TerrainStoreError> {
tiles.sort_by_key(|tile| (tile.lat_index, tile.lon_index));
for pair in tiles.windows(2) {
if (pair[0].lat_index, pair[0].lon_index) == (pair[1].lat_index, pair[1].lon_index) {
return Err(TerrainStoreError::DuplicateTile {
lat_index: pair[0].lat_index,
lon_index: pair[0].lon_index,
});
}
}
let index_end = STORE_HEADER_LEN
.checked_add(
tiles
.len()
.checked_mul(STORE_INDEX_RECORD_LEN)
.ok_or_else(|| TerrainStoreError::Parse {
reason: "tile index length overflows usize".to_string(),
})?,
)
.ok_or_else(|| TerrainStoreError::Parse {
reason: "tile index end overflows usize".to_string(),
})?;
let data_offset = align_up(index_end, STORE_ALIGNMENT)?;
let mut offsets = Vec::with_capacity(tiles.len());
let mut cursor = data_offset;
for tile in &tiles {
cursor = align_up(cursor, STORE_ALIGNMENT)?;
offsets.push(cursor);
cursor = cursor
.checked_add(tile.data.len())
.ok_or_else(|| TerrainStoreError::Parse {
reason: "store length overflows usize".to_string(),
})?;
}
let mut out = vec![0u8; cursor];
out[..STORE_MAGIC.len()].copy_from_slice(STORE_MAGIC);
write_u16(&mut out, HEADER_VERSION_OFFSET, STORE_VERSION);
out[HEADER_DATUM_OFFSET] = VerticalDatum::Egm96MslOrthometric.tag();
write_u32(
&mut out,
HEADER_TILE_COUNT_OFFSET,
u32::try_from(tiles.len()).map_err(|_| TerrainStoreError::Parse {
reason: "tile count exceeds u32".to_string(),
})?,
);
write_u64(
&mut out,
HEADER_INDEX_OFFSET_OFFSET,
STORE_HEADER_LEN as u64,
);
write_u64(&mut out, HEADER_DATA_OFFSET_OFFSET, data_offset as u64);
write_u64(&mut out, HEADER_TOTAL_LEN_OFFSET, cursor as u64);
for (idx, tile) in tiles.iter().enumerate() {
let record_offset = STORE_HEADER_LEN + idx * STORE_INDEX_RECORD_LEN;
let offset = offsets[idx];
let data_len = tile.data.len();
let expected_len = (tile.lon_count as usize)
.checked_mul(tile.lat_count as usize)
.and_then(|count| count.checked_mul(2))
.ok_or_else(|| TerrainStoreError::Parse {
reason: format!(
"tile ({},{}) data length overflows usize",
tile.lat_index, tile.lon_index
),
})?;
if data_len != expected_len {
return Err(TerrainStoreError::Parse {
reason: format!(
"tile ({},{}) data length must be {expected_len}, got {data_len}",
tile.lat_index, tile.lon_index
),
});
}
let record = &mut out[record_offset..record_offset + STORE_INDEX_RECORD_LEN];
write_i32(record, INDEX_LAT_OFFSET, tile.lat_index);
write_i32(record, INDEX_LON_OFFSET, tile.lon_index);
write_u32(record, INDEX_LON_COUNT_OFFSET, tile.lon_count);
write_u32(record, INDEX_LAT_COUNT_OFFSET, tile.lat_count);
write_u64(record, INDEX_DATA_OFFSET_OFFSET, offset as u64);
write_u64(record, INDEX_DATA_LEN_OFFSET, data_len as u64);
write_u64(record, INDEX_CHECKSUM_OFFSET, fnv1a64(&tile.data));
write_f64(record, INDEX_MIN_LAT_OFFSET, tile.min_latitude_deg);
write_f64(record, INDEX_MIN_LON_OFFSET, tile.min_longitude_deg);
write_f64(record, INDEX_MAX_LAT_OFFSET, tile.max_latitude_deg);
write_f64(record, INDEX_MAX_LON_OFFSET, tile.max_longitude_deg);
record[INDEX_DATUM_OFFSET] = tile.vertical_datum.tag();
out[offset..offset + data_len].copy_from_slice(&tile.data);
}
Ok(out)
}
fn align_up(value: usize, alignment: usize) -> core::result::Result<usize, TerrainStoreError> {
let rem = value % alignment;
if rem == 0 {
Ok(value)
} else {
value
.checked_add(alignment - rem)
.ok_or_else(|| TerrainStoreError::Parse {
reason: "aligned offset overflows usize".to_string(),
})
}
}
fn ensure_zero(
bytes: &[u8],
start: usize,
end: usize,
context: &str,
) -> core::result::Result<(), TerrainStoreError> {
if start > end || end > bytes.len() {
return Err(TerrainStoreError::Parse {
reason: format!("{context} range is out of bounds"),
});
}
if bytes[start..end].iter().any(|&byte| byte != 0) {
return Err(TerrainStoreError::Parse {
reason: format!("{context} must be zero-filled"),
});
}
Ok(())
}
fn fnv1a64(bytes: &[u8]) -> u64 {
bytes.iter().fold(FNV_OFFSET_BASIS, |hash, byte| {
(hash ^ u64::from(*byte)).wrapping_mul(FNV_PRIME)
})
}
fn read_u16(bytes: &[u8], offset: usize) -> core::result::Result<u16, TerrainStoreError> {
Ok(u16::from_le_bytes(read_array(bytes, offset)?))
}
fn read_u32(bytes: &[u8], offset: usize) -> core::result::Result<u32, TerrainStoreError> {
Ok(u32::from_le_bytes(read_array(bytes, offset)?))
}
fn read_i32(bytes: &[u8], offset: usize) -> core::result::Result<i32, TerrainStoreError> {
Ok(i32::from_le_bytes(read_array(bytes, offset)?))
}
fn read_u64(bytes: &[u8], offset: usize) -> core::result::Result<u64, TerrainStoreError> {
Ok(u64::from_le_bytes(read_array(bytes, offset)?))
}
fn read_f64(bytes: &[u8], offset: usize) -> core::result::Result<f64, TerrainStoreError> {
Ok(f64::from_le_bytes(read_array(bytes, offset)?))
}
fn read_array<const N: usize>(
bytes: &[u8],
offset: usize,
) -> core::result::Result<[u8; N], TerrainStoreError> {
let end = offset
.checked_add(N)
.ok_or_else(|| TerrainStoreError::Parse {
reason: "numeric field offset overflows usize".to_string(),
})?;
let slice = bytes
.get(offset..end)
.ok_or_else(|| TerrainStoreError::Parse {
reason: "numeric field extends past record".to_string(),
})?;
slice.try_into().map_err(|_| TerrainStoreError::Parse {
reason: "numeric field has wrong length".to_string(),
})
}
fn write_u16(bytes: &mut [u8], offset: usize, value: u16) {
bytes[offset..offset + 2].copy_from_slice(&value.to_le_bytes());
}
fn write_u32(bytes: &mut [u8], offset: usize, value: u32) {
bytes[offset..offset + 4].copy_from_slice(&value.to_le_bytes());
}
fn write_i32(bytes: &mut [u8], offset: usize, value: i32) {
bytes[offset..offset + 4].copy_from_slice(&value.to_le_bytes());
}
fn write_u64(bytes: &mut [u8], offset: usize, value: u64) {
bytes[offset..offset + 8].copy_from_slice(&value.to_le_bytes());
}
fn write_f64(bytes: &mut [u8], offset: usize, value: f64) {
bytes[offset..offset + 8].copy_from_slice(&value.to_le_bytes());
}