gistools/readers/tile/
raster.rs

1use super::{
2    GetRasterTileValue, S2TileMetadata, TileFetcher, TileMetadata, TileReader, WMTileMetadata,
3};
4use crate::{
5    geometry::{Source, merc_to_ll, xyz_to_bbox},
6    parsers::FeatureReader,
7};
8use alloc::{format, string::String, vec, vec::Vec};
9use core::marker::PhantomData;
10use image::RgbaImage;
11use s2_tilejson::{Metadata, Scheme, UnknownMetadata};
12use s2json::{
13    BBox, BBox3D, Face, VectorFeature, VectorFeatureType, VectorGeometry, VectorMultiPoint,
14    VectorPoint,
15};
16use std::{
17    fs,
18    path::{Path, PathBuf},
19};
20
21/// # Raster Tiles Fetcher
22///
23/// ## Description
24/// Read an entire archive of raster tiles, where the max zoom data is iterated upon
25///
26/// Supports reading either RGB(A) data, RGB(A) encoded elevation data, or build your own structure.
27///
28/// Implements the [`FeatureReader`] and [`TileFetcher`] traits
29///
30/// ## Usage
31///
32/// The methods you have access to:
33/// - [`RasterTileFetcher::new`]: Create a new RasterTileFetcher
34/// - [`RasterTileFetcher::get_metadata`]: Get the metadata of the tileset
35/// - [`RasterTileFetcher::has_tile_wm`]: Check if it is WM tile
36/// - [`RasterTileFetcher::has_tile_s2`]: Check if it is S2 tile
37/// - [`RasterTileFetcher::get_tile_wm`]: Get an WM tile
38/// - [`RasterTileFetcher::get_tile_s2`]: Get an S2 tile
39/// - [`RasterTileFetcher::iter`]: Iterate over the tiles
40///
41/// ```rust
42/// use gistools::{parsers::{RGBA, FeatureReader}, readers::{TileFetcher, RasterTileFetcher}};
43/// use std::path::PathBuf;
44///
45/// let mut path = PathBuf::from(env!("CARGO_MANIFEST_DIR"));
46/// path.push("tests/readers/tile/fixtures/wm/satellite");
47///
48/// // read the RGBA data of each tile. Each pixel is stored as a VectorPoint
49/// let reader = RasterTileFetcher::<RGBA>::new(path, Some(1));
50/// let tiles: Vec<_> = reader.iter().collect();
51/// assert_eq!(tiles.len(), 4);
52/// ```
53///
54/// ## Links
55/// - <https://satakagi.github.io/mapsForWebWS2020-docs/QuadTreeCompositeTilingAndVectorTileStandard.html>
56/// - <https://cesium.com/blog/2015/04/07/quadtree-cheatseet/>
57#[derive(Debug, Clone)]
58pub struct RasterTileFetcher<D: Clone + Default + GetRasterTileValue> {
59    path: PathBuf,
60    threshold: Option<u8>,
61    metadata: Metadata,
62    phantom: PhantomData<D>,
63}
64impl<D: Clone + Default + GetRasterTileValue> TileFetcher<D, D, RasterTileReader<D>>
65    for RasterTileFetcher<D>
66{
67    fn new<P: AsRef<Path>>(path: P, threshold: Option<u8>) -> RasterTileFetcher<D> {
68        let path = path.as_ref().to_path_buf();
69        let metadata_path = path.join("metadata.json");
70        let meta_string: String = fs::read_to_string(&metadata_path).unwrap();
71        let unknown_meta: UnknownMetadata = serde_json::from_str(&meta_string).unwrap();
72
73        RasterTileFetcher {
74            path,
75            threshold,
76            metadata: unknown_meta.to_metadata(),
77            phantom: PhantomData,
78        }
79    }
80    fn get_metadata(&self) -> &Metadata {
81        &self.metadata
82    }
83
84    fn has_tile_wm(&self, zoom: u8, x: u32, y: u32) -> bool {
85        let Metadata { extension, .. } = self.get_metadata();
86        let tile_path = self.path.join(format!("{zoom}/{x}/{y}.{extension}"));
87
88        tile_path.exists()
89    }
90
91    fn has_tile_s2(&self, face: Face, zoom: u8, x: u32, y: u32) -> bool {
92        let Metadata { extension, .. } = self.get_metadata();
93        let tile_path =
94            self.path.join(format!("{}/{}/{}/{}.{}", u8::from(face), zoom, x, y, extension));
95
96        tile_path.exists()
97    }
98
99    fn get_tile_wm(&self, zoom: u8, x: u32, y: u32) -> RasterTileReader<D> {
100        RasterTileReader::new(self.path.clone(), self.get_metadata(), 0.into(), zoom, x, y, false)
101    }
102
103    fn get_tile_s2(&self, face: Face, zoom: u8, x: u32, y: u32) -> RasterTileReader<D> {
104        RasterTileReader::new(self.path.clone(), self.get_metadata(), face, zoom, x, y, true)
105    }
106}
107impl<D: Clone + Default + GetRasterTileValue> FeatureReader<TileMetadata, D, D>
108    for RasterTileFetcher<D>
109{
110    type FeatureIterator<'a>
111        = RasterIterator<'a, D>
112    where
113        Self: 'a;
114
115    fn iter(&self) -> Self::FeatureIterator<'_> {
116        let Metadata { minzoom, maxzoom, .. } = self.get_metadata();
117        let mut stack = vec![(0.into(), 0, 0, 0)];
118        if self.is_s2() {
119            for face in [Face::Face1, Face::Face2, Face::Face3, Face::Face4, Face::Face5] {
120                stack.push((face, 0, 0, 0));
121            }
122        }
123        let threshold = self.threshold.unwrap_or(*maxzoom);
124        RasterIterator {
125            container: self,
126            stack,
127            minzoom: *minzoom,
128            threshold,
129            pool_size: 1,
130            thread_id: 0,
131            index: 0,
132        }
133    }
134
135    #[cfg(feature = "std")]
136    fn par_iter(&self, pool_size: usize, thread_id: usize) -> Self::FeatureIterator<'_> {
137        let Metadata { minzoom, maxzoom, .. } = self.get_metadata();
138        let mut stack = vec![(0.into(), 0, 0, 0)];
139        if self.is_s2() {
140            for face in [Face::Face1, Face::Face2, Face::Face3, Face::Face4, Face::Face5] {
141                stack.push((face, 0, 0, 0));
142            }
143        }
144        let threshold = self.threshold.unwrap_or(*maxzoom);
145        RasterIterator {
146            container: self,
147            stack,
148            minzoom: *minzoom,
149            threshold,
150            pool_size: pool_size as u64,
151            thread_id: thread_id as u64,
152            index: 0,
153        }
154    }
155}
156
157/// Iterator for the S2 Raster Tile Fetcher
158#[derive(Debug)]
159pub struct RasterIterator<'a, D: Clone + Default + GetRasterTileValue> {
160    container: &'a RasterTileFetcher<D>,
161    stack: Vec<(Face, u8, u32, u32)>,
162    minzoom: u8,
163    threshold: u8,
164    pool_size: u64,
165    thread_id: u64,
166    index: u64,
167}
168impl<D: Clone + Default + GetRasterTileValue> Iterator for RasterIterator<'_, D> {
169    type Item = VectorFeature<TileMetadata, D, D>;
170    fn next(&mut self) -> Option<Self::Item> {
171        let is_s2 = self.container.is_s2();
172        while let Some((face, zoom, x, y)) = self.stack.pop() {
173            // if zoom not reached yet, push children and continue
174            let has_tile = if is_s2 {
175                self.container.has_tile_s2(face, zoom, x, y)
176            } else {
177                self.container.has_tile_wm(zoom, x, y)
178            };
179            if zoom < self.minzoom || (zoom != self.threshold && has_tile) {
180                self.stack.extend(vec![
181                    (face, zoom + 1, x * 2, y * 2),
182                    (face, zoom + 1, x * 2 + 1, y * 2),
183                    (face, zoom + 1, x * 2, y * 2 + 1),
184                    (face, zoom + 1, x * 2 + 1, y * 2 + 1),
185                ]);
186                continue;
187            } else if zoom == self.threshold && has_tile {
188                let idx = self.index;
189                self.index += 1;
190                if self.pool_size > 1 && idx % self.pool_size != self.thread_id {
191                    continue; // skip, belongs to another thread
192                }
193                let tile = if is_s2 {
194                    self.container.get_tile_s2(face, zoom, x, y)
195                } else {
196                    self.container.get_tile_wm(zoom, x, y)
197                };
198                return Some(tile.build_feature());
199            }
200        }
201        None
202    }
203}
204
205/// Raster Tile Reader
206#[derive(Debug)]
207pub struct RasterTileReader<D: Clone + Default + GetRasterTileValue> {
208    /// Tile Metadata
209    pub metadata: TileMetadata,
210    /// Tile Image
211    pub image: RgbaImage,
212    /// If the tile is S2 or WM
213    pub is_s2: bool,
214    /// Tile's scheme is TMS
215    pub tms_style: bool,
216    /// Help D pass
217    _phantom: PhantomData<D>,
218}
219impl<D: Clone + Default + GetRasterTileValue> TileReader<D, D> for RasterTileReader<D> {
220    fn new(
221        path: PathBuf,
222        metadata: &Metadata,
223        face: Face,
224        zoom: u8,
225        x: u32,
226        y: u32,
227        is_s2: bool,
228    ) -> Self {
229        let Metadata { extension, scheme, .. } = metadata;
230        let tile_path = if is_s2 {
231            path.join(format!("{}/{}/{}/{}.{}", u8::from(face), zoom, x, y, extension))
232        } else {
233            path.join(format!("{zoom}/{x}/{y}.{extension}"))
234        };
235        let image = image::open(&tile_path).unwrap().to_rgba8();
236        let metadata = if is_s2 {
237            TileMetadata::S2(S2TileMetadata { face, zoom, x, y })
238        } else {
239            TileMetadata::WM(super::WMTileMetadata { zoom, x, y })
240        };
241
242        RasterTileReader {
243            metadata,
244            image,
245            is_s2,
246            tms_style: *scheme == Scheme::Tms,
247            _phantom: PhantomData,
248        }
249    }
250
251    fn build_feature(&self) -> VectorFeature<TileMetadata, D, D> {
252        match &self.metadata {
253            TileMetadata::S2(s2) => self.build_feature_s2(s2),
254            TileMetadata::WM(wm) => self.build_feature_wm(wm),
255        }
256    }
257}
258impl<D: Clone + Default + GetRasterTileValue> RasterTileReader<D> {
259    fn build_feature_s2(&self, metadata: &S2TileMetadata) -> VectorFeature<TileMetadata, D, D> {
260        let S2TileMetadata { x, y, zoom, face } = metadata;
261        let tile_size = self.image.width() as usize;
262        let mut bbox = BBox3D::default();
263        // Get the bounding box of the tile in lon-lat
264        let BBox { left: min_s, bottom: min_t, right: max_s, top: max_t } =
265            BBox::from_st_zoom(*x as f64, *y as f64, *zoom);
266        let s_step = (max_s - min_s) / (tile_size as f64);
267        let t_step = (max_t - min_t) / (tile_size as f64);
268        let mut coordinates: VectorMultiPoint<D> = vec![];
269
270        for py in 0..tile_size {
271            let y = min_s + ((py as f64) + 0.5) * t_step; // Center of the row
272            for px in 0..tile_size {
273                let x = min_t + ((px as f64) + 0.5) * s_step; // Center of the column
274                let pixel = self.image.get_pixel(px as u32, py as u32);
275                let m_value =
276                    D::get_raster_tile_value(pixel.0[0], pixel.0[1], pixel.0[2], Some(pixel.0[3]));
277
278                let vp = VectorPoint::new_xy(x, y, Some(m_value));
279                bbox.extend_from_point(&vp);
280                coordinates.push(vp);
281            }
282        }
283
284        VectorFeature {
285            _type: VectorFeatureType::S2Feature,
286            face: *face,
287            metadata: Some(self.metadata.clone()),
288            geometry: VectorGeometry::new_multipoint(coordinates, Some(bbox)),
289            ..Default::default()
290        }
291    }
292
293    fn build_feature_wm(&self, metadata: &WMTileMetadata) -> VectorFeature<TileMetadata, D, D> {
294        let WMTileMetadata { x, y, zoom } = metadata;
295        let tile_size = self.image.width() as usize;
296        let mut bbox = BBox3D::default();
297        // Get the bounding box of the tile in lon-lat
298        let (west, south, east, north) = xyz_to_bbox(
299            *x,
300            *y,
301            *zoom,
302            Some(self.tms_style),
303            Some(Source::Google),
304            Some(tile_size as u16),
305        );
306        let x_step = (east - west) / (tile_size as f64);
307        let y_step = (north - south) / (tile_size as f64);
308        let mut coordinates: VectorMultiPoint<D> = vec![];
309
310        for py in 0..tile_size {
311            let py_f64 = py as f64;
312            let y_pos = north - (py_f64 + 0.5) * y_step; // Center of the row
313            for px in 0..tile_size {
314                let px_f64 = px as f64;
315                let x_pos = west + (px_f64 + 0.5) * x_step; // Center of the column
316                let (lon, lat) = merc_to_ll((x_pos, y_pos));
317                let pixel = self.image.get_pixel(px as u32, py as u32);
318                let m_value =
319                    D::get_raster_tile_value(pixel.0[0], pixel.0[1], pixel.0[2], Some(pixel.0[3]));
320
321                let vp = VectorPoint::new_xy(lon, lat, Some(m_value));
322                bbox.extend_from_point(&vp);
323                coordinates.push(vp);
324            }
325        }
326
327        VectorFeature {
328            _type: VectorFeatureType::VectorFeature,
329            metadata: Some(self.metadata.clone()),
330            geometry: VectorGeometry::new_multipoint(coordinates, Some(bbox)),
331            ..Default::default()
332        }
333    }
334}
gistools/readers/tile/raster.rs

gistools/readers/tile/
raster.rs
