grib/grid/
lambert.rs

1#[cfg(feature = "gridpoints-proj")]
2use crate::LatLons;
3use crate::{
4    GridPointIndex,
5    def::grib2::template::{Template3_30, param_set},
6    error::GribError,
7    grid::AngleUnit,
8};
9
10impl crate::GridShortName for Template3_30 {
11    fn short_name(&self) -> &'static str {
12        "lambert"
13    }
14}
15
16impl GridPointIndex for Template3_30 {
17    fn grid_shape(&self) -> (usize, usize) {
18        (self.ni as usize, self.nj as usize)
19    }
20
21    fn scanning_mode(&self) -> &param_set::ScanningMode {
22        &self.scanning_mode
23    }
24}
25
26#[cfg(feature = "gridpoints-proj")]
27#[cfg_attr(docsrs, doc(cfg(feature = "gridpoints-proj")))]
28impl LatLons for Template3_30 {
29    type Iter<'a> = std::vec::IntoIter<(f32, f32)>;
30
31    fn latlons_unchecked<'a>(&'a self) -> Result<Self::Iter<'a>, GribError> {
32        let angle_units = self.angle_unit();
33        let lad = self.lad as f64 * angle_units;
34        let lov = self.lov as f64 * angle_units;
35        let latin1 = self.latin1 as f64 * angle_units;
36        let latin2 = self.latin2 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        let proj_def = format!(
44            "+a={a} +b={b} +proj=lcc +lat_0={lad} +lon_0={lov} +lat_1={latin1} +lat_2={latin2}"
45        );
46
47        let dx = self.dx as f64 * 1e-3;
48        let dy = self.dy as f64 * 1e-3;
49        let dx = if !self.scanning_mode.scans_positively_for_i() && dx > 0. {
50            -dx
51        } else {
52            dx
53        };
54        let dy = if !self.scanning_mode.scans_positively_for_j() && dy > 0. {
55            -dy
56        } else {
57            dy
58        };
59
60        super::helpers::latlons_from_projection_definition_and_first_point(
61            &proj_def,
62            (
63                self.first_point_lat as f64 * angle_units,
64                self.first_point_lon as f64 * angle_units,
65            ),
66            (dx, dy),
67            self.ij()?,
68        )
69    }
70}
71
72impl AngleUnit for Template3_30 {
73    fn angle_unit(&self) -> f64 {
74        1e-6
75    }
76}
77
78#[cfg(test)]
79mod tests {
80    use super::*;
81
82    #[cfg(feature = "gridpoints-proj")]
83    #[test]
84    fn lambert_grid_latlon_computation() -> Result<(), Box<dyn std::error::Error>> {
85        use crate::grid::helpers::test_helpers::assert_coord_almost_eq;
86        // grid point definition extracted from testdata/ds.critfireo.bin.xz
87        let grid_def = Template3_30 {
88            earth_shape: param_set::EarthShape {
89                shape: 1,
90                spherical_earth_radius_scale_factor: 0,
91                spherical_earth_radius_scaled_value: 6371200,
92                major_axis_scale_factor: 0,
93                major_axis_scaled_value: 0,
94                minor_axis_scale_factor: 0,
95                minor_axis_scaled_value: 0,
96            },
97            ni: 2145,
98            nj: 1377,
99            first_point_lat: 20190000,
100            first_point_lon: 238449996,
101            resolution_and_component_flags: param_set::ResolutionAndComponentFlags(0b00000000),
102            lad: 25000000,
103            lov: 265000000,
104            dx: 2539703,
105            dy: 2539703,
106            projection_centre: param_set::ProjectionCentreFlag(0b00000000),
107            scanning_mode: param_set::ScanningMode(0b01010000),
108            latin1: 25000000,
109            latin2: 25000000,
110            south_pole_lat: -90000000,
111            south_pole_lon: 0,
112        };
113        let latlons = grid_def.latlons()?.collect::<Vec<_>>();
114
115        // Following lat/lon values are taken from the calculation results using pygrib.
116        let delta = 1e-4;
117        assert_coord_almost_eq(latlons[0], (20.19, -121.550004), delta);
118        assert_coord_almost_eq(latlons[1], (20.19442682, -121.52621665), delta);
119        assert_coord_almost_eq(
120            latlons[latlons.len() - 2],
121            (50.10756403, -60.91298217),
122            delta,
123        );
124        assert_coord_almost_eq(
125            latlons[latlons.len() - 1],
126            (50.1024611, -60.88202274),
127            delta,
128        );
129
130        Ok(())
131    }
132}
grib/grid/lambert.rs

grib/grid/
lambert.rs
