use std::collections::HashMap;
use std::fs;
use std::path::{Path, PathBuf};
const UHL_SIZE: usize = 80;
const DSI_SIZE: usize = 648;
const ACC_SIZE: usize = 2700;
const DATA_OFFSET: usize = UHL_SIZE + DSI_SIZE + ACC_SIZE;
const DATA_SENTINEL: u8 = 0xAA;
const DTED_SUFFIX: &str = concat!("_1arc_v3.d", "t", "2");
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum DtedInterpolation {
NearestPosting,
Bilinear,
}
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub struct DtedLookupOptions {
pub interpolation: DtedInterpolation,
}
impl Default for DtedLookupOptions {
fn default() -> Self {
Self {
interpolation: DtedInterpolation::Bilinear,
}
}
}
#[derive(Debug)]
pub struct DtedTerrain {
root: PathBuf,
tiles: HashMap<(i32, i32), DtedTile>,
}
impl DtedTerrain {
pub fn new(root: impl Into<PathBuf>) -> Self {
Self {
root: root.into(),
tiles: HashMap::new(),
}
}
pub fn height_m(&mut self, longitude_deg: f64, latitude_deg: f64) -> Result<f64, String> {
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, String> {
let Some(tile) = self.load_tile(longitude_deg, latitude_deg)? else {
return Ok(0.0);
};
if options.interpolation == DtedInterpolation::NearestPosting {
return tile
.get_elevation(longitude_deg, latitude_deg)
.map(|v| v as f64);
}
let postings_per_deg_lon = tile.lon_count - 1;
let postings_per_deg_lat = tile.lat_count - 1;
let lon_idx = (longitude_deg - tile.origin_longitude) * postings_per_deg_lon as f64;
let lat_idx = (latitude_deg - tile.origin_latitude) * 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.origin_longitude + (lon_lo + di) as f64 / postings_per_deg_lon as f64;
let posting_lat =
tile.origin_latitude + (lat_lo + dj) as f64 / postings_per_deg_lat as f64;
z += w * tile.get_elevation(posting_lon, posting_lat)? as f64;
}
}
Ok(z)
}
fn load_tile(&mut self, longitude: f64, latitude: f64) -> Result<Option<&DtedTile>, String> {
let mut selected = None;
for grid_idx in terrain_grid_candidates(longitude, latitude) {
if !self.tiles.contains_key(&grid_idx) {
let Some(path) = self.terrain_path_for_grid(grid_idx.0, grid_idx.1) else {
continue;
};
if !path.is_file() {
continue;
}
let tile = DtedTile::from_path(path)?;
self.tiles.insert(grid_idx, tile);
}
if let Some(tile) = self.tiles.get(&grid_idx) {
if tile.contains(longitude, latitude) {
selected = Some(grid_idx);
break;
}
}
}
Ok(selected.and_then(|grid_idx| self.tiles.get(&grid_idx)))
}
fn terrain_path_for_grid(&self, latitude_index: i32, longitude_index: i32) -> Option<PathBuf> {
let tile_name = format!(
"{}_{}{}",
format_lat(latitude_index),
format_lon(longitude_index),
DTED_SUFFIX
);
let direct = self.root.join(&tile_name);
if direct.is_file() {
return Some(direct);
}
let block_dir = terrain_block_dir(latitude_index, longitude_index);
let nested = self.root.join(&block_dir).join(&tile_name);
if nested.is_file() {
return Some(nested);
}
let sibling = self.root.parent()?.join(&block_dir).join(&tile_name);
sibling.is_file().then_some(sibling)
}
}
#[derive(Debug)]
pub struct DtedTile {
origin_latitude: f64,
origin_longitude: f64,
lon_count: usize,
lat_count: usize,
data_block_length: usize,
bytes: Vec<u8>,
}
impl DtedTile {
pub fn from_path(path: impl AsRef<Path>) -> Result<Self, String> {
let bytes =
fs::read(path.as_ref()).map_err(|e| format!("{}: {e}", path.as_ref().display()))?;
if bytes.len() < DATA_OFFSET {
return Err(format!(
"{} is too short for DTED headers",
path.as_ref().display()
));
}
if &bytes[0..4] != b"UHL1" {
return Err(format!("{} missing UHL1 header", path.as_ref().display()));
}
let origin_longitude =
parse_dted_coord(std::str::from_utf8(&bytes[4..12]).map_err(|e| e.to_string())?)?;
let origin_latitude =
parse_dted_coord(std::str::from_utf8(&bytes[12..20]).map_err(|e| e.to_string())?)?;
let lon_count = parse_ascii_usize(&bytes[47..51])?;
let lat_count = parse_ascii_usize(&bytes[51..55])?;
let data_block_length = 12 + 2 * lat_count;
let expected_len = DATA_OFFSET + lon_count * data_block_length;
if bytes.len() < expected_len {
return Err(format!(
"{} has {} bytes but expected at least {}",
path.as_ref().display(),
bytes.len(),
expected_len
));
}
Ok(Self {
origin_latitude,
origin_longitude,
lon_count,
lat_count,
data_block_length,
bytes,
})
}
pub fn get_elevation(&self, longitude: f64, latitude: f64) -> Result<i16, String> {
if !self.contains(longitude, latitude) {
return Err(format!(
"point ({longitude},{latitude}) is outside DTED tile ({},{})",
self.origin_longitude, self.origin_latitude
));
}
let latitude_index =
py_round_to_usize((latitude - self.origin_latitude) * (self.lat_count - 1) as f64)?;
let longitude_index =
py_round_to_usize((longitude - self.origin_longitude) * (self.lon_count - 1) as f64)?;
if latitude_index >= self.lat_count || longitude_index >= self.lon_count {
return Err(format!(
"posting index out of bounds lon={} lat={}",
longitude_index, latitude_index
));
}
let block_start = DATA_OFFSET + longitude_index * self.data_block_length;
let block_end = block_start + self.data_block_length;
let block = &self.bytes[block_start..block_end];
if block[0] != DATA_SENTINEL {
return Err(format!(
"DTED block {longitude_index} missing data sentinel"
));
}
let checksum = i32::from_be_bytes([
block[block.len() - 4],
block[block.len() - 3],
block[block.len() - 2],
block[block.len() - 1],
]);
let sum = block[..block.len() - 4]
.iter()
.fold(0i32, |acc, b| acc + i32::from(*b));
if sum != checksum {
return Err(format!(
"DTED checksum failed for block {longitude_index}: expected {checksum}, found {sum}"
));
}
let sample_start = 8 + latitude_index * 2;
let raw = i16::from_be_bytes([block[sample_start], block[sample_start + 1]]);
Ok(convert_signed_magnitude(raw))
}
fn contains(&self, longitude: f64, latitude: f64) -> bool {
latitude >= self.origin_latitude
&& latitude <= self.origin_latitude + 1.0
&& longitude >= self.origin_longitude
&& longitude <= self.origin_longitude + 1.0
}
}
fn terrain_grid(longitude: f64, latitude: f64) -> (i32, i32) {
(latitude.floor() as i32, longitude.floor() as i32)
}
fn terrain_grid_candidates(longitude: f64, latitude: f64) -> Vec<(i32, i32)> {
let (lat, lon) = terrain_grid(longitude, latitude);
let mut out = vec![(lat, lon)];
let on_lat_edge = latitude == latitude.floor();
let on_lon_edge = longitude == longitude.floor();
if on_lat_edge {
out.push((lat - 1, lon));
}
if on_lon_edge {
out.push((lat, lon - 1));
}
if on_lat_edge && on_lon_edge {
out.push((lat - 1, lon - 1));
}
out
}
fn format_lat(latitude_index: i32) -> String {
if latitude_index >= 0 {
format!("n{:02}", latitude_index)
} else {
format!("s{:02}", -latitude_index)
}
}
fn format_lon(longitude_index: i32) -> String {
if longitude_index >= 0 {
format!("e{:03}", longitude_index)
} else {
format!("w{:03}", -longitude_index)
}
}
fn terrain_block_dir(latitude_index: i32, longitude_index: i32) -> String {
format!(
"{}_{}",
format_lat(block_origin(latitude_index)),
format_lon(block_origin(longitude_index))
)
}
fn block_origin(index: i32) -> i32 {
(index / 10) * 10
}
fn parse_ascii_usize(bytes: &[u8]) -> Result<usize, String> {
std::str::from_utf8(bytes)
.map_err(|e| e.to_string())?
.trim()
.parse::<usize>()
.map_err(|e| e.to_string())
}
fn parse_dted_coord(input: &str) -> Result<f64, String> {
let hemi = input
.chars()
.last()
.ok_or_else(|| "empty DTED coordinate".to_string())?;
let sign = match hemi {
'S' | 'W' => -1.0,
'N' | 'E' => 1.0,
_ => return Err(format!("invalid DTED hemisphere {hemi}")),
};
let coord = &input[..input.len() - 1];
let seconds_index = if coord.as_bytes().get(coord.len().saturating_sub(2)) == Some(&b'.') {
coord.len() - 4
} else {
coord.len() - 2
};
let minutes_index = seconds_index - 2;
let degree = coord[..minutes_index]
.parse::<i32>()
.map_err(|e| e.to_string())?;
let minute = coord[minutes_index..seconds_index]
.parse::<i32>()
.map_err(|e| e.to_string())?;
let second = coord[seconds_index..]
.parse::<f64>()
.map_err(|e| e.to_string())?;
Ok(sign * (degree as f64 + ((minute as f64 + second / 60.0) / 60.0)))
}
fn py_round_to_usize(value: f64) -> Result<usize, String> {
if value < 0.0 {
return Err(format!("cannot round negative posting index {value}"));
}
let lo = value.floor();
let frac = value - lo;
let rounded = if frac < 0.5 {
lo
} else if frac > 0.5 {
lo + 1.0
} else {
let lo_i = lo as u64;
if lo_i.is_multiple_of(2) {
lo
} else {
lo + 1.0
}
};
Ok(rounded as usize)
}
fn convert_signed_magnitude(raw: i16) -> i16 {
if raw < 0 {
(-32768i32 - i32::from(raw)) as i16
} else {
raw
}
}
#[cfg(test)]
mod tests {
use std::path::PathBuf;
use bitexact::parity::f64_from_hex;
use serde_json::Value;
use super::{terrain_block_dir, DtedInterpolation, DtedLookupOptions, DtedTerrain};
#[test]
fn terrain_block_dir_matches_reference_bucket_names() {
assert_eq!(terrain_block_dir(36, -107), "n30_w100");
assert_eq!(terrain_block_dir(32, -117), "n30_w110");
assert_eq!(terrain_block_dir(43, -112), "n40_w110");
assert_eq!(terrain_block_dir(20, -103), "n20_w100");
}
fn fixture_path(name: &str) -> PathBuf {
PathBuf::from(env!("CARGO_MANIFEST_DIR"))
.join("tests")
.join("fixtures")
.join("dted")
.join(name)
}
fn bits(v: &Value) -> f64 {
f64_from_hex(v.as_str().expect("hex-bit string")).expect("valid f64 bits")
}
#[test]
fn dted_lookup_matches_generated_fixture_bits() {
let raw =
std::fs::read_to_string(fixture_path("dted_points.json")).expect("read dted fixture");
let doc: Value = serde_json::from_str(&raw).expect("parse dted fixture");
assert_eq!(doc["schema"], "gnss-dted-points-v1");
let mut terrain = DtedTerrain::new(fixture_path("tiles"));
let nearest = DtedLookupOptions {
interpolation: DtedInterpolation::NearestPosting,
};
let bilinear = DtedLookupOptions {
interpolation: DtedInterpolation::Bilinear,
};
let mut checked = 0usize;
for case in doc["nearest_cases"].as_array().expect("nearest_cases") {
let lon = bits(&case["longitude_bits"]);
let lat = bits(&case["latitude_bits"]);
let got = terrain
.height_m_with_options(lon, lat, nearest)
.expect("nearest DTED height");
let want = bits(&case["elevation_bits"]);
assert_eq!(
got.to_bits(),
want.to_bits(),
"nearest DTED {},{}",
lon,
lat
);
checked += 1;
}
for case in doc["bilinear_cases"].as_array().expect("bilinear_cases") {
let lon = bits(&case["longitude_bits"]);
let lat = bits(&case["latitude_bits"]);
let got = terrain
.height_m_with_options(lon, lat, bilinear)
.expect("bilinear DTED height");
let want = bits(&case["elevation_bits"]);
assert_eq!(
got.to_bits(),
want.to_bits(),
"bilinear DTED {},{}",
lon,
lat
);
checked += 1;
}
assert!(checked > 0, "empty DTED fixture");
}
}