grib/grid/
earth.rs

1use crate::def::grib2::template::param_set::EarthShape;
2
3impl EarthShape {
4    pub fn radii(&self) -> Option<(f64, f64)> {
5        let radii = match self.shape {
6            0 => (6367470.0, 6367470.0),
7            1 => {
8                let radius = f64::from(self.spherical_earth_radius_scaled_value)
9                    * f64::powf(10., f64::from(self.spherical_earth_radius_scale_factor));
10                (radius, radius)
11            }
12            2 => (6378160.0, 6356775.0),
13            3 => {
14                let (major, minor) = self.radii_defined();
15                (major * 1000., minor * 1000.)
16            }
17            4 => (6378137.0, 6356752.314),
18            5 => (6378137.0, 6356752.3142), // WGS84
19            6 => (6371229.0, 6371229.0),
20            7 => self.radii_defined(),
21            8 => (6371200.0, 6371200.0),
22            9.. => return None,
23        };
24        Some(radii)
25    }
26
27    fn radii_defined(&self) -> (f64, f64) {
28        let major = f64::from(self.major_axis_scaled_value)
29            * f64::powf(10., f64::from(self.major_axis_scale_factor));
30        let minor = f64::from(self.minor_axis_scaled_value)
31            * f64::powf(10., f64::from(self.minor_axis_scale_factor));
32        (major, minor)
33    }
34}
35
36#[cfg(test)]
37mod tests {
38    use std::{
39        fs::File,
40        io::{BufReader, Read},
41    };
42
43    use grib_template_helpers::TryFromSlice;
44
45    use super::*;
46
47    fn get_uncompressed<P>(file_path: P) -> Result<Vec<u8>, std::io::Error>
48    where
49        P: AsRef<std::path::Path>,
50    {
51        let mut buf = Vec::new();
52
53        let f = File::open(&file_path)?;
54        let mut f = BufReader::new(f);
55        match file_path.as_ref().extension().map(|s| s.as_encoded_bytes()) {
56            Some(b"gz") => {
57                let mut f = flate2::read::GzDecoder::new(f);
58                f.read_to_end(&mut buf)?;
59            }
60            Some(b"xz") => {
61                let mut f = xz2::bufread::XzDecoder::new(f);
62                f.read_to_end(&mut buf)?;
63            }
64            _ => {
65                f.read_to_end(&mut buf)?;
66            }
67        };
68
69        Ok(buf)
70    }
71
72    #[test]
73    fn radii_for_shape_1() -> Result<(), Box<dyn std::error::Error>> {
74        let buf = get_uncompressed("testdata/ds.critfireo.bin.xz")?;
75        let mut pos = 0x83;
76        let earth_actual = EarthShape::try_from_slice(&buf, &mut pos)?;
77        let earth_expected = EarthShape {
78            shape: 1,
79            spherical_earth_radius_scale_factor: 0,
80            spherical_earth_radius_scaled_value: 6371200,
81            major_axis_scale_factor: 0,
82            major_axis_scaled_value: 0,
83            minor_axis_scale_factor: 0,
84            minor_axis_scaled_value: 0,
85        };
86        assert_eq!(earth_actual, earth_expected);
87        assert_eq!(earth_actual.radii(), Some((6_371_200., 6_371_200.)));
88
89        Ok(())
90    }
91
92    #[test]
93    fn radii_for_shape_2() -> Result<(), Box<dyn std::error::Error>> {
94        let buf = get_uncompressed(
95            "testdata/MRMS_ReflectivityAtLowestAltitude_00.50_20230406-120039.grib2.gz",
96        )?;
97        let mut pos = 0x33;
98        let earth_actual = EarthShape::try_from_slice(&buf, &mut pos)?;
99        let earth_expected = EarthShape {
100            shape: 2,
101            spherical_earth_radius_scale_factor: 1,
102            spherical_earth_radius_scaled_value: 6367470,
103            major_axis_scale_factor: 1,
104            major_axis_scaled_value: 6378160,
105            minor_axis_scale_factor: 1,
106            minor_axis_scaled_value: 6356775,
107        };
108        assert_eq!(earth_actual, earth_expected);
109        assert_eq!(earth_actual.radii(), Some((6_378_160.0, 6_356_775.0)));
110
111        Ok(())
112    }
113
114    #[test]
115    fn radii_for_shape_4() -> Result<(), Box<dyn std::error::Error>> {
116        let buf = get_uncompressed(
117            "testdata/Z__C_RJTD_20160822020000_NOWC_GPV_Ggis10km_Pphw10_FH0000-0100_grib2.bin",
118        )?;
119        let mut pos = 0x33;
120        let earth_actual = EarthShape::try_from_slice(&buf, &mut pos)?;
121        let earth_expected = EarthShape {
122            shape: 4,
123            spherical_earth_radius_scale_factor: 0xff,
124            spherical_earth_radius_scaled_value: 0xffffffff,
125            major_axis_scale_factor: 1,
126            major_axis_scaled_value: 63781370,
127            minor_axis_scale_factor: 1,
128            minor_axis_scaled_value: 63567523,
129        };
130        assert_eq!(earth_actual, earth_expected);
131        assert_eq!(earth_actual.radii(), Some((6_378_137.0, 6_356_752.314)));
132
133        Ok(())
134    }
135
136    #[test]
137    fn radii_for_shape_6() -> Result<(), Box<dyn std::error::Error>> {
138        let buf = get_uncompressed("testdata/gdas.t12z.pgrb2.0p25.f000.0-10.xz")?;
139        let mut pos = 0x33;
140        let earth_actual = EarthShape::try_from_slice(&buf, &mut pos)?;
141        let earth_expected = EarthShape {
142            shape: 6,
143            spherical_earth_radius_scale_factor: 0,
144            spherical_earth_radius_scaled_value: 0,
145            major_axis_scale_factor: 0,
146            major_axis_scaled_value: 0,
147            minor_axis_scale_factor: 0,
148            minor_axis_scaled_value: 0,
149        };
150        assert_eq!(earth_actual, earth_expected);
151        assert_eq!(earth_actual.radii(), Some((6_371_229.0, 6_371_229.0)));
152
153        Ok(())
154    }
155}
grib/grid/earth.rs

grib/grid/
earth.rs
