Publish the Gaussian latitudes computation function as a utility function

src/datatypes/sections.rs

@@ -9,7 +9,7 @@ use crate::{
     grid::{
         GaussianGridDefinition, GridPointIterator, LambertGridDefinition, LatLonGridDefinition,
     },
-    utils::{read_as, GribInt},
+    helpers::{read_as, GribInt},
     GridPointIndexIterator, PolarStereographicGridDefinition,
 };
 
@@ -235,7 +235,8 @@ impl GridDefinitionTemplateValues {
         }
     }
 
-    /// Returns an iterator over latitudes and longitudes of grid points.
+    /// Returns an iterator over latitudes and longitudes of grid points in
+    /// degrees.
     ///
     /// Note that this is a low-level API and it is not checked that the number
     /// of iterator iterations is consistent with the number of grid points

src/decoder/complex.rs

@@ -18,7 +18,7 @@ use crate::{
         DecodeError, Grib2SubmessageDecoder,
     },
     error::*,
-    utils::{read_as, GribInt},
+    helpers::{read_as, GribInt},
 };
 
 #[derive(Debug, Clone, PartialEq, Eq, Hash)]

src/decoder/complex/diff.rs

@@ -2,7 +2,7 @@ use super::{
     missing::DecodedValue::{self, Normal},
     ComplexPackingDecodeError,
 };
-use crate::{decoder::DecodeError, error::GribError, utils::grib_int_from_bytes};
+use crate::{decoder::DecodeError, error::GribError, helpers::grib_int_from_bytes};
 
 pub(crate) struct SpatialDifferencingExtraDescriptors<'a> {
     slice: &'a [u8],

src/decoder/jpeg2000.rs

@@ -6,7 +6,7 @@ use crate::{
         Grib2SubmessageDecoder,
     },
     error::*,
-    utils::read_as,
+    helpers::read_as,
 };
 
 mod ext;

src/decoder/param.rs

@@ -1,4 +1,4 @@
-use crate::utils::{read_as, GribInt};
+use crate::helpers::{read_as, GribInt};
 
 pub(crate) struct SimplePackingParam {
     pub(crate) ref_val: f32,

src/decoder/png.rs

@@ -4,7 +4,7 @@ use crate::{
         simple::{SimplePackingDecodeError, SimplePackingDecodeIterator},
         stream::NBitwiseIterator,
     },
-    utils::read_as,
+    helpers::read_as,
     DecodeError, Grib2SubmessageDecoder, GribError,
 };
 

src/decoder/run_length.rs

@@ -1,7 +1,7 @@
 use crate::{
     decoder::{stream::NBitwiseIterator, DecodeError, Grib2SubmessageDecoder},
     error::*,
-    utils::read_as,
+    helpers::read_as,
 };
 
 #[derive(Debug, Clone, PartialEq, Eq, Hash)]

src/decoder/simple.rs

@@ -7,7 +7,7 @@ use crate::{
         DecodeError, Grib2SubmessageDecoder,
     },
     error::*,
