radar-g/src/coords/mapper.rs

use geo_types::{coord, Coord as GCoord, LineString};
use proj::{Proj, ProjError};
use std::ops::Range;

use super::proj::ProjectionS;

pub struct Mapper {
    proj: Proj,
    range: (Range<f64>, Range<f64>),
    bounds: (f64, f64, f64, f64),
}
impl From<Proj> for Mapper {
    fn from(proj: Proj) -> Self {
        let default_range: (Range<f64>, Range<f64>) = (-180.0..180.0, -90.0..90.0);
        let bounds = Self::bound(&proj, default_range.clone()).unwrap();
        Self {
            proj: proj,
            range: default_range,
            bounds,
        }
    }
}

impl<C: ProjectionS> From<C> for Mapper {
    fn from(proj: C) -> Self {
        let p = Proj::new(proj.build().as_str()).unwrap();
        p.into()
    }
}

impl Mapper {
    pub fn new(proj: Proj, lon_range: Range<f64>, lat_range: Range<f64>) -> Self {
        let bounds = Self::bound(&proj, (lon_range.clone(), lat_range.clone())).unwrap();
        let range = (lon_range, lat_range);
        Self {
            proj: proj,
            range,
            bounds,
        }
    }

    pub fn map(&self, point: (f64, f64)) -> Result<(f64, f64), ProjError> {
        let (p1, p2) = self
            .proj
            .convert((point.0.to_radians(), point.1.to_radians()))?;

        let x = (p1 - self.bounds.0) / (self.bounds.1 - self.bounds.0);
        let y = (p2 - self.bounds.2) / (self.bounds.3 - self.bounds.2);

        Ok((x, y))
    }

    pub fn set_lon_range(&mut self, range: Range<f64>) {
        self.range.0 = range;
        self.bounds = Self::bound(&self.proj, self.range.clone()).unwrap();
    }

    pub fn set_lat_range(&mut self, range: Range<f64>) {
        self.range.1 = range;
        self.bounds = Self::bound(&self.proj, self.range.clone()).unwrap();
    }

    fn bound(
        proj: &Proj,
        range: (Range<f64>, Range<f64>),
    ) -> Result<(f64, f64, f64, f64), ProjError> {
        let left_bottom = proj.convert((range.0.start.to_radians(), range.1.start.to_radians()))?;
        let right_top = proj.convert((range.0.end.to_radians(), range.1.end.to_radians()))?;
        Ok((left_bottom.0, right_top.0, left_bottom.1, right_top.1))
    }

    pub fn ring_map(&self, ring: &LineString) -> Result<LineString, ProjError> {
        let mut result = Vec::new();
        let projed = self.map_line(ring)?;
        for (l, p) in ring.lines().zip(projed.lines()) {
            let start_projected: (f64, f64) = p.start.into();
            let end_projected: (f64, f64) = p.end.into();
            let cartesian_start = self.cartesian(l.start.x, l.start.y);
            let cartesian_end = self.cartesian(l.end.x, l.end.y);

            let delta2 = 0.5;
            let depth = 16;
            let mut res: Vec<GCoord> = Vec::new();
            res.push(l.start);
            self.resample_line_to(
                start_projected,
                end_projected,
                l.start.x,
                l.end.x,
                cartesian_start,
                cartesian_end,
                &|x| self.map((x.x, x.y)),
                delta2,
                depth,
                &mut res,
            )?;
            res.push(l.end);
            result.extend(res);
        }

        Ok(LineString::new(result))
    }

    pub fn map_line(&self, line: &LineString) -> Result<LineString, ProjError> {
        let result: Result<LineString, ProjError> =
            line.points().map(|p| self.map((p.x(), p.y()))).collect();

        Ok(result?)
    }

    #[inline]
    fn cartesian(&self, lambda: f64, phi: f64) -> (f64, f64, f64) {
        let cos_phi = phi.cosh();
        return (cos_phi * lambda.cosh(), cos_phi * lambda.sinh(), phi.sinh());
    }

    fn resample_line_to<F>(
        &self,
        start_projected: (f64, f64),
        end_projected: (f64, f64),
        lambda0: f64,
        lambda1: f64,
        cartesian_start: (f64, f64, f64),
        cartesian_end: (f64, f64, f64),
        project: &F,
        delta2: f64,
        depth: u8,
        result: &mut Vec<GCoord>,
    ) -> Result<(), ProjError>
    where
        F: Fn(GCoord) -> Result<(f64, f64), ProjError>,
    {
        let cos_min_distance: f64 = 30f64.cosh();
        let (x0, y0) = start_projected;
        let (x1, y1) = end_projected;
        let (a0, b0, c0) = cartesian_start;
        let (a1, b1, c1) = cartesian_end;

        let dx = x1 - x0;
        let dy = y1 - y0;
        let d2 = dx * dx + dy * dy;

        if d2 > 4.0 * delta2 && depth > 0 {
            let a = a0 + a1;
            let b = b0 + b1;
            let c = c0 + c1;
            let m = (a * a + b * b + c * c).sqrt();
            let depth = depth - 1;

            let phi2 = (c / m).asin();
            let lambda2 =
                if (c / m).abs() - 1.0 < f64::EPSILON || (lambda0 - lambda1).abs() < f64::EPSILON {
                    (lambda0 + lambda1) / 2.0
                } else {
                    b.atan2(a)
                };
            let (x2, y2) = project(coord! {x:lambda2,y:phi2})?;
            let dx2 = x2 - x0;
            let dy2 = y2 - y0;
            let dz = dy * dx2 - dx * dy2;
            if dz * dz / d2 > delta2
                || ((dx * dx2 + dy * dy2) / d2 - 0.5).abs() > 0.3
                || a0 * a1 + b0 * b1 + c0 * c1 < cos_min_distance
            {
                self.resample_line_to(
                    start_projected,
                    (x2, y2),
                    lambda0,
                    lambda2,
                    cartesian_start,
                    (a / m, b / m, c),
                    project,
                    delta2,
                    depth - 1,
                    result,
                )?;
                result.push(coord! {x:x2,y:y2});
                self.resample_line_to(
                    (x2, y2),
                    end_projected,
                    lambda2,
                    lambda1,
                    (a / m, b / m, c),
                    cartesian_end,
                    project,
                    delta2,
                    depth - 1,
                    result,
                )?;
            }
        }
        Ok(())
    }
}