radar-g/src/coords/mod.rs

use num_traits::AsPrimitive;
use num_traits::Num;
pub mod mapper;
pub mod proj;
pub mod wgs84;

pub use mapper::Mapper;
// pub use wgs84::LatLonCoord;

pub type ScreenCoord = (f64, f64);
type Lat = f64;
type Lon = f64;

pub trait Coord<T: Num> {
    fn map(&self, axis_1: T, axis_2: T) -> ScreenCoord;
    fn dim1_range(&self) -> (T, T);
    fn dim2_range(&self) -> (T, T);
}

#[derive(Clone, Copy)]
pub struct Range(pub f64, pub f64);

impl Range {
    pub fn key_points(&self, max_points: usize) -> Vec<f64> {
        if max_points == 0 {
            return vec![];
        }

        let range = (self.0.min(self.1) as f64, self.1.max(self.0) as f64);
        assert!(!(range.0.is_nan() || range.1.is_nan()));
        if (range.0 - range.1).abs() < std::f64::EPSILON {
            return vec![range.0 as f64];
        }

        let mut scale = (10f64).powf((range.1 - range.0).log(10.0).floor());
        // The value granularity controls how we round the values.
        // To avoid generating key points like 1.00000000001, we round to the nearest multiple of the
        // value granularity.
        // By default, we make the granularity as the 1/10 of the scale.
        let mut value_granularity = scale / 10.0;
        fn rem_euclid(a: f64, b: f64) -> f64 {
            let ret = if b > 0.0 {
                a - (a / b).floor() * b
            } else {
                a - (a / b).ceil() * b
            };
            if (ret - b).abs() < std::f64::EPSILON {
                0.0
            } else {
                ret
            }
        }

        // At this point we need to make sure that the loop invariant:
        // The scale must yield number of points than requested
        if 1 + ((range.1 - range.0) / scale).floor() as usize > max_points {
            scale *= 10.0;
            value_granularity *= 10.0;
        }

        'outer: loop {
            let old_scale = scale;
            for nxt in [2.0, 5.0, 10.0].iter() {
                let mut new_left = range.0 - rem_euclid(range.0, old_scale / nxt);
                if new_left < range.0 {
                    new_left += old_scale / nxt;
                }
                let new_right = range.1 - rem_euclid(range.1, old_scale / nxt);

                let npoints = 1.0 + ((new_right - new_left) / old_scale * nxt);

                if npoints.round() as usize > max_points {
                    break 'outer;
                }

                scale = old_scale / nxt;
            }
            scale = old_scale / 10.0;
            value_granularity /= 10.0;
        }

        let mut ret = vec![];
        // In some extreme cases, left might be too big, so that (left + scale) - left == 0 due to
        // floating point error.
        // In this case, we may loop forever. To avoid this, we need to use two variables to store
        // the current left value. So we need keep a left_base and a left_relative.
        let left = {
            let mut value = range.0 - rem_euclid(range.0, scale);
            if value < range.0 {
                value += scale;
            }
            value
        };
        let left_base = (left / value_granularity).floor() * value_granularity;
        let mut left_relative = left - left_base;
        let right = range.1 - rem_euclid(range.1, scale);
        while (right - left_relative - left_base) >= -std::f64::EPSILON {
            let new_left_relative = (left_relative / value_granularity).round() * value_granularity;
            if new_left_relative < 0.0 {
                left_relative += value_granularity;
            }
            ret.push((left_relative + left_base) as f64);
            left_relative += scale;
        }
        return ret;
    }
}

impl<T: AsPrimitive<f64> + Num> From<(T, T)> for Range {
    fn from(value: (T, T)) -> Self {
        let value = (value.0.as_(), value.1.as_());
        let (_min, _max) = (value.0.min(value.1), value.0.max(value.1));
        Self(_min, _max)
    }
}

impl<T: Num + AsPrimitive<f64>> From<std::ops::Range<T>> for Range {
    fn from(value: std::ops::Range<T>) -> Self {
        let value = (value.start.as_(), value.end.as_());
        let (_min, _max) = (value.0.min(value.1), value.0.max(value.1));
        Self(_min, _max)
    }
}

// impl<T, Raw> RadarData2d<T, Raw>
// where
//     T: Num + Clone + PartialEq + PartialOrd,
//     Raw: ndarray::Data<Elem = T> + Clone + ndarray::RawDataClone,
// {
//     pub fn mapped(&self, coord: impl Borrow<Mapper>) -> Result<AfterMapping2d<T>, ProjError> {
//         let mapper: &Mapper = coord.borrow();
//         self.map_by_fn(|x| mapper.map(x))
//     }

//     pub fn map_by_fn<F>(&self, f: F) -> Result<AfterMapping2d<T>, ProjError>
//     where
//         F: Fn((f64, f64)) -> Result<(f64, f64), ProjError>,
//     {
//         let mesh_dim1_len = self.dim1.len();
//         let mesh_dim2_len = self.dim2.len();

//         let mut d1 = Array2::<f64>::zeros((mesh_dim2_len, mesh_dim1_len));
//         let mut d2 = Array2::<f64>::zeros((mesh_dim2_len, mesh_dim1_len));

//         for (i, v) in self.dim1.iter().enumerate() {
//             for (j, u) in self.dim2.iter().enumerate() {
//                 let (x, y) = f((*v, *u))?;
//                 d1[[j, i]] = x;
//                 d2[[j, i]] = y;
//             }
//         }

//         Ok(AfterMapping2d {
//             dim1: d1,
//             dim2: d2,
//             data: self.data.view(),
//             coord_type: self.coord_type.clone(),
//         })
//     }
// }