grib/grid/
earth.rs

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