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    #[cfg_attr(docsrs, doc(cfg(feature = "gridpoints-proj")))]
120    pub fn latlons(&self) -> Result<std::vec::IntoIter<(f32, f32)>, GribError> {
121        let lad = self.lad as f64 * 1e-6;
122        let lov = self.lov as f64 * 1e-6;
123        let (a, b) = self.earth_shape.radii().ok_or_else(|| {
124            GribError::NotSupported(format!(
125                "unknown value of Code Table 3.2 (shape of the Earth): {}",
126                self.earth_shape.shape_of_the_earth
127            ))
128        })?;
129
130        if self.projection_centre.has_unsupported_flags() {
131            let ProjectionCentreFlag(flag) = self.projection_centre;
132            return Err(GribError::NotSupported(format!("projection centre {flag}")));
133        }
134        let lat_origin = if self
135            .projection_centre
136            .contains_north_pole_on_projection_plane()
137        {
138            90.
139        } else {
140            -90.
141        };
142
143        let proj_def =
144            format!("+a={a} +b={b} +proj=stere +lat_ts={lad} +lat_0={lat_origin} +lon_0={lov}");
145
146        let dx = self.dx as f64 * 1e-3;
147        let dy = self.dy as f64 * 1e-3;
148        let dx = if !self.scanning_mode.scans_positively_for_i() && dx > 0. {
149            -dx
150        } else {
151            dx
152        };
153        let dy = if !self.scanning_mode.scans_positively_for_j() && dy > 0. {
154            -dy
155        } else {
156            dy
157        };
158
159        super::helpers::latlons_from_projection_definition_and_first_point(
160            &proj_def,
161            (
162                self.first_point_lat as f64 * 1e-6,
163                self.first_point_lon as f64 * 1e-6,
164            ),
165            (dx, dy),
166            self.ij()?,
167        )
168    }
169
170    pub(crate) fn from_buf(buf: &[u8]) -> Self {
171        let earth_shape = EarthShapeDefinition::from_buf(buf);
172        let ni = read_as!(u32, buf, 16);
173        let nj = read_as!(u32, buf, 20);
174        let first_point_lat = read_as!(u32, buf, 24).as_grib_int();
175        let first_point_lon = read_as!(u32, buf, 28).as_grib_int();
176        let lad = read_as!(u32, buf, 33).as_grib_int();
177        let lov = read_as!(u32, buf, 37).as_grib_int();
178        let dx = read_as!(u32, buf, 41);
179        let dy = read_as!(u32, buf, 45);
180        let projection_centre = read_as!(u8, buf, 49);
181        let scanning_mode = read_as!(u8, buf, 50);
182        Self {
183            earth_shape,
184            ni,
185            nj,
186            first_point_lat,
187            first_point_lon,
188            lad,
189            lov,
190            dx,
191            dy,
192            projection_centre: ProjectionCentreFlag(projection_centre),
193            scanning_mode: ScanningMode(scanning_mode),
194        }
195    }
196}
197
198#[cfg(test)]
199mod tests {
200    use std::io::{BufReader, Read};
201
202    use super::*;
203
204    #[test]
205    fn polar_stereographic_grid_definition_from_buf() -> Result<(), Box<dyn std::error::Error>> {
206        let mut buf = Vec::new();
207
208        let f = std::fs::File::open(
209            "testdata/CMC_RDPA_APCP-024-0100cutoff_SFC_0_ps10km_2023121806_000.grib2.xz",
210        )?;
211        let f = BufReader::new(f);
212        let mut f = xz2::bufread::XzDecoder::new(f);
213        f.read_to_end(&mut buf)?;
214
215        let actual = PolarStereographicGridDefinition::from_buf(&buf[0x33..]);
216        let expected = PolarStereographicGridDefinition {
217            earth_shape: EarthShapeDefinition {
218                shape_of_the_earth: 6,
219                scale_factor_of_radius_of_spherical_earth: 0xff,
220                scaled_value_of_radius_of_spherical_earth: 0xffffffff,
221                scale_factor_of_earth_major_axis: 0xff,
222                scaled_value_of_earth_major_axis: 0xffffffff,
223                scale_factor_of_earth_minor_axis: 0xff,
224                scaled_value_of_earth_minor_axis: 0xffffffff,
225            },
226            ni: 935,
227            nj: 824,
228            first_point_lat: 18145030,
229            first_point_lon: 217107456,
230            lad: 60000000,
231            lov: 249000000,
232            dx: 10000000,
233            dy: 10000000,
234            projection_centre: ProjectionCentreFlag(0b00000000),
235            scanning_mode: ScanningMode(0b01000000),
236        };
237        assert_eq!(actual, expected);
238
239        Ok(())
240    }
241
242    #[cfg(feature = "gridpoints-proj")]
243    #[test]
244    fn polar_stereographic_grid_latlon_computation() -> Result<(), Box<dyn std::error::Error>> {
245        use crate::grid::helpers::test_helpers::assert_coord_almost_eq;
246        let grid_def = PolarStereographicGridDefinition {
247            earth_shape: EarthShapeDefinition {
248                shape_of_the_earth: 6,
249                scale_factor_of_radius_of_spherical_earth: 0xff,
250                scaled_value_of_radius_of_spherical_earth: 0xffffffff,
251                scale_factor_of_earth_major_axis: 0xff,
252                scaled_value_of_earth_major_axis: 0xffffffff,
253                scale_factor_of_earth_minor_axis: 0xff,
254                scaled_value_of_earth_minor_axis: 0xffffffff,
255            },
256            ni: 935,
257            nj: 824,
258            first_point_lat: 18145030,
259            first_point_lon: 217107456,
260            lad: 60000000,
261            lov: 249000000,
262            dx: 10000000,
263            dy: 10000000,
264            projection_centre: ProjectionCentreFlag(0b00000000),
265            scanning_mode: ScanningMode(0b01000000),
266        };
267        let latlons = grid_def.latlons()?.collect::<Vec<_>>();
268
269        // Following lat/lon values are taken from the calculation results using pygrib.
270        let delta = 1e-10;
271        assert_coord_almost_eq(latlons[0], (18.14503, -142.892544), delta);
272        assert_coord_almost_eq(latlons[1], (18.17840149, -142.83604096), delta);
273        assert_coord_almost_eq(
274            latlons[latlons.len() - 2],
275            (45.4865147, -10.15230394),
276            delta,
277        );
278        assert_coord_almost_eq(
279            latlons[latlons.len() - 1],
280            (45.40545211, -10.17442147),
281            delta,
282        );
283
284        Ok(())
285    }
286}
grib/grid/polar_stereographic.rs

grib/grid/
polar_stereographic.rs