-    utils::read_as,
+    helpers::read_as,
 };
 
 pub(crate) enum SimplePackingDecodeIteratorWrapper<I> {

src/grid.rs

@@ -1,8 +1,11 @@
 use helpers::RegularGridIterator;
 
 pub use self::{
-    earth::EarthShapeDefinition, gaussian::GaussianGridDefinition, lambert::LambertGridDefinition,
-    latlon::LatLonGridDefinition, polar_stereographic::PolarStereographicGridDefinition,
+    earth::EarthShapeDefinition,
+    gaussian::{compute_gaussian_latitudes, GaussianGridDefinition},
+    lambert::LambertGridDefinition,
+    latlon::LatLonGridDefinition,
+    polar_stereographic::PolarStereographicGridDefinition,
 };
 
 /// An iterator over latitudes and longitudes of grid points in a submessage.

src/grid/earth.rs

@@ -1,4 +1,4 @@
-use crate::utils::read_as;
+use crate::helpers::read_as;
 
 #[derive(Debug, PartialEq, Eq)]
 pub struct EarthShapeDefinition {

src/grid/gaussian.rs

@@ -4,7 +4,7 @@ use super::{
 };
 use crate::{
     error::GribError,
-    utils::{read_as, GribInt},
+    helpers::{read_as, GribInt},
 };
 
 #[derive(Debug, PartialEq, Eq)]
@@ -50,7 +50,8 @@ impl GaussianGridDefinition {
         Ok(iter)
     }
 
-    /// Returns an iterator over latitudes and longitudes of grid points.
+    /// Returns an iterator over latitudes and longitudes of grid points in
+    /// degrees.
     ///
     /// Note that this is a low-level API and it is not checked that the number
     /// of iterator iterations is consistent with the number of grid points
@@ -61,7 +62,7 @@ impl GaussianGridDefinition {
         }
 
         let ij = self.ij()?;
-        let mut lat = compute_gaussian_latitudes(self.nj as usize)
+        let mut lat = compute_gaussian_latitudes_in_degrees(self.nj as usize)
             .map_err(|e| GribError::Unknown(e.to_owned()))?;
         if self.scanning_mode.scans_positively_for_j() {
             lat.reverse()
@@ -108,13 +109,39 @@ impl GaussianGridDefinition {
     }
 }
 
-fn compute_gaussian_latitudes(div: usize) -> Result<Vec<f64>, &'static str> {
-    let lat: Option<Vec<_>> = legendre_roots_iterator(div)
-        .map(|x| x.map(|i| i.asin().to_degrees()))
+fn compute_gaussian_latitudes_in_degrees(div: usize) -> Result<Vec<f64>, &'static str> {
+    let lat: Option<Vec<_>> = compute_gaussian_latitudes(div)
+        .map(|x| x.map(|i| i.to_degrees()))
         .collect();
     lat.ok_or("finding root for Legendre polynomial failed")
 }
 
+/// Computes Gaussian latitudes in radians.
+///
+/// The Newton-Raphson method is used for the computation.
+/// If the computation does not converge and no solution is obtained after the
+/// predefined number of iterations (10), the solution will have the value
+/// `None`.
+///
+/// # Examples
+///
+/// ```
+/// let mut iter = grib::utils::compute_gaussian_latitudes(0);
+/// assert_eq!(iter.next(), None);
+///
+/// let mut iter = grib::utils::compute_gaussian_latitudes(1);
+/// assert_eq!(iter.next(), Some(Some(0.0)));
+/// assert_eq!(iter.next(), None);
+///
+/// let mut iter = grib::utils::compute_gaussian_latitudes(2);
+/// assert!((iter.next().unwrap().unwrap() - (1.0 / 3.0_f64.sqrt()).asin()).abs() < 1e-15);
+/// assert!((iter.next().unwrap().unwrap() - (-1.0 / 3.0_f64.sqrt()).asin()).abs() < 1e-15);
+/// assert_eq!(iter.next(), None);
+/// ```
+pub fn compute_gaussian_latitudes(div: usize) -> impl Iterator<Item = Option<f64>> {
+    legendre_roots_iterator(div).map(|x| x.map(|i| i.asin()))
+}
+
 // Finds roots (zero points) of the Legendre polynomial using Newton–Raphson
 // method.
 //
@@ -330,7 +357,7 @@ mod tests {
     #[test]
     fn gaussian_latitudes_computation_compared_with_numerical_solutions() {
         let n = 160;
-        let result = compute_gaussian_latitudes(n);
+        let result = compute_gaussian_latitudes_in_degrees(n);
         assert!(result.is_ok());
 
         let actual = result.unwrap().into_iter().take(n / 2);

src/grid/lambert.rs

@@ -1,7 +1,7 @@
 use super::{earth::EarthShapeDefinition, GridPointIndexIterator, ScanningMode};
 use crate::{
     error::GribError,
-    utils::{read_as, GribInt},
+    helpers::{read_as, GribInt},
 };
 
 #[derive(Debug, PartialEq, Eq)]
@@ -111,7 +111,8 @@ impl LambertGridDefinition {
         Ok(iter)
     }
 
-    /// Returns an iterator over latitudes and longitudes of grid points.
+    /// Returns an iterator over latitudes and longitudes of grid points in
+    /// degrees.
     ///
     /// Note that this is a low-level API and it is not checked that the number
     /// of iterator iterations is consistent with the number of grid points

src/grid/latlon.rs

@@ -4,7 +4,7 @@ use super::{
 };
 use crate::{
     error::GribError,
-    utils::{read_as, GribInt},
+    helpers::{read_as, GribInt},
 };
 
 #[derive(Debug, PartialEq, Eq)]
@@ -83,7 +83,8 @@ impl LatLonGridDefinition {
         Ok(iter)
     }
 
-    /// Returns an iterator over latitudes and longitudes of grid points.
+    /// Returns an iterator over latitudes and longitudes of grid points in
+    /// degrees.
     ///
     /// Note that this is a low-level API and it is not checked that the number
     /// of iterator iterations is consistent with the number of grid points

src/grid/polar_stereographic.rs

@@ -1,7 +1,7 @@
 use super::{earth::EarthShapeDefinition, GridPointIndexIterator, ScanningMode};
 use crate::{
     error::GribError,
-    utils::{read_as, GribInt},
+    helpers::{read_as, GribInt},
     ProjectionCentreFlag,
 };
 
@@ -109,7 +109,8 @@ impl PolarStereographicGridDefinition {
         Ok(iter)
     }
 
-    /// Returns an iterator over latitudes and longitudes of grid points.
+    /// Returns an iterator over latitudes and longitudes of grid points in
+    /// degrees.
     ///
     /// Note that this is a low-level API and it is not checked that the number
     /// of iterator iterations is consistent with the number of grid points

src/helpers.rs

@@ -0,0 +1,158 @@
+pub(crate) trait GribInt<I> {
+    fn as_grib_int(&self) -> I;
+}
+
+macro_rules! add_impl_for_ints {
+    ($(($ty_src:ty, $ty_dst:ty),)*) => ($(
+        impl GribInt<$ty_dst> for $ty_src {
+            fn as_grib_int(&self) -> $ty_dst {
+                if self.leading_zeros() == 0 {
+                    let abs = (self << 1 >> 1) as $ty_dst;
+                    -abs
+                } else {
+                    *self as $ty_dst
+                }
+            }
+        }
+    )*);
+}
+
+add_impl_for_ints! {
+    (u8, i8),
+    (u16, i16),
+    (u32, i32),
+    (u64, i64),
+}
+
+macro_rules! read_as {
+    ($ty:ty, $buf:ident, $start:expr) => {{
+        let end = $start + std::mem::size_of::<$ty>();
+        <$ty>::from_be_bytes($buf[$start..end].try_into().unwrap())
+    }};
+}
+pub(crate) use read_as;
+
+pub(crate) fn grib_int_from_bytes(bytes: &[u8]) -> i32 {
+    let len = bytes.len();
+    // Although there is logic that can be used to generalize, not so many patterns
+    // exist that generalization is necessary.
+    match len {
+        1 => i32::from(read_as!(u8, bytes, 0).as_grib_int()),
+        2 => i32::from(read_as!(u16, bytes, 0).as_grib_int()),
+        3 => {
+            let first = read_as!(u8, bytes, 0);
+            let positive = first.leading_zeros() != 0;
+            let rest = i32::from(read_as!(u16, bytes, 1));
+            let abs = i32::from(first << 1 >> 1) * 0x10000 + rest;
+            if positive {
+                abs
+            } else {
+                -abs
+            }
+        }
+        4 => read_as!(u32, bytes, 0).as_grib_int(),
+        _ => unimplemented!(),
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn into_grib_i8() {
+        let input: Vec<u8> = vec![0b01000000, 0b00000001, 0b10000001, 0b11000000];
+        let output: Vec<i8> = vec![64, 1, -1, -64];
+
+        let mut actual = Vec::new();
+        let mut pos = 0;
+        while pos < input.len() {
+            let val = u8::from_be_bytes(input[pos..pos + 1].try_into().unwrap());
+            pos += 1;
+            let val = val.as_grib_int();
+            actual.push(val);
+        }
+
+        assert_eq!(actual, output);
+    }
+
+    #[test]
+    fn into_grib_i16() {
+        let input: Vec<u8> = vec![
+            0b00000000, 0b01000000, 0b00000000, 0b00000001, 0b10000000, 0b00000001, 0b10000000,
+            0b01000000,
+        ];
+        let output: Vec<i16> = vec![64, 1, -1, -64];
+
+        let mut actual = Vec::new();
+        let mut pos = 0;
+        while pos < input.len() {
+            let val = u16::from_be_bytes(input[pos..pos + 2].try_into().unwrap());
+            pos += 2;
+            let val = val.as_grib_int();
+            actual.push(val);
+        }
+
+        assert_eq!(actual, output);
+    }
+
+    macro_rules! test_conversion_from_bytes_to_grib_int {
+        ($(($name:ident, $input:expr, $expected:expr),)*) => ($(
+            #[test]
+            fn $name() {
+                let bytes = $input;
+                let actual = grib_int_from_bytes(&bytes);
+                let expected = $expected;
+                assert_eq!(actual, expected)
+            }
+        )*);
+    }
+
+    test_conversion_from_bytes_to_grib_int! {
+        (
+            conversion_from_bytes_to_grib_int_for_1_byte_positive,
+            vec![0b01010101],
+            0b01010101
+        ),
+        (
+            conversion_from_bytes_to_grib_int_for_1_byte_negative,
+            vec![0b11010101],
+            -0b01010101
+        ),
+        (
+            conversion_from_bytes_to_grib_int_for_2_bytes_positive,
+            vec![0b01010101, 0b10101010],
+            0b0101_0101_1010_1010
+        ),
+        (
+            conversion_from_bytes_to_grib_int_for_2_bytes_negative,
+            vec![0b11010101, 0b10101010],
+            -0b0101_0101_1010_1010
+        ),
+        (
+            conversion_from_bytes_to_grib_int_for_3_bytes_positive,
+            vec![0b01010101, 0b10101010, 0b10101010],
+            0b0101_0101_1010_1010_1010_1010
+        ),
+        (
+            conversion_from_bytes_to_grib_int_for_3_bytes_negative,
+            vec![0b11010101, 0b10101010, 0b10101010],
+            -0b0101_0101_1010_1010_1010_1010
+        ),
+        (
+            conversion_from_bytes_to_grib_int_for_3_bytes_negative_starting_from_0x80,
+            vec![0b10000000, 0b10101010, 0b10101010],
+            -0b0000_0000_1010_1010_1010_1010
+        ),
+        (
+            conversion_from_bytes_to_grib_int_for_4_bytes_positive,
+            vec![0b01010101, 0b10101010, 0b10101010, 0b10101010],
+            0b0101_0101_1010_1010_1010_1010_1010_1010
+        ),
+        (
+            conversion_from_bytes_to_grib_int_for_4_bytes_negative,
+            vec![0b11010101, 0b10101010, 0b10101010, 0b10101010],
+            -0b0101_0101_1010_1010_1010_1010_1010_1010
+        ),
+    }
+}