diff --git a/Cargo.lock b/Cargo.lock index cf27cd6fc..54776df3f 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -4509,6 +4509,7 @@ dependencies = [ "crs-definitions", "csv", "flatgeobuf", + "geo", "geozero", "image", "laz", diff --git a/ROADMAP.md b/ROADMAP.md index 8ec5a6132..732506709 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -10,8 +10,8 @@ Many features are first versions: there are no full 3D scenes or enterprise port - Data: files by their web address; GeoJSON, CSV/TSV (x/y in any EPSG system, or WKT), Excel and OpenDocument spreadsheets, GeoRSS, LandXML, encoded polylines, Shapefile (zipped too, with the coordinate system in its `.prj`), GeoPackage, GeoParquet, FlatGeobuf, KML/KMZ (styles included), GPX, GML, MapInfo MIF, OSM PBF, TopoJSON, DXF and LiDAR (LAS/LAZ/COPC) in, all of them - in the web app too (large files read in place, a view at a time); the vector formats GDAL - reads when it is installed, PostGIS, MySQL and MariaDB tables read live (desktop), and SQL queries + in the web app too (large files read in place, a view at a time); Esri File Geodatabases without GDAL (desktop; 3D multipatches + as footprints with their heights); the vector formats GDAL reads when it is installed, PostGIS, MySQL and MariaDB tables read live (desktop), and SQL queries and DuckDB databases (DuckDB Spatial) when DuckDB is; GeoJSON, CSV, GeoPackage, GeoParquet, Shapefile, FlatGeobuf, KML/KMZ, GPX and WKT out; standalone tables; geotagged photos and KML ground overlays; geocoding (Find Place, address tables); a Data Sources pane (favourites, recent, saved services, a library) and an SQL pane diff --git a/crates/data/Cargo.toml b/crates/data/Cargo.toml index cad9f3d08..ddb1ec459 100644 --- a/crates/data/Cargo.toml +++ b/crates/data/Cargo.toml @@ -28,6 +28,8 @@ parquet = { workspace = true, features = ["zstd"] } url = { workspace = true } bytes = { workspace = true } flatgeobuf = { workspace = true } +# Polygon union: a multipatch's footprint (File Geodatabase). +geo = { workspace = true } geozero = { workspace = true } calamine = { workspace = true } # LAZ (compressed LiDAR) point clouds. diff --git a/crates/data/src/filegdb.rs b/crates/data/src/filegdb.rs new file mode 100644 index 000000000..2f2102b08 --- /dev/null +++ b/crates/data/src/filegdb.rs @@ -0,0 +1,1168 @@ +//! Esri File Geodatabases (a `.gdb` folder), read without GDAL: every feature class and table, with +//! their fields and coordinate system. +//! +//! The format is ESRI's own and undocumented; this follows the description published by the GDAL +//! project (`dump_gdbtable`'s FGDB-Spec). A geodatabase is a folder of tables, `aXXXXXXXX.gdbtable` +//! (fields and rows) with `aXXXXXXXX.gdbtablx` (where each row starts), numbered by their row in +//! the system catalog, `a00000001.gdbtable`. +//! +//! Geometry is read in 2D. A multipatch (3D surfaces, such as the buildings of swissBUILDINGS3D) +//! becomes its footprint: the union of its faces seen from above, with its lowest and highest Z in +//! the fields `z_min` and `z_max`, and the difference in `height` (an extrusion stands it up by it). +//! In the browser, where a panic in the union can't be caught, a multipatch is the hull of its +//! points. Curves are drawn as straight segments between their ends. +//! Compressed tables, raster fields and the newer sparse row layouts aren't read: they give +//! [`DataError::Unsupported`], and the caller can read the file through GDAL instead. +//! +//! Files are hostile input: every count and offset is checked, rows that can't be read are skipped +//! and counted. + +use crate::raster::{Crs, Transformer}; +use crate::{DataError, Dataset, Field, FieldType, Geometry, GeometryType, ReadReport, Result, Value}; +use rustgis_geom::Coord; +use std::path::Path; + +/// The largest table read (its `.gdbtable` is read whole). +pub const MAX_TABLE_BYTES: usize = if cfg!(target_arch = "wasm32") { 512 << 20 } else { 2 << 30 }; +/// Most tables read from one geodatabase. +const MAX_TABLES: usize = 512; +/// Most fields a table may have. +const MAX_FIELDS: usize = 4096; +/// Most parts and points one geometry may have (its blob's length bounds them first). +const MAX_PARTS: u64 = 1 << 24; +/// Most faces of one multipatch merged into its footprint; beyond, its outline is its points' hull. +const MAX_FACES: usize = 200_000; +/// Most faces of one multipatch merged by a union; past it, the hull too (a union's time grows +/// faster than its faces, so one crafted row mustn't stall loading). The browser merges none: a +/// panic in the union library can't be caught there, and hostile input mustn't end the page, so +/// multipatches are always the hulls of their points (see [`footprint`]). On the desktop the catch +/// relies on the workspace's default `panic = "unwind"`. +const MAX_UNION_FACES: usize = 50_000; +/// Most faces merged for one table, whatever its rows ask: past it, multipatches are hulls, so a +/// crafted table can't make loading take hours. +const MAX_TABLE_FACES: usize = if cfg!(target_arch = "wasm32") { 500_000 } else { 4_000_000 }; + +/// Whether a folder is a File Geodatabase: named `.gdb`, holding a system catalog. +pub fn is_geodatabase(dir: &Path) -> bool { + crate::source::is_gdb_name(&dir.to_string_lossy()) && dir.join("a00000001.gdbtable").is_file() +} + +/// What a geodatabase gave: its layers, and the tables that couldn't be read (by name, with why), +/// which GDAL may read instead. +#[derive(Debug, Default)] +pub struct Read { + pub layers: Vec<(Dataset, ReadReport)>, + pub unread: Vec<(String, String)>, +} + +/// Reads the feature classes and tables of a geodatabase folder; `layer` picks one by name +/// (ignoring case). +pub fn read_dir(dir: &Path, layer: Option<&str>) -> Result { + let name = dir.file_stem().and_then(|s| s.to_str()).unwrap_or("geodatabase").to_string(); + let mut out = read(&name, &|file: &str| read_file(dir, file), layer)?; + for (ds, _) in &mut out.layers { + ds.source = Some(dir.display().to_string()); + } + Ok(out) +} + +/// The names of a geodatabase folder's feature classes and tables, from its catalog alone. +pub fn layers(dir: &Path) -> Result> { + let name = dir.file_stem().and_then(|s| s.to_str()).unwrap_or("geodatabase"); + let read = |file: &str| -> Result> { + read_file(dir, file)?.ok_or_else(|| DataError::Format(format!("{name}: not a File Geodatabase (no {file})"))) + }; + Ok(own_tables(catalog(name, &read("a00000001.gdbtable")?, &read("a00000001.gdbtablx")?)?).map(|(_, n, _)| n).collect()) +} + +/// One file of a geodatabase folder, `None` when it isn't there; one past [`MAX_TABLE_BYTES`] is +/// refused before it is read. +fn read_file(dir: &Path, file: &str) -> Result>> { + let path = dir.join(file); + match std::fs::metadata(&path) { + Ok(m) if m.len() > MAX_TABLE_BYTES as u64 => Err(DataError::Unsupported(format!("{file} is too large to read here"))), + Ok(_) => Ok(Some(std::fs::read(&path)?)), + Err(e) if e.kind() == std::io::ErrorKind::NotFound => Ok(None), + Err(e) => Err(e.into()), + } +} + +/// The tables the catalog lists that hold the geodatabase's data, not its own bookkeeping +/// (`GDB_…`) nor nothing at all: their number, name and whether they are compressed. What +/// [`layers`] lists and [`read`] reads (or says it can't). +fn own_tables(listed: Vec<(u64, String, Option)>) -> impl Iterator { + listed + .into_iter() + .filter(|(_, n, _)| !n.starts_with("GDB_")) + .filter_map(|(row, n, format)| match format { + None | Some(0.0) => Some((row, n, false)), + Some(1.0) => Some((row, n, true)), + _ => None, + }) + .take(MAX_TABLES) +} + +/// Each table the system catalog lists: its number, name and storage (0: a table of its own; +/// 1: compressed; 2: none). +fn catalog(name: &str, t: &[u8], x: &[u8]) -> Result)>> { + let table = Table::parse(t)?; + let rows = row_offsets(x)?; + let format_at = table.field("FileFormat"); + let Some(name_at) = table.field("Name") else { + return Err(DataError::Format(format!("{name}: the system catalog has no Name field"))); + }; + Ok(rows + .into_iter() + .filter_map(|(row, offset)| { + let values = table.row(t, offset)?; + let Some(Value::Text(name)) = values.get(name_at).map(|v| v.value.clone()) else { return None }; + Some((row, name, format_at.and_then(|i| values.get(i)).and_then(|v| v.value.as_f64()))) + }) + .collect()) +} + +/// Reads the feature classes and tables of a geodatabase whose files `open` gives by name +/// (`a00000001.gdbtable`; `None` for a file it doesn't have). +pub fn read(name: &str, open: &dyn Fn(&str) -> Result>>, layer: Option<&str>) -> Result { + let t = open("a00000001.gdbtable")?.ok_or_else(|| DataError::Format(format!("{name}: not a File Geodatabase (no system catalog)")))?; + let x = open("a00000001.gdbtablx")?.ok_or_else(|| DataError::Format(format!("{name}: the system catalog has no row index")))?; + let listed = catalog(name, &t, &x)?; + // The coordinate system codes the item definitions give, by table name. + let wkids = listed + .iter() + .find(|(_, n, _)| n == "GDB_Items") + .and_then(|(row, _, _)| { + let file = format!("a{row:08x}"); + Some(wkids(&open(&format!("{file}.gdbtable")).ok()??, &open(&format!("{file}.gdbtablx")).ok()??)) + }) + .unwrap_or_default(); + let mut out = Read::default(); + for (row, table_name, compressed) in own_tables(listed) { + if layer.is_some_and(|l| !l.eq_ignore_ascii_case(&table_name)) { + continue; + } + if compressed { + out.unread.push((table_name, "it is compressed".into())); + continue; + } + // A table that can't be opened (too large, not readable) is said; the others are read. + let file = format!("a{row:08x}"); + let files = open(&format!("{file}.gdbtable")).and_then(|t| Ok((t, open(&format!("{file}.gdbtablx"))?))); + let (t, x) = match files { + Ok((Some(t), Some(x))) => (t, x), + // Listed, but its files aren't there: said by name, as `layers` lists it. + Ok((t, x)) => { + let missing: Vec = [(t.is_none(), "gdbtable"), (x.is_none(), "gdbtablx")] + .iter() + .filter(|(gone, _)| *gone) + .map(|(_, e)| format!("{file}.{e}")) + .collect(); + let why = match missing.as_slice() { + [one] => format!("its file {one} is missing"), + _ => format!("its files {} are missing", missing.join(" and ")), + }; + out.unread.push((table_name, why)); + continue; + } + Err(e) => { + out.unread.push((table_name, e.to_string())); + continue; + } + }; + // Items that share the table's name but disagree on the code give none: the WKT decides. + let mut codes = wkids.iter().filter(|(n, _)| n.eq_ignore_ascii_case(&table_name)).map(|&(_, c)| c); + let wkid = codes.next().filter(|&c| codes.all(|d| d == c)); + match read_table(&table_name, &t, &x, wkid) { + Ok(Some(read)) => out.layers.push(read), + Ok(None) => {} + Err(e) => out.unread.push((table_name, e.to_string())), + } + } + if out.layers.is_empty() { + return Err(match (out.unread.first(), layer) { + (Some((table, why)), _) => DataError::Unsupported(format!("{name}: {table}: {why}")), + (None, Some(l)) => DataError::Format(format!("{name}: no layer '{l}'")), + (None, None) => DataError::Format(format!("{name}: no feature classes or tables")), + }); + } + // A note several layers share is said once. + let mut said: Vec = Vec::new(); + for (_, report) in &mut out.layers { + report.messages.retain(|m| !said.contains(m)); + said.extend(report.messages.iter().cloned()); + } + Ok(out) +} + +/// The coordinate system code (`LatestWKID`, else `WKID`) of each item `GDB_Items` defines, by +/// name: the WKT a table holds is Esri's, which names some projections differently. +fn wkids(t: &[u8], x: &[u8]) -> Vec<(String, u32)> { + let (Ok(table), Ok(rows)) = (Table::parse(t), row_offsets(x)) else { return Vec::new() }; + let (Some(name_at), Some(def_at)) = (table.field("Name"), table.field("Definition")) else { return Vec::new() }; + let tag = |xml: &str, tag: &str| -> Option { + let open = format!("<{tag}>"); + let start = xml.find(&open)? + open.len(); + xml.get(start..)?.split('<').next()?.trim().parse().ok() + }; + rows.into_iter() + .filter_map(|(_, offset)| { + let values = table.row(t, offset)?; + let Some(Value::Text(name)) = values.get(name_at).map(|v| v.value.clone()) else { return None }; + let Some(Value::Text(def)) = values.get(def_at).map(|v| &v.value) else { return None }; + Some((name, tag(def, "LatestWKID").or_else(|| tag(def, "WKID"))?)) + }) + .collect() +} + +/// A coordinate system by its EPSG or Esri code. +fn crs_of_code(code: u32) -> Option { + Crs::from_epsg(code).or_else(|| { + let proj4 = crate::raster::crs::esri_proj4(code)?; + Crs::custom(crate::raster::crs::CustomCrs { + proj4: proj4.to_string(), + label: format!("ESRI:{code}"), + geokeys: Vec::new(), + doubles: Vec::new(), + ascii: String::new(), + }) + .ok() + }) +} + +/// One table as a dataset: `None` for a table with no rows and no fields worth showing. `wkid` is +/// its coordinate system's code, when the item definitions give one. +fn read_table(name: &str, t: &[u8], x: &[u8], wkid: Option) -> Result> { + let table = Table::parse(t)?; + let rows = row_offsets(x)?; + let mut report = ReadReport::default(); + let layer_type = match table.shape { + _ if table.geometry.is_none() => GeometryType::None, + Shape::Point | Shape::Multipoint => GeometryType::Point, + Shape::Polyline => GeometryType::Polyline, + Shape::Polygon | Shape::Multipatch => GeometryType::Polygon, + Shape::None => GeometryType::None, + Shape::Other(code) => return Err(DataError::Unsupported(format!("geometry type {code}"))), + }; + // Its coordinate system: by its code, else from the geometry field's WKT. + let by_code = wkid.filter(|_| table.geometry.is_some()).and_then(crs_of_code).and_then(|c| match c.transformer() { + Ok(t) if c != Crs::LonLat => Some((Some(c), Some(t))), + Ok(_) => Some((None, None)), + Err(_) => None, + }); + let wkt = table.geometry.as_ref().map(|g| g.wkt.as_str()).filter(|w| !w.is_empty() && !w.starts_with('{')); + let (crs, to_lonlat) = match (by_code, wkt) { + (Some(by_code), _) => by_code, + (None, Some(wkt)) => match crate::raster::wkt::read_wkt(wkt) { + Ok((c, note)) => { + if let Some(n) = note { + report.messages.push(format!("coordinate system: {n}")); + } + match c.transformer() { + Ok(t) if c != Crs::LonLat => (Some(c), Some(t)), + Ok(_) => (None, None), + Err(e) => { + report.messages.push(format!("coordinate system: {e}; coordinates are read as longitude/latitude")); + (None, None) + } + } + } + Err(e) => { + report.messages.push(format!("coordinate system: {e}; coordinates are read as longitude/latitude")); + (None, None) + } + }, + (None, None) => (None, None), + }; + // The attribute fields: every field but the geometry, and the object ID first. + let mut fields: Vec = Vec::new(); + let mut columns: Vec = Vec::new(); + for (i, f) in table.fields.iter().enumerate() { + let ty = match f.kind { + Kind::Geometry | Kind::Binary | Kind::Raster => continue, + Kind::Int16 | Kind::Int32 | Kind::Int64 | Kind::ObjectId => FieldType::Integer, + Kind::Float32 | Kind::Float64 => FieldType::Double, + Kind::Text | Kind::Xml | Kind::Guid | Kind::GlobalId | Kind::Date | Kind::DateOnly | Kind::TimeOnly | Kind::DateTimeOffset => { + FieldType::Text + } + }; + let alias = (!f.alias.is_empty() && f.alias != f.name).then(|| f.alias.clone()); + fields.push(Field { name: f.name.clone(), ty, alias }); + columns.push(i); + } + let patches = table.shape == Shape::Multipatch; + // A multipatch's lowest and highest Z, and the difference (what an extrusion stands it up by), + // in three fields after the table's own (named apart from them). + let z_at = patches.then(|| { + let free = |base: &str| -> String { + let taken = |n: &str| fields.iter().any(|f| f.name.eq_ignore_ascii_case(n)); + (0..100) + .map(|k| if k == 0 { base.to_string() } else { format!("{base}_{k}") }) + .find(|n| !taken(n)) + .unwrap_or_else(|| base.to_string()) + }; + let names = [free("z_min"), free("z_max"), free("height")]; + let at = fields.len(); + fields.extend(names.into_iter().map(|name| Field { name, ty: FieldType::Double, alias: None })); + at + }); + if rows.is_empty() && fields.is_empty() { + return Ok(None); + } + let mut ds = Dataset::new(name, layer_type, fields); + ds.crs = crs; + let geometry_at = table.fields.iter().position(|f| f.kind == Kind::Geometry); + let (mut curves, mut wall_hulls, mut unmerged, mut unprojected) = (0usize, 0usize, 0usize, 0usize); + let mut budget = MAX_TABLE_FACES; + for (row, offset) in rows { + let Some(mut values) = table.row(t, offset) else { + report.skipped += 1; + continue; + }; + let mut attrs: Vec = Vec::with_capacity(ds.fields.len()); + for &c in &columns { + attrs.push(match table.fields.get(c).map(|f| f.kind) { + Some(Kind::ObjectId) => Value::Int(i64::try_from(row).unwrap_or(i64::MAX)), + // Taken, not cloned: each column is read once, and the geometry's blob stays. + _ => values.get_mut(c).map(|v| std::mem::replace(&mut v.value, Value::Null)).unwrap_or(Value::Null), + }); + } + let mut geometry = None; + if let (Some(g), Some(def)) = (geometry_at.and_then(|i| values.get(i)).and_then(|v| v.blob.as_ref()), table.geometry.as_ref()) { + match decode_geometry(g, def, &mut budget) { + Ok(Decoded::Geometry(geom, curved)) => { + curves += curved; + geometry = Some(geom); + } + Ok(Decoded::Patch(patch)) => { + if let (Some((lo, hi)), Some(at)) = (patch.z, z_at) { + attrs.resize(at, Value::Null); + attrs.extend([Value::Float(lo), Value::Float(hi), Value::Float(hi - lo)]); + } + wall_hulls += usize::from(patch.hull == Some(Hull::NoFaces)); + unmerged += usize::from(patch.hull == Some(Hull::Unmerged)); + geometry = patch.footprint; + } + Ok(Decoded::Empty) => {} + Err(()) => report.skipped += 1, + } + } + let geometry = match (geometry, &to_lonlat) { + (Some(g), Some(t)) => { + let g = to_lonlat_geometry(g, t); + unprojected += usize::from(g.is_none()); + g + } + (g, _) => g, + }; + ds.push_next(geometry, attrs); + } + if curves > 0 { + report.messages.push(format!("{curves} curved segments are drawn as straight lines")); + } + if unprojected > 0 { + report + .messages + .push(format!("{unprojected} features have coordinates outside their coordinate system: they are kept without geometry")); + } + if patches { + report.messages.push("3D multipatches are drawn as their footprints, with their lowest and highest Z in z_min and z_max, and the difference in height".into()); + } + if wall_hulls > 0 { + report.messages.push(format!("{wall_hulls} multipatches with no faces seen from above are drawn as the outline of their points")); + } + if unmerged > 0 { + report.messages.push(if cfg!(target_arch = "wasm32") { + format!("{unmerged} multipatches are drawn as the outline of their points (the desktop app draws their footprints)") + } else { + format!("{unmerged} multipatches with too many faces to merge are drawn as the outline of their points") + }); + } + Ok(Some((ds, report))) +} + +/// Projects a geometry to longitude/latitude; `None` when a point can't be. +fn to_lonlat_geometry(g: Geometry, t: &Transformer) -> Option { + let mut ok = true; + let out = g.map_coords(&|c| t.to_lonlat(c).unwrap_or(Coord::new(f64::NAN, f64::NAN))); + out.for_each_coord(|c: Coord| ok &= c.is_finite()); + ok.then_some(out) +} + +// --------------------------------------------------------------------------------------------- +// Bytes + +/// A reader over a byte slice that never reads past its end. +struct Cur<'a> { + b: &'a [u8], + i: usize, +} + +impl<'a> Cur<'a> { + fn new(b: &'a [u8], i: usize) -> Self { + Cur { b, i } + } + fn bytes(&mut self, n: usize) -> Option<&'a [u8]> { + let out = self.b.get(self.i..self.i.checked_add(n)?)?; + self.i += n; + Some(out) + } + fn skip(&mut self, n: usize) -> Option<()> { + self.bytes(n).map(|_| ()) + } + fn u8(&mut self) -> Option { + self.bytes(1).and_then(|b| b.first().copied()) + } + fn arr(&mut self) -> Option<[u8; N]> { + self.bytes(N).and_then(|b| b.try_into().ok()) + } + fn i16(&mut self) -> Option { + self.arr().map(i16::from_le_bytes) + } + fn i32(&mut self) -> Option { + self.arr().map(i32::from_le_bytes) + } + fn u32(&mut self) -> Option { + self.arr().map(u32::from_le_bytes) + } + fn i64(&mut self) -> Option { + self.arr().map(i64::from_le_bytes) + } + fn f32(&mut self) -> Option { + self.arr().map(f32::from_le_bytes) + } + fn f64(&mut self) -> Option { + self.arr().map(f64::from_le_bytes) + } + /// An unsigned little-endian base-128 number (at most 64 bits). + fn varuint(&mut self) -> Option { + let mut v: u64 = 0; + for shift in (0..64).step_by(7) { + let part = u64::from(self.u8()?); + // Bits past the 64th are refused, not dropped. + if shift > 57 && (part & 0x7F) >> (64 - shift) != 0 { + return None; + } + v |= (part & 0x7F) << shift; + if part & 0x80 == 0 { + return Some(v); + } + } + None + } + /// A signed base-128 number: the first byte's 0x40 bit is the sign. + fn varint(&mut self) -> Option { + let c = self.u8()?; + let negative = c & 0x40 != 0; + let mut v: u64 = u64::from(c & 0x3F); + let mut more = c & 0x80 != 0; + let mut shift = 6u32; + while more { + let part = u64::from(self.u8()?); + // Bits past the 64th are refused, not dropped. + if shift >= 64 || (shift > 57 && (part & 0x7F) >> (64 - shift) != 0) { + return None; + } + v |= (part & 0x7F) << shift; + shift += 7; + more = part & 0x80 != 0; + } + // A magnitude past `i64::MAX` (which `-i64::MIN` would overflow) is refused. + let v = i64::try_from(v).ok()?; + Some(if negative { -v } else { v }) + } + /// A count, bounded by `max` and by the bytes left (each item takes at least `min_bytes`). + fn count(&mut self, max: u64, min_bytes: usize) -> Option { + let n = self.varuint()?; + let left = (self.b.len().saturating_sub(self.i) / min_bytes.max(1)) as u64; + (n <= max && n <= left.max(1)).then(|| usize::try_from(n).ok()).flatten() + } +} + +fn utf16(b: &[u8]) -> String { + let units: Vec = b.as_chunks::<2>().0.iter().map(|c| u16::from_le_bytes(*c)).collect(); + String::from_utf16_lossy(&units) +} + +// --------------------------------------------------------------------------------------------- +// Tables + +/// A field's type. +#[derive(Clone, Copy, Debug, PartialEq)] +enum Kind { + Int16, + Int32, + Float32, + Float64, + Text, + Date, + ObjectId, + Geometry, + Binary, + Raster, + Guid, + GlobalId, + Xml, + Int64, + DateOnly, + TimeOnly, + DateTimeOffset, +} + +impl Kind { + fn of(code: u8) -> Option { + Some(match code { + 0 => Kind::Int16, + 1 => Kind::Int32, + 2 => Kind::Float32, + 3 => Kind::Float64, + 4 => Kind::Text, + 5 => Kind::Date, + 6 => Kind::ObjectId, + 7 => Kind::Geometry, + 8 => Kind::Binary, + 9 => Kind::Raster, + 10 => Kind::Guid, + 11 => Kind::GlobalId, + 12 => Kind::Xml, + 13 => Kind::Int64, + 14 => Kind::DateOnly, + 15 => Kind::TimeOnly, + 16 => Kind::DateTimeOffset, + _ => return None, + }) + } +} + +#[derive(Clone, Debug)] +struct FieldDef { + name: String, + alias: String, + kind: Kind, + /// Its bit in each row's null flags, for a field that can be null. + null_bit: Option, +} + +/// How a geometry field's coordinates are stored. +#[derive(Clone, Debug)] +struct GeomDef { + wkt: String, + x0: f64, + y0: f64, + xy_scale: f64, + z0: f64, + z_scale: f64, +} + +/// A table's kind of geometry. +#[derive(Clone, Copy, Debug, PartialEq)] +enum Shape { + None, + Point, + Multipoint, + Polyline, + Polygon, + Multipatch, + Other(u8), +} + +struct Table { + fields: Vec, + geometry: Option, + shape: Shape, + /// Strings are UTF-8 (else UTF-16). + utf8: bool, + /// How many fields have a bit in each row's null flags. + nullable: usize, +} + +/// One field's value in a row: a value, or a geometry's bytes. +struct Cell<'a> { + value: Value, + blob: Option<&'a [u8]>, +} + +impl Table { + fn parse(t: &[u8]) -> Result { + let bad = |what: &str| DataError::Format(format!("unreadable table ({what})")); + let mut c = Cur::new(t, 32); + let at = c.i64().and_then(|o| usize::try_from(o).ok()).ok_or_else(|| bad("header"))?; + let mut c = Cur::new(t, at); + c.skip(8).ok_or_else(|| bad("field header"))?; + let flags = c.u32().ok_or_else(|| bad("field header"))?; + let count = c.i16().and_then(|n| usize::try_from(n).ok()).filter(|&n| n <= MAX_FIELDS).ok_or_else(|| bad("field count"))?; + let has_z = flags & 0x8000_0000 != 0; + let has_m = flags & 0x4000_0000 != 0; + let shape = match flags & 0xFF { + 0 => Shape::None, + 1 => Shape::Point, + 2 => Shape::Multipoint, + 3 => Shape::Polyline, + 4 => Shape::Polygon, + 9 => Shape::Multipatch, + other => Shape::Other(other as u8), + }; + let mut fields = Vec::with_capacity(count); + let mut geometry = None; + for _ in 0..count { + let field = (|| -> Option<(FieldDef, Option)> { + let n = usize::from(c.u8()?); + let name = utf16(c.bytes(n * 2)?); + let n = usize::from(c.u8()?); + let alias = utf16(c.bytes(n * 2)?); + let code = c.u8()?; + let kind = Kind::of(code)?; + let mut geom = None; + let flag = match kind { + Kind::Text => { + c.skip(4)?; + let flag = c.u8()?; + let default = usize::try_from(c.varuint()?).ok()?; + if flag & 4 != 0 { + c.skip(default)?; + } + flag + } + Kind::ObjectId => { + c.skip(1)?; + c.u8()? + } + Kind::Geometry => { + c.skip(1)?; + let flag = c.u8()?; + let n = usize::try_from(c.i16()?).ok()?; + let wkt = utf16(c.bytes(n)?); + let g = c.u8()?; + let (z, m) = (g & 2 != 0, g & 4 != 0); + let (x0, y0, xy_scale) = (c.f64()?, c.f64()?, c.f64()?); + if m { + c.skip(16)?; + } + let (z0, z_scale) = if z { (c.f64()?, c.f64()?) } else { (0.0, 1.0) }; + // Tolerances, then the extent (and its Z and M ranges). + c.skip(8 + if m { 8 } else { 0 } + if z { 8 } else { 0 } + 32)?; + if has_z { + c.skip(16)?; + } + if has_m { + c.skip(16)?; + } + // The spatial index's grid sizes. + c.skip(1)?; + let grids = usize::try_from(c.u32()?).ok().filter(|&n| n <= 3)?; + c.skip(grids * 8)?; + let usable = |s: f64| s.is_finite() && s > 0.0; + if !usable(xy_scale) || !x0.is_finite() || !y0.is_finite() || (z && (!usable(z_scale) || !z0.is_finite())) { + return None; + } + geom = Some(GeomDef { wkt, x0, y0, xy_scale, z0, z_scale }); + flag + } + Kind::Binary | Kind::Guid | Kind::GlobalId | Kind::Xml => { + c.skip(1)?; + c.u8()? + } + // Raster fields have a layout of their own; such tables are left to GDAL. + Kind::Raster => return None, + _ => { + c.skip(1)?; + let flag = c.u8()?; + let default = usize::from(c.u8()?); + if flag & 4 != 0 { + c.skip(default)?; + } + flag + } + }; + let nullable = kind != Kind::ObjectId && flag & 1 != 0; + Some((FieldDef { name, alias, kind, null_bit: nullable.then_some(0) }, geom)) + })(); + let Some((f, g)) = field else { + return Err(DataError::Unsupported("a field this reader doesn't know (a raster field?)".into())); + }; + if g.is_some() && geometry.is_none() { + geometry = g; + } + fields.push(f); + } + // The null flags have a bit for each field that can be null, in order. + let mut nullable = 0usize; + for f in &mut fields { + if f.null_bit.is_some() { + f.null_bit = Some(nullable); + nullable += 1; + } + } + Ok(Table { fields, geometry, shape, utf8: flags & 0x100 != 0, nullable }) + } + + fn field(&self, name: &str) -> Option { + self.fields.iter().position(|f| f.name.eq_ignore_ascii_case(name)) + } + + /// A row's values, by field (the object ID's is `Null`: it isn't stored). `None` when the row + /// can't be read. + fn row<'a>(&self, t: &'a [u8], offset: usize) -> Option>> { + let mut c = Cur::new(t, offset); + let size = usize::try_from(c.i32()?).ok()?; + let body = c.bytes(size)?; + let mut c = Cur::new(body, 0); + let nulls = c.bytes(self.nullable.div_ceil(8))?; + let mut out = Vec::with_capacity(self.fields.len()); + for f in &self.fields { + let null = f.null_bit.is_some_and(|bit| nulls.get(bit / 8).is_some_and(|b| b >> (bit % 8) & 1 == 1)); + if null || f.kind == Kind::ObjectId { + out.push(Cell { value: Value::Null, blob: None }); + continue; + } + let text = |b: &[u8]| if self.utf8 { String::from_utf8_lossy(b).into_owned() } else { utf16(b) }; + let cell = |value: Value| Cell { value, blob: None }; + out.push(match f.kind { + Kind::Int16 => cell(Value::Int(i64::from(c.i16()?))), + Kind::Int32 => cell(Value::Int(i64::from(c.i32()?))), + Kind::Int64 => cell(Value::Int(c.i64()?)), + Kind::Float32 => cell(float(f64::from(c.f32()?))), + Kind::Float64 => cell(float(c.f64()?)), + Kind::Text | Kind::Xml => { + let n = usize::try_from(c.varuint()?).ok()?; + cell(Value::Text(text(c.bytes(n)?))) + } + Kind::Date => cell(date(c.f64()?, "YYYY-MM-DDTHH:mm:ss")), + Kind::DateOnly => cell(date(c.f64()?, "YYYY-MM-DD")), + Kind::TimeOnly => cell(date(c.f64()?, "HH:mm:ss")), + Kind::DateTimeOffset => { + let days = c.f64()?; + let minutes = c.i16()?; + cell(match date(days, "YYYY-MM-DDTHH:mm:ss") { + Value::Text(s) => Value::Text(format!( + "{s}{}{:02}:{:02}", + if minutes < 0 { '-' } else { '+' }, + minutes.unsigned_abs() / 60, + minutes.unsigned_abs() % 60 + )), + v => v, + }) + } + Kind::Guid | Kind::GlobalId => cell(Value::Text(guid(c.bytes(16)?))), + Kind::Geometry => { + let n = usize::try_from(c.varuint()?).ok()?; + Cell { value: Value::Null, blob: Some(c.bytes(n)?) } + } + Kind::Binary => { + let n = usize::try_from(c.varuint()?).ok()?; + c.skip(n)?; + cell(Value::Null) + } + Kind::ObjectId | Kind::Raster => cell(Value::Null), + }); + } + Some(out) + } +} + +fn float(v: f64) -> Value { + if v.is_finite() { Value::Float(v) } else { Value::Null } +} + +/// Days since 1899-12-30 as text in `fmt`. +fn date(days: f64, fmt: &str) -> Value { + // 1970-01-01 is day 25,569. + let ms = (days - 25_569.0) * 86_400_000.0; + if ms.is_finite() && ms.abs() < 1e16 { Value::Text(crate::time::format(ms, fmt)) } else { Value::Null } +} + +/// A GUID as `{XXXXXXXX-XXXX-XXXX-XXXX-XXXXXXXXXXXX}` (its first three groups little-endian). +fn guid(b: &[u8]) -> String { + let h = |r: &[usize]| r.iter().filter_map(|&i| b.get(i)).map(|x| format!("{x:02X}")).collect::(); + format!("{{{}-{}-{}-{}-{}}}", h(&[3, 2, 1, 0]), h(&[5, 4]), h(&[7, 6]), h(&[8, 9]), h(&[10, 11, 12, 13, 14, 15])) +} + +/// Where each row starts in its table, from the `.gdbtablx`: `(object ID, offset)` for every row +/// not deleted. +fn row_offsets(x: &[u8]) -> Result> { + let bad = |what: &str| DataError::Format(format!("unreadable row index ({what})")); + let mut c = Cur::new(x, 0); + let version = c.i32().ok_or_else(|| bad("header"))?; + let blocks = c.i32().and_then(|n| usize::try_from(n).ok()).ok_or_else(|| bad("header"))?; + let rows_v3 = c.i32().ok_or_else(|| bad("header"))?; + let width = c.i32().and_then(|n| usize::try_from(n).ok()).filter(|w| (4..=6).contains(w)).ok_or_else(|| bad("offset size"))?; + let stored = blocks.checked_mul(1024).ok_or_else(|| bad("block count"))?; + let trailer = stored.checked_mul(width).and_then(|n| n.checked_add(16)).filter(|&n| n <= x.len()).ok_or_else(|| bad("block count"))?; + let mut c = Cur::new(x, trailer); + // Which 1024-row blocks are stored, when some are left out; and how many rows there are. + let (bitmap, total): (Option<&[u8]>, usize) = match version { + 3 => { + let words = c.i32().and_then(|n| usize::try_from(n).ok()).unwrap_or(0); + let total = usize::try_from(rows_v3).map_err(|_| bad("row count"))?; + if words == 0 { + (None, total) + } else { + c.skip(12).ok_or_else(|| bad("bitmap"))?; + (Some(c.bytes(words.checked_mul(4).ok_or_else(|| bad("bitmap"))?).ok_or_else(|| bad("bitmap"))?), total) + } + } + 4 => { + let total = c.i64().and_then(|n| usize::try_from(n).ok()).ok_or_else(|| bad("row count"))?; + let varying = c.i32().ok_or_else(|| bad("trailer"))?; + if varying != 0 { + return Err(DataError::Unsupported("a row index with missing blocks".into())); + } + (None, total) + } + _ => return Err(DataError::Unsupported(format!("row index version {version}"))), + }; + // Block by block (1024 rows each): a block the bitmap leaves out is passed over whole, and no + // more blocks are looked at than are stored or the bitmap has bits for, whatever the row count + // in the header says. + let in_blocks = match bitmap { + None => total.min(stored).div_ceil(1024), + Some(bits) => total.div_ceil(1024).min(bits.len().saturating_mul(8)), + }; + let mut out = Vec::new(); + let mut stored_before = 0usize; + for block in 0..in_blocks { + let present = bitmap.is_none_or(|bits| bits.get(block / 8).is_some_and(|b| b >> (block % 8) & 1 == 1)); + if !present { + continue; + } + if stored_before >= blocks { + break; + } + let first = stored_before * 1024; + stored_before += 1; + for k in 0..1024 { + let row = block * 1024 + k; + if row >= total { + break; + } + let at = 16 + (first + k) * width; + let Some(bytes) = x.get(at..at + width) else { break }; + let offset = bytes.iter().rev().fold(0u64, |v, &b| v << 8 | u64::from(b)); + if offset == 0 { + continue; + } + if let Ok(offset) = usize::try_from(offset) { + out.push((row as u64 + 1, offset)); + } + } + } + Ok(out) +} + +// --------------------------------------------------------------------------------------------- +// Geometry + +enum Decoded { + /// A geometry, and how many curved segments it draws straight. + Geometry(Geometry, usize), + Patch(Patch), + Empty, +} + +/// A multipatch seen from above. +struct Patch { + footprint: Option, + /// Its lowest and highest Z. + z: Option<(f64, f64)>, + /// Why the footprint is the hull of its points, when it is. + hull: Option, +} + +/// Why a multipatch's footprint is the hull of its points. +#[derive(Clone, Copy, Debug, PartialEq)] +enum Hull { + /// No face shows from above (only walls). + NoFaces, + /// Its faces weren't merged: too many, past the table's budget, or the union failed. + Unmerged, +} + +const HAS_Z: u64 = 0x8000_0000; +const HAS_CURVES: u64 = 0x2000_0000; + +/// `budget` is how many more faces the table's multipatches may merge (see [`MAX_TABLE_FACES`]). +fn decode_geometry(b: &[u8], def: &GeomDef, budget: &mut usize) -> std::result::Result { + let mut c = Cur::new(b, 0); + let ty = c.varuint().ok_or(())?; + let x = |v: f64| v / def.xy_scale + def.x0; + let y = |v: f64| v / def.xy_scale + def.y0; + match ty & 0xFF { + 0 => Ok(Decoded::Empty), + // Points (2D, Z, M, ZM, and the general type). + 1 | 9 | 11 | 21 | 52 => { + let (vx, vy) = (c.varuint().ok_or(())?, c.varuint().ok_or(())?); + if vx == 0 || vy == 0 { + return Ok(Decoded::Empty); + } + Ok(Decoded::Geometry(Geometry::Point(Coord::new(x((vx - 1) as f64), y((vy - 1) as f64))), 0)) + } + // Multipoints. + 8 | 18 | 20 | 28 | 53 => { + let n = c.count(MAX_PARTS, 2).ok_or(())?; + for _ in 0..4 { + c.varuint().ok_or(())?; + } + let pts = deltas(&mut c, n, def)?; + if pts.is_empty() { Ok(Decoded::Empty) } else { Ok(Decoded::Geometry(Geometry::MultiPoint(pts), 0)) } + } + // Lines and polygons (and the general types, with curves). + code @ (3 | 10 | 13 | 23 | 50 | 5 | 15 | 19 | 25 | 51) => { + let general = code == 50 || code == 51; + let polygon = matches!(code, 5 | 15 | 19 | 25 | 51); + let n = c.count(MAX_PARTS, 2).ok_or(())?; + let parts = c.count(MAX_PARTS, 1).ok_or(())?; + let curves = if general && ty & HAS_CURVES != 0 { c.count(MAX_PARTS, 1).ok_or(())? } else { 0 }; + for _ in 0..4 { + c.varuint().ok_or(())?; + } + let sizes = part_sizes(&mut c, n, parts)?; + let pts = deltas(&mut c, n, def)?; + let rings = split(&pts, &sizes); + if rings.is_empty() { + return Ok(Decoded::Empty); + } + let geom = if polygon { + // Clockwise rings are outer; each counterclockwise one is a hole in the one before. + // A ring with no area (its points on a line) is neither, and is left out. + let rings: Vec<(Vec, bool)> = rings + .into_iter() + .filter_map(|r| { + let a = signed_area(&r); + (a != 0.0 && a.is_finite()).then_some((r, a < 0.0)) + }) + .collect(); + if rings.is_empty() { + return Ok(Decoded::Empty); + } + crate::io::group_rings(rings) + } else { + Geometry::Line(rings) + }; + Ok(Decoded::Geometry(geom, curves)) + } + // Multipatches: 3D surfaces (the general type, or the Shapefile ones GDAL writes, which + // always have Z). + code @ (31 | 32 | 54) => { + let has_z = code != 54 || ty & HAS_Z != 0; + let n = c.count(MAX_PARTS, 2).ok_or(())?; + c.varuint().ok_or(())?; + let parts = c.count(MAX_PARTS, 1).ok_or(())?; + for _ in 0..4 { + c.varuint().ok_or(())?; + } + let sizes = part_sizes(&mut c, n, parts)?; + let mut kinds = Vec::with_capacity(parts); + for _ in 0..parts { + kinds.push(c.varuint().ok_or(())? & 0xF); + } + let pts = deltas(&mut c, n, def)?; + // Heights, when stored (a multipatch without them lies flat). + let mut zs = Vec::with_capacity(if has_z { n } else { 0 }); + if has_z { + let mut dz: i64 = 0; + for _ in 0..n { + dz = dz.checked_add(c.varint().ok_or(())?).ok_or(())?; + zs.push(dz as f64 / def.z_scale + def.z0); + } + } + let z = zs.iter().fold(None, |r: Option<(f64, f64)>, &v| match r { + _ if !v.is_finite() => r, + None => Some((v, v)), + Some((lo, hi)) => Some((lo.min(v), hi.max(v))), + }); + let pts: Vec = pts.iter().enumerate().map(|(i, p)| [p.x, p.y, zs.get(i).copied().unwrap_or(0.0)]).collect(); + let (footprint, hull) = footprint(&split(&pts, &sizes), &kinds, budget); + Ok(Decoded::Patch(Patch { footprint, z, hull })) + } + _ => Err(()), + } +} + +/// The point counts of `parts` parts holding `n` points: all but the last are stored. +fn part_sizes(c: &mut Cur, n: usize, parts: usize) -> std::result::Result, ()> { + if parts == 0 { + return Ok(Vec::new()); + } + let mut sizes = Vec::with_capacity(parts); + let mut sum = 0usize; + for _ in 1..parts { + let s = c.count(MAX_PARTS, 1).ok_or(())?; + sum = sum.checked_add(s).ok_or(())?; + sizes.push(s); + } + sizes.push(n.checked_sub(sum).ok_or(())?); + Ok(sizes) +} + +/// `n` points stored as deltas from the one before (across parts). +fn deltas(c: &mut Cur, n: usize, def: &GeomDef) -> std::result::Result, ()> { + let (mut dx, mut dy) = (0i64, 0i64); + let mut pts = Vec::with_capacity(n); + for _ in 0..n { + dx = dx.checked_add(c.varint().ok_or(())?).ok_or(())?; + dy = dy.checked_add(c.varint().ok_or(())?).ok_or(())?; + pts.push(Coord::new(dx as f64 / def.xy_scale + def.x0, dy as f64 / def.xy_scale + def.y0)); + } + Ok(pts) +} + +/// Points cut into parts of the given sizes (sizes past the points are cut short). +fn split(pts: &[T], sizes: &[usize]) -> Vec> { + let mut out = Vec::with_capacity(sizes.len()); + let mut at = 0usize; + for &s in sizes { + let end = at.saturating_add(s).min(pts.len()); + if let Some(part) = pts.get(at..end).filter(|p| !p.is_empty()) { + out.push(part.to_vec()); + } + at = end; + } + out +} + +/// Twice the signed area of a ring (positive counterclockwise). +fn signed_area(ring: &[Coord]) -> f64 { + // From the first point: products of projected coordinates in the millions (or more) would round + // a small ring's area away. + let o = ring.first().copied().unwrap_or(Coord::new(0.0, 0.0)); + let local = |c: &Coord| (c.x - o.x, c.y - o.y); + ring.iter() + .zip(ring.iter().cycle().skip(1)) + .map(|(a, b)| { + let ((ax, ay), (bx, by)) = (local(a), local(b)); + ax * by - bx * ay + }) + .sum() +} + +/// A point with its height. +type P3 = [f64; 3]; + +/// A multipatch's faces, each its outer ring and then its holes: triangles from strips, fans and +/// triangle lists, and rings. At most [`MAX_FACES`]; the flag says some were left out. +fn faces(parts: &[Vec], kinds: &[u64]) -> (Vec>>, bool) { + let mut faces: Vec>> = Vec::new(); + let mut truncated = false; + for (part, &kind) in parts.iter().zip(kinds) { + // Faces past the cap are left out, within a part too (one strip can hold millions). + let room = MAX_FACES.saturating_sub(faces.len()); + let has = match kind { + 0 | 1 => part.len().saturating_sub(2), + 6 => part.len() / 3, + 3 => 0, + _ => 1, + }; + if has > room { + truncated = true; + } + if room == 0 { + break; + } + match kind { + // Triangle strip, fan, triangles. + 0 => faces.extend(part.windows(3).take(room).map(|w| vec![w.to_vec()])), + 1 => { + if let Some((first, rest)) = part.split_first() { + faces.extend( + rest.windows(2) + .filter_map(|w| match w { + [a, b] => Some(vec![vec![*first, *a, *b]]), + _ => None, + }) + .take(room), + ); + } + } + 6 => faces.extend(part.as_chunks::<3>().0.iter().take(room).map(|w| vec![w.to_vec()])), + // An inner ring is a hole in the face before. + 3 => match faces.last_mut() { + Some(face) => face.push(part.clone()), + None => faces.push(vec![part.clone()]), + }, + // Outer, first and plain rings. + _ => faces.push(vec![part.clone()]), + } + } + (faces, truncated) +} + +/// A multipatch's footprint: the union of its faces seen from above. Faces standing on edge, such +/// as walls, show nothing and are left out. The footprint is the hull of its points when no face +/// shows, or when its faces aren't merged: more than [`MAX_FACES`], more seen from above than +/// [`MAX_UNION_FACES`] or than the table's `budget` has left (which they then take from it). +fn footprint(parts: &[Vec], kinds: &[u64], budget: &mut usize) -> (Option, Option) { + use geo::{ConvexHull, Orient, orient::Direction}; + // Where no union runs (the browser), the faces aren't built at all: the hull stands in. + let unions = !cfg!(target_arch = "wasm32"); + let (faces, truncated) = if unions { faces(parts, kinds) } else { (Vec::new(), false) }; + let all = || parts.iter().flatten(); + let (x0, x1) = all().fold((f64::INFINITY, f64::NEG_INFINITY), |(a, b), p| (a.min(p[0]), b.max(p[0]))); + let (y0, y1) = all().fold((f64::INFINITY, f64::NEG_INFINITY), |(a, b), p| (a.min(p[1]), b.max(p[1]))); + let size = (x1 - x0).max(y1 - y0); + if !size.is_finite() { + return (None, None); + } + // A face shows from above by the cosine of its tilt (its projected area over its area): walls, + // even leaning ones, show as slivers, which cover nothing and which the union can't always + // take. Faces narrower than a ten-thousandth of the patch are left out too (a flat patch). + let shows = |ring: &[P3]| { + let (mut nx, mut ny, mut nz) = (0.0, 0.0, 0.0); + let (mut longest, mut area2) = (0.0f64, 0.0); + // From the ring's first point: projected coordinates run to millions, whose products would + // lose the millimetres a small face is measured in. + let o = ring.first().copied().unwrap_or([0.0; 3]); + let local: Vec = ring.iter().map(|p| [p[0] - o[0], p[1] - o[1], p[2] - o[2]]).collect(); + for (a, b) in local.iter().zip(local.iter().cycle().skip(1)) { + nx += (a[1] - b[1]) * (a[2] + b[2]); + ny += (a[2] - b[2]) * (a[0] + b[0]); + nz += (a[0] - b[0]) * (a[1] + b[1]); + longest = longest.max((b[0] - a[0]).hypot(b[1] - a[1])); + area2 += a[0] * b[1] - b[0] * a[1]; + } + let n = (nx * nx + ny * ny + nz * nz).sqrt(); + n > 0.0 && nz.abs() / n > 0.05 && longest > 0.0 && area2.abs() / longest > size * 1e-4 + }; + let ring2 = |r: &[P3]| geo::LineString::from(r.iter().map(|p| (p[0], p[1])).collect::>()); + // Past the cap, a partial union would be wrong: the hull below stands in. + let seen: Vec> = faces + .iter() + .filter(|_| !truncated) + .filter_map(|face| { + let (outer, holes) = face.split_first()?; + (outer.len() >= 3 && shows(outer)) + .then(|| geo::Polygon::new(ring2(outer), holes.iter().map(|h| ring2(h)).collect()).orient(Direction::Default)) + }) + .collect(); + let merge = unions && !truncated && seen.len() <= MAX_UNION_FACES && seen.len() <= *budget; + // The union library may still trip over faces it can't place: the hull stands in then. (In the + // browser, where a panic can't be caught, nothing is merged: see `MAX_UNION_FACES`.) + let union = if seen.is_empty() || !merge { + None + } else { + // Spent before the union runs, so faces whose union panics (and become the hull) count too: + // the budget bounds the work tried, not the footprints made. + *budget -= seen.len(); + std::panic::catch_unwind(std::panic::AssertUnwindSafe(|| geo::unary_union(seen.iter()))).ok().filter(|u| !u.0.is_empty()) + }; + let (shape, hull) = match union { + Some(u) => (u, None), + None => { + let points: Vec> = all().map(|p| geo::Point::new(p[0], p[1])).collect(); + if points.len() < 3 { + return (None, None); + } + let why = if unions && seen.is_empty() && !truncated { Hull::NoFaces } else { Hull::Unmerged }; + (geo::MultiPolygon::new(vec![geo::MultiPoint::from(points).convex_hull()]), Some(why)) + } + }; + let ring = |ls: &geo::LineString| ls.0.iter().map(|c| Coord::new(c.x, c.y)).collect::>(); + let polys: Vec>> = shape + .0 + .iter() + .filter(|p| p.exterior().0.len() >= 4) + .map(|p| std::iter::once(ring(p.exterior())).chain(p.interiors().iter().map(ring)).collect()) + .collect(); + if polys.is_empty() { (None, None) } else { (Some(Geometry::Polygon(polys)), hull) } +} + +#[cfg(test)] +mod tests; diff --git a/crates/data/src/filegdb/tests.rs b/crates/data/src/filegdb/tests.rs new file mode 100644 index 000000000..f285b190e --- /dev/null +++ b/crates/data/src/filegdb/tests.rs @@ -0,0 +1,625 @@ +use super::*; + +fn varuint(mut v: u64) -> Vec { + let mut out = Vec::new(); + loop { + let b = (v & 0x7F) as u8; + v >>= 7; + if v == 0 { + out.push(b); + return out; + } + out.push(b | 0x80); + } +} + +fn varint(v: i64) -> Vec { + let (neg, mut m) = (v < 0, v.unsigned_abs()); + let mut first = (m & 0x3F) as u8 | if neg { 0x40 } else { 0 }; + m >>= 6; + if m == 0 { + return vec![first]; + } + first |= 0x80; + let mut out = vec![first]; + out.extend(varuint(m)); + out +} + +fn utf16_bytes(s: &str) -> Vec { + s.encode_utf16().flat_map(u16::to_le_bytes).collect() +} + +/// A field for the test writer. +enum F { + Oid, + Text(&'static str), + Int(&'static str), + Shape { z: bool, wkt: &'static str }, +} + +/// One row's stored values (the object ID stores nothing): `None` is null. +enum V { + Text(Option<&'static str>), + Int(Option), + Blob(Vec), +} + +/// A minimal table (version 3) and its row index, as ArcGIS lays them out. +fn table(layer_flags: u32, fields: &[F], rows: &[Vec]) -> (Vec, Vec) { + let mut f = Vec::new(); + f.extend(4i32.to_le_bytes()); + f.extend((layer_flags | 0x100).to_le_bytes()); + f.extend((fields.len() as i16).to_le_bytes()); + let name = |f: &mut Vec, n: &str| { + f.push(n.encode_utf16().count() as u8); + f.extend(utf16_bytes(n)); + f.push(0); + }; + for field in fields { + match field { + F::Oid => { + name(&mut f, "OBJECTID"); + f.extend([6, 4, 2]); + } + F::Text(n) => { + name(&mut f, n); + f.push(4); + f.extend(255i32.to_le_bytes()); + f.extend([5, 0]); + } + F::Int(n) => { + name(&mut f, n); + f.extend([1, 4, 5, 0]); + } + F::Shape { z, wkt } => { + name(&mut f, "SHAPE"); + f.extend([7, 0, 7]); + let w = utf16_bytes(wkt); + f.extend((w.len() as i16).to_le_bytes()); + f.extend(w); + f.push(if *z { 3 } else { 1 }); + for v in [0.0f64, 0.0, 1000.0] { + f.extend(v.to_le_bytes()); + } + if *z { + f.extend(0.0f64.to_le_bytes()); + f.extend(1000.0f64.to_le_bytes()); + } + f.extend(0.001f64.to_le_bytes()); + if *z { + f.extend(0.001f64.to_le_bytes()); + } + for v in [0.0f64, 0.0, 100.0, 100.0] { + f.extend(v.to_le_bytes()); + } + if layer_flags & 0x8000_0000 != 0 { + f.extend([0.0f64, 100.0].iter().flat_map(|v| v.to_le_bytes())); + } + f.push(0); + f.extend(1u32.to_le_bytes()); + f.extend(10.0f64.to_le_bytes()); + } + } + } + // Every field but the object ID can be null here (the shape's flag 7 says so), so each has a + // bit in the null flags. + let nullable = fields.iter().filter(|f| !matches!(f, F::Oid)).count(); + let mut t = vec![0u8; 40]; + t.extend((f.len() as i32).to_le_bytes()); + t.extend(&f); + let mut offsets = Vec::new(); + for row in rows { + let mut nulls = vec![0u8; nullable.div_ceil(8)]; + let mut body = Vec::new(); + for (bit, v) in row.iter().enumerate() { + match v { + V::Text(None) | V::Int(None) => { + if let Some(b) = nulls.get_mut(bit / 8) { + *b |= 1 << (bit % 8); + } + } + V::Text(Some(s)) => { + body.extend(varuint(s.len() as u64)); + body.extend(s.as_bytes()); + } + V::Int(Some(i)) => body.extend(i.to_le_bytes()), + V::Blob(b) => { + body.extend(varuint(b.len() as u64)); + body.extend(b); + } + } + } + offsets.push(t.len() as u32); + t.extend(((nulls.len() + body.len()) as i32).to_le_bytes()); + t.extend(nulls); + t.extend(body); + } + let size = t.len() as i64; + t[0..4].copy_from_slice(&3i32.to_le_bytes()); + t[4..8].copy_from_slice(&(rows.len() as i32).to_le_bytes()); + t[12..16].copy_from_slice(&5i32.to_le_bytes()); + t[24..32].copy_from_slice(&size.to_le_bytes()); + t[32..40].copy_from_slice(&40i64.to_le_bytes()); + let mut x = Vec::new(); + for v in [3i32, 1, rows.len() as i32, 4] { + x.extend(v.to_le_bytes()); + } + for i in 0..1024 { + x.extend(offsets.get(i).copied().unwrap_or(0).to_le_bytes()); + } + x.extend([0u8; 16]); + (t, x) +} + +/// A polygon blob (type 5): one ring, scale 1000 from 0. +fn polygon(ring: &[(f64, f64)]) -> Vec { + let mut b = varuint(5); + b.extend(varuint(ring.len() as u64)); + b.extend(varuint(1)); + b.extend([0u8; 4]); + let (mut px, mut py) = (0i64, 0i64); + for &(x, y) in ring { + let (ix, iy) = ((x * 1000.0).round() as i64, (y * 1000.0).round() as i64); + b.extend(varint(ix - px)); + b.extend(varint(iy - py)); + (px, py) = (ix, iy); + } + b +} + +/// A multipatch blob (type 54 with Z): its parts' types and points. +fn multipatch(parts: &[(u64, Vec<(f64, f64, f64)>)]) -> Vec { + let n: usize = parts.iter().map(|(_, p)| p.len()).sum(); + let mut b = varuint(54 | HAS_Z); + b.extend(varuint(n as u64)); + b.extend(varuint(0)); + b.extend(varuint(parts.len() as u64)); + b.extend([0u8; 4]); + for (_, p) in parts.iter().take(parts.len().saturating_sub(1)) { + b.extend(varuint(p.len() as u64)); + } + for (kind, _) in parts { + b.extend(varuint(*kind)); + } + let (mut px, mut py) = (0i64, 0i64); + for (_, p) in parts { + for &(x, y, _) in p { + let (ix, iy) = ((x * 1000.0).round() as i64, (y * 1000.0).round() as i64); + b.extend(varint(ix - px)); + b.extend(varint(iy - py)); + (px, py) = (ix, iy); + } + } + let mut pz = 0i64; + for (_, p) in parts { + for &(_, _, z) in p { + let iz = (z * 1000.0).round() as i64; + b.extend(varint(iz - pz)); + pz = iz; + } + } + b +} + +const LONLAT: &str = r#"GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]"#; + +/// A geodatabase's files: a catalog naming the given tables, which follow as a00000002…. +fn geodatabase(tables: Vec<(&'static str, (Vec, Vec))>) -> Vec<(String, Vec)> { + let rows: Vec> = std::iter::once("GDB_SystemCatalog") + .chain(tables.iter().map(|(n, _)| *n)) + .map(|n| vec![V::Text(Some(n)), V::Int(Some(0))]) + .collect(); + let (t, x) = table(0, &[F::Oid, F::Text("Name"), F::Int("FileFormat")], &rows); + let mut files = vec![("a00000001.gdbtable".to_string(), t), ("a00000001.gdbtablx".to_string(), x)]; + for (i, (_, (t, x))) in tables.into_iter().enumerate() { + files.push((format!("a{:08x}.gdbtable", i + 2), t)); + files.push((format!("a{:08x}.gdbtablx", i + 2), x)); + } + files +} + +fn read_files(files: &[(String, Vec)]) -> Result> { + read("test", &|name: &str| Ok(files.iter().find(|(n, _)| n == name).map(|(_, b)| b.clone())), None).map(|r| r.layers) +} + +/// A square (clockwise, as ArcGIS stores outer rings) from (x, y) with side `s`. +fn square(x: f64, y: f64, s: f64) -> Vec<(f64, f64)> { + vec![(x, y), (x, y + s), (x + s, y + s), (x + s, y), (x, y)] +} + +#[test] +fn numbers_decode() { + for v in [0u64, 1, 63, 64, 127, 128, 300, 1 << 40, u64::MAX >> 1] { + assert_eq!(Cur::new(&varuint(v), 0).varuint(), Some(v), "{v}"); + } + for v in [0i64, 1, -1, 63, -63, 64, -64, -65, 1_000_000, -1_000_000, 1 << 50] { + assert_eq!(Cur::new(&varint(v), 0).varint(), Some(v), "{v}"); + } + // Endless continuation bits and cut-short numbers are refused. + assert_eq!(Cur::new(&[0xFF; 12], 0).varuint(), None); + assert_eq!(Cur::new(&[0xFF; 12], 0).varint(), None); + assert_eq!(Cur::new(&[0x80], 0).varuint(), None); + assert_eq!(guid(&[0x67, 0x45, 0x23, 0x01, 0xAB, 0x89, 0xEF, 0xCD, 1, 2, 3, 4, 5, 6, 7, 8]), "{01234567-89AB-CDEF-0102-030405060708}"); + // 2000-01-01 is day 36,526 after 1899-12-30. + assert_eq!(date(36_526.5, "YYYY-MM-DDTHH:mm:ss"), Value::Text("2000-01-01T12:00:00".into())); +} + +#[test] +fn reads_a_feature_class_with_fields_nulls_and_holes() { + let outer = square(0.0, 0.0, 10.0); + // A hole: counterclockwise, in the same blob as a second part. + let hole: Vec<(f64, f64)> = square(2.0, 2.0, 2.0).into_iter().rev().collect(); + let mut blob = varuint(5); + blob.extend(varuint((outer.len() + hole.len()) as u64)); + blob.extend(varuint(2)); + blob.extend([0u8; 4]); + blob.extend(varuint(outer.len() as u64)); + let (mut px, mut py) = (0i64, 0i64); + for &(x, y) in outer.iter().chain(&hole) { + let (ix, iy) = ((x * 1000.0) as i64, (y * 1000.0) as i64); + blob.extend(varint(ix - px)); + blob.extend(varint(iy - py)); + (px, py) = (ix, iy); + } + let fields = [F::Oid, F::Shape { z: false, wkt: LONLAT }, F::Text("name"), F::Int("count")]; + let rows = vec![ + vec![V::Blob(blob), V::Text(Some("Gebäude")), V::Int(Some(7))], + vec![V::Blob(polygon(&square(20.0, 20.0, 5.0))), V::Text(None), V::Int(None)], + ]; + let files = geodatabase(vec![("Parcels", table(4, &fields, &rows))]); + let layers = read_files(&files).unwrap(); + assert_eq!(layers.len(), 1); + let (ds, report) = &layers[0]; + assert_eq!(ds.name, "Parcels"); + assert_eq!(ds.geometry_type, GeometryType::Polygon); + assert_eq!(report.skipped, 0, "{report:?}"); + let names: Vec<&str> = ds.fields.iter().map(|f| f.name.as_str()).collect(); + assert_eq!(names, ["OBJECTID", "name", "count"]); + let rows: Vec<&crate::Feature> = ds.features.iter().collect(); + assert_eq!(rows[0].attrs, vec![Value::Int(1), Value::Text("Gebäude".into()), Value::Int(7)]); + assert_eq!(rows[1].attrs, vec![Value::Int(2), Value::Null, Value::Null]); + match &rows[0].geometry { + Some(Geometry::Polygon(p)) => { + assert_eq!(p.len(), 1, "one polygon"); + assert_eq!(p[0].len(), 2, "with its hole"); + assert_eq!(p[0][0][2], Coord::new(10.0, 10.0)); + } + g => panic!("{g:?}"), + } +} + +#[test] +fn a_multipatch_is_its_footprint_with_its_heights() { + // A 10 m box from z 400 to 410: its roof as two triangles (a strip), a floor ring, and walls, + // which show nothing from above. + let roof = (0, vec![(0.0, 0.0, 410.0), (0.0, 10.0, 410.0), (10.0, 0.0, 410.0), (10.0, 10.0, 410.0)]); + let floor = (2, vec![(0.0, 0.0, 400.0), (0.0, 10.0, 400.0), (10.0, 10.0, 400.0), (10.0, 0.0, 400.0), (0.0, 0.0, 400.0)]); + let wall = (2, vec![(0.0, 0.0, 400.0), (0.0, 0.0, 410.0), (10.0, 0.0, 410.0), (10.0, 0.0, 400.0), (0.0, 0.0, 400.0)]); + // A leaning wall: a sliver from above, which is left out too. + let leaning = (6, vec![(0.0, 10.0, 400.0), (5.0, 10.003, 410.0), (10.0, 10.0, 400.0)]); + // Only walls: the outline of its points stands in. + let walls_only = vec![ + (2, vec![(30.0, 30.0, 0.0), (30.0, 30.0, 5.0), (40.0, 30.0, 5.0), (40.0, 30.0, 0.0), (30.0, 30.0, 0.0)]), + (2, vec![(40.0, 30.0, 0.0), (40.0, 30.0, 5.0), (40.0, 40.0, 5.0), (40.0, 40.0, 0.0), (40.0, 30.0, 0.0)]), + ]; + let fields = [F::Oid, F::Shape { z: true, wkt: LONLAT }, F::Text("name")]; + let rows = vec![ + vec![V::Blob(multipatch(&[roof, floor, wall, leaning])), V::Text(Some("box"))], + vec![V::Blob(multipatch(&walls_only)), V::Text(Some("walls"))], + ]; + let files = geodatabase(vec![("Buildings", table(9 | 0x8000_0000, &fields, &rows))]); + let layers = read_files(&files).unwrap(); + let (ds, report) = &layers[0]; + assert_eq!(ds.geometry_type, GeometryType::Polygon); + assert_eq!(ds.fields.iter().map(|f| f.name.as_str()).collect::>(), ["OBJECTID", "name", "z_min", "z_max", "height"]); + let rows: Vec<&crate::Feature> = ds.features.iter().collect(); + assert_eq!(rows[0].attrs[2..], [Value::Float(400.0), Value::Float(410.0), Value::Float(10.0)]); + let area = |g: &Option| match g { + Some(Geometry::Polygon(p)) => p.iter().map(|poly| poly.first().map_or(0.0, |r| signed_area(r).abs() / 2.0)).sum::(), + _ => 0.0, + }; + assert!((area(&rows[0].geometry) - 100.0).abs() < 1e-6, "{:?}", rows[0].geometry); + assert!((area(&rows[1].geometry) - 50.0).abs() < 1e-6, "the hull of an L of walls: {:?}", rows[1].geometry); + assert!(report.messages.iter().any(|m| m.contains("footprints")), "{report:?}"); + assert!(report.messages.iter().any(|m| m.contains("1 multipatches with no faces")), "{report:?}"); +} + +#[test] +fn row_index_with_left_out_blocks() { + // Rows 1025 and 1026 stored in the one block present (the second of two). + let mut x = Vec::new(); + for v in [3i32, 1, 2048, 4] { + x.extend(v.to_le_bytes()); + } + for i in 0..1024u32 { + x.extend( + match i { + 0 => 100u32, + 1 => 0, + 2 => 200, + _ => 0, + } + .to_le_bytes(), + ); + } + for v in [1i32, 2, 1, 1] { + x.extend(v.to_le_bytes()); + } + x.extend(0b10u32.to_le_bytes()); + assert_eq!(row_offsets(&x).unwrap(), vec![(1025, 100), (1027, 200)]); + // A version 4 index with missing blocks isn't read here. + let mut x = Vec::new(); + for v in [4i32, 0, 0, 4] { + x.extend(v.to_le_bytes()); + } + x.extend(10i64.to_le_bytes()); + x.extend(30i32.to_le_bytes()); + assert!(matches!(row_offsets(&x), Err(DataError::Unsupported(_)))); +} + +/// A row index claiming billions of rows in blocks its bitmap leaves out is passed over block by +/// block, not row by row (it hung before). +#[test] +fn a_huge_row_count_with_absent_blocks_is_quick() { + let mut x = Vec::new(); + for v in [3i32, 0, i32::MAX, 4] { + x.extend(v.to_le_bytes()); + } + for v in [1i32, (i32::MAX / 1024) + 1, 0, 1] { + x.extend(v.to_le_bytes()); + } + x.extend(0u32.to_le_bytes()); + let start = std::time::Instant::now(); + assert_eq!(row_offsets(&x).unwrap(), Vec::new()); + assert!(start.elapsed() < std::time::Duration::from_secs(1), "{:?}", start.elapsed()); +} + +/// One strip of a million points makes no more faces than the cap (each is an allocation). +#[test] +fn a_huge_strip_stops_at_the_face_cap() { + let strip: Vec = (0..1_000_000).map(|i| [f64::from(i), f64::from(i % 2), 0.0]).collect(); + let n = strip.len(); + for kind in [0, 1, 6] { + let (f, truncated) = faces(std::slice::from_ref(&strip), &[kind]); + assert_eq!((f.len(), truncated), (MAX_FACES, true), "kind {kind}"); + } + // Later parts get only the room left. + assert_eq!(faces(&[strip[..10].to_vec(), strip.clone()], &[0, 0]).0.len(), MAX_FACES); + assert!(n / 3 > MAX_FACES); + // Within the cap, nothing is left out. + assert_eq!(faces(&[strip[..10].to_vec()], &[0]), (strip[..10].windows(3).map(|w| vec![w.to_vec()]).collect(), false)); + // A truncated patch is drawn as its points' hull, not as part of its faces. + let mut budget = MAX_TABLE_FACES; + let (shape, hull) = footprint(std::slice::from_ref(&strip), &[0], &mut budget); + assert!(hull == Some(Hull::Unmerged) && shape.is_some()); + assert_eq!(budget, MAX_TABLE_FACES, "nothing was merged"); +} + +/// A ring with no area (its points on a line) is neither an outer ring nor a hole: no sliver. +#[test] +fn a_ring_with_no_area_is_left_out() { + let fields = [F::Oid, F::Shape { z: false, wkt: LONLAT }]; + let line = [(0.0, 0.0), (5.0, 0.0), (10.0, 0.0), (0.0, 0.0)]; + let rows = vec![vec![V::Blob(polygon(&line))], vec![V::Blob(polygon(&square(0.0, 0.0, 10.0)))]]; + let files = geodatabase(vec![("Flat", table(4, &fields, &rows))]); + let layers = read_files(&files).unwrap(); + let (ds, report) = &layers[0]; + let rows: Vec<&crate::Feature> = ds.features.iter().collect(); + assert!(!matches!(rows.first().and_then(|f| f.geometry.as_ref()), Some(Geometry::Polygon(_))), "{report:?}"); + assert!(matches!(rows.get(1).and_then(|f| f.geometry.as_ref()), Some(Geometry::Polygon(p)) if p.len() == 1)); +} + +/// A small ring far from the origin keeps its area: products of its coordinates would round it away. +#[test] +fn a_small_ring_far_from_the_origin_has_an_area() { + let ring: Vec = + [(0.0, 0.0), (0.0, 1.0), (1.0, 1.0), (1.0, 0.0), (0.0, 0.0)].iter().map(|&(x, y)| Coord::new(1e9 + x, 1e9 + y)).collect(); + assert_eq!(signed_area(&ring), -2.0, "clockwise, twice its area"); +} + +/// A table's multipatches merge no more faces than its budget: past it, they are hulls. +#[test] +fn the_face_budget_bounds_a_tables_unions() { + // A 10 m square roof as a strip of two triangles, and a wall. + let parts = vec![ + vec![[0.0, 0.0, 5.0], [0.0, 10.0, 5.0], [10.0, 0.0, 5.0], [10.0, 10.0, 5.0]], + vec![[0.0, 0.0, 0.0], [0.0, 0.0, 5.0], [10.0, 0.0, 5.0], [10.0, 0.0, 0.0], [0.0, 0.0, 0.0]], + ]; + let kinds = [0, 2]; + let mut budget = 3; + let (shape, hull) = footprint(&parts, &kinds, &mut budget); + assert_eq!((hull, budget), (None, 1), "two faces show and are merged: {shape:?}"); + let (shape, hull) = footprint(&parts, &kinds, &mut budget); + assert_eq!((hull, budget), (Some(Hull::Unmerged), 1), "past the budget: {shape:?}"); + assert!(shape.is_some()); +} + +/// Numbers whose bits run past 64 are refused: one decoding to `i64::MIN` panicked when negated. +#[test] +fn numbers_past_64_bits_are_refused() { + // Sign set, magnitude 2^63: eight continuation bytes reach bit 62, then 0x02. + let mut b = vec![0xC0]; + b.extend([0x80; 8]); + b.push(0x02); + assert_eq!(Cur::new(&b, 0).varint(), None); + // The largest magnitude still reads. + let mut b = vec![0xFF]; + b.extend([0xFF; 8]); + b.push(0x01); + assert_eq!(Cur::new(&b, 0).varint(), Some(-i64::MAX)); + // Unsigned: bit 63 reads, a bit past it doesn't. + let mut b = vec![0x80; 9]; + b.push(0x01); + assert_eq!(Cur::new(&b, 0).varuint(), Some(1 << 63)); + let last = b.len() - 1; + b[last] = 0x02; + assert_eq!(Cur::new(&b, 0).varuint(), None); +} + +/// A table that can't be opened (too large, not readable) is said by name; the others are read. +#[test] +fn an_unreadable_table_leaves_the_others() { + let fields = [F::Oid, F::Text("name")]; + let rows = vec![vec![V::Text(Some("a"))]]; + let files = geodatabase(vec![("Good", table(0, &fields, &rows)), ("Huge", table(0, &fields, &rows))]); + let open = |name: &str| -> Result>> { + if name == "a00000003.gdbtable" { + return Err(DataError::Unsupported(format!("{name} is too large to read here"))); + } + Ok(files.iter().find(|(n, _)| n == name).map(|(_, b)| b.clone())) + }; + let read = read("test", &open, None).unwrap(); + assert_eq!(read.layers.iter().map(|(d, _)| d.name.as_str()).collect::>(), ["Good"]); + assert_eq!(read.unread.len(), 1); + assert_eq!(read.unread[0].0, "Huge"); + assert!(read.unread[0].1.contains("too large"), "{:?}", read.unread); +} + +/// A table the catalog lists whose files are missing is said by name, not passed over. +#[test] +fn a_table_with_missing_files_is_said() { + let fields = [F::Oid, F::Text("name")]; + let rows = vec![vec![V::Text(Some("a"))]]; + let mut files = geodatabase(vec![("Good", table(0, &fields, &rows)), ("Gone", table(0, &fields, &rows))]); + files.retain(|(n, _)| n != "a00000003.gdbtablx"); + let open = |name: &str| -> Result>> { Ok(files.iter().find(|(n, _)| n == name).map(|(_, b)| b.clone())) }; + let read = read("test", &open, None).unwrap(); + assert_eq!(read.layers.len(), 1); + assert_eq!(read.unread.len(), 1, "{:?}", read.unread); + assert_eq!(read.unread[0].0, "Gone"); + assert_eq!(read.unread[0].1, "its file a00000003.gdbtablx is missing"); +} + +/// A feature whose coordinates its coordinate system can't place keeps its row, and is said. +#[test] +fn coordinates_that_cant_be_converted_are_said() { + const LAEA: &str = r#"PROJCS["ETRS89 / LAEA Europe",GEOGCS["ETRS89",DATUM["European_Terrestrial_Reference_System_1989",SPHEROID["GRS 1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Lambert_Azimuthal_Equal_Area"],PARAMETER["latitude_of_center",52],PARAMETER["longitude_of_center",10],PARAMETER["false_easting",4321000],PARAMETER["false_northing",3210000],UNIT["metre",1],AUTHORITY["EPSG","3035"]]"#; + let fields = [F::Oid, F::Shape { z: false, wkt: LAEA }, F::Text("name")]; + let rows = vec![ + vec![V::Blob(polygon(&square(4_321_000.0, 3_210_000.0, 10.0))), V::Text(Some("in"))], + // Far past the far side of the Earth. + vec![V::Blob(polygon(&square(9e8, 9e8, 10.0))), V::Text(Some("out"))], + ]; + let files = geodatabase(vec![("T", table(4, &fields, &rows))]); + let layers = read_files(&files).unwrap(); + let (ds, report) = &layers[0]; + let rows: Vec<&crate::Feature> = ds.features.iter().collect(); + assert_eq!(rows.len(), 2); + assert!(rows[0].geometry.is_some() && rows[1].geometry.is_none(), "{:?}", rows[1].geometry); + assert!(report.messages.iter().any(|m| m.starts_with("1 features have coordinates outside")), "{report:?}"); +} + +#[test] +fn compressed_tables_and_unknown_files_are_said() { + let (t, x) = table( + 0, + &[F::Oid, F::Text("Name"), F::Int("FileFormat")], + &[vec![V::Text(Some("GDB_SystemCatalog")), V::Int(Some(0))], vec![V::Text(Some("Packed")), V::Int(Some(1))]], + ); + let files = vec![("a00000001.gdbtable".to_string(), t), ("a00000001.gdbtablx".to_string(), x)]; + assert!(matches!(read_files(&files), Err(DataError::Unsupported(m)) if m.contains("compressed"))); + assert!(matches!(read_files(&[]), Err(DataError::Format(m)) if m.contains("not a File Geodatabase"))); +} + +/// Damaged files are refused or read in part, never a panic. +#[test] +fn hostile_bytes_never_panic() { + let fields = [F::Oid, F::Shape { z: true, wkt: LONLAT }, F::Text("name"), F::Int("count")]; + let rows = vec![ + vec![V::Blob(multipatch(&[(0, vec![(0.0, 0.0, 1.0), (0.0, 5.0, 1.0), (5.0, 0.0, 2.0)])])), V::Text(Some("a")), V::Int(Some(1))], + vec![V::Blob(polygon(&square(1.0, 1.0, 3.0))), V::Text(None), V::Int(Some(2))], + ]; + let files = geodatabase(vec![("T", table(9 | 0x8000_0000, &fields, &rows))]); + assert!(read_files(&files).is_ok()); + let mut seed = 0x2545_F491_4F6C_DD1Du64; + let mut next = || { + seed ^= seed << 13; + seed ^= seed >> 7; + seed ^= seed << 17; + seed + }; + for _ in 0..3000 { + let mut damaged = files.clone(); + let which = (next() % damaged.len() as u64) as usize; + let bytes = &mut damaged[which].1; + match next() % 3 { + // Flipped bytes. + 0 => { + for _ in 0..1 + next() % 8 { + let at = (next() % bytes.len().max(1) as u64) as usize; + if let Some(b) = bytes.get_mut(at) { + *b = next() as u8; + } + } + } + // Cut short. + 1 => bytes.truncate((next() % (bytes.len() as u64 + 1)) as usize), + // Counts and sizes at their extremes. + _ => { + let at = (next() % bytes.len().max(1) as u64) as usize; + for (k, b) in [0xFFu8, 0xFF, 0xFF, 0x7F].into_iter().enumerate() { + if let Some(x) = bytes.get_mut(at + k) { + *x = b; + } + } + } + } + let _ = read_files(&damaged); + } +} + +/// A geodatabase GDAL writes (when GDAL is installed), read back without it: points with typed +/// fields, and a TIN, which GDAL stores as a multipatch (a triangle fan). +#[test] +fn reads_what_gdal_writes_when_gdal_is_installed() { + if !crate::source::gdal::available() { + return; + } + let dir = std::env::temp_dir().join(format!("rustgis-filegdb-{}", std::process::id())); + let _ = std::fs::remove_dir_all(&dir); + std::fs::create_dir_all(&dir).unwrap(); + let gdb = dir.join("written.gdb"); + let run = |args: &[&str]| { + let ok = std::process::Command::new("ogr2ogr").args(args).status().is_ok_and(|s| s.success()); + assert!(ok, "ogr2ogr {args:?}"); + }; + let points = dir.join("points.csv"); + std::fs::write(&points, "WKT,name,count,height,day\n\"POINT (8.5 47.4)\",Zürich,3,1.5,2024-05-06\n\"POINT (7.4 46.9)\",,,,\n").unwrap(); + let tins = dir.join("tins.csv"); + std::fs::write( + &tins, + "WKT,name\n\"TIN Z (((0 0 400,0 10 400,10 10 400,0 0 400)),((0 0 400,10 10 400,10 0 400,0 0 400)),((0 0 400,10 0 400,10 0 410,0 0 400)))\",block\n", + ) + .unwrap(); + let (gdb_s, points_s, tins_s) = (gdb.display().to_string(), points.display().to_string(), tins.display().to_string()); + run(&["-f", "OpenFileGDB", "-a_srs", "EPSG:4326", "-nlt", "POINT", "-nln", "places", "-oo", "AUTODETECT_TYPE=YES", &gdb_s, &points_s]); + run(&["-f", "OpenFileGDB", "-update", "-a_srs", "EPSG:2056", "-nlt", "TINZ", "-nln", "blocks", &gdb_s, &tins_s]); + assert!(is_geodatabase(&gdb)); + assert_eq!(layers(&gdb).unwrap(), ["places", "blocks"]); + let read = read_dir(&gdb, None).unwrap().layers; + let (places, report) = &read[0]; + assert_eq!(report.skipped, 0, "{report:?}"); + let rows: Vec<&crate::Feature> = places.features.iter().collect(); + assert_eq!(rows.len(), 2); + assert!( + matches!(rows[0].geometry, Some(Geometry::Point(c)) if (c.x - 8.5).abs() < 1e-7 && (c.y - 47.4).abs() < 1e-7), + "{:?}", + rows[0].geometry + ); + let at = |n: &str| places.field_index(n).unwrap(); + assert_eq!(rows[0].attrs[at("name")], Value::Text("Zürich".into())); + assert_eq!(rows[0].attrs[at("count")], Value::Int(3)); + assert_eq!(rows[0].attrs[at("height")], Value::Float(1.5)); + assert!(matches!(&rows[0].attrs[at("day")], Value::Text(d) if d.starts_with("2024-05-06")), "{:?}", rows[0].attrs); + assert_eq!(rows[1].attrs[at("count")], Value::Null); + let (blocks, _) = &read[1]; + assert_eq!(blocks.crs, Some(Crs::Epsg(2056)), "the coordinate system, by its code"); + let block = blocks.features.iter().next().unwrap(); + let z = |n: &str| block.attrs[blocks.field_index(n).unwrap()].clone(); + assert_eq!((z("z_min"), z("z_max"), z("height")), (Value::Float(400.0), Value::Float(410.0), Value::Float(10.0))); + assert!(matches!(&block.geometry, Some(Geometry::Polygon(p)) if p.len() == 1), "{:?}", block.geometry); + // One layer by name, ignoring case. + assert_eq!(read_dir(&gdb, Some("BLOCKS")).unwrap().layers.len(), 1); + let _ = std::fs::remove_dir_all(&dir); +} diff --git a/crates/data/src/lib.rs b/crates/data/src/lib.rs index 9e1b69b34..6f06c19d0 100644 --- a/crates/data/src/lib.rs +++ b/crates/data/src/lib.rs @@ -6,6 +6,7 @@ pub mod dirs; pub mod dxf; +pub mod filegdb; pub mod geometry; pub mod georss; pub mod gml; diff --git a/crates/engine/src/cmd/layer.rs b/crates/engine/src/cmd/layer.rs index 427f5a959..8d8c6c969 100644 --- a/crates/engine/src/cmd/layer.rs +++ b/crates/engine/src/cmd/layer.rs @@ -480,6 +480,59 @@ fn gdb_path(path: &str) -> std::borrow::Cow<'_, str> { } } +/// Adds a File Geodatabase's feature classes and tables (`"layer"` picks one), read without GDAL; +/// through GDAL, when it is installed, if they can't be (compressed tables, raster fields). +#[cfg(not(target_arch = "wasm32"))] +fn add_geodatabase(s: &mut Session, p: &Value, path: &str) -> CmdResult { + let layer = ps(p, "layer").map(str::trim).filter(|l| !l.is_empty()); + let read = match rustgis_data::filegdb::read_dir(Path::new(path), layer) { + Ok(read) => read, + Err(e) if rustgis_data::source::gdal::available() => { + return add_gdal(s, p, path).map_err(|g| format!("{path}: {e}; through GDAL: {g}")); + } + Err(e) => return Err(format!("{path}: {e}")), + }; + let reads = + read.layers.into_iter().map(|(dataset, report)| rustgis_data::io::LayerRead { dataset, report, styles: Vec::new() }).collect(); + // The tables read, then those it couldn't read (compressed, raster fields), through GDAL when it + // is installed; else, or if GDAL can't either, said. The second go on whether the first added. + let mut results: Vec = Vec::new(); + let mut failed: Vec = Vec::new(); + match add_reads(s, p, path, reads, Vec::new()) { + Ok(r) => results.push(r), + Err(e) => failed.push(e), + } + let gdal = rustgis_data::source::gdal::available(); + for (table, why) in read.unread { + let mut q = p.clone(); + if let Some(o) = q.as_object_mut() { + o.insert("layer".into(), json!(table)); + o.remove("name"); + } + match gdal.then(|| add_gdal(s, &q, path)) { + Some(Ok(r)) => results.push(r), + Some(Err(g)) => s.notify(format!("{path}: {table} not read: {why}; through GDAL: {g}"), true), + None => s.notify(format!("{path}: {table} not read: {why} (GDAL, when installed, reads it)"), true), + } + } + if results.is_empty() { + return Err(if failed.is_empty() { format!("{path}: nothing could be added") } else { failed.join("; ") }); + } + // Something was added (through GDAL): the layers read natively that couldn't be are still said. + for e in failed { + s.notify(e, true); + } + // One result, as adding a single file gives it: every layer, the first, and the features. + let layers: Vec = results.iter().filter_map(|r| r.get("layers").and_then(Value::as_array)).flatten().cloned().collect(); + let tables: Vec = results.iter().filter_map(|r| r.get("tables").and_then(Value::as_array)).flatten().cloned().collect(); + let features: u64 = results.iter().filter_map(|r| r.get("features").and_then(Value::as_u64)).sum(); + let mut out = json!({"layers": layers, "tables": tables, "features": features}); + if let (Some(o), Some(first)) = (out.as_object_mut(), layers.first()) { + o.insert("layer".into(), first.clone()); + } + Ok(out) +} + fn add_data(s: &mut Session, p: &Value) -> CmdResult { let path = s_req(p, "path")?; let path = &*gdb_path(path); @@ -568,6 +621,10 @@ fn add_data(s: &mut Session, p: &Value) -> CmdResult { return super::tiles::add(s, &q); } #[cfg(not(target_arch = "wasm32"))] + if rustgis_data::filegdb::is_geodatabase(Path::new(path)) { + return add_geodatabase(s, p, path); + } + #[cfg(not(target_arch = "wasm32"))] if rustgis_data::source::gdal::is_gdal_name(path) { return add_gdal(s, p, path); } @@ -666,6 +723,18 @@ pub(crate) fn add_layers(s: &mut Session, p: &Value, file_name: &str, bytes: &[u return Err(format!("{file_name}: no layer '{want}' (layers: {})", names.join(", "))); } } + add_reads(s, p, file_name, reads, overlays) +} + +/// Adds the layers a file gave and its ground overlays: a layer that can't be added is skipped +/// and said, the others are still added. A name given goes to the first layer added. +fn add_reads( + s: &mut Session, + p: &Value, + file_name: &str, + reads: Vec, + overlays: Vec<(String, rustgis_data::raster::Raster)>, +) -> CmdResult { let mut added = Vec::new(); let mut failed = Vec::new(); let mut q = p.clone(); @@ -828,6 +897,15 @@ fn file_layers(_s: &mut Session, p: &Value) -> CmdResult { let path = s_req(p, "path")?; let path = &*gdb_path(path); let err = |e: rustgis_data::DataError| format!("{path}: {e}"); + // A geodatabase read without GDAL; through GDAL if it can't be and GDAL is installed. + #[cfg(not(target_arch = "wasm32"))] + if rustgis_data::filegdb::is_geodatabase(Path::new(path)) { + match rustgis_data::filegdb::layers(Path::new(path)) { + Ok(layers) => return Ok(json!({"count": layers.len(), "layers": layers})), + Err(_) if rustgis_data::source::gdal::available() => {} + Err(e) => return Err(err(e)), + } + } #[cfg(not(target_arch = "wasm32"))] if rustgis_data::source::gdal::is_gdal_name(path) { let layers = rustgis_data::source::gdal::layers(Path::new(path)).map_err(err)?; diff --git a/docs/user-guide/data-formats.md b/docs/user-guide/data-formats.md index 4b2bc89d8..a285daa80 100644 --- a/docs/user-guide/data-formats.md +++ b/docs/user-guide/data-formats.md @@ -31,6 +31,7 @@ web app reads GeoPackages without SQLite, with a small SQLite file reader writte | GeoPackage | `.gpkg` | Add Data, or upload in the web app | Every feature table becomes a layer, up to 32 at once (the result counts the rest in `skippedCount` and names up to 100 of them in `skipped`); `"layer"` picks one, and `info.fileLayers` lists them all. Read in place by its R-tree; any EPSG system is reprojected. The web app reads the file with its own SQLite reader (tables scanned, the envelope in each geometry header tested first; no R-tree); UTF-16 databases aren't read there. | | GeoParquet | `.parquet`, `.geoparquet` | Add Data, or upload in the web app | WKB geometry, with the CRS from the file's `geo` metadata. Row groups are pruned by the bbox covering. Snappy, gzip, LZ4 and zstd compression are read, in the browser too; Brotli isn't. | | FlatGeobuf | `.fgb` | Add Data, or upload | Read in place by its packed R-tree. | +| Esri File Geodatabase | `.gdb` folder (desktop) | A file dialog can't pick a folder: in Add Data pick any `.gdbtable` file inside it; or double-click the `.gdb` in the Browser, or drop the folder on the window | Read without GDAL: every feature class becomes a layer and every table a standalone table, with their fields (dates as text) and coordinate system (by its code, EPSG or Esri; else its WKT). `"layer"` picks one, and `info.fileLayers` lists them. 3D multipatches (such as swissBUILDINGS3D's buildings) are drawn as their footprints, the union of their faces seen from above, with their lowest and highest Z in `z_min` and `z_max` and the difference in `height`, which an extrusion (Style › Extrusion) can stand them up by; one with no face seen from above (only walls) is drawn as the outline of its points. Curves are drawn as straight segments between their ends. Compressed tables, raster fields and some newer row indexes aren't read: with GDAL installed, such a geodatabase opens through GDAL instead (see below). Loaded whole, up to 2 GB a table. | | Standalone tables | CSV, or any layer without geometry | Add Data | Listed under the map in Layers; join or relate them to layers. | Files with more than 100,000 features (GeoPackage, GeoParquet, FlatGeobuf, and Shapefile on the @@ -47,7 +48,6 @@ GDAL](data.md#formats-through-gdal). | Format | Extensions or source | Notes | | --- | --- | --- | -| Esri File Geodatabase | `.gdb` folder | Every feature class becomes a layer. A file dialog can't pick a folder: pick any `.gdbtable` file inside it, double-click the `.gdb` in the Browser, or drop the folder on the window. | | AutoCAD DWG | `.dwg` | Through GDAL's CAD driver, when it is built in. Give `"crs"` (such as `EPSG:2263`) when the drawing doesn't say its coordinate system. | | MapInfo tables | `.tab` | Keep the companion files together. | | OpenStreetMap XML | `.osm` | | diff --git a/docs/user-guide/data.md b/docs/user-guide/data.md index bf1af26a3..bee7f5738 100644 --- a/docs/user-guide/data.md +++ b/docs/user-guide/data.md @@ -18,7 +18,8 @@ Formats](data-formats.md). | GML, MapInfo MIF/MID, OpenStreetMap PBF, TopoJSON, AutoCAD DXF (ASCII) | Yes | No | | LiDAR point clouds: LAS, LAZ, COPC (an evenly spread sample) | Yes | No | | Geotagged photos (JPEG with EXIF GPS): points with their thumbnails | Yes | No | -| Through GDAL: File Geodatabase (`.gdb` folder), AutoCAD DWG, MapInfo TAB, OpenStreetMap XML… | Yes (desktop, when GDAL is installed) | No | +| Esri File Geodatabase (`.gdb` folder; 3D multipatches as footprints) | Yes (desktop) | No | +| Through GDAL: AutoCAD DWG, MapInfo TAB, OpenStreetMap XML, File Geodatabases RustGIS can't read itself… | Yes (desktop, when GDAL is installed) | No | | PostgreSQL/PostGIS databases, read live (large tables a view at a time) | Yes (desktop) | No | | MySQL and MariaDB databases, read live (large tables a view at a time) | Yes (desktop) | No | | DuckDB databases (`.duckdb`, `.ddb`) and SQL queries | Yes (desktop, when DuckDB is installed) | No | @@ -59,7 +60,8 @@ reprojected on the fly as they draw. When [GDAL](https://gdal.org) is installed (its `ogr2ogr` and `ogrinfo` programs on the `PATH`, or the folder holding them in `RUSTGIS_GDAL_BIN`; GDAL 3.7 or later), the desktop app also opens -the vector formats GDAL reads: File Geodatabases, DWG, MapInfo TAB, OpenStreetMap XML +the vector formats GDAL reads: DWG, MapInfo TAB, OpenStreetMap XML, File Geodatabases RustGIS +can't read itself (compressed ones, or ones with raster fields) and more. RustGIS doesn't link to GDAL: each layer is converted once, by `ogr2ogr`, into a FlatGeobuf with a spatial index in WGS 84 (GDAL reprojects it), kept in the cache folder (`~/.cache/rustgis/gdal` on Linux, `~/Library/Caches/RustGIS/gdal` on macOS,