grib/grid/
mercator.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_10, param_set},
8    error::GribError,
9    grid::AngleUnit,
10    projection,
11};
12
13impl crate::GridShortName for Template3_10 {
14    fn short_name(&self) -> &'static str {
15        "mercator"
16    }
17}
18
19impl GridPointIndex for Template3_10 {
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_10 {
30    type Iter<'a> = std::vec::IntoIter<(f32, f32)>;
31
32    fn latlons_unchecked<'a>(&'a self) -> Result<Self::Iter<'a>, GribError> {
33        if self.orientation != 0 {
34            return Err(GribError::NotSupported(format!(
35                "Mercator grid orientation {}",
36                self.orientation
37            )));
38        }
39
40        if !self.is_consistent_for_j() {
41            return Err(GribError::InvalidValueError(
42                "Latitudes for first/last grid points are not consistent with scanning mode"
43                    .to_owned(),
44            ));
45        }
46
47        let angle_units = self.angle_unit();
48        let first_point_lon = self.first_point_lon as f64 * angle_units;
49        let last_point_lon = self.last_point_lon as f64 * angle_units;
50        let lon_diff = last_point_lon - first_point_lon;
51        let (first_point_lon, last_point_lon) =
52            if self.scanning_mode.scans_positively_for_i() && lon_diff < 0. {
53                (first_point_lon, last_point_lon + 360.)
54            } else if !self.scanning_mode.scans_positively_for_i() && lon_diff > 0. {
55                (first_point_lon + 360., last_point_lon)
56            } else {
57                (first_point_lon, last_point_lon)
58            };
59        let first_point = (self.first_point_lat as f64 * angle_units, first_point_lon);
60        let last_point = (self.last_point_lat as f64 * angle_units, last_point_lon);
61        // ensure that all points are within the interval `[lon_0 - 180., lon_0
62        // + 180.]` so that the delta can be calculated correctly.
63        let lon_0 = (first_point_lon + last_point_lon) / 2.;
64
65        let (a, b) = self.earth_shape.radii().ok_or_else(|| {
66            GribError::NotSupported(format!(
67                "unknown value of Code Table 3.2 (shape of the Earth): {}",
68                self.earth_shape.shape
69            ))
70        })?;
71        let params = projection::MercParams {
72            ellipsoid: projection::Ellipsoid::from_a_and_b(a, b),
73            lat_ts: self.lad as f64 * angle_units,
74            lon_0,
75        };
76
77        #[cfg(feature = "gridpoints-proj")]
78        {
79            super::helpers::latlons_from_projection_with_first_point_and_last_point(
80                &params.proj_args(),
81                first_point,
82                last_point,
83                self.grid_shape(),
84                self.ij()?,
85            )
86        }
87
88        #[cfg(not(feature = "gridpoints-proj"))]
89        {
90            let projection = projection::Merc::new(&params)?;
91            let (first_point_lat, first_point_lon) = first_point;
92            let (first_corner_x, first_corner_y) = projection.project(
93                &(first_point_lon.to_radians(), first_point_lat.to_radians()),
94                false,
95            )?;
96            let (last_point_lat, last_point_lon) = last_point;
97            let (last_corner_x, last_corner_y) = projection.project(
98                &(last_point_lon.to_radians(), last_point_lat.to_radians()),
99                false,
100            )?;
101
102            let dx = (last_corner_x - first_corner_x) / (self.ni - 1) as f64;
103            let dy = (last_corner_y - first_corner_y) / (self.nj - 1) as f64;
104            let latlons = self
105                .ij()?
106                .map(|(i, j)| {
107                    projection
108                        .project(
109                            &(
110                                first_corner_x + dx * i as f64,
111                                first_corner_y + dy * j as f64,
112                            ),
113                            true,
114                        )
115                        .map(|(lon, lat)| (lat.to_degrees() as f32, lon.to_degrees() as f32))
116                })
117                .collect::<Result<Vec<_>, _>>()?;
118            Ok(latlons.into_iter())
119        }
120    }
121}
122
123impl AngleUnit for Template3_10 {
124    fn angle_unit(&self) -> f64 {
125        1e-6
126    }
127}
128
129impl Template3_10 {
130    pub(crate) fn is_consistent_for_j(&self) -> bool {
131        let lat_diff = self.last_point_lat - self.first_point_lat;
132        !((lat_diff > 0) ^ self.scanning_mode.scans_positively_for_j())
133    }
134}
135
136#[cfg(test)]
137mod tests {
138    use super::*;
139    use crate::grid::helpers::test_helpers::assert_coord_almost_eq;
140
141    fn grid_definition() -> Template3_10 {
142        // grid point definition extracted from testdata/ds.wwa.bin.xz
143        Template3_10 {
144            earth_shape: param_set::EarthShape {
145                shape: 1,
146                spherical_earth_radius: param_set::ScaledValue {
147                    scale_factor: 0,
148                    scaled_value: 6371200,
149                },
150                major_axis: param_set::ScaledValue {
151                    scale_factor: 0,
152                    scaled_value: 0,
153                },
154                minor_axis: param_set::ScaledValue {
155                    scale_factor: 0,
156                    scaled_value: 0,
157                },
158            },
159            ni: 2517,
160            nj: 1793,
161            first_point_lat: -30419200,
162            first_point_lon: 129906005,
163            resolution_and_component_flags: param_set::ResolutionAndComponentFlags(0b00000000),
164            lad: 20000000,
165            last_point_lat: 80010000,
166            last_point_lon: 10710000,
167            scanning_mode: param_set::ScanningMode(0b01010000),
168            orientation: 0,
169            di: 10000000,
170            dj: 10000000,
171        }
172    }
173
174    #[test]
175    fn mercator_grid_latlon_computation() -> Result<(), Box<dyn std::error::Error>> {
176        let grid_def = grid_definition();
177        let latlons = grid_def.latlons()?.collect::<Vec<_>>();
178
179        // Following lat/lon values are taken from the calculation results using
180        // pygrib.
181        let num_points = latlons.len();
182        let ni = grid_def.ni as usize;
183        let delta = 3e-5;
184        // lat[0], lon[0]
185        assert_coord_almost_eq(latlons[0], (-30.4192, 129.906005), delta);
186        // lat[0], lon[1]
187        assert_coord_almost_eq(latlons[1], (-30.4192, 130.00171406), delta);
188        // lat[1], lon[-1]
189        assert_coord_almost_eq(latlons[ni], (-30.33658686, 10.71), delta);
190        // lat[-2], lon[0]
191        assert_coord_almost_eq(
192            latlons[num_points - ni - 1],
193            (79.9933742, 129.906005),
194            delta,
195        );
196        // lat[-1], lon[-2]
197        assert_coord_almost_eq(latlons[num_points - 2], (80.01, 10.61429094), delta);
198        // lat[-1], lon[-1]
199        assert_coord_almost_eq(latlons[num_points - 1], (80.01, 10.71), delta);
200
201        Ok(())
202    }
203}
grib/grid/mercator.rs

grib/grid/
mercator.rs
