grib/grid/
lambert.rs

1#[cfg(feature = "gridpoints-proj")]
2use crate::projection::OsgeoProj;
3#[cfg(not(feature = "gridpoints-proj"))]
4use crate::projection::Project;
5use crate::{
6    GridPointIndex, LatLons,
7    def::grib2::template::{Template3_30, param_set},
8    error::GribError,
9    grid::AngleUnit,
10    projection,
11};
12
13impl crate::GridShortName for Template3_30 {
14    fn short_name(&self) -> &'static str {
15        "lambert"
16    }
17}
18
19impl GridPointIndex for Template3_30 {
20    fn grid_shape(&self) -> (usize, usize) {
21        (self.ni as usize, self.nj as usize)
22    }
23
24    fn scanning_mode(&self) -> &param_set::ScanningMode {
25        &self.scanning_mode
26    }
27}
28
29impl LatLons for Template3_30 {
30    #[cfg(feature = "gridpoints-proj")]
31    type Iter<'a> = super::helpers::ProjectionLatLonIterator;
32    #[cfg(not(feature = "gridpoints-proj"))]
33    type Iter<'a> = std::vec::IntoIter<(f32, f32)>;
34
35    fn latlons_unchecked<'a>(&'a self) -> Result<Self::Iter<'a>, GribError> {
36        let angle_units = self.angle_unit();
37        let lad = self.lad as f64 * angle_units;
38        let lov = self.lov as f64 * angle_units;
39        let latin1 = self.latin1 as f64 * angle_units;
40        let latin2 = self.latin2 as f64 * angle_units;
41        let (a, b) = self.earth_shape.radii().ok_or_else(|| {
42            GribError::NotSupported(format!(
43                "unknown value of Code Table 3.2 (shape of the Earth): {}",
44                self.earth_shape.shape
45            ))
46        })?;
47        let params = projection::LccParams {
48            ellipsoid: projection::Ellipsoid::from_a_and_b(a, b),
49            lat_0: lad,
50            lon_0: lov,
51            lat_1: latin1,
52            lat_2: latin2,
53        };
54
55        let dx = self.dx as f64 * 1e-3;
56        let dy = self.dy as f64 * 1e-3;
57        let dx = if !self.scanning_mode.scans_positively_for_i() && dx > 0. {
58            -dx
59        } else {
60            dx
61        };
62        let dy = if !self.scanning_mode.scans_positively_for_j() && dy > 0. {
63            -dy
64        } else {
65            dy
66        };
67
68        let first_point = (
69            self.first_point_lat as f64 * angle_units,
70            self.first_point_lon as f64 * angle_units,
71        );
72
73        #[cfg(feature = "gridpoints-proj")]
74        {
75            super::helpers::latlons_from_projection_with_first_point_and_delta(
76                &params.proj_args(),
77                first_point,
78                (dx, dy),
79                self.ij()?,
80            )
81        }
82
83        #[cfg(not(feature = "gridpoints-proj"))]
84        {
85            let projection = projection::Lcc::new(&params)?;
86            let (first_point_lat, first_point_lon) = first_point;
87            let (first_corner_x, first_corner_y) = projection.project(
88                &(first_point_lon.to_radians(), first_point_lat.to_radians()),
89                false,
90            )?;
91
92            let latlons = self
93                .ij()?
94                .map(|(i, j)| {
95                    projection
96                        .project(
97                            &(
98                                first_corner_x + dx * i as f64,
99                                first_corner_y + dy * j as f64,
100                            ),
101                            true,
102                        )
103                        .map(|(lon, lat)| (lat.to_degrees() as f32, lon.to_degrees() as f32))
104                })
105                .collect::<Result<Vec<_>, _>>()?;
106            Ok(latlons.into_iter())
107        }
108    }
109}
110
111impl AngleUnit for Template3_30 {
112    fn angle_unit(&self) -> f64 {
113        1e-6
114    }
115}
116
117#[cfg(test)]
118mod tests {
119    use super::*;
120
121    #[test]
122    fn lambert_grid_latlon_computation() -> Result<(), Box<dyn std::error::Error>> {
123        use crate::grid::helpers::test_helpers::assert_coord_almost_eq;
124        // grid point definition extracted from testdata/ds.critfireo.bin.xz
125        let grid_def = Template3_30 {
126            earth_shape: param_set::EarthShape {
127                shape: 1,
128                spherical_earth_radius: param_set::ScaledValue {
129                    scale_factor: 0,
130                    scaled_value: 6371200,
131                },
132                major_axis: param_set::ScaledValue {
133                    scale_factor: 0,
134                    scaled_value: 0,
135                },
136                minor_axis: param_set::ScaledValue {
137                    scale_factor: 0,
138                    scaled_value: 0,
139                },
140            },
141            ni: 2145,
142            nj: 1377,
143            first_point_lat: 20190000,
144            first_point_lon: 238449996,
145            resolution_and_component_flags: param_set::ResolutionAndComponentFlags(0b00000000),
146            lad: 25000000,
147            lov: 265000000,
148            dx: 2539703,
149            dy: 2539703,
150            projection_centre: param_set::ProjectionCentreFlag(0b00000000),
151            scanning_mode: param_set::ScanningMode(0b01010000),
152            latin1: 25000000,
153            latin2: 25000000,
154            south_pole_lat: -90000000,
155            south_pole_lon: 0,
156        };
157        let latlons = grid_def.latlons()?.collect::<Vec<_>>();
158
159        // Following lat/lon values are taken from the calculation results using
160        // pygrib.
161        let delta = 1e-4;
162        assert_coord_almost_eq(latlons[0], (20.19, -121.550004), delta);
163        assert_coord_almost_eq(latlons[1], (20.19442682, -121.52621665), delta);
164        assert_coord_almost_eq(
165            latlons[latlons.len() - 2],
166            (50.10756403, -60.91298217),
167            delta,
168        );
169        assert_coord_almost_eq(
170            latlons[latlons.len() - 1],
171            (50.1024611, -60.88202274),
172            delta,
173        );
174
175        Ok(())
176    }
177}
grib/grid/lambert.rs

grib/grid/
lambert.rs
