grib/grid/
polar_stereographic.rs

1use super::{earth::EarthShapeDefinition, GridPointIndexIterator, ScanningMode};
2use crate::{
3    error::GribError,
4    helpers::{read_as, GribInt},
5    ProjectionCentreFlag,
6};
7
8#[derive(Debug, PartialEq, Eq)]
9pub struct PolarStereographicGridDefinition {
10    pub earth_shape: EarthShapeDefinition,
11    pub ni: u32,
12    pub nj: u32,
13    pub first_point_lat: i32,
14    pub first_point_lon: i32,
15    pub lad: i32,
16    pub lov: i32,
17    pub dx: u32,
18    pub dy: u32,
19    pub projection_centre: ProjectionCentreFlag,
20    pub scanning_mode: ScanningMode,
21}
22
23impl PolarStereographicGridDefinition {
24    /// Returns the shape of the grid, i.e. a tuple of the number of grids in
25    /// the i and j directions.
26    ///
27    /// Examples
28    ///
29    /// ```
30    /// let def = grib::PolarStereographicGridDefinition {
31    ///     earth_shape: grib::EarthShapeDefinition {
32    ///         shape_of_the_earth: 6,
33    ///         scale_factor_of_radius_of_spherical_earth: 0xff,
34    ///         scaled_value_of_radius_of_spherical_earth: 0xffffffff,
35    ///         scale_factor_of_earth_major_axis: 0xff,
36    ///         scaled_value_of_earth_major_axis: 0xffffffff,
37    ///         scale_factor_of_earth_minor_axis: 0xff,
38    ///         scaled_value_of_earth_minor_axis: 0xffffffff,
39    ///     },
40    ///     ni: 2,
41    ///     nj: 3,
42    ///     first_point_lat: 0,
43    ///     first_point_lon: 0,
44    ///     lad: 0,
45    ///     lov: 0,
46    ///     dx: 1000,
47    ///     dy: 1000,
48    ///     projection_centre: grib::ProjectionCentreFlag(0b00000000),
49    ///     scanning_mode: grib::ScanningMode(0b01000000),
50    /// };
51    /// let shape = def.grid_shape();
52    /// assert_eq!(shape, (2, 3));
53    /// ```
54    pub fn grid_shape(&self) -> (usize, usize) {
55        (self.ni as usize, self.nj as usize)
56    }
57
58    /// Returns the grid type.
59    pub fn short_name(&self) -> &'static str {
60        "polar_stereographic"
61    }
62
63    /// Returns an iterator over `(i, j)` of grid points.
64    ///
65    /// Note that this is a low-level API and it is not checked that the number
66    /// of iterator iterations is consistent with the number of grid points
67    /// defined in the data.
68    ///
69    /// Examples
70    ///
71    /// ```
72    /// let def = grib::PolarStereographicGridDefinition {
73    ///     earth_shape: grib::EarthShapeDefinition {
74    ///         shape_of_the_earth: 6,
75    ///         scale_factor_of_radius_of_spherical_earth: 0xff,
76    ///         scaled_value_of_radius_of_spherical_earth: 0xffffffff,
77    ///         scale_factor_of_earth_major_axis: 0xff,
78    ///         scaled_value_of_earth_major_axis: 0xffffffff,
79    ///         scale_factor_of_earth_minor_axis: 0xff,
80    ///         scaled_value_of_earth_minor_axis: 0xffffffff,
81    ///     },
82    ///     ni: 2,
83    ///     nj: 3,
84    ///     first_point_lat: 0,
85    ///     first_point_lon: 0,
86    ///     lad: 0,
87    ///     lov: 0,
88    ///     dx: 1000,
89    ///     dy: 1000,
90    ///     projection_centre: grib::ProjectionCentreFlag(0b00000000),
91    ///     scanning_mode: grib::ScanningMode(0b01000000),
92    /// };
93    /// let ij = def.ij();
94    /// assert!(ij.is_ok());
95    ///
96    /// let mut ij = ij.unwrap();
97    /// assert_eq!(ij.next(), Some((0, 0)));
98    /// assert_eq!(ij.next(), Some((1, 0)));
99    /// assert_eq!(ij.next(), Some((0, 1)));
100    /// ```
101    pub fn ij(&self) -> Result<GridPointIndexIterator, GribError> {
102        if self.scanning_mode.has_unsupported_flags() {
103            let ScanningMode(mode) = self.scanning_mode;
104            return Err(GribError::NotSupported(format!("scanning mode {mode}")));
105        }
106
107        let iter =
108            GridPointIndexIterator::new(self.ni as usize, self.nj as usize, self.scanning_mode);
109        Ok(iter)
110    }
111
112    /// Returns an iterator over latitudes and longitudes of grid points in
113    /// degrees.
114    ///
115    /// Note that this is a low-level API and it is not checked that the number
116    /// of iterator iterations is consistent with the number of grid points
117    /// defined in the data.
118    #[cfg(feature = "gridpoints-proj")]
119    pub fn latlons(&self) -> Result<std::vec::IntoIter<(f32, f32)>, GribError> {
120        let lad = self.lad as f64 * 1e-6;
121        let lov = self.lov as f64 * 1e-6;
122        let (a, b) = self.earth_shape.radii().ok_or_else(|| {
123            GribError::NotSupported(format!(
124                "unknown value of Code Table 3.2 (shape of the Earth): {}",
125                self.earth_shape.shape_of_the_earth
126            ))
127        })?;
128
129        if self.projection_centre.has_unsupported_flags() {
130            let ProjectionCentreFlag(flag) = self.projection_centre;
131            return Err(GribError::NotSupported(format!("projection centre {flag}")));
132        }
133        let lat_origin = if self
134            .projection_centre
135            .contains_north_pole_on_projection_plane()
136        {
137            90.
138        } else {
139            -90.
140        };
141
142        let proj_def =
143            format!("+a={a} +b={b} +proj=stere +lat_ts={lad} +lat_0={lat_origin} +lon_0={lov}");
144
145        let dx = self.dx as f64 * 1e-3;
146        let dy = self.dy as f64 * 1e-3;
147        let dx = if !self.scanning_mode.scans_positively_for_i() && dx > 0. {
148            -dx
149        } else {
150            dx
151        };
152        let dy = if !self.scanning_mode.scans_positively_for_j() && dy > 0. {
153            -dy
154        } else {
155            dy
156        };
157
158        super::helpers::latlons_from_projection_definition_and_first_point(
159            &proj_def,
160            (
161                self.first_point_lat as f64 * 1e-6,
162                self.first_point_lon as f64 * 1e-6,
163            ),
164            (dx, dy),
165            self.ij()?,
166        )
167    }
168
169    pub(crate) fn from_buf(buf: &[u8]) -> Self {
170        let earth_shape = EarthShapeDefinition::from_buf(buf);
171        let ni = read_as!(u32, buf, 16);
172        let nj = read_as!(u32, buf, 20);
173        let first_point_lat = read_as!(u32, buf, 24).as_grib_int();
174        let first_point_lon = read_as!(u32, buf, 28).as_grib_int();
175        let lad = read_as!(u32, buf, 33).as_grib_int();
176        let lov = read_as!(u32, buf, 37).as_grib_int();
177        let dx = read_as!(u32, buf, 41);
178        let dy = read_as!(u32, buf, 45);
179        let projection_centre = read_as!(u8, buf, 49);
180        let scanning_mode = read_as!(u8, buf, 50);
181        Self {
182            earth_shape,
183            ni,
184            nj,
185            first_point_lat,
186            first_point_lon,
187            lad,
188            lov,
189            dx,
190            dy,
191            projection_centre: ProjectionCentreFlag(projection_centre),
192            scanning_mode: ScanningMode(scanning_mode),
193        }
194    }
195}
196
197#[cfg(test)]
198mod tests {
199    use std::io::{BufReader, Read};
200
201    use super::*;
202
203    #[test]
204    fn polar_stereographic_grid_definition_from_buf() -> Result<(), Box<dyn std::error::Error>> {
205        let mut buf = Vec::new();
206
207        let f = std::fs::File::open(
208            "testdata/CMC_RDPA_APCP-024-0100cutoff_SFC_0_ps10km_2023121806_000.grib2.xz",
209        )?;
210        let f = BufReader::new(f);
211        let mut f = xz2::bufread::XzDecoder::new(f);
212        f.read_to_end(&mut buf)?;
213
214        let actual = PolarStereographicGridDefinition::from_buf(&buf[0x33..]);
215        let expected = PolarStereographicGridDefinition {
216            earth_shape: EarthShapeDefinition {
217                shape_of_the_earth: 6,
218                scale_factor_of_radius_of_spherical_earth: 0xff,
219                scaled_value_of_radius_of_spherical_earth: 0xffffffff,
220                scale_factor_of_earth_major_axis: 0xff,
221                scaled_value_of_earth_major_axis: 0xffffffff,
222                scale_factor_of_earth_minor_axis: 0xff,
223                scaled_value_of_earth_minor_axis: 0xffffffff,
224            },
225            ni: 935,
226            nj: 824,
227            first_point_lat: 18145030,
228            first_point_lon: 217107456,
229            lad: 60000000,
230            lov: 249000000,
231            dx: 10000000,
232            dy: 10000000,
233            projection_centre: ProjectionCentreFlag(0b00000000),
234            scanning_mode: ScanningMode(0b01000000),
235        };
236        assert_eq!(actual, expected);
237
238        Ok(())
239    }
240
241    #[cfg(feature = "gridpoints-proj")]
242    #[test]
243    fn polar_stereographic_grid_latlon_computation() -> Result<(), Box<dyn std::error::Error>> {
244        use crate::grid::helpers::test_helpers::assert_coord_almost_eq;
245        let grid_def = PolarStereographicGridDefinition {
246            earth_shape: EarthShapeDefinition {
247                shape_of_the_earth: 6,
248                scale_factor_of_radius_of_spherical_earth: 0xff,
249                scaled_value_of_radius_of_spherical_earth: 0xffffffff,
250                scale_factor_of_earth_major_axis: 0xff,
251                scaled_value_of_earth_major_axis: 0xffffffff,
252                scale_factor_of_earth_minor_axis: 0xff,
253                scaled_value_of_earth_minor_axis: 0xffffffff,
254            },
255            ni: 935,
256            nj: 824,
257            first_point_lat: 18145030,
258            first_point_lon: 217107456,
259            lad: 60000000,
260            lov: 249000000,
261            dx: 10000000,
262            dy: 10000000,
263            projection_centre: ProjectionCentreFlag(0b00000000),
264            scanning_mode: ScanningMode(0b01000000),
265        };
266        let latlons = grid_def.latlons()?.collect::<Vec<_>>();
267
268        // Following lat/lon values are taken from the calculation results using pygrib.
269        let delta = 1e-10;
270        assert_coord_almost_eq(latlons[0], (18.14503, -142.892544), delta);
271        assert_coord_almost_eq(latlons[1], (18.17840149, -142.83604096), delta);
272        assert_coord_almost_eq(
273            latlons[latlons.len() - 2],
274            (45.4865147, -10.15230394),
275            delta,
276        );
277        assert_coord_almost_eq(
278            latlons[latlons.len() - 1],
279            (45.40545211, -10.17442147),
280            delta,
281        );
282
283        Ok(())
284    }
285}
grib/grid/polar_stereographic.rs

grib/grid/
polar_stereographic.rs
