grib/grid/
rotated_ll.rs

1use super::GridPointIndexIterator;
2use crate::{
3    GridPointIndex, LatLons,
4    def::grib2::template::{Template3_1, param_set::Rotation},
5    error::GribError,
6    grid::helpers::RegularGridIterator,
7};
8
9impl crate::GridShortName for Template3_1 {
10    fn short_name(&self) -> &'static str {
11        "rotated_ll"
12    }
13}
14
15impl GridPointIndex for Template3_1 {
16    fn grid_shape(&self) -> (usize, usize) {
17        self.rotated.grid_shape()
18    }
19
20    fn scanning_mode(&self) -> &crate::def::grib2::template::param_set::ScanningMode {
21        self.rotated.scanning_mode()
22    }
23
24    fn ij(&self) -> Result<GridPointIndexIterator, GribError> {
25        self.rotated.ij()
26    }
27}
28
29impl LatLons for Template3_1 {
30    type Iter<'a>
31        = Unrotate<RegularGridIterator>
32    where
33        Self: 'a;
34
35    fn latlons<'a>(&'a self) -> Result<Self::Iter<'a>, GribError> {
36        let iter = Unrotate::new(self.rotated.latlons()?, &self.rotation);
37        Ok(iter)
38    }
39}
40
41#[derive(Clone)]
42pub struct Unrotate<I> {
43    latlons: I,
44    sinφp: f32,
45    cosφp: f32,
46    λp: f32,
47    gamma: f32,
48}
49
50impl<I> Unrotate<I> {
51    fn new(latlons: I, rot: &Rotation) -> Self {
52        let φp = (rot.south_pole_lat as f32 * 1e-6).to_radians();
53        let λp = (rot.south_pole_lon as f32 * 1e-6).to_radians();
54        let gamma = rot.rot_angle.to_radians();
55
56        // south pole to north pole
57        let φp = -φp;
58        let λp = λp + std::f32::consts::PI;
59
60        let (sinφp, cosφp) = φp.sin_cos();
61        Self {
62            latlons,
63            sinφp,
64            cosφp,
65            λp,
66            gamma,
67        }
68    }
69}
70
71impl<I> Iterator for Unrotate<I>
72where
73    I: Iterator<Item = (f32, f32)>,
74{
75    type Item = (f32, f32);
76
77    fn next(&mut self) -> Option<Self::Item> {
78        let (lat, lon) = self.latlons.next()?;
79        let λr = lon.to_radians();
80        let φr = lat.to_radians();
81
82        let λr = λr - self.gamma;
83
84        let (sinφr, cosφr) = φr.sin_cos();
85        let (sinλr, cosλr) = λr.sin_cos();
86
87        let sinφ = self.sinφp * sinφr + self.cosφp * cosφr * cosλr;
88        let φ = sinφ.asin();
89
90        let y = cosφr * sinλr;
91        let x = self.cosφp * sinφr - self.sinφp * cosφr * cosλr;
92        let λ = self.λp - y.atan2(x);
93
94        let latlon = (φ.to_degrees(), λ.to_degrees());
95        Some(latlon)
96    }
97
98    fn size_hint(&self) -> (usize, Option<usize>) {
99        self.latlons.size_hint()
100    }
101}
102
103#[cfg(test)]
104mod tests {
105
106    use super::*;
107
108    macro_rules! test_rotation{
109        ($(($name:ident, $rot:expr, $input:expr, $expected:expr),)*) => ($(
110            #[test]
111            fn $name() {
112                let rot = $rot;
113                let latlons = vec![$input];
114                let mut iter = Unrotate::new(latlons.into_iter(), &rot);
115                let actual = iter.next().unwrap();
116                let expected = $expected;
117                dbg!(actual);
118                dbg!(expected);
119                assert!((actual.0 - expected.0).abs() < 0.00001);
120                assert!((actual.1 - expected.1).abs() < 0.001);
121            }
122        )*);
123    }
124
125    test_rotation! {
126        (
127            no_rotation,
128            Rotation {
129                south_pole_lat: -90000000,
130                south_pole_lon: 0,
131                rot_angle: 0.,
132            },
133            (-12.302501_f32, 345.178780_f32),
134            (-12.302501_f32, 345.178780_f32)
135        ),
136        (
137            // grid point definition extracted from
138            // testdata/20260219T00Z_MSC_HRDPS_CAPE_Sfc_RLatLon0.0225_PT000H.grib2
139            rotation_for_first_point,
140            Rotation {
141                south_pole_lat: -36088520,
142                south_pole_lon: 245305142,
143                rot_angle: 0.,
144            },
145            (-12.302501_f32, 345.178780_f32),
146            // taken from results from pygrib
147            (39.626032, -133.62952 + 720.)
148        ),
149    }
150}
grib/grid/rotated_ll.rs

grib/grid/
rotated_ll.rs
