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::{AngleUnit, 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_unchecked<'a>(&'a self) -> Result<Self::Iter<'a>, GribError> {
36        let iter = Unrotate::new(
37            self.rotated.latlons_unchecked()?,
38            &self.rotation,
39            self.angle_unit() as f32,
40        );
41        Ok(iter)
42    }
43}
44
45impl AngleUnit for Template3_1 {
46    fn angle_unit(&self) -> f64 {
47        self.rotated.grid.angle_unit()
48    }
49}
50
51#[derive(Clone)]
52pub struct Unrotate<I> {
53    latlons: I,
54    sinφp: f32,
55    cosφp: f32,
56    λp: f32,
57    gamma: f32,
58}
59
60impl<I> Unrotate<I> {
61    fn new(latlons: I, rot: &Rotation, angle_units: f32) -> Self {
62        let φp = (rot.south_pole_lat as f32 * angle_units).to_radians();
63        let λp = (rot.south_pole_lon as f32 * angle_units).to_radians();
64        let gamma = (rot.rot_angle * angle_units).to_radians();
65
66        // south pole to north pole
67        let φp = -φp;
68        let λp = λp + std::f32::consts::PI;
69
70        let (sinφp, cosφp) = φp.sin_cos();
71        Self {
72            latlons,
73            sinφp,
74            cosφp,
75            λp,
76            gamma,
77        }
78    }
79}
80
81impl<I> Iterator for Unrotate<I>
82where
83    I: Iterator<Item = (f32, f32)>,
84{
85    type Item = (f32, f32);
86
87    fn next(&mut self) -> Option<Self::Item> {
88        let (lat, lon) = self.latlons.next()?;
89        let λr = lon.to_radians();
90        let φr = lat.to_radians();
91
92        let λr = λr - self.gamma;
93
94        let (sinφr, cosφr) = φr.sin_cos();
95        let (sinλr, cosλr) = λr.sin_cos();
96
97        let sinφ = self.sinφp * sinφr + self.cosφp * cosφr * cosλr;
98        let φ = sinφ.asin();
99
100        let y = cosφr * sinλr;
101        let x = self.cosφp * sinφr - self.sinφp * cosφr * cosλr;
102        let λ = self.λp - y.atan2(x);
103
104        let latlon = (φ.to_degrees(), λ.to_degrees());
105        Some(latlon)
106    }
107
108    fn size_hint(&self) -> (usize, Option<usize>) {
109        self.latlons.size_hint()
110    }
111}
112
113#[cfg(test)]
114mod tests {
115
116    use super::*;
117
118    macro_rules! test_rotation{
119        ($(($name:ident, $rot:expr, $input:expr, $expected:expr),)*) => ($(
120            #[test]
121            fn $name() {
122                let rot = $rot;
123                let latlons = vec![$input];
124                let mut iter = Unrotate::new(latlons.into_iter(), &rot, 1e-6);
125                let actual = iter.next().unwrap();
126                let expected = $expected;
127                dbg!(actual);
128                dbg!(expected);
129                assert!((actual.0 - expected.0).abs() < 0.00001);
130                assert!((actual.1 - expected.1).abs() < 0.001);
131            }
132        )*);
133    }
134
135    test_rotation! {
136        (
137            no_rotation,
138            Rotation {
139                south_pole_lat: -90000000,
140                south_pole_lon: 0,
141                rot_angle: 0.,
142            },
143            (-12.302501_f32, 345.178780_f32),
144            (-12.302501_f32, 345.178780_f32)
145        ),
146        (
147            // grid point definition extracted from
148            // testdata/20260219T00Z_MSC_HRDPS_CAPE_Sfc_RLatLon0.0225_PT000H.grib2
149            rotation_for_first_point,
150            Rotation {
151                south_pole_lat: -36088520,
152                south_pole_lon: 245305142,
153                rot_angle: 0.,
154            },
155            (-12.302501_f32, 345.178780_f32),
156            // taken from results from pygrib
157            (39.626032, -133.62952 + 720.)
158        ),
159    }
160}
grib/grid/rotated_ll.rs

grib/grid/
rotated_ll.rs
