grib/grid/
polar_stereographic.rs

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

grib/grid/
polar_stereographic.rs
