1#[cfg(feature = "gridpoints-proj")]
2use proj::Proj;
3
4#[allow(unused_imports)]
5use crate::GribError;
6use crate::{GridPointIndexIterator, def::grib2::template::param_set::ScanningMode};
7
8pub(crate) fn evenly_spaced_longitudes(
9 start_microdegree: u32,
10 end_microdegree: u32,
11 div: usize,
12 angle_units: f32,
13 scanning_mode: ScanningMode,
14) -> Vec<f32> {
15 let is_consistent =
16 !((end_microdegree > start_microdegree) ^ scanning_mode.scans_positively_for_i());
17
18 let (start, end) = (start_microdegree as f32, end_microdegree as f32);
19 let (start, end) = if is_consistent {
20 (start, end)
21 } else if start_microdegree > end_microdegree {
22 (start, end + 360_000_000_f32)
23 } else {
24 (start + 360_000_000_f32, end)
25 };
26
27 let lons = evenly_spaced_degrees(start, end, div, angle_units);
28
29 if is_consistent {
30 lons
31 } else {
32 lons.into_iter()
33 .map(|x| if x < 360.0 { x } else { x - 360.0 })
34 .collect()
35 }
36}
37
38pub(crate) fn evenly_spaced_degrees(
39 start_microdegree: f32,
40 end_microdegree: f32,
41 div: usize,
42 angle_units: f32,
43) -> Vec<f32> {
44 let delta = (end_microdegree - start_microdegree) / div as f32;
45 (0..=div)
46 .map(move |x| (start_microdegree + x as f32 * delta) * angle_units)
47 .collect()
48}
49
50#[derive(Clone)]
52pub struct RegularGridIterator {
53 lat: Vec<f32>,
54 lon: Vec<f32>,
55 ij: GridPointIndexIterator,
56}
57
58impl RegularGridIterator {
59 pub(crate) fn new(lat: Vec<f32>, lon: Vec<f32>, ij: GridPointIndexIterator) -> Self {
60 Self { lat, lon, ij }
61 }
62}
63
64impl Iterator for RegularGridIterator {
65 type Item = (f32, f32);
66
67 fn next(&mut self) -> Option<Self::Item> {
68 let (i, j) = self.ij.next()?;
69 Some((self.lat[j], self.lon[i]))
70 }
71
72 fn size_hint(&self) -> (usize, Option<usize>) {
73 self.ij.size_hint()
74 }
75}
76
77#[derive(Clone)]
78#[cfg(feature = "gridpoints-proj")]
79pub struct ProjectionLatLonIterator {
80 xy: std::vec::IntoIter<(f64, f64)>,
81}
82
83#[cfg(feature = "gridpoints-proj")]
84impl Iterator for ProjectionLatLonIterator {
85 type Item = (f32, f32);
86
87 fn next(&mut self) -> Option<Self::Item> {
88 self.xy
89 .next()
90 .map(|(lon, lat)| (lat.to_degrees() as f32, lon.to_degrees() as f32))
91 }
92
93 fn size_hint(&self) -> (usize, Option<usize>) {
94 self.xy.size_hint()
95 }
96}
97
98#[cfg(feature = "gridpoints-proj")]
99pub(crate) fn latlons_from_projection_with_first_point_and_delta(
100 proj_def: &str,
101 first_point_latlon_in_degrees: (f64, f64),
102 delta_in_meters: (f64, f64),
103 indices: GridPointIndexIterator,
104) -> Result<ProjectionLatLonIterator, GribError> {
105 let projection = Proj::new(proj_def)?;
106 let (first_point_lat, first_point_lon) = first_point_latlon_in_degrees;
107 let (first_corner_x, first_corner_y) = projection.project(
108 (first_point_lon.to_radians(), first_point_lat.to_radians()),
109 false,
110 )?;
111
112 let (dx, dy) = delta_in_meters;
113 let mut xy = indices
114 .map(|(i, j)| {
115 (
116 first_corner_x + dx * i as f64,
117 first_corner_y + dy * j as f64,
118 )
119 })
120 .collect::<Vec<_>>();
121
122 projection.project_array(&mut xy, true)?;
123
124 Ok(ProjectionLatLonIterator { xy: xy.into_iter() })
125}
126
127#[cfg(feature = "gridpoints-proj")]
128pub(crate) fn latlons_from_projection_with_first_point_and_last_point(
129 proj_def: &str,
130 first_point_latlon_in_degrees: (f64, f64),
131 last_point_latlon_in_degrees: (f64, f64),
132 (ni, nj): (usize, usize),
133 indices: GridPointIndexIterator,
134) -> Result<std::vec::IntoIter<(f32, f32)>, GribError> {
135 let projection = Proj::new(proj_def)?;
136 let (first_point_lat, first_point_lon) = first_point_latlon_in_degrees;
137 let (first_corner_x, first_corner_y) = projection.project(
138 (first_point_lon.to_radians(), first_point_lat.to_radians()),
139 false,
140 )?;
141 let (last_point_lat, last_point_lon) = last_point_latlon_in_degrees;
142 let (last_corner_x, last_corner_y) = projection.project(
143 (last_point_lon.to_radians(), last_point_lat.to_radians()),
144 false,
145 )?;
146
147 let dx = (last_corner_x - first_corner_x) / (ni - 1) as f64;
148 let dy = (last_corner_y - first_corner_y) / (nj - 1) as f64;
149 let mut xy = indices
150 .map(|(i, j)| {
151 (
152 first_corner_x + dx * i as f64,
153 first_corner_y + dy * j as f64,
154 )
155 })
156 .collect::<Vec<_>>();
157
158 let lonlat = projection.project_array(&mut xy, true)?;
159 let latlon = lonlat
160 .iter_mut()
161 .map(|(lon, lat)| (lat.to_degrees() as f32, lon.to_degrees() as f32))
162 .collect::<Vec<_>>();
163
164 Ok(latlon.into_iter())
165}
166
167pub(crate) fn normalize_latlon((lat, lon): (f32, f32)) -> (f32, f32) {
168 let lon = (lon + 540.) % 360. - 180.;
169 (lat, lon)
170}
171
172#[cfg(feature = "gridpoints-proj")]
173impl From<proj::ProjCreateError> for GribError {
174 fn from(e: proj::ProjCreateError) -> Self {
175 Self::Unknown(e.to_string())
176 }
177}
178
179#[cfg(feature = "gridpoints-proj")]
180impl From<proj::ProjError> for GribError {
181 fn from(e: proj::ProjError) -> Self {
182 Self::Unknown(e.to_string())
183 }
184}
185
186#[cfg(test)]
187pub(crate) mod test_helpers {
188 macro_rules! assert_almost_eq {
189 ($a1:expr, $a2:expr, $d:expr) => {
190 if $a1 - $a2 > $d {
191 panic!(
192 "assertion a1 - a2 <= delta failed\n a1 - a2: {} - {}\n delta: {}",
193 $a1, $a2, $d
194 );
195 } else if $a2 - $a1 > $d {
196 panic!(
197 "assertion a2 - a1 <= delta failed\n a2 - a1: {} - {}\n delta: {}",
198 $a2, $a1, $d
199 );
200 }
201 };
202 }
203 pub(crate) use assert_almost_eq;
204
205 macro_rules! test_assert_almost_eq_do_not_panic {
206 ($((
207 $name:ident,
208 $a1:expr,
209 $a2:expr,
210 $d:expr
211 ),)*) => ($(
212 #[test]
213 fn $name() {
214 assert_almost_eq!($a1, $a2, $d)
215 }
216 )*);
217 }
218
219 test_assert_almost_eq_do_not_panic! {
220 (assert_almost_eq_does_not_panic_for_positive_lt_positive, 1.01, 1.02, 0.1),
221 (assert_almost_eq_does_not_panic_for_positive_gt_positive, 1.02, 1.01, 0.1),
222 (assert_almost_eq_does_not_panic_for_negative_lt_negative, -1.02, -1.01, 0.1),
223 (assert_almost_eq_does_not_panic_for_negative_gt_negative, -1.01, -1.02, 0.1),
224 (assert_almost_eq_does_not_panic_for_positive_negative, 0.01, -0.01, 0.1),
225 (assert_almost_eq_does_not_panic_for_negative_positive, -0.01, 0.01, 0.1),
226 }
227
228 macro_rules! test_assert_almost_eq_panic {
229 ($((
230 $name:ident,
231 $a1:expr,
232 $a2:expr,
233 $d:expr,
234 $message:expr
235 ),)*) => ($(
236 #[test]
237 #[should_panic(expected = $message)]
238 fn $name() {
239 assert_almost_eq!($a1, $a2, $d)
240 }
241 )*);
242 }
243
244 test_assert_almost_eq_panic! {
245 (
246 assert_almost_eq_panics_for_positive_lt_positive, 1.01, 1.02, 0.001,
247 " a2 - a1: 1.02 - 1.01\n delta: 0.001"
248 ),
249 (
250 assert_almost_eq_panics_for_positive_gt_positive, 1.02, 1.01, 0.001,
251 " a1 - a2: 1.02 - 1.01\n delta: 0.001"
252 ),
253 (
254 assert_almost_eq_panics_for_negative_lt_negative, -1.02, -1.01, 0.001,
255 " a2 - a1: -1.01 - -1.02\n delta: 0.001"
256 ),
257 (
258 assert_almost_eq_panics_for_negative_gt_negative, -1.01, -1.02, 0.001,
259 " a1 - a2: -1.01 - -1.02\n delta: 0.001"
260 ),
261 (
262 assert_almost_eq_panics_for_positive_negative, 0.01, -0.01, 0.001,
263 " a1 - a2: 0.01 - -0.01\n delta: 0.001"
264 ),
265 (
266 assert_almost_eq_panics_for_negative_positive, -0.01, 0.01, 0.001,
267 " a2 - a1: 0.01 - -0.01\n delta: 0.001"
268 ),
269 (
270 assert_almost_eq_panic_message_containing_trailing_zeros, -0.0100, 0.0100, 0.0010,
271 " a2 - a1: 0.01 - -0.01\n delta: 0.001"
272 ),
273 }
274
275 #[allow(dead_code)]
276 pub(crate) fn assert_coord_almost_eq((x1, y1): (f32, f32), (x2, y2): (f32, f32), delta: f32) {
277 assert_almost_eq!(x1, x2, delta);
278 assert_almost_eq!(y1, y2, delta);
279 }
280}
281
282#[cfg(test)]
283mod tests {
284 use super::*;
285
286 macro_rules! test_lat_lon_grid_iter {
287 ($(($name:ident, $scanning_mode:expr, $expected:expr),)*) => ($(
288 #[test]
289 fn $name() {
290 let lat = (0..3).into_iter().map(|i| i as f32).collect::<Vec<_>>();
291 let lon = (10..12).into_iter().map(|i| i as f32).collect::<Vec<_>>();
292 let scanning_mode = ScanningMode($scanning_mode);
293 let ij = GridPointIndexIterator::new((lon.len(), lat.len()), scanning_mode).unwrap();
294 let actual = RegularGridIterator::new(lat, lon, ij).collect::<Vec<_>>();
295 assert_eq!(actual, $expected);
296 }
297 )*);
298 }
299
300 test_lat_lon_grid_iter! {
301 (
302 lat_lon_grid_iter_with_scanning_mode_0b00000000,
303 0b00000000,
304 vec![
305 (0., 10.),
306 (0., 11.),
307 (1., 10.),
308 (1., 11.),
309 (2., 10.),
310 (2., 11.),
311 ]
312 ),
313 (
314 lat_lon_grid_iter_with_scanning_mode_0b00100000,
315 0b00100000,
316 vec![
317 (0., 10.),
318 (1., 10.),
319 (2., 10.),
320 (0., 11.),
321 (1., 11.),
322 (2., 11.),
323 ]
324 ),
325 (
326 lat_lon_grid_iter_with_scanning_mode_0b00010000,
327 0b00010000,
328 vec![
329 (0., 10.),
330 (0., 11.),
331 (1., 11.),
332 (1., 10.),
333 (2., 10.),
334 (2., 11.),
335 ]
336 ),
337 (
338 lat_lon_grid_iter_with_scanning_mode_0b00110000,
339 0b00110000,
340 vec![
341 (0., 10.),
342 (1., 10.),
343 (2., 10.),
344 (2., 11.),
345 (1., 11.),
346 (0., 11.),
347 ]
348 ),
349 }
350
351 #[test]
352 fn lat_lon_grid_iterator_size_hint() {
353 let lat = (0..3).map(|i| i as f32).collect::<Vec<_>>();
354 let lon = (10..12).map(|i| i as f32).collect::<Vec<_>>();
355 let scanning_mode = ScanningMode(0b00000000);
356 let ij = GridPointIndexIterator::new((lon.len(), lat.len()), scanning_mode).unwrap();
357 let mut iter = RegularGridIterator::new(lat, lon, ij);
358
359 assert_eq!(iter.size_hint(), (6, Some(6)));
360 let _ = iter.next();
361 assert_eq!(iter.size_hint(), (5, Some(5)));
362 }
363
364 #[test]
365 fn latlon_normalization() {
366 assert_eq!(normalize_latlon((90., -180.)), (90., -180.));
367 assert_eq!(normalize_latlon((90., 0.)), (90., 0.));
368 assert_eq!(normalize_latlon((90., 179.)), (90., 179.));
369 assert_eq!(normalize_latlon((90., 180.)), (90., -180.));
370 assert_eq!(normalize_latlon((90., 360.)), (90., 0.));
371 assert_eq!(normalize_latlon((90., 540.)), (90., -180.));
372 }
373}