use geo::CoordsIter as _;
use geo_types::{Coord, Geometry};
use crate::codecs::hilbert::{hilbert_curve_params_from_bounds, hilbert_sort_key};
use crate::codecs::morton::morton_sort_key;
use crate::encoder::model::CurveParams;
use crate::tile::{TileFeature, TileLayer};
#[derive(Debug, Clone, Copy, Default, PartialEq, Eq, strum::EnumIter, strum::EnumCount)]
pub enum SortStrategy {
#[default]
Unsorted,
SpatialMorton,
SpatialHilbert,
Id,
}
impl TileLayer {
#[hotpath::measure]
pub fn sort(&mut self, strategy: SortStrategy, params: CurveParams) {
match strategy {
SortStrategy::SpatialMorton | SortStrategy::SpatialHilbert => {
let curve_key = if let SortStrategy::SpatialMorton = strategy {
morton_sort_key
} else {
hilbert_sort_key
};
self.features_mut().sort_by_cached_key(|f| {
first_vertex(f.geometry()).map_or(u64::MAX, |c| u64::from(curve_key(c, params)))
});
}
SortStrategy::Id => {
self.features_mut().sort_by_cached_key(TileFeature::id);
}
SortStrategy::Unsorted => {
}
}
}
#[hotpath::measure]
#[must_use]
pub fn curve_params(&self) -> CurveParams {
let (min_val, max_val) = self
.features()
.iter()
.flat_map(|f| f.geometry().coords_iter())
.fold((i32::MAX, i32::MIN), |(min, max), c| {
(min.min(c.x).min(c.y), max.max(c.x).max(c.y))
});
hilbert_curve_params_from_bounds(min_val, max_val)
}
}
fn first_vertex(geom: &Geometry<i32>) -> Option<Coord<i32>> {
match geom {
Geometry::<i32>::Point(p) => Some(p.0),
Geometry::<i32>::Line(l) => Some(l.start),
Geometry::<i32>::LineString(ls) => ls.0.first().copied(),
Geometry::<i32>::Polygon(p) => p.exterior().0.first().copied(),
Geometry::<i32>::MultiPoint(mp) => mp.0.first().map(|p| p.0),
Geometry::<i32>::MultiLineString(mls) => mls.0.first().and_then(|ls| ls.0.first().copied()),
Geometry::<i32>::MultiPolygon(mp) => {
mp.0.first().and_then(|p| p.exterior().0.first().copied())
}
Geometry::<i32>::Triangle(t) => Some(t.v1()),
Geometry::<i32>::Rect(r) => Some(r.min()),
Geometry::<i32>::GeometryCollection(gc) => gc.0.first().and_then(first_vertex),
}
}
pub(crate) fn spatial_sort_likely_to_help(layer: &TileLayer) -> bool {
const SPATIAL_HELP_COVERAGE: f64 = 0.8;
let extent = f64::from(layer.extent().get());
if extent <= 0.0 || layer.features().is_empty() {
return true;
}
let (min_x, max_x, min_y, max_y) = layer
.features()
.iter()
.filter_map(|f| first_vertex(f.geometry()))
.fold(
(i32::MAX, i32::MIN, i32::MAX, i32::MIN),
|(min_x, max_x, min_y, max_y), Coord::<i32> { x, y }| {
(min_x.min(x), max_x.max(x), min_y.min(y), max_y.max(y))
},
);
if min_x > max_x || min_y > max_y {
return true;
}
let range_x = f64::from(max_x) - f64::from(min_x);
let range_y = f64::from(max_y) - f64::from(min_y);
let spread_x = range_x / extent;
let spread_y = range_y / extent;
!(spread_x > SPATIAL_HELP_COVERAGE && spread_y > SPATIAL_HELP_COVERAGE)
}
#[cfg(test)]
mod tests {
use geo_types::{
Coord, Geometry as GeoGeom, Geometry, GeometryCollection, Line, LineString,
MultiLineString, MultiPoint, MultiPolygon, Point, Polygon, Rect, Triangle,
};
use rstest::rstest;
use super::spatial_sort_likely_to_help;
use crate::decoder::{GeometryType, GeometryValues, RawGeometry};
use crate::encoder::model::CurveParams;
use crate::encoder::{
Codecs, Encoder, EncoderConfig, ExplicitEncoder, IntEncoder, SortStrategy, stage_tile,
};
use crate::test_helpers::{assert_empty, dec, into_layer01, parser};
use crate::tile::{TileFeature, TileLayer};
use crate::{Layer, LazyParsed};
fn pt(x: i32, y: i32) -> Geometry<i32> {
GeoGeom::Point(Point::new(x, y))
}
fn ls(coords: &[(i32, i32)]) -> Geometry<i32> {
GeoGeom::LineString(LineString::new(
coords.iter().map(|&(x, y)| Coord { x, y }).collect(),
))
}
fn poly_square(x0: i32, y0: i32, side: i32) -> Geometry<i32> {
let ring = LineString::new(vec![
Coord { x: x0, y: y0 },
Coord {
x: x0 + side,
y: y0,
},
Coord {
x: x0 + side,
y: y0 + side,
},
Coord {
x: x0,
y: y0 + side,
},
Coord { x: x0, y: y0 },
]);
GeoGeom::Polygon(Polygon::new(ring, vec![]))
}
fn roundtrip_geom(decoded: &GeometryValues) -> GeometryValues {
let mut enc = Encoder::default();
let mut codecs = Codecs::default();
decoded
.clone()
.write_to(&mut enc, &mut codecs)
.expect("encode failed");
let buf = enc.data().to_vec();
let parsed = assert_empty(RawGeometry::from_bytes(&buf, &mut parser()));
let mut d = dec();
let result = LazyParsed::Raw(parsed)
.into_parsed(&mut d)
.expect("decode failed");
assert!(
d.consumed() > 0,
"decoder should consume bytes after decode"
);
result
}
fn canonical(geoms: &[Geometry<i32>]) -> GeometryValues {
let mut decoded = GeometryValues::default();
for g in geoms {
decoded.push_geom(g);
}
roundtrip_geom(&decoded)
}
fn layer_after_sort(geoms: &[Geometry<i32>], ids: &[u64], strategy: SortStrategy) -> TileLayer {
let features: Vec<TileFeature> = geoms
.iter()
.zip(ids.iter())
.map(|(g, &id)| TileFeature {
id: Some(id),
geometry: g.clone(),
properties: vec![],
#[cfg(feature = "unstable-v2")]
m_values: vec![],
#[cfg(feature = "unstable-v2")]
nested: vec![],
#[cfg(feature = "unstable-v2")]
z: Vec::new(),
})
.collect();
let mut layer = TileLayer::from_parts("test", 4096, vec![], features).unwrap();
let params = layer.curve_params();
layer.sort(strategy, params);
layer
}
fn assert_sort_roundtrip(
geoms: &[Geometry<i32>],
ids: &[u64],
strategy: SortStrategy,
expected: &[Geometry<i32>],
) {
let layer = layer_after_sort(geoms, ids, strategy);
let mut sorted_decoded = GeometryValues::default();
for f in layer.features() {
sorted_decoded.push_geom(f.geometry());
}
let after_roundtrip = roundtrip_geom(&sorted_decoded);
let expected_canonical = canonical(expected);
assert_eq!(
after_roundtrip, expected_canonical,
"\nsorted geometry did not match expected after encode->decode round-trip\
\nvector_types after sort: {:?}\
\nvector_types expected: {:?}",
sorted_decoded.vector_types, expected_canonical.vector_types,
);
}
#[test]
fn pure_points_id_sort_roundtrip() {
assert_sort_roundtrip(
&[pt(0, 0), pt(1, 1), pt(2, 2)],
&[3, 2, 1],
SortStrategy::Id,
&[pt(2, 2), pt(1, 1), pt(0, 0)],
);
}
#[test]
fn pure_linestrings_id_sort_roundtrip() {
assert_sort_roundtrip(
&[ls(&[(0, 0), (0, 10)]), ls(&[(5, 5), (10, 10)])],
&[2, 1],
SortStrategy::Id,
&[ls(&[(5, 5), (10, 10)]), ls(&[(0, 0), (0, 10)])],
);
}
#[test]
fn point_line_point_id_sort_to_line_point_point_roundtrip() {
assert_sort_roundtrip(
&[pt(0, 0), ls(&[(1, 0), (1, 5)]), pt(5, 5)],
&[3, 1, 2],
SortStrategy::Id,
&[ls(&[(1, 0), (1, 5)]), pt(5, 5), pt(0, 0)],
);
}
#[test]
fn point_line_point_id_sort_to_point_point_line_roundtrip() {
assert_sort_roundtrip(
&[pt(0, 0), ls(&[(1, 0), (1, 5)]), pt(5, 5)],
&[1, 3, 2],
SortStrategy::Id,
&[pt(0, 0), pt(5, 5), ls(&[(1, 0), (1, 5)])],
);
}
#[test]
fn point_polygon_point_id_sort_roundtrip() {
assert_sort_roundtrip(
&[pt(0, 0), poly_square(10, 10, 5), pt(5, 5)],
&[2, 1, 3],
SortStrategy::Id,
&[poly_square(10, 10, 5), pt(0, 0), pt(5, 5)],
);
}
#[test]
fn point_line_point_morton_sort_roundtrip() {
assert_sort_roundtrip(
&[pt(2, 0), ls(&[(0, 0), (0, 5)]), pt(1, 0)],
&[1, 2, 3],
SortStrategy::SpatialMorton,
&[ls(&[(0, 0), (0, 5)]), pt(1, 0), pt(2, 0)],
);
}
#[test]
fn id_sort_already_sorted_is_identity_roundtrip() {
let geoms = &[pt(0, 0), ls(&[(1, 0), (1, 5)]), pt(5, 5)];
assert_sort_roundtrip(geoms, &[1, 2, 3], SortStrategy::Id, geoms);
}
#[test]
fn id_column_co_permuted_with_geometry() {
let layer = layer_after_sort(
&[pt(0, 0), ls(&[(1, 0), (1, 5)]), pt(5, 5)],
&[3, 1, 2],
SortStrategy::Id,
);
let ids: Vec<Option<u64>> = layer.features().iter().map(TileFeature::id).collect();
assert_eq!(ids, vec![Some(1u64), Some(2), Some(3)]);
let geom_types: Vec<&str> = layer
.features()
.iter()
.map(|f| GeometryType::try_from(f.geometry()).unwrap().into())
.collect();
assert_eq!(geom_types, vec!["LineString", "Point", "Point"]);
}
fn build_tile_layer(geoms: &[Geometry<i32>], ids: &[Option<u64>]) -> TileLayer {
assert_eq!(geoms.len(), ids.len());
TileLayer::from_parts(
"test",
4096,
vec![],
geoms
.iter()
.zip(ids.iter())
.map(|(g, &id)| TileFeature {
id,
geometry: g.clone(),
properties: vec![],
#[cfg(feature = "unstable-v2")]
m_values: vec![],
#[cfg(feature = "unstable-v2")]
nested: vec![],
#[cfg(feature = "unstable-v2")]
z: Vec::new(),
})
.collect(),
)
.unwrap()
}
fn sort_encode_decode(tile: TileLayer, sort: SortStrategy) -> TileLayer {
let enc_cfg = EncoderConfig::default();
let enc = Encoder::with_explicit(enc_cfg, ExplicitEncoder::for_id(IntEncoder::varint()));
let mut codecs = Codecs::default();
let enc = stage_tile(tile, sort, false, enc_cfg.tessellate())
.encode_into(enc, &mut codecs)
.expect("encode failed");
let buf = enc.into_layer_bytes().expect("into_layer_bytes failed");
let mut p = parser();
let layer_back = assert_empty(Layer::from_bytes(&buf, &mut p));
assert!(p.reserved() > 0, "parser should reserve bytes after parse");
let layer01 = into_layer01(layer_back);
let mut d = dec();
let tile = layer01.into_tile(&mut d).expect("decode after sort failed");
assert!(
d.consumed() > 0,
"decoder should consume bytes after decode"
);
tile
}
fn vertices_from_source(source: &TileLayer) -> Vec<i32> {
let mut geom = GeometryValues::default();
for f in source.features() {
geom.push_geom(f.geometry());
}
geom.vertices().unwrap_or_default().to_vec()
}
#[test]
fn test_shared_morton_shift() {
let tile = build_tile_layer(&[pt(0, -10), pt(-10, 0)], &[Some(1), Some(2)]);
let source = sort_encode_decode(tile, SortStrategy::SpatialMorton);
let verts = vertices_from_source(&source);
assert_eq!(verts, vec![0, -10, -10, 0]);
}
#[test]
fn test_id_sort_nulls_first() {
let tile = build_tile_layer(&[pt(2, 2), pt(1, 1), pt(0, 0)], &[Some(10), None, Some(5)]);
let source = sort_encode_decode(tile, SortStrategy::Id);
let ids: Vec<Option<u64>> = source.features().iter().map(TileFeature::id).collect();
assert_eq!(ids, vec![None, Some(5), Some(10)]);
let verts = vertices_from_source(&source);
assert_eq!(verts, vec![1, 1, 0, 0, 2, 2]);
}
#[test]
fn test_mixed_geometry_morton_sort() {
let tile = build_tile_layer(
&[pt(2, 0), ls(&[(0, 0), (0, 5)]), pt(1, 0)],
&[Some(1), Some(2), Some(3)],
);
let source = sort_encode_decode(tile, SortStrategy::SpatialMorton);
let types: Vec<_> = source
.features()
.iter()
.map(|f| GeometryType::try_from(f.geometry()).unwrap())
.collect();
assert_eq!(
types,
vec![
GeometryType::LineString,
GeometryType::Point,
GeometryType::Point
]
);
let verts = vertices_from_source(&source);
assert_eq!(verts, vec![0, 0, 0, 5, 1, 0, 2, 0]);
}
#[test]
fn id_sort_separates_the_two_largest_ids() {
let layer = layer_after_sort(
&[pt(0, 0), pt(1, 1)],
&[u64::MAX, u64::MAX - 1],
SortStrategy::Id,
);
let ids: Vec<Option<u64>> = layer.features().iter().map(TileFeature::id).collect();
assert_eq!(ids, vec![Some(u64::MAX - 1), Some(u64::MAX)]);
}
#[test]
fn id_sort_puts_the_missing_id_before_zero() {
let layer = build_tile_layer(&[pt(0, 0), pt(1, 1)], &[Some(0), None]);
let mut layer = layer;
layer.sort(SortStrategy::Id, CurveParams { shift: 0, bits: 1 });
let ids: Vec<Option<u64>> = layer.features().iter().map(TileFeature::id).collect();
assert_eq!(ids, vec![None, Some(0)]);
}
#[test]
fn spatial_sort_reads_the_first_vertex_of_every_geometry_kind() {
let at = |k: i32| Coord { x: k, y: k };
let geoms = [
GeoGeom::Rect(Rect::new(at(0), at(20))),
GeoGeom::Triangle(Triangle::new(at(1), at(11), at(21))),
GeoGeom::MultiPolygon(MultiPolygon(vec![Polygon::new(
LineString(vec![at(2), at(12), at(22), at(2)]),
vec![],
)])),
GeoGeom::MultiLineString(MultiLineString(vec![LineString(vec![at(3), at(13)])])),
GeoGeom::MultiPoint(MultiPoint(vec![Point::from(at(4)), Point::from(at(14))])),
GeoGeom::Polygon(Polygon::new(
LineString(vec![at(5), at(15), at(25), at(5)]),
vec![],
)),
GeoGeom::LineString(LineString(vec![at(6), at(16)])),
GeoGeom::Line(Line::new(at(7), at(17))),
GeoGeom::Point(Point::from(at(8))),
GeoGeom::GeometryCollection(GeometryCollection(vec![GeoGeom::Point(Point::from(at(
9,
)))])),
GeoGeom::LineString(LineString(vec![])),
];
let ids: Vec<Option<u64>> = (0..u64::try_from(geoms.len()).unwrap()).map(Some).collect();
let mut layer = build_tile_layer(&geoms, &ids);
let params = layer.curve_params();
layer.sort(SortStrategy::SpatialMorton, params);
let sorted: Vec<Option<u64>> = layer.features().iter().map(TileFeature::id).collect();
assert_eq!(sorted, ids);
}
fn points_help_spatial_sort(coords: &[(i32, i32)]) -> bool {
let geoms: Vec<Geometry<i32>> = coords.iter().map(|&(x, y)| pt(x, y)).collect();
let ids = vec![None; geoms.len()];
spatial_sort_likely_to_help(&build_tile_layer(&geoms, &ids))
}
#[rstest]
#[case::no_features(&[], true)]
#[case::one_feature(&[(100, 100)], true)]
#[case::clustered(&[(0, 0), (100, 100)], true)]
#[case::spread_on_x_only(&[(0, 0), (4000, 100)], true)]
#[case::spread_on_both_axes(&[(0, 0), (4000, 4000)], false)]
#[case::spread_on_y_only(&[(0, 0), (100, 4000)], true)]
#[case::exactly_at_the_coverage_limit(&[(0, 0), (3276, 3276)], true)]
#[case::just_past_the_coverage_limit(&[(0, 0), (3277, 3277)], false)]
#[case::wider_than_the_extent(&[(-9000, -9000), (9000, 9000)], false)]
#[case::wider_than_i32(&[(i32::MIN, i32::MIN), (i32::MAX, i32::MAX)], false)]
fn spatial_sort_help_heuristic(#[case] coords: &[(i32, i32)], #[case] expected: bool) {
assert_eq!(points_help_spatial_sort(coords), expected);
}
#[test]
fn vertexless_features_keep_the_spatial_sort() {
let geoms = [GeoGeom::LineString(LineString(vec![]))];
let layer = build_tile_layer(&geoms, &[None]);
assert!(spatial_sort_likely_to_help(&layer));
}
#[test]
fn curve_params_of_a_vertexless_layer_is_the_degenerate_grid() {
let layer = build_tile_layer(&[GeoGeom::LineString(LineString(vec![]))], &[None]);
assert_eq!(layer.curve_params(), CurveParams { shift: 0, bits: 1 });
}
#[rstest]
#[case::morton(SortStrategy::SpatialMorton, &[0, 1, 2, 3])]
#[case::hilbert(SortStrategy::SpatialHilbert, &[0, 2, 3, 1])]
fn spatial_sort_of_the_unit_square_follows_the_chosen_curve(
#[case] strategy: SortStrategy,
#[case] expected: &[u64],
) {
let geoms = [pt(0, 0), pt(1, 0), pt(0, 1), pt(1, 1)];
let ids: Vec<Option<u64>> = (0..4).map(Some).collect();
let mut layer = build_tile_layer(&geoms, &ids);
let params = layer.curve_params();
assert_eq!(params, CurveParams { shift: 0, bits: 1 });
layer.sort(strategy, params);
let sorted: Vec<Option<u64>> = layer.features().iter().map(TileFeature::id).collect();
assert_eq!(
sorted,
expected.iter().copied().map(Some).collect::<Vec<_>>()
);
}
#[test]
fn unsorted_leaves_the_feature_order_alone() {
let geoms = [pt(9, 9), pt(0, 0), pt(5, 5)];
let ids = [Some(3u64), Some(1), Some(2)];
let mut layer = build_tile_layer(&geoms, &ids);
let params = layer.curve_params();
layer.sort(SortStrategy::Unsorted, params);
let sorted: Vec<Option<u64>> = layer.features().iter().map(TileFeature::id).collect();
assert_eq!(sorted, ids.to_vec());
}
#[test]
fn every_vertexless_geometry_kind_sorts_after_the_one_feature_with_a_vertex() {
let empty_ring = || Polygon::new(LineString(vec![]), vec![]);
let geoms = [
GeoGeom::LineString(LineString(vec![])),
GeoGeom::Polygon(empty_ring()),
GeoGeom::MultiPoint(MultiPoint(vec![])),
GeoGeom::MultiLineString(MultiLineString(vec![])),
GeoGeom::MultiLineString(MultiLineString(vec![LineString(vec![])])),
GeoGeom::MultiPolygon(MultiPolygon(vec![])),
GeoGeom::MultiPolygon(MultiPolygon(vec![empty_ring()])),
GeoGeom::GeometryCollection(GeometryCollection(vec![])),
GeoGeom::GeometryCollection(GeometryCollection(vec![GeoGeom::LineString(LineString(
vec![],
))])),
pt(1, 1),
];
let ids: Vec<Option<u64>> = (0..u64::try_from(geoms.len()).unwrap()).map(Some).collect();
let mut layer = build_tile_layer(&geoms, &ids);
let params = layer.curve_params();
layer.sort(SortStrategy::SpatialMorton, params);
let sorted: Vec<Option<u64>> = layer.features().iter().map(TileFeature::id).collect();
assert_eq!(
sorted,
vec![
Some(9),
Some(0),
Some(1),
Some(2),
Some(3),
Some(4),
Some(5),
Some(6),
Some(7),
Some(8)
]
);
}
#[test]
fn the_spatial_sort_heuristic_ignores_features_without_a_vertex() {
let geoms = [
pt(0, 0),
GeoGeom::LineString(LineString(vec![])),
pt(4000, 4000),
];
let layer = build_tile_layer(&geoms, &[None, None, None]);
assert!(!spatial_sort_likely_to_help(&layer));
}
}