gistools/data_structures/
point_grid.rs

1use super::{IndexPoint, PointIndex};
2use crate::{
3    data_store::{KV, KVStore, Vector, VectorStore},
4    geometry::{LonLat, S2CellId, S2Point},
5    parsers::{FeatureReader, RGBA},
6    util::{
7        GetInterpolateValue, Interpolatable, InterpolationFunction, InterpolationMethod,
8        get_interpolation,
9    },
10};
11use alloc::{collections::BTreeSet, fmt::Debug, string::String, vec, vec::Vec};
12use libm::{floor, fmin, log2};
13use s2json::{BBox, Face, GetXY, JSONCollection, Projection};
14use serde::{Deserialize, Serialize, de::DeserializeOwned};
15
16/// Options for grid clustering
17#[allow(unpredictable_function_pointer_comparisons)]
18#[derive(Debug, Default, Clone, Serialize, Deserialize, PartialEq)]
19pub struct GridOptions<V: Clone + Debug> {
20    /// projection to use
21    pub projection: Option<Projection>,
22    /// Name of the layer to build when requesting a tile
23    pub layer_name: Option<String>,
24    /// min zoom to generate clusters on
25    pub minzoom: Option<u8>,
26    /// max zoom level to cluster the points on
27    pub maxzoom: Option<u8>,
28    /// Used by cell search to specify the type of interpolation to use.
29    /// The recommendation is IDW as you want to prioritize closest data points. [default: 'idw']
30    #[serde(skip)]
31    pub maxzoom_interpolation: Option<InterpolationMethod>,
32    /// Used by cell search to specify the type of interpolation to use.
33    /// From experimentation, lanczos is a fast algorithm that maintains the quality of the data
34    /// [default: 'lanczos']
35    #[serde(skip)]
36    pub interpolation: Option<InterpolationMethod>,
37    /// Used by cell search to specify the interpolation function to use [default: the Default of the Value type]
38    #[serde(skip)]
39    pub get_value: Option<GetInterpolateValue<IndexPoint<V>, V>>,
40    /// Grid size, assumed pixel ratio.
41    pub grid_size: Option<u32>,
42    /// Used by the cell search to specify the tile buffer size in pixels. [default: 0]
43    pub buffer_size: Option<u32>,
44}
45
46/// An export of the data as a grid
47#[derive(Debug)]
48pub struct TileGrid<'a, V> {
49    /// name of the layer
50    pub name: String,
51    /// size of the grid including the buffer
52    pub size: u32,
53    /// flattened array of number or RGBA.
54    /// The size of the array is grid_size * grid_size
55    /// Access the position as `grid_size * y + x`
56    pub data: &'a Vec<V>,
57}
58
59/// A Point Grid that uses locally allocated data to store the data
60pub type LocalPointGrid<V = RGBA> =
61    PointGrid<V, Vector<S2CellId, IndexPoint<V>>, KV<S2CellId, Vec<V>>>;
62
63/// # Grid Cluster
64///
65/// ## Description
66/// A cluster store to build grid data of grid_size x grid_size. The resultant tiles are filled.
67/// Useful for building raster tiles or other grid like data (temperature, precipitation, wind, etc).
68#[derive(Debug)]
69pub struct PointGrid<
70    V: Interpolatable + Serialize + DeserializeOwned + Debug = RGBA,
71    S: VectorStore<S2CellId, IndexPoint<V>> = Vector<S2CellId, IndexPoint<V>>,
72    G: KVStore<S2CellId, Vec<V>> = KV<S2CellId, Vec<V>>,
73> {
74    projection: Projection,
75    layer_name: String,
76    minzoom: u8,
77    maxzoom: u8,
78    buffer_size: u32,
79    maxzoom_interpolation: InterpolationFunction<S2Point, IndexPoint<V>, V>,
80    interpolation: InterpolationFunction<S2Point, IndexPoint<V>, V>,
81    get_value: GetInterpolateValue<IndexPoint<V>, V>,
82    grid_size: u32, // a default is a 512x512 pixel tile
83    point_index: PointIndex<V, S>,
84    grid_tile_store: G,
85}
86impl<
87    V: Interpolatable + Serialize + DeserializeOwned + Debug,
88    S: VectorStore<S2CellId, IndexPoint<V>>,
89    G: KVStore<S2CellId, Vec<V>>,
90> PointGrid<V, S, G>
91{
92    /// Create a new point grid
93    pub fn new(options: Option<GridOptions<V>>) -> Self {
94        let options = options.unwrap_or_default();
95        let maxzoom_interpolation =
96            options.maxzoom_interpolation.unwrap_or(InterpolationMethod::IDW);
97        let interpolation = options.interpolation.unwrap_or(InterpolationMethod::Lanczos);
98        let projection = options.projection.unwrap_or(Projection::S2);
99        PointGrid {
100            projection,
101            maxzoom_interpolation: get_interpolation(maxzoom_interpolation),
102            interpolation: get_interpolation(interpolation),
103            get_value: options.get_value.unwrap_or(|v| v.data),
104            buffer_size: options.buffer_size.unwrap_or(0),
105            minzoom: options.minzoom.unwrap_or(0),
106            maxzoom: options.maxzoom.unwrap_or(16),
107            layer_name: options.layer_name.unwrap_or(String::from("default")),
108            grid_size: options.grid_size.unwrap_or(512),
109            point_index: PointIndex::new(None, Some(projection)),
110            grid_tile_store: G::new(None),
111        }
112    }
113
114    // /// Insert a cell with the point and its corresponding data to the index
115    // pub fn insert(&mut self, cell: S2CellId, point: S2Point, data: Option<V>) {
116    //     self.point_index.insert(cell, point, data);
117    // }
118
119    /// Add a lon-lat pair as with any shape
120    pub fn insert_point<T: GetXY>(&mut self, point: T, data: Option<V>) {
121        self.point_index.insert_point(point, data);
122    }
123
124    /// Add a lon-lat pair to the cluster
125    pub fn insert_lon_lat(&mut self, ll: LonLat<V>) {
126        self.point_index.insert_lon_lat(ll);
127    }
128
129    /// Insert an STPoint to the index
130    pub fn insert_face_st(&mut self, face: Face, s: f64, t: f64, data: V) {
131        self.point_index.insert_face_st(face, s, t, data);
132    }
133
134    /// Add all points from a reader. It will try to use the M-value first, but if it doesn't exist
135    /// it will use the feature properties data
136    pub fn insert_reader<L: Clone, F: FeatureReader<L, V, V>>(&mut self, reader: &F) {
137        self.point_index.insert_reader::<L, F>(reader);
138    }
139
140    /// Add a vector feature. It will try to use the M-value first, but if it doesn't exist
141    /// it will use the feature properties data
142    pub fn insert_feature<L: Clone>(&mut self, data: JSONCollection<L, V, V>) {
143        self.point_index.insert_feature::<L>(data);
144    }
145
146    ///  Get the point data as a grid of a tile
147    pub fn get_tile(&self, id: S2CellId) -> Option<TileGrid<'_, V>> {
148        let Self { layer_name, grid_size, buffer_size, .. } = self;
149        let data = self.grid_tile_store.get(id);
150
151        if let Some(data) = data {
152            return Some(TileGrid {
153                name: layer_name.clone(),
154                size: grid_size + buffer_size * 2,
155                data,
156            });
157        }
158
159        None
160    }
161
162    /// Build the grid cluster tiles
163    pub fn build_clusters(&mut self) {
164        // build tiles at maxzoom
165        let mut parents = self.cluster_maxzoom();
166        // work upwards, take the 4 children and merge them
167        for zoom in (self.minzoom..self.maxzoom).rev() {
168            parents = self.custer_zoom(zoom, &parents);
169        }
170    }
171
172    /// Using the point index, build grids at maxzoom by doing searches for each gridpoint.
173    /// return the parent cells
174    fn cluster_maxzoom(&mut self) -> BTreeSet<S2CellId> {
175        let mut parents = BTreeSet::<S2CellId>::new();
176        let grid_length = (self.grid_size + self.buffer_size * 2) as usize;
177        let buffer_size_f64 = self.buffer_size as f64;
178        let grid_size_f64: f64 = self.grid_size as f64;
179        // if the grid is 512 x 512, log2 is 9, meaning the quadtree must split 9 times to analyze
180        // each individual pixel. Make sure we don't dive past 30 levels as that's the limit of the spec.
181        let zoom_grid_level =
182            fmin(self.maxzoom as f64 + floor(log2(grid_size_f64)) - 1., 30.) as u8;
183
184        self.point_index.sort();
185        for (cell, _) in self.point_index.iter() {
186            let maxzoom_id = cell.parent(Some(self.maxzoom)); // idParent(cell, maxzoom);
187            // if maxzoom_id grid tile already exists, skip
188            if self.grid_tile_store.has(maxzoom_id) {
189                continue;
190            }
191            // prep variables and grid result
192            let face = cell.face();
193            let BBox { left: s_min, bottom: t_min, right: s_max, top: t_max } =
194                maxzoom_id.bounds_st(Some(self.maxzoom));
195            let s_pixel = (s_max - s_min) / grid_size_f64;
196            let t_pixel = (t_max - t_min) / grid_size_f64;
197            let s_start = s_min - s_pixel * buffer_size_f64;
198            let t_start = t_min - t_pixel * buffer_size_f64;
199            let mut grid: Vec<V> = vec![V::default(); grid_length * grid_length];
200            // iterate through the grid and do searches for each position. Interpolate the data to the
201            // position and store the result in the grid.
202            for y in 0..grid_length {
203                for x in 0..grid_length {
204                    let t = t_start + (y as f64) * t_pixel;
205                    let mut s = s_start + (x as f64) * s_pixel;
206                    if self.projection == Projection::WG {
207                        s = (s + 1.) % 1.;
208                    } // ensure within 0-1 range via wrapping to the other side
209                    // search for points within a reasonable cell size
210                    let mut grid_level_search = zoom_grid_level;
211                    let mut point_shapes: Vec<IndexPoint<V>>;
212                    let fst_point = S2CellId::from_face_st(face, s, t);
213                    let mut st_cell = fst_point.parent(Some(zoom_grid_level));
214                    loop {
215                        point_shapes = self
216                            .point_index
217                            .search_range(st_cell, None, None)
218                            .into_iter()
219                            .map(|(_, point)| point.clone())
220                            .collect();
221                        if !point_shapes.is_empty()
222                            || grid_level_search <= zoom_grid_level - 3
223                            || grid_level_search == 0
224                        {
225                            break;
226                        }
227                        grid_level_search -= 1;
228                        st_cell = st_cell.parent(Some(grid_level_search));
229                    }
230                    if point_shapes.is_empty() {
231                        continue;
232                    }
233                    let point: S2Point = fst_point.to_point();
234                    grid[y * grid_length + x] =
235                        (self.maxzoom_interpolation)(&point, &point_shapes, self.get_value);
236                }
237            }
238            // store the grid and add the parent cell for future upscaling
239            self.grid_tile_store.set(maxzoom_id, grid);
240            if self.maxzoom != 0 {
241                parents.insert(maxzoom_id.parent(Some(self.maxzoom - 1)));
242            }
243        }
244
245        parents
246    }
247
248    /// Build the parent cells. We simply search for the children of the cell and merge/downsample.
249    ///
250    /// ## Parameters
251    /// - `zoom`: the current zoom we are upscaling to
252    /// - `cells`: the cells to build grids for
253    ///
254    /// ## Returns
255    /// returns the parent cells for the next round of upscaling
256    fn custer_zoom(&mut self, zoom: u8, cells: &BTreeSet<S2CellId>) -> BTreeSet<S2CellId> {
257        let mut parents = BTreeSet::<S2CellId>::new();
258        let grid_length = (self.grid_size + self.buffer_size * 2) as usize;
259        let half_grid_length = grid_length / 2;
260
261        for cell in cells {
262            let mut grid: Vec<V> = vec![V::default(); grid_length * grid_length];
263            let (face, cell_zoom, i, j) = cell.to_face_ij();
264            let [bl_id, br_id, tl_id, tr_id] = S2CellId::children_ij(face, cell_zoom, i, j);
265            // for each child, downsample into the result grid
266            self.downsample_grid(bl_id, &mut grid, 0, 0);
267            self.downsample_grid(br_id, &mut grid, half_grid_length, 0);
268            self.downsample_grid(tl_id, &mut grid, 0, half_grid_length);
269            self.downsample_grid(tr_id, &mut grid, half_grid_length, half_grid_length);
270            // store the grid and add the parent cell for future upscaling
271            self.grid_tile_store.set(*cell, grid);
272            if zoom != 0 {
273                parents.insert(cell.parent(Some(zoom - 1)));
274            }
275        }
276
277        parents
278    }
279
280    /// Upscale a grid into the target grid at x,y position
281    fn downsample_grid(&self, cell_id: S2CellId, target: &mut [V], x: usize, y: usize) {
282        let grid = if let Some(grid) = self.grid_tile_store.get(cell_id) {
283            grid
284        } else {
285            return;
286        };
287
288        let Self { grid_size, buffer_size, .. } = self;
289        let grid_length = (grid_size + buffer_size * 2) as usize;
290        let half_grid_length = grid_length / 2;
291        let mid_point = S2Point::new(0.5, 0.5, 0.);
292        let null_val = V::default();
293
294        for j in 0..half_grid_length {
295            for i in 0..half_grid_length {
296                let base_j = j * 2;
297                let base_i = i * 2;
298
299                let source_points: Vec<IndexPoint<V>> = [
300                    ((0.0, 0.0), base_j * grid_length + base_i),
301                    ((1.0, 0.0), base_j * grid_length + base_i + 1),
302                    ((0.0, 1.0), (base_j + 1) * grid_length + base_i),
303                    ((1.0, 1.0), (base_j + 1) * grid_length + base_i + 1),
304                ]
305                .into_iter()
306                .filter_map(|((x, y), idx)| {
307                    let value = grid[idx];
308                    if value != null_val {
309                        Some(IndexPoint::new(S2Point::new(x, y, 0.0), value))
310                    } else {
311                        None
312                    }
313                })
314                .collect();
315
316                if source_points.is_empty() {
317                    continue;
318                }
319
320                target[(j + y) * grid_length + (i + x)] =
321                    (self.interpolation)(&mid_point, &source_points, self.get_value);
322            }
323        }
324    }
325}
gistools/data_structures/point_grid.rs

gistools/data_structures/
point_grid.rs
