1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
use geo_types::{
Coordinate,
Rect,
};
use h3ron::index::Index;
use crate::{error::Error, sphere::{
area_linearring,
area_rect,
}, transform::Transform, AxisOrder};
use h3ron::{H3_MIN_RESOLUTION, H3_MAX_RESOLUTION};
pub enum ResolutionSearchMode {
MinDiff,
SmallerThanPixel,
}
pub fn nearest_h3_resolution(shape: &[usize], transform: &Transform, axis_order: &AxisOrder, search_mode: ResolutionSearchMode) -> Result<u8, Error> {
if shape.len() != 2 {
return Err(Error::UnsupportedArrayShape);
}
if shape[0] == 0 || shape[1] == 0 {
return Err(Error::EmptyArray);
}
let bbox_array = Rect::new(
transform * &Coordinate::from((0.0_f64, 0.0_f64)),
transform * &Coordinate::from((
(shape[axis_order.x_axis()] - 1) as f64,
(shape[axis_order.y_axis()] - 1) as f64
)),
);
let area_pixel = area_rect(&bbox_array) / (shape[axis_order.x_axis()] * shape[axis_order.y_axis()]) as f64;
let center_of_array = bbox_array.center();
let mut nearest_h3_res = 0;
let mut area_difference = None;
for h3_res in H3_MIN_RESOLUTION..=H3_MAX_RESOLUTION {
let area_h3_index = area_linearring(Index::from_coordinate(¢er_of_array, h3_res)
.polygon()
.exterior());
match search_mode {
ResolutionSearchMode::SmallerThanPixel => if area_h3_index <= area_pixel {
nearest_h3_res = h3_res;
break;
}
ResolutionSearchMode::MinDiff => {
let new_area_difference = if area_h3_index > area_pixel {
area_h3_index - area_pixel
} else {
area_pixel - area_h3_index
};
if let Some(old_area_difference) = area_difference {
if old_area_difference < new_area_difference {
nearest_h3_res = h3_res - 1;
break;
} else {
area_difference = Some(new_area_difference);
}
} else {
area_difference = Some(new_area_difference);
}
}
}
}
Ok(nearest_h3_res)
}
#[cfg(test)]
mod tests {
use crate::resolution::{nearest_h3_resolution, ResolutionSearchMode};
use crate::transform::Transform;
use crate::AxisOrder;
#[test]
fn test_nearest_h3_resolution() {
let gt = Transform::from_rasterio(&[
0.0011965049999999992, 0.0, 8.11377, 0.0, -0.001215135, 49.40792
]);
let h3_res1 = nearest_h3_resolution(&[2000_usize, 2000_usize], >, &AxisOrder::YX, ResolutionSearchMode::MinDiff).unwrap();
assert_eq!(h3_res1, 10);
let h3_res2 = nearest_h3_resolution(&[2000_usize, 2000_usize], >, &AxisOrder::YX, ResolutionSearchMode::SmallerThanPixel).unwrap();
assert_eq!(h3_res2, 11);
}
}