grib/grid/
polar_stereographic.rs

1#[cfg(feature = "gridpoints-proj")]
2use crate::error::GribError;
3use crate::{
4    GridPointIndex, LatLons,
5    def::grib2::template::{Template3_20, param_set},
6    grid::AngleUnit,
7};
8
9impl crate::GridShortName for Template3_20 {
10    fn short_name(&self) -> &'static str {
11        "polar_stereographic"
12    }
13}
14
15impl GridPointIndex for Template3_20 {
16    fn grid_shape(&self) -> (usize, usize) {
17        (self.ni as usize, self.nj as usize)
18    }
19
20    fn scanning_mode(&self) -> &param_set::ScanningMode {
21        &self.scanning_mode
22    }
23}
24
25#[cfg(feature = "gridpoints-proj")]
26#[cfg_attr(docsrs, doc(cfg(feature = "gridpoints-proj")))]
27impl LatLons for Template3_20 {
28    type Iter<'a>
29        = std::vec::IntoIter<(f32, f32)>
30    where
31        Self: 'a;
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
41            ))
42        })?;
43
44        if self.projection_centre.has_unsupported_flags() {
45            let param_set::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 proj_def =
58            format!("+a={a} +b={b} +proj=stere +lat_ts={lad} +lat_0={lat_origin} +lon_0={lov}");
59
60        let dx = self.dx as f64 * 1e-3;
61        let dy = self.dy as f64 * 1e-3;
62        let dx = if !self.scanning_mode.scans_positively_for_i() && dx > 0. {
63            -dx
64        } else {
65            dx
66        };
67        let dy = if !self.scanning_mode.scans_positively_for_j() && dy > 0. {
68            -dy
69        } else {
70            dy
71        };
72
73        super::helpers::latlons_from_projection_definition_and_first_point(
74            &proj_def,
75            (
76                self.first_point_lat as f64 * angle_units,
77                self.first_point_lon as f64 * angle_units,
78            ),
79            (dx, dy),
80            self.ij()?,
81        )
82    }
83}
84
85impl AngleUnit for Template3_20 {
86    fn angle_unit(&self) -> f64 {
87        1e-6
88    }
89}
90
91#[cfg(test)]
92mod tests {
93    use super::*;
94
95    #[cfg(feature = "gridpoints-proj")]
96    #[test]
97    fn polar_stereographic_grid_latlon_computation() -> Result<(), Box<dyn std::error::Error>> {
98        use crate::grid::helpers::test_helpers::assert_coord_almost_eq;
99        // grid point definition extracted from
100        // testdata/CMC_RDPA_APCP-024-0100cutoff_SFC_0_ps10km_2023121806_000.grib2.xz
101        let grid_def = Template3_20 {
102            earth_shape: param_set::EarthShape {
103                shape: 6,
104                spherical_earth_radius_scale_factor: 0xff,
105                spherical_earth_radius_scaled_value: 0xffffffff,
106                major_axis_scale_factor: 0xff,
107                major_axis_scaled_value: 0xffffffff,
108                minor_axis_scale_factor: 0xff,
109                minor_axis_scaled_value: 0xffffffff,
110            },
111            ni: 935,
112            nj: 824,
113            first_point_lat: 18145030,
114            first_point_lon: 217107456,
115            resolution_and_component_flags: param_set::ResolutionAndComponentFlags(0b00001000),
116            lad: 60000000,
117            lov: 249000000,
118            dx: 10000000,
119            dy: 10000000,
120            projection_centre: param_set::ProjectionCentreFlag(0b00000000),
121            scanning_mode: param_set::ScanningMode(0b01000000),
122        };
123        let latlons = grid_def.latlons()?.collect::<Vec<_>>();
124
125        // Following lat/lon values are taken from the calculation results using pygrib.
126        let delta = 1e-4;
127        assert_coord_almost_eq(latlons[0], (18.14503, -142.892544), delta);
128        assert_coord_almost_eq(latlons[1], (18.17840149, -142.83604096), delta);
129        assert_coord_almost_eq(
130            latlons[latlons.len() - 2],
131            (45.4865147, -10.15230394),
132            delta,
133        );
134        assert_coord_almost_eq(
135            latlons[latlons.len() - 1],
136            (45.40545211, -10.17442147),
137            delta,
138        );
139
140        Ok(())
141    }
142}
grib/grid/polar_stereographic.rs

grib/grid/
polar_stereographic.rs
