utiles_cover/
cover_geotypes.rs

1//! Tile cover for geojson object(s) based on mapbox's tile-cover alg
2#![expect(
3    clippy::cast_possible_truncation,
4    clippy::cast_possible_wrap,
5    clippy::cast_sign_loss
6)]
7use crate::Result;
8use std::collections::{BTreeMap, HashSet};
9use utiles_core::{Tile, lnglat2tile_frac, simplify, tile, utile};
10
11#[derive(Debug, Clone, Copy, PartialEq, Eq)]
12pub struct GeoTypesCoverOptions {
13    pub zoom: u8,
14    pub minzoom: Option<u8>,
15}
16
17impl From<u8> for GeoTypesCoverOptions {
18    fn from(zoom: u8) -> Self {
19        Self {
20            zoom,
21            minzoom: None,
22        }
23    }
24}
25
26#[expect(clippy::cast_precision_loss)]
27fn line_string_cover(
28    tiles_set: &mut HashSet<Tile>,
29    ls: &geo_types::LineString<f64>,
30    maxzoom: u8,
31    mut ring: Option<&mut Vec<(u32, u32)>>,
32) {
33    if ls.0.len() < 2 {
34        return;
35    }
36    let mut prev_x: Option<i64> = None;
37    let mut prev_y: Option<i64> = None;
38    let mut y_value: Option<i64> = None;
39    let minxy = (1u32 << maxzoom) - 1; // Maximum valid tile index at this zoom level
40
41    // for segment in coords.windows(2) {
42    for segment in ls.lines() {
43        // let start_coord = segment
44        let start_coord = segment.start;
45        let stop_coord = segment.end;
46
47        let (x0f, y0f, _) = lnglat2tile_frac(start_coord.x, start_coord.y, maxzoom);
48        let (x1f, y1f, _) = lnglat2tile_frac(stop_coord.x, stop_coord.y, maxzoom);
49
50        let dx = x1f - x0f;
51        let dy = y1f - y0f;
52
53        // Directly check for zero movement
54        if dx == 0.0 && dy == 0.0 {
55            continue;
56        }
57        let sx = dx.signum() as i64;
58        let sy = dy.signum() as i64;
59
60        let mut x = x0f.floor() as i64;
61        let mut y = y0f.floor() as i64;
62        y_value = Some(y);
63
64        let tdx = if dx == 0.0 {
65            f64::INFINITY
66        } else {
67            (sx as f64 / dx).abs()
68        };
69        let tdy = if dy == 0.0 {
70            f64::INFINITY
71        } else {
72            (sy as f64 / dy).abs()
73        };
74
75        let mut t_max_x = if dx == 0.0 {
76            f64::INFINITY
77        } else {
78            ((if dx > 0.0 { 1.0 } else { 0.0 } + x as f64 - x0f) / dx).abs()
79        };
80        let mut t_max_y = if dy == 0.0 {
81            f64::INFINITY
82        } else {
83            ((if dy > 0.0 { 1.0 } else { 0.0 } + y as f64 - y0f) / dy).abs()
84        };
85
86        // Add initial tile
87        if prev_x != Some(x) || prev_y != Some(y) {
88            let tile = utile!(x as u32, y as u32, maxzoom);
89            tiles_set.insert(tile);
90            if let Some(ring) = &mut ring
91                && prev_y != Some(y)
92            {
93                ring.push((x as u32, y as u32));
94            }
95            prev_x = Some(x);
96            prev_y = Some(y);
97        }
98
99        // the MAX number of tiles to check to dx.abs() + dy.abs()
100        let mut max_it = (dx.abs() + dy.abs()) as i64;
101        while (t_max_x < 1.0 || t_max_y < 1.0) && max_it >= 0 {
102            if t_max_x < t_max_y {
103                t_max_x += tdx;
104                x += sx;
105            } else {
106                t_max_y += tdy;
107                y += sy;
108            }
109
110            // Sanity check if x or y is outside valid tile ranges
111            if x > i64::from(minxy) + 1 || y > i64::from(minxy) + 1 {
112                break;
113            }
114
115            if prev_x != Some(x) || prev_y != Some(y) {
116                let tile = utile!(x as u32, y as u32, maxzoom);
117                tiles_set.insert(tile);
118                if let Some(ring) = &mut ring
119                    && prev_y != Some(y)
120                {
121                    ring.push((x as u32, y as u32));
122                }
123                prev_x = Some(x);
124                prev_y = Some(y);
125            }
126
127            max_it -= 1; // Decrement the number of steps remaining
128        }
129    }
130
131    // adjust the ring if needed
132    if let Some(ring) = &mut ring
133        && let (Some(first_ring), Some(y_value)) = (ring.first(), y_value)
134        && y_value == i64::from(first_ring.1)
135    {
136        ring.pop();
137    }
138}
139
140fn polygon_cover(
141    tiles_set: &mut HashSet<Tile>,
142    geom: &geo_types::Polygon<f64>,
143    zoom: u8,
144) {
145    // Collect all x-intersections per scanline (y)
146    let mut scanlines: BTreeMap<u32, Vec<u32>> = BTreeMap::new();
147
148    let rings = vec![geom.exterior()] // Start with the exterior ring
149        .into_iter()
150        .chain(geom.interiors().iter()); // Add any interior rings
151
152    for ring_pts in rings {
153        // First collect the polygon boundary into `boundary` as tile-edge samples
154        let mut boundary: Vec<(u32, u32)> = Vec::new();
155        line_string_cover(tiles_set, ring_pts, zoom, Some(&mut boundary));
156        if boundary.is_empty() {
157            continue;
158        }
159
160        // Iterate each edge, including the closing edge from last→first
161        let iter = boundary
162            .windows(2)
163            .map(|w| (w[0], w[1]))
164            .chain(std::iter::once((boundary[boundary.len() - 1], boundary[0])));
165
166        for ((x0f, y0f), (x1f, y1f)) in iter {
167            let x0 = x0f as i32;
168            let y0 = y0f as i32;
169            let x1 = x1f as i32;
170            let y1 = y1f as i32;
171
172            // skip fully horizontal edges
173            if y0 == y1 {
174                continue;
175            }
176
177            // y ranges over the scanlines this edge crosses
178            let (ymin, ymax) = (y0.min(y1), y0.max(y1));
179            let dx = x1 - x0;
180            let dy = y1 - y0;
181
182            for y in ymin..ymax {
183                // parametric t along the edge at integer y
184                let t = f64::from(y - y0) / f64::from(dy);
185                let x = t.mul_add(f64::from(dx), f64::from(x0)).floor() as u32;
186                scanlines.entry(y as u32).or_default().push(x);
187            }
188        }
189    }
190
191    // Fill in between pairs of intersections on each scanline
192    for (y, mut xs) in scanlines {
193        xs.sort_unstable();
194        // take (x_start, x_end) from each pair of sorted intersections
195        for pair in xs.chunks(2).filter(|c| c.len() == 2) {
196            let x_start = pair[0];
197            let x_end = pair[1];
198            tiles_set.extend((x_start..=x_end).map(|x| Tile::new(x, y, zoom)));
199        }
200    }
201}
202
203fn gt_geometry2tiles(
204    geom: &geo_types::Geometry,
205    opts: GeoTypesCoverOptions,
206) -> Result<HashSet<Tile>> {
207    let mut tilescoverage: HashSet<Tile> = HashSet::new();
208
209    match geom {
210        geo_types::Geometry::Point(pt) => {
211            let tile = tile(pt.x(), pt.y(), opts.zoom, None)?;
212            tilescoverage.insert(tile);
213        }
214        geo_types::Geometry::MultiPoint(pts) => {
215            let it = pts
216                .iter()
217                .filter_map(|pt| tile(pt.x(), pt.y(), opts.zoom, None).ok());
218            tilescoverage.extend(it);
219        }
220        geo_types::Geometry::Line(ln) => {
221            let ls = geo_types::LineString::from(ln);
222            line_string_cover(&mut tilescoverage, &ls, opts.zoom, None);
223        }
224        geo_types::Geometry::LineString(ls) => {
225            line_string_cover(&mut tilescoverage, ls, opts.zoom, None);
226        }
227        geo_types::Geometry::MultiLineString(mls) => {
228            mls.iter().for_each(|ls| {
229                line_string_cover(&mut tilescoverage, ls, opts.zoom, None);
230            });
231        }
232        geo_types::Geometry::Polygon(poly) => {
233            polygon_cover(&mut tilescoverage, poly, opts.zoom);
234        }
235        geo_types::Geometry::MultiPolygon(mpoly) => {
236            for poly in mpoly.iter() {
237                polygon_cover(&mut tilescoverage, poly, opts.zoom);
238            }
239        }
240        geo_types::Geometry::GeometryCollection(gjcoll) => {
241            for g in gjcoll {
242                tilescoverage.extend(gt_geometry2tiles(g, opts)?);
243            }
244        }
245        geo_types::Geometry::Rect(r) => {
246            let poly = geo_types::Polygon::from(*r);
247            polygon_cover(&mut tilescoverage, &poly, opts.zoom);
248        }
249        geo_types::Geometry::Triangle(t) => {
250            let poly = geo_types::Polygon::from(*t);
251            polygon_cover(&mut tilescoverage, &poly, opts.zoom);
252        }
253    }
254
255    match opts.minzoom {
256        Some(z) => {
257            let cov = simplify(&tilescoverage, Some(z));
258            Ok(cov)
259        }
260        None => Ok(tilescoverage),
261    }
262}
263
264/// Convert a `geo_types::Geometry` to a set of tiles at the specified zoom level.
265///
266/// # Errors
267///
268/// If the `GeoJSON` object is invalid or if the conversion fails due to
269/// projecting coordinate issues.
270pub fn geometry2tiles<T>(geom: &geo_types::Geometry, opts: T) -> Result<HashSet<Tile>>
271where
272    T: Into<GeoTypesCoverOptions>,
273{
274    let opts = opts.into();
275    gt_geometry2tiles(geom, opts)
276}
utiles_cover/cover_geotypes.rs

utiles_cover/
cover_geotypes.rs
