1use geo::CoordsIter as _;
4use geo_types::{Coord, Geometry};
5
6use crate::codecs::hilbert::{hilbert_curve_params_from_bounds, hilbert_sort_key};
7use crate::codecs::morton::morton_sort_key;
8use crate::encoder::model::CurveParams;
9use crate::tile::{TileFeature, TileLayer};
10
11#[derive(Debug, Clone, Copy, Default, PartialEq, Eq, strum::EnumIter, strum::EnumCount)]
17pub enum SortStrategy {
18 #[default]
22 Unsorted,
23
24 SpatialMorton,
31
32 SpatialHilbert,
36
37 Id,
39}
40
41impl TileLayer {
42 #[hotpath::measure]
52 pub fn sort(&mut self, strategy: SortStrategy, params: CurveParams) {
53 match strategy {
54 SortStrategy::SpatialMorton | SortStrategy::SpatialHilbert => {
55 let curve_key = if let SortStrategy::SpatialMorton = strategy {
56 morton_sort_key
57 } else {
58 hilbert_sort_key
59 };
60 self.features_mut().sort_by_cached_key(|f| {
61 first_vertex(f.geometry()).map_or(u64::MAX, |c| u64::from(curve_key(c, params)))
62 });
63 }
64 SortStrategy::Id => {
65 self.features_mut().sort_by_cached_key(TileFeature::id);
66 }
67 SortStrategy::Unsorted => {
68 }
70 }
71 }
72
73 #[hotpath::measure]
77 #[must_use]
78 pub fn curve_params(&self) -> CurveParams {
79 let (min_val, max_val) = self
80 .features()
81 .iter()
82 .flat_map(|f| f.geometry().coords_iter())
83 .fold((i32::MAX, i32::MIN), |(min, max), c| {
84 (min.min(c.x).min(c.y), max.max(c.x).max(c.y))
85 });
86 hilbert_curve_params_from_bounds(min_val, max_val)
87 }
88}
89
90fn first_vertex(geom: &Geometry<i32>) -> Option<Coord<i32>> {
92 match geom {
93 Geometry::<i32>::Point(p) => Some(p.0),
94 Geometry::<i32>::Line(l) => Some(l.start),
95 Geometry::<i32>::LineString(ls) => ls.0.first().copied(),
96 Geometry::<i32>::Polygon(p) => p.exterior().0.first().copied(),
97 Geometry::<i32>::MultiPoint(mp) => mp.0.first().map(|p| p.0),
98 Geometry::<i32>::MultiLineString(mls) => mls.0.first().and_then(|ls| ls.0.first().copied()),
99 Geometry::<i32>::MultiPolygon(mp) => {
100 mp.0.first().and_then(|p| p.exterior().0.first().copied())
101 }
102 Geometry::<i32>::Triangle(t) => Some(t.v1()),
103 Geometry::<i32>::Rect(r) => Some(r.min()),
104 Geometry::<i32>::GeometryCollection(gc) => gc.0.first().and_then(first_vertex),
105 }
106}
107
108pub(crate) fn spatial_sort_likely_to_help(layer: &TileLayer) -> bool {
115 const SPATIAL_HELP_COVERAGE: f64 = 0.8;
116
117 let extent = f64::from(layer.extent().get());
118 if extent <= 0.0 || layer.features().is_empty() {
119 return true;
120 }
121
122 let (min_x, max_x, min_y, max_y) = layer
123 .features()
124 .iter()
125 .filter_map(|f| first_vertex(f.geometry()))
126 .fold(
127 (i32::MAX, i32::MIN, i32::MAX, i32::MIN),
128 |(min_x, max_x, min_y, max_y), Coord::<i32> { x, y }| {
129 (min_x.min(x), max_x.max(x), min_y.min(y), max_y.max(y))
130 },
131 );
132
133 if min_x > max_x || min_y > max_y {
134 return true;
135 }
136
137 let range_x = f64::from(max_x) - f64::from(min_x);
138 let range_y = f64::from(max_y) - f64::from(min_y);
139
140 let spread_x = range_x / extent;
141 let spread_y = range_y / extent;
142
143 !(spread_x > SPATIAL_HELP_COVERAGE && spread_y > SPATIAL_HELP_COVERAGE)
144}
145
146#[cfg(test)]
147mod tests {
148 use geo_types::{
149 Coord, Geometry as GeoGeom, Geometry, GeometryCollection, Line, LineString,
150 MultiLineString, MultiPoint, MultiPolygon, Point, Polygon, Rect, Triangle,
151 };
152 use rstest::rstest;
153
154 use super::spatial_sort_likely_to_help;
155 use crate::decoder::{GeometryType, GeometryValues, RawGeometry};
156 use crate::encoder::model::CurveParams;
157 use crate::encoder::{
158 Codecs, Encoder, EncoderConfig, ExplicitEncoder, IntEncoder, SortStrategy, stage_tile,
159 };
160 use crate::test_helpers::{assert_empty, dec, into_layer01, parser};
161 use crate::tile::{TileFeature, TileLayer};
162 use crate::{Layer, LazyParsed};
163
164 fn pt(x: i32, y: i32) -> Geometry<i32> {
165 GeoGeom::Point(Point::new(x, y))
166 }
167
168 fn ls(coords: &[(i32, i32)]) -> Geometry<i32> {
169 GeoGeom::LineString(LineString::new(
170 coords.iter().map(|&(x, y)| Coord { x, y }).collect(),
171 ))
172 }
173
174 fn poly_square(x0: i32, y0: i32, side: i32) -> Geometry<i32> {
175 let ring = LineString::new(vec![
176 Coord { x: x0, y: y0 },
177 Coord {
178 x: x0 + side,
179 y: y0,
180 },
181 Coord {
182 x: x0 + side,
183 y: y0 + side,
184 },
185 Coord {
186 x: x0,
187 y: y0 + side,
188 },
189 Coord { x: x0, y: y0 },
190 ]);
191 GeoGeom::Polygon(Polygon::new(ring, vec![]))
192 }
193
194 fn roundtrip_geom(decoded: &GeometryValues) -> GeometryValues {
196 let mut enc = Encoder::default();
197 let mut codecs = Codecs::default();
198 decoded
199 .clone()
200 .write_to(&mut enc, &mut codecs)
201 .expect("encode failed");
202 let buf = enc.data().to_vec();
203
204 let parsed = assert_empty(RawGeometry::from_bytes(&buf, &mut parser()));
205 let mut d = dec();
206 let result = LazyParsed::Raw(parsed)
207 .into_parsed(&mut d)
208 .expect("decode failed");
209 assert!(
210 d.consumed() > 0,
211 "decoder should consume bytes after decode"
212 );
213 result
214 }
215
216 fn canonical(geoms: &[Geometry<i32>]) -> GeometryValues {
218 let mut decoded = GeometryValues::default();
219 for g in geoms {
220 decoded.push_geom(g);
221 }
222 roundtrip_geom(&decoded)
223 }
224
225 fn layer_after_sort(geoms: &[Geometry<i32>], ids: &[u64], strategy: SortStrategy) -> TileLayer {
228 let features: Vec<TileFeature> = geoms
229 .iter()
230 .zip(ids.iter())
231 .map(|(g, &id)| TileFeature {
232 id: Some(id),
233 geometry: g.clone(),
234 properties: vec![],
235 #[cfg(feature = "unstable-v2")]
236 m_values: vec![],
237 #[cfg(feature = "unstable-v2")]
238 nested: vec![],
239 })
240 .collect();
241
242 let mut layer = TileLayer::from_parts("test", 4096, vec![], features).unwrap();
243
244 let params = layer.curve_params();
245 layer.sort(strategy, params);
246 layer
247 }
248
249 fn assert_sort_roundtrip(
251 geoms: &[Geometry<i32>],
252 ids: &[u64],
253 strategy: SortStrategy,
254 expected: &[Geometry<i32>],
255 ) {
256 let layer = layer_after_sort(geoms, ids, strategy);
257
258 let mut sorted_decoded = GeometryValues::default();
259 for f in layer.features() {
260 sorted_decoded.push_geom(f.geometry());
261 }
262
263 let after_roundtrip = roundtrip_geom(&sorted_decoded);
264 let expected_canonical = canonical(expected);
265
266 assert_eq!(
267 after_roundtrip, expected_canonical,
268 "\nsorted geometry did not match expected after encode->decode round-trip\
269 \nvector_types after sort: {:?}\
270 \nvector_types expected: {:?}",
271 sorted_decoded.vector_types, expected_canonical.vector_types,
272 );
273 }
274
275 #[test]
278 fn pure_points_id_sort_roundtrip() {
279 assert_sort_roundtrip(
280 &[pt(0, 0), pt(1, 1), pt(2, 2)],
281 &[3, 2, 1],
282 SortStrategy::Id,
283 &[pt(2, 2), pt(1, 1), pt(0, 0)],
284 );
285 }
286
287 #[test]
290 fn pure_linestrings_id_sort_roundtrip() {
291 assert_sort_roundtrip(
292 &[ls(&[(0, 0), (0, 10)]), ls(&[(5, 5), (10, 10)])],
293 &[2, 1],
294 SortStrategy::Id,
295 &[ls(&[(5, 5), (10, 10)]), ls(&[(0, 0), (0, 10)])],
296 );
297 }
298
299 #[test]
302 fn point_line_point_id_sort_to_line_point_point_roundtrip() {
303 assert_sort_roundtrip(
304 &[pt(0, 0), ls(&[(1, 0), (1, 5)]), pt(5, 5)],
305 &[3, 1, 2],
306 SortStrategy::Id,
307 &[ls(&[(1, 0), (1, 5)]), pt(5, 5), pt(0, 0)],
308 );
309 }
310
311 #[test]
312 fn point_line_point_id_sort_to_point_point_line_roundtrip() {
313 assert_sort_roundtrip(
314 &[pt(0, 0), ls(&[(1, 0), (1, 5)]), pt(5, 5)],
315 &[1, 3, 2],
316 SortStrategy::Id,
317 &[pt(0, 0), pt(5, 5), ls(&[(1, 0), (1, 5)])],
318 );
319 }
320
321 #[test]
324 fn point_polygon_point_id_sort_roundtrip() {
325 assert_sort_roundtrip(
326 &[pt(0, 0), poly_square(10, 10, 5), pt(5, 5)],
327 &[2, 1, 3],
328 SortStrategy::Id,
329 &[poly_square(10, 10, 5), pt(0, 0), pt(5, 5)],
330 );
331 }
332
333 #[test]
336 fn point_line_point_morton_sort_roundtrip() {
337 assert_sort_roundtrip(
338 &[pt(2, 0), ls(&[(0, 0), (0, 5)]), pt(1, 0)],
339 &[1, 2, 3],
340 SortStrategy::SpatialMorton,
341 &[ls(&[(0, 0), (0, 5)]), pt(1, 0), pt(2, 0)],
342 );
343 }
344
345 #[test]
348 fn id_sort_already_sorted_is_identity_roundtrip() {
349 let geoms = &[pt(0, 0), ls(&[(1, 0), (1, 5)]), pt(5, 5)];
350 assert_sort_roundtrip(geoms, &[1, 2, 3], SortStrategy::Id, geoms);
351 }
352
353 #[test]
356 fn id_column_co_permuted_with_geometry() {
357 let layer = layer_after_sort(
358 &[pt(0, 0), ls(&[(1, 0), (1, 5)]), pt(5, 5)],
359 &[3, 1, 2],
360 SortStrategy::Id,
361 );
362
363 let ids: Vec<Option<u64>> = layer.features().iter().map(TileFeature::id).collect();
364 assert_eq!(ids, vec![Some(1u64), Some(2), Some(3)]);
365
366 let geom_types: Vec<&str> = layer
368 .features()
369 .iter()
370 .map(|f| GeometryType::try_from(f.geometry()).unwrap().into())
371 .collect();
372 assert_eq!(geom_types, vec!["LineString", "Point", "Point"]);
373 }
374
375 fn build_tile_layer(geoms: &[Geometry<i32>], ids: &[Option<u64>]) -> TileLayer {
377 assert_eq!(geoms.len(), ids.len());
378 TileLayer::from_parts(
379 "test",
380 4096,
381 vec![],
382 geoms
383 .iter()
384 .zip(ids.iter())
385 .map(|(g, &id)| TileFeature {
386 id,
387 geometry: g.clone(),
388 properties: vec![],
389 #[cfg(feature = "unstable-v2")]
390 m_values: vec![],
391 #[cfg(feature = "unstable-v2")]
392 nested: vec![],
393 })
394 .collect(),
395 )
396 .unwrap()
397 }
398
399 fn sort_encode_decode(tile: TileLayer, sort: SortStrategy) -> TileLayer {
402 let enc_cfg = EncoderConfig::default();
403 let enc = Encoder::with_explicit(enc_cfg, ExplicitEncoder::for_id(IntEncoder::varint()));
404 let mut codecs = Codecs::default();
405 let enc = stage_tile(tile, sort, false, enc_cfg.tessellate())
406 .encode_into(enc, &mut codecs)
407 .expect("encode failed");
408
409 let buf = enc.into_layer_bytes().expect("into_layer_bytes failed");
411
412 let mut p = parser();
413 let layer_back = assert_empty(Layer::from_bytes(&buf, &mut p));
414 assert!(p.reserved() > 0, "parser should reserve bytes after parse");
415
416 let layer01 = into_layer01(layer_back);
417
418 let mut d = dec();
419 let tile = layer01.into_tile(&mut d).expect("decode after sort failed");
420 assert!(
421 d.consumed() > 0,
422 "decoder should consume bytes after decode"
423 );
424 tile
425 }
426
427 fn vertices_from_source(source: &TileLayer) -> Vec<i32> {
429 let mut geom = GeometryValues::default();
430 for f in source.features() {
431 geom.push_geom(f.geometry());
432 }
433 geom.vertices().unwrap_or_default().to_vec()
434 }
435
436 #[test]
437 fn test_shared_morton_shift() {
438 let tile = build_tile_layer(&[pt(0, -10), pt(-10, 0)], &[Some(1), Some(2)]);
445 let source = sort_encode_decode(tile, SortStrategy::SpatialMorton);
446
447 let verts = vertices_from_source(&source);
448 assert_eq!(verts, vec![0, -10, -10, 0]);
449 }
450
451 #[test]
452 fn test_id_sort_nulls_first() {
453 let tile = build_tile_layer(&[pt(2, 2), pt(1, 1), pt(0, 0)], &[Some(10), None, Some(5)]);
454 let source = sort_encode_decode(tile, SortStrategy::Id);
455
456 let ids: Vec<Option<u64>> = source.features().iter().map(TileFeature::id).collect();
457 assert_eq!(ids, vec![None, Some(5), Some(10)]);
459
460 let verts = vertices_from_source(&source);
461 assert_eq!(verts, vec![1, 1, 0, 0, 2, 2]);
463 }
464
465 #[test]
466 fn test_mixed_geometry_morton_sort() {
467 let tile = build_tile_layer(
475 &[pt(2, 0), ls(&[(0, 0), (0, 5)]), pt(1, 0)],
476 &[Some(1), Some(2), Some(3)],
477 );
478 let source = sort_encode_decode(tile, SortStrategy::SpatialMorton);
479
480 let types: Vec<_> = source
481 .features()
482 .iter()
483 .map(|f| GeometryType::try_from(f.geometry()).unwrap())
484 .collect();
485
486 assert_eq!(
487 types,
488 vec![
489 GeometryType::LineString,
490 GeometryType::Point,
491 GeometryType::Point
492 ]
493 );
494
495 let verts = vertices_from_source(&source);
496 assert_eq!(verts, vec![0, 0, 0, 5, 1, 0, 2, 0]);
498 }
499
500 #[test]
501 fn id_sort_separates_the_two_largest_ids() {
502 let layer = layer_after_sort(
503 &[pt(0, 0), pt(1, 1)],
504 &[u64::MAX, u64::MAX - 1],
505 SortStrategy::Id,
506 );
507
508 let ids: Vec<Option<u64>> = layer.features().iter().map(TileFeature::id).collect();
509 assert_eq!(ids, vec![Some(u64::MAX - 1), Some(u64::MAX)]);
510 }
511
512 #[test]
513 fn id_sort_puts_the_missing_id_before_zero() {
514 let layer = build_tile_layer(&[pt(0, 0), pt(1, 1)], &[Some(0), None]);
515 let mut layer = layer;
516 layer.sort(SortStrategy::Id, CurveParams { shift: 0, bits: 1 });
517
518 let ids: Vec<Option<u64>> = layer.features().iter().map(TileFeature::id).collect();
519 assert_eq!(ids, vec![None, Some(0)]);
520 }
521
522 #[test]
523 fn spatial_sort_reads_the_first_vertex_of_every_geometry_kind() {
524 let at = |k: i32| Coord { x: k, y: k };
525 let geoms = [
526 GeoGeom::Rect(Rect::new(at(0), at(20))),
527 GeoGeom::Triangle(Triangle::new(at(1), at(11), at(21))),
528 GeoGeom::MultiPolygon(MultiPolygon(vec![Polygon::new(
529 LineString(vec![at(2), at(12), at(22), at(2)]),
530 vec![],
531 )])),
532 GeoGeom::MultiLineString(MultiLineString(vec![LineString(vec![at(3), at(13)])])),
533 GeoGeom::MultiPoint(MultiPoint(vec![Point::from(at(4)), Point::from(at(14))])),
534 GeoGeom::Polygon(Polygon::new(
535 LineString(vec![at(5), at(15), at(25), at(5)]),
536 vec![],
537 )),
538 GeoGeom::LineString(LineString(vec![at(6), at(16)])),
539 GeoGeom::Line(Line::new(at(7), at(17))),
540 GeoGeom::Point(Point::from(at(8))),
541 GeoGeom::GeometryCollection(GeometryCollection(vec![GeoGeom::Point(Point::from(at(
542 9,
543 )))])),
544 GeoGeom::LineString(LineString(vec![])),
545 ];
546 let ids: Vec<Option<u64>> = (0..u64::try_from(geoms.len()).unwrap()).map(Some).collect();
547
548 let mut layer = build_tile_layer(&geoms, &ids);
549 let params = layer.curve_params();
550 layer.sort(SortStrategy::SpatialMorton, params);
553
554 let sorted: Vec<Option<u64>> = layer.features().iter().map(TileFeature::id).collect();
555 assert_eq!(sorted, ids);
556 }
557
558 fn points_help_spatial_sort(coords: &[(i32, i32)]) -> bool {
559 let geoms: Vec<Geometry<i32>> = coords.iter().map(|&(x, y)| pt(x, y)).collect();
560 let ids = vec![None; geoms.len()];
561 spatial_sort_likely_to_help(&build_tile_layer(&geoms, &ids))
562 }
563
564 #[rstest]
565 #[case::no_features(&[], true)]
566 #[case::one_feature(&[(100, 100)], true)]
567 #[case::clustered(&[(0, 0), (100, 100)], true)]
568 #[case::spread_on_x_only(&[(0, 0), (4000, 100)], true)]
569 #[case::spread_on_both_axes(&[(0, 0), (4000, 4000)], false)]
570 #[case::spread_on_y_only(&[(0, 0), (100, 4000)], true)]
571 #[case::exactly_at_the_coverage_limit(&[(0, 0), (3276, 3276)], true)]
572 #[case::just_past_the_coverage_limit(&[(0, 0), (3277, 3277)], false)]
573 #[case::wider_than_the_extent(&[(-9000, -9000), (9000, 9000)], false)]
574 #[case::wider_than_i32(&[(i32::MIN, i32::MIN), (i32::MAX, i32::MAX)], false)]
575 fn spatial_sort_help_heuristic(#[case] coords: &[(i32, i32)], #[case] expected: bool) {
576 assert_eq!(points_help_spatial_sort(coords), expected);
577 }
578
579 #[test]
580 fn vertexless_features_keep_the_spatial_sort() {
581 let geoms = [GeoGeom::LineString(LineString(vec![]))];
582 let layer = build_tile_layer(&geoms, &[None]);
583 assert!(spatial_sort_likely_to_help(&layer));
584 }
585
586 #[test]
587 fn curve_params_of_a_vertexless_layer_is_the_degenerate_grid() {
588 let layer = build_tile_layer(&[GeoGeom::LineString(LineString(vec![]))], &[None]);
589 assert_eq!(layer.curve_params(), CurveParams { shift: 0, bits: 1 });
590 }
591
592 #[rstest]
593 #[case::morton(SortStrategy::SpatialMorton, &[0, 1, 2, 3])]
594 #[case::hilbert(SortStrategy::SpatialHilbert, &[0, 2, 3, 1])]
595 fn spatial_sort_of_the_unit_square_follows_the_chosen_curve(
596 #[case] strategy: SortStrategy,
597 #[case] expected: &[u64],
598 ) {
599 let geoms = [pt(0, 0), pt(1, 0), pt(0, 1), pt(1, 1)];
600 let ids: Vec<Option<u64>> = (0..4).map(Some).collect();
601
602 let mut layer = build_tile_layer(&geoms, &ids);
603 let params = layer.curve_params();
604 assert_eq!(params, CurveParams { shift: 0, bits: 1 });
605 layer.sort(strategy, params);
606
607 let sorted: Vec<Option<u64>> = layer.features().iter().map(TileFeature::id).collect();
608 assert_eq!(
609 sorted,
610 expected.iter().copied().map(Some).collect::<Vec<_>>()
611 );
612 }
613
614 #[test]
615 fn unsorted_leaves_the_feature_order_alone() {
616 let geoms = [pt(9, 9), pt(0, 0), pt(5, 5)];
617 let ids = [Some(3u64), Some(1), Some(2)];
618
619 let mut layer = build_tile_layer(&geoms, &ids);
620 let params = layer.curve_params();
621 layer.sort(SortStrategy::Unsorted, params);
622
623 let sorted: Vec<Option<u64>> = layer.features().iter().map(TileFeature::id).collect();
624 assert_eq!(sorted, ids.to_vec());
625 }
626
627 #[test]
628 fn every_vertexless_geometry_kind_sorts_after_the_one_feature_with_a_vertex() {
629 let empty_ring = || Polygon::new(LineString(vec![]), vec![]);
630 let geoms = [
631 GeoGeom::LineString(LineString(vec![])),
632 GeoGeom::Polygon(empty_ring()),
633 GeoGeom::MultiPoint(MultiPoint(vec![])),
634 GeoGeom::MultiLineString(MultiLineString(vec![])),
635 GeoGeom::MultiLineString(MultiLineString(vec![LineString(vec![])])),
636 GeoGeom::MultiPolygon(MultiPolygon(vec![])),
637 GeoGeom::MultiPolygon(MultiPolygon(vec![empty_ring()])),
638 GeoGeom::GeometryCollection(GeometryCollection(vec![])),
639 GeoGeom::GeometryCollection(GeometryCollection(vec![GeoGeom::LineString(LineString(
640 vec![],
641 ))])),
642 pt(1, 1),
643 ];
644 let ids: Vec<Option<u64>> = (0..u64::try_from(geoms.len()).unwrap()).map(Some).collect();
645
646 let mut layer = build_tile_layer(&geoms, &ids);
647 let params = layer.curve_params();
648 layer.sort(SortStrategy::SpatialMorton, params);
649
650 let sorted: Vec<Option<u64>> = layer.features().iter().map(TileFeature::id).collect();
651 assert_eq!(
652 sorted,
653 vec![
654 Some(9),
655 Some(0),
656 Some(1),
657 Some(2),
658 Some(3),
659 Some(4),
660 Some(5),
661 Some(6),
662 Some(7),
663 Some(8)
664 ]
665 );
666 }
667
668 #[test]
669 fn the_spatial_sort_heuristic_ignores_features_without_a_vertex() {
670 let geoms = [
671 pt(0, 0),
672 GeoGeom::LineString(LineString(vec![])),
673 pt(4000, 4000),
674 ];
675 let layer = build_tile_layer(&geoms, &[None, None, None]);
676 assert!(!spatial_sort_likely_to_help(&layer));
677 }
678}