osm-katana/tests/metric_dp.rs
//! **Row 29a — the simplifier measures in DEGREES, so a shore is thinned by its
//! COMPASS BEARING.**
//!
//! # The defect
//!
//! `tiles::douglas_peucker_into` computes its perpendicular distance straight out of
//! `(lon, lat)` pairs:
//!
//! ```text
//! cross = (px - ax) * dy - (py - ay) * dx // dx in degrees LONGITUDE
//! // dy in degrees LATITUDE
//! ```
//!
//! A degree of latitude is ~111 km everywhere. A degree of longitude is 111 km at the
//! equator and **55.6 km at 60°N** — the middle of Sweden. So the same physical
//! deviation, measured across a NORTH–SOUTH feature (where the deviation is in
//! longitude), reads about **twice as large** as it does across an EAST–WEST one, and
//! `eps` admits it where it would have rejected its east–west twin.
//!
//! The consequence is not "slightly wrong vertex counts": it is that the generaliser's
//! aggressiveness is a function of which way a coastline happens to run. Two identical
//! bays, one facing north and one facing east, come back at different levels of detail
//! on the same rung.
//!
//! # What is asserted, and why it is APPLIED OUTPUT
//!
//! Everything below goes through a real `spatial::pack_batches` and reads the vertex
//! counts back out of the written `geometry_z9` column. The kernel is `pub(crate)`, but
//! that is not why: a simplifier tested through its own function proves the arithmetic
//! and not the pipeline, and the rung columns are what the map actually reads.
//!
//! Run: `cargo test -p osm-katana --test metric_dp -- --nocapture`
use std::sync::Arc;
use arrow::array::{Array, ArrayRef, BinaryArray, BinaryBuilder, Int32Builder, Int64Builder, StringBuilder};
use arrow::datatypes::Schema;
use arrow::record_batch::RecordBatch;
use osm_katana::spatial::{PackOptions, pack_batches};
use parquet::arrow::arrow_reader::ParquetRecordBatchReaderBuilder;
/// Metres per degree of latitude — the one conversion this test uses, written out.
const M_PER_DEG: f64 = 111_320.0;
/// The latitude the comparison is made at. 60°N is the middle of Sweden and
/// `cos 60° = 0.5` exactly, which is what makes the expected factor a round 2.
const LAT: f64 = 60.0;
/// **One shore, in METRES**: an along-track run with a cross-track wobble whose
/// amplitude RAMPS through the simplifier's epsilon.
///
/// The ramp is the point. A wobble far above `eps` is kept whichever units the kernel
/// works in, and one far below is dropped either way; only amplitudes *near* the
/// threshold change their verdict when the threshold is effectively doubled. A fixture
/// with a single amplitude would come back either all-kept or all-dropped in both
/// orientations and would prove nothing.
fn shore_metres(n: usize, length_m: f64, max_amp_m: f64) -> Vec<(f64, f64)> {
(0..n)
.map(|i| {
let t = i as f64 / (n - 1) as f64;
let s = t * length_m;
// Amplitude ramps 0 → max across the run; wavelength stays put.
let amp = t * max_amp_m;
let w = amp * (s / 400.0).sin();
(s, w)
})
.collect()
}
/// Place a metric shore running EAST–WEST at [`LAT`]: along-track is longitude, so the
/// cross-track deviation lands in LATITUDE.
fn place_ew(shore: &[(f64, f64)], lon0: f64, lat0: f64) -> Vec<(f64, f64)> {
let k = (LAT.to_radians()).cos();
shore
.iter()
.map(|&(s, w)| (lon0 + s / (M_PER_DEG * k), lat0 + w / M_PER_DEG))
.collect()
}
/// The same metric shore running NORTH–SOUTH: along-track is latitude, so the
/// cross-track deviation lands in LONGITUDE — where a degree is only 55.6 km.
fn place_ns(shore: &[(f64, f64)], lon0: f64, lat0: f64) -> Vec<(f64, f64)> {
let k = (LAT.to_radians()).cos();
shore
.iter()
.map(|&(s, w)| (lon0 + w / (M_PER_DEG * k), lat0 + s / M_PER_DEG))
.collect()
}
/// The number of points a WKB `LineString` declares.
///
/// Read straight out of the header — byte order, then the `u32` type, then the `u32`
/// point count — rather than through a decoder, so this test cannot pass because a
/// decoder and an encoder agree with each other about something wrong.
fn wkb_line_points(b: &[u8]) -> u32 {
assert!(b.len() >= 9, "not a WKB LineString: {} bytes", b.len());
assert_eq!(b[0], 1, "expected little-endian WKB");
let ty = u32::from_le_bytes([b[1], b[2], b[3], b[4]]);
assert_eq!(ty, 2, "expected WKB type 2 (LineString), got {ty}");
u32::from_le_bytes([b[5], b[6], b[7], b[8]])
}
/// Pack the two shores and hand back `(ew_vertices, ns_vertices)` from `geometry_z9`,
/// read off disk.
fn packed_z9_vertices(rows: &[Vec<(f64, f64)>]) -> Vec<u32> {
let mut meta = std::collections::HashMap::new();
meta.insert(
String::from(osm_katana::metadata::GEO_KEY),
osm_katana::metadata::geo_metadata(&["LineString"], None),
);
let schema = Arc::new(Schema::new_with_metadata(
osm_katana::shared::schema::ways_schema(false).fields().clone(),
meta,
));
let (mut id, mut g, mut tags, mut ver) = (
Int64Builder::new(),
BinaryBuilder::new(),
StringBuilder::new(),
Int32Builder::new(),
);
let mut cls = osm_katana::shared::schema::ClassColumns::new();
for (i, pts) in rows.iter().enumerate() {
id.append_value(i as i64);
g.append_value(osm_katana::shared::geometry::encode_wkb_line(pts));
// `natural=coastline` puts the row on rung 6, so it is admitted at every rung
// and `geometry_z9` is non-null for all of them.
let t = r#"{"natural":"coastline"}"#;
tags.append_value(t);
ver.append_value(1);
cls.push(t);
}
let mut cols: Vec<ArrayRef> = vec![
Arc::new(id.finish()),
Arc::new(g.finish()),
Arc::new(tags.finish()),
Arc::new(ver.finish()),
];
cols.extend(cls.finish());
let batch = RecordBatch::try_new(schema.clone(), cols).unwrap();
let dir = tempfile::tempdir().unwrap();
let dst = dir.path().join("ways.parquet");
pack_batches(vec![batch], &schema, &dst, &PackOptions::default(), &[]).unwrap();
// Read `geometry_z9` back, keyed by `id` so the pack's Hilbert reordering cannot
// silently swap the two rows under us — which it would, and which would make this
// test compare each shore with the other.
let rd = ParquetRecordBatchReaderBuilder::try_new(std::fs::File::open(&dst).unwrap())
.unwrap()
.build()
.unwrap();
let mut by_id: std::collections::BTreeMap<i64, u32> = std::collections::BTreeMap::new();
for b in rd {
let b = b.unwrap();
let ids = b
.column_by_name("id")
.unwrap()
.as_any()
.downcast_ref::<arrow::array::Int64Array>()
.unwrap();
let gz = b
.column_by_name("geometry_z9")
.expect("the packed table carries no geometry_z9 column")
.as_any()
.downcast_ref::<BinaryArray>()
.unwrap();
for i in 0..b.num_rows() {
assert!(
!gz.is_null(i),
"row {} has a NULL geometry_z9 — it is not admitted at rung 9, so this \
test would be measuring nothing",
ids.value(i)
);
by_id.insert(ids.value(i), wkb_line_points(gz.value(i)));
}
}
by_id.into_values().collect()
}
/// ★ **RED-when-broken: the same shore is thinned differently depending on which way
/// it runs.**
///
/// Two metrically IDENTICAL shores at 60°N, one east–west and one north–south, through
/// a real pack, counted off disk. They must come back within **±5 %** of each other.
///
/// **RED ARM, seen 2026-09-05** with the degrees kernel: the east–west shore keeps
/// **111** vertices and the north–south one **171** — a 54 % spread, because a
/// cross-track deviation expressed in longitude at 60°N reads twice its true size and
/// clears an epsilon its east–west twin does not.
///
/// The non-vacuity arms matter as much: both shores must keep a real fraction of their
/// vertices (a simplifier that dropped everything, or kept everything, would sail
/// through the ±5 %), and the fixture must actually straddle the threshold.
#[test]
fn a_shore_is_thinned_by_its_shape_and_not_by_its_bearing() {
let eps_deg = osm_katana::shared::generalise::generalisation_eps(9);
let eps_m = eps_deg * M_PER_DEG;
// 20 km of shore, wobble ramping to 4x epsilon: the run crosses the threshold about
// a quarter of the way along, so roughly three quarters of the wobbles survive.
let shore = shore_metres(16_384, 200_000.0, 4.0 * eps_m);
let ew = place_ew(&shore, 17.0, LAT);
let ns = place_ns(&shore, 17.0, LAT);
let v = packed_z9_vertices(&[ew.clone(), ns.clone()]);
let (kept_ew, kept_ns) = (v[0], v[1]);
println!("z9 eps {eps_deg:.3e}° ≈ {eps_m:.1} m · E-W kept {kept_ew} · N-S kept {kept_ns}");
// NON-VACUITY first, so a collapsed or a no-op simplifier cannot pass the ratio.
for (what, kept) in [("E-W", kept_ew), ("N-S", kept_ns)] {
assert!(
kept > 8,
"{what} came back with {kept} vertices — the simplifier collapsed the shore, \
so the comparison below is between two straight lines"
);
assert!(
(kept as usize) < ew.len() / 2,
"{what} kept {kept} of {} vertices — the simplifier did almost nothing, so \
the comparison below is between two unsimplified copies",
ew.len()
);
}
let spread = (kept_ew as f64 - kept_ns as f64).abs() / (kept_ew.max(kept_ns) as f64);
assert!(
spread <= 0.05,
"the same shore kept {kept_ew} vertices running EAST-WEST and {kept_ns} running \
NORTH-SOUTH — a {:.1} % spread at {LAT}°N. The Douglas-Peucker kernel measures \
perpendicular distance in raw degrees, so a deviation expressed in LONGITUDE \
reads 1/cos(lat) times its true size (2x here) and clears an epsilon its \
east-west twin does not. The generaliser's aggressiveness must not be a \
function of a feature's compass bearing.",
spread * 100.0
);
}
/// ★ **RED-when-broken: the fix is a constant, not a projection.**
///
/// The anti-identity arm. Scaling longitude by `cos(lat)` must have **no effect at the
/// equator** — `cos 0 = 1` — and a large one at 60°N. A "fix" that scaled by a constant,
/// or that scaled both axes, would pass the bearing test above and fail here.
#[test]
fn the_correction_is_the_latitude_and_vanishes_at_the_equator() {
let eps_m = osm_katana::shared::generalise::generalisation_eps(9) * M_PER_DEG;
let shore = shore_metres(16_384, 200_000.0, 4.0 * eps_m);
// At the EQUATOR a degree of longitude is a degree of latitude, so the two
// orientations must already agree — before and after the fix alike.
let k_eq = 1.0_f64;
let ew_eq: Vec<(f64, f64)> = shore
.iter()
.map(|&(s, w)| (17.0 + s / (M_PER_DEG * k_eq), 0.0 + w / M_PER_DEG))
.collect();
let ns_eq: Vec<(f64, f64)> = shore
.iter()
.map(|&(s, w)| (17.0 + w / (M_PER_DEG * k_eq), 0.0 + s / M_PER_DEG))
.collect();
let v = packed_z9_vertices(&[ew_eq, ns_eq]);
let spread = (v[0] as f64 - v[1] as f64).abs() / (v[0].max(v[1]) as f64);
println!("equator: E-W {} · N-S {} · spread {:.2} %", v[0], v[1], spread * 100.0);
assert!(
spread <= 0.05,
"at the equator the two orientations already agree in degrees, so a correction \
that changes anything here is not cos(lat): E-W {} vs N-S {}",
v[0],
v[1]
);
// And the correction must BITE at 60°N: a north–south shore there keeps a
// measurably different count than the same shore at the equator would if the kernel
// were still working in raw degrees. Compare the 60°N north–south shore against the
// equator's north–south shore — with the fix they agree, without it they do not.
let ns_60 = place_ns(&shore, 17.0, LAT);
let n60 = packed_z9_vertices(&[ns_60])[0];
let spread_lat = (n60 as f64 - v[1] as f64).abs() / (n60.max(v[1]) as f64);
println!("N-S at 0°N {} vs at {LAT}°N {} · spread {:.2} %", v[1], n60, spread_lat * 100.0);
assert!(
spread_lat <= 0.05,
"the same metric shore kept {} vertices at the equator and {n60} at {LAT}°N — \
the simplifier is still latitude-dependent, so the correction is not being \
applied (or is applied to the wrong axis)",
v[1]
);
}