rasterh3/
resolution.rs

1use geo::{AffineOps, AffineTransform};
2use geo_types::Rect;
3use h3o::{LatLng, Resolution};
4
5use crate::{error::Error, sphere::AreaOnSphere, AxisOrder};
6
7#[derive(Copy, Clone)]
8pub enum ResolutionSearchMode {
9    /// Chose the H3 resolution where the difference in the area of a pixel and the h3index is
10    /// as small as possible.
11    MinDiff,
12
13    /// Chose the H3 resolution where the area of the h3index is smaller than the area of a pixel.
14    SmallerThanPixel,
15}
16
17impl ResolutionSearchMode {
18    /// Find the H3 resolution closed to the size of a pixel in an array
19    /// of the given shape with the given transform.
20    pub fn nearest_h3_resolution(
21        &self,
22        shape: [usize; 2],
23        transform: &AffineTransform<f64>,
24        axis_order: &AxisOrder,
25    ) -> Result<Resolution, Error> {
26        if shape[0] == 0 || shape[1] == 0 {
27            return Err(Error::EmptyArray);
28        }
29        let bbox_array = Rect::new(
30            (0.0_f64, 0.0_f64),
31            (
32                (shape[axis_order.x_axis()] - 1) as f64,
33                (shape[axis_order.y_axis()] - 1) as f64,
34            ),
35        )
36        .affine_transform(transform);
37        let area_pixel = bbox_array.area_on_sphere_m2()
38            / (shape[axis_order.x_axis()] * shape[axis_order.y_axis()]) as f64;
39        let center_of_array: LatLng = bbox_array.center().try_into()?;
40
41        let mut nearest_h3_res = Resolution::Zero;
42        let mut area_difference = None;
43        for h3_res in Resolution::range(Resolution::Zero, Resolution::Fifteen) {
44            let area_h3_index = center_of_array.to_cell(h3_res).area_m2();
45
46            match self {
47                Self::SmallerThanPixel => {
48                    if area_h3_index <= area_pixel {
49                        nearest_h3_res = h3_res;
50                        break;
51                    }
52                }
53
54                Self::MinDiff => {
55                    let new_area_difference = if area_h3_index > area_pixel {
56                        area_h3_index - area_pixel
57                    } else {
58                        area_pixel - area_h3_index
59                    };
60                    if let Some(old_area_difference) = area_difference {
61                        if old_area_difference < new_area_difference {
62                            nearest_h3_res = h3_res.pred().unwrap_or(Resolution::Zero);
63                            break;
64                        } else {
65                            area_difference = Some(new_area_difference);
66                        }
67                    } else {
68                        area_difference = Some(new_area_difference);
69                    }
70                }
71            }
72        }
73
74        Ok(nearest_h3_res)
75    }
76}
77
78#[cfg(test)]
79mod tests {
80    use h3o::Resolution;
81
82    use crate::resolution::ResolutionSearchMode;
83    use crate::AxisOrder;
84
85    #[test]
86    fn test_nearest_h3_resolution() {
87        // transform of the included r.tiff
88        let gt = crate::transform::from_rasterio(&[
89            0.0011965049999999992,
90            0.0,
91            8.11377,
92            0.0,
93            -0.001215135,
94            49.40792,
95        ]);
96        let h3_res1 = ResolutionSearchMode::MinDiff
97            .nearest_h3_resolution([2000_usize, 2000_usize], &gt, &AxisOrder::YX)
98            .unwrap();
99        assert_eq!(h3_res1, Resolution::Ten); // TODO: validate
100
101        let h3_res2 = ResolutionSearchMode::SmallerThanPixel
102            .nearest_h3_resolution([2000_usize, 2000_usize], &gt, &AxisOrder::YX)
103            .unwrap();
104        assert_eq!(h3_res2, Resolution::Eleven); // TODO: validate
105    }
106}
rasterh3/resolution.rs

rasterh3/
resolution.rs
