grib/grid/
polar_stereographic.rs

1#[cfg(feature = "gridpoints-proj")]
2use crate::projection::OsgeoProj;
3use crate::{
4    GridPointIndex, LatLons,
5    def::grib2::{table, template::Template3_20},
6    error::GribError,
7    grid::AngleUnit,
8    projection,
9};
10
11impl crate::GridShortName for Template3_20 {
12    fn short_name(&self) -> &'static str {
13        "polar_stereographic"
14    }
15}
16
17impl GridPointIndex for Template3_20 {
18    fn grid_shape(&self) -> (usize, usize) {
19        (self.ni as usize, self.nj as usize)
20    }
21
22    fn scanning_mode(&self) -> &table::ScanningMode {
23        &self.scanning_mode
24    }
25}
26
27impl LatLons for Template3_20 {
28    #[cfg(feature = "gridpoints-proj")]
29    type Iter<'a> = super::helpers::ProjectionLatLonIterator;
30    #[cfg(not(feature = "gridpoints-proj"))]
31    type Iter<'a> = std::vec::IntoIter<(f64, f64)>;
32
33    fn latlons_unchecked<'a>(&'a self) -> Result<Self::Iter<'a>, GribError> {
34        let angle_units = self.angle_unit();
35        let lad = self.lad as f64 * angle_units;
36        let lov = self.lov as f64 * angle_units;
37        let (a, b) = self.earth_shape.radii().ok_or_else(|| {
38            GribError::NotSupported(format!(
39                "unknown value of Code Table 3.2 (shape of the Earth): {}",
40                self.earth_shape.shape.0
41            ))
42        })?;
43
44        if self.projection_centre.has_unsupported_flags() {
45            let table::ProjectionCentreFlag(flag) = self.projection_centre;
46            return Err(GribError::NotSupported(format!("projection centre {flag}")));
47        }
48        let lat_origin = if self
49            .projection_centre
50            .contains_north_pole_on_projection_plane()
51        {
52            90.
53        } else {
54            -90.
55        };
56
57        let params = projection::StereParams {
58            ellipsoid: projection::Ellipsoid::from_a_and_b(a, b),
59            lat_ts: lad,
60            lat_0: lat_origin,
61            lon_0: lov,
62            k_0: 1.0,
63        };
64
65        let dx = self.dx as f64 * 1e-3;
66        let dy = self.dy as f64 * 1e-3;
67        let dx = if !self.scanning_mode.scans_positively_for_i() && dx > 0. {
68            -dx
69        } else {
70            dx
71        };
72        let dy = if !self.scanning_mode.scans_positively_for_j() && dy > 0. {
73            -dy
74        } else {
75            dy
76        };
77
78        let first_point = (
79            self.first_point_lat as f64 * angle_units,
80            self.first_point_lon as f64 * angle_units,
81        );
82
83        #[cfg(feature = "gridpoints-proj")]
84        {
85            super::helpers::latlons_from_projection_with_first_point_and_delta(
86                &params.proj_args(),
87                first_point,
88                (dx, dy),
89                self.ij()?,
90            )
91        }
92
93        #[cfg(not(feature = "gridpoints-proj"))]
94        {
95            let projection = projection::Stere::new(&params)?;
96            super::helpers::latlons_from_projection_with_first_point_and_delta(
97                &projection,
98                first_point,
99                (dx, dy),
100                self.ij()?,
101            )
102        }
103    }
104}
105
106impl AngleUnit for Template3_20 {
107    fn angle_unit(&self) -> f64 {
108        1e-6
109    }
110}
111
112#[cfg(test)]
113mod tests {
114    use super::*;
115    use crate::def::grib2::template::param_set;
116
117    #[test]
118    fn polar_stereographic_grid_latlon_computation() -> Result<(), Box<dyn std::error::Error>> {
119        use crate::grid::helpers::test_helpers::assert_coord_almost_eq;
120        // grid point definition extracted from
121        // testdata/CMC_RDPA_APCP-024-0100cutoff_SFC_0_ps10km_2023121806_000.
122        // grib2.xz
123        let grid_def = Template3_20 {
124            earth_shape: param_set::EarthShape {
125                shape: table::EarthShapeType(6),
126                spherical_earth_radius: param_set::ScaledValue {
127                    scale_factor: 0xff,
128                    scaled_value: 0xffffffff,
129                },
130                major_axis: param_set::ScaledValue {
131                    scale_factor: 0xff,
132                    scaled_value: 0xffffffff,
133                },
134                minor_axis: param_set::ScaledValue {
135                    scale_factor: 0xff,
136                    scaled_value: 0xffffffff,
137                },
138            },
139            ni: 935,
140            nj: 824,
141            first_point_lat: 18145030,
142            first_point_lon: 217107456,
143            resolution_and_component_flags: table::ResolutionAndComponentFlags(0b00001000),
144            lad: 60000000,
145            lov: 249000000,
146            dx: 10000000,
147            dy: 10000000,
148            projection_centre: table::ProjectionCentreFlag(0b00000000),
149            scanning_mode: table::ScanningMode(0b01000000),
150        };
151        let latlons = grid_def.latlons()?.collect::<Vec<_>>();
152
153        // Following lat/lon values are taken from the calculation results using
154        // pygrib.
155        let delta = 1e-4;
156        assert_coord_almost_eq(latlons[0], (18.14503, -142.892544), delta);
157        assert_coord_almost_eq(latlons[1], (18.17840149, -142.83604096), delta);
158        assert_coord_almost_eq(
159            latlons[latlons.len() - 2],
160            (45.4865147, -10.15230394),
161            delta,
162        );
163        assert_coord_almost_eq(
164            latlons[latlons.len() - 1],
165            (45.40545211, -10.17442147),
166            delta,
167        );
168
169        Ok(())
170    }
171}
grib/grid/polar_stereographic.rs

grib/grid/
polar_stereographic.rs
