Fapra-OSM-2/src/coordinates.rs

use crate::ACCURACY_BOUNDARY;
use geojson::{Position, Value};
use std::convert::From;
use std::f64::consts::{FRAC_PI_2, PI, TAU};
use serde::{Serialize, Deserialize};


/// Spherical coordinates in radians.
#[derive(Clone, Copy, Debug, PartialEq, Default, Serialize, Deserialize)]
pub struct RadianCoordinate {
    pub lat: f64,
    pub lon: f64,
}

/// Spherical coordinates in degrees.
#[derive(Clone, Copy, Debug, PartialEq, Default, Serialize, Deserialize)]
pub struct DegreeCoordinate {
    pub lat: f64,
    pub lon: f64,
}

impl From<DegreeCoordinate> for RadianCoordinate {
    fn from(other: DegreeCoordinate) -> Self {
        RadianCoordinate::from_degrees(other.lat, other.lon)
    }
}

impl From<DegreeCoordinate> for geojson::Value {
    fn from(coordinate: DegreeCoordinate) -> Self {
        Value::Point(vec![coordinate.lon, coordinate.lat])
    }
}

impl From<RadianCoordinate> for Position {
    fn from(coordinate: RadianCoordinate) -> Self {
        let coordinate = DegreeCoordinate::from(coordinate);
        vec![coordinate.lon, coordinate.lat]
    }
}

impl From<RadianCoordinate> for geojson::Value {
    fn from(coordinate: RadianCoordinate) -> Self {
        let degrees: DegreeCoordinate = coordinate.into();
        degrees.into()
    }
}

impl From<RadianCoordinate> for DegreeCoordinate {
    fn from(other: RadianCoordinate) -> Self {
        DegreeCoordinate::from_radians(other.lat, other.lon)
    }
}


impl RadianCoordinate {
    /// Builds a RadianCoordinate from latitude and longitude given in
    /// degrees.
    pub fn from_degrees(lat: f64, lon: f64) -> Self {
        let lat = lat / 90.0 * FRAC_PI_2;
        let lon = normalize_lon(lon / 180.0 * PI);

        RadianCoordinate { lat, lon }
    }

    /// gives a normalizes version of the Coordinate
    pub fn normalize(&self) -> RadianCoordinate {
        RadianCoordinate {
            lat: self.lat,
            lon: normalize_lon(self.lon),
        }
    }

    /// returns the longitude of a point when the coordinate system is
    /// transformed, such that `north_pole` is the north pole of that system.
    pub fn get_transformed_longitude(&self, north_pole: &RadianCoordinate) -> f64 {
        if self.lat == FRAC_PI_2 {
            return self.lon;
        };

        let top = (self.lon - north_pole.lon).sin() * self.lat.cos();
        let bottom = (self.lat.sin() * north_pole.lat.cos())
            - (self.lat.cos() * north_pole.lat.sin() * (self.lon - north_pole.lon).cos());

        bottom.atan2(top)
    }

    /// returns `true` when the point is on a great circle with point `a` and `b`
    /// or numerically very close to that
    pub fn on_great_circle(&self, a: &RadianCoordinate, b: &RadianCoordinate) -> bool {
        let lat_a = a.get_transformed_longitude(&self);
        let lat_b = b.get_transformed_longitude(&self);

        (lat_a - lat_b).abs().rem_euclid(PI) < ACCURACY_BOUNDARY
    }

    /// calculates whether `other` is antipodal to this point.
    pub fn antipodal(&self, other: &RadianCoordinate) -> bool {
        // if the distance between both points is very close to PI, they are
        // antipodal

        let c = self.distance_to(other);

        let diff = (c - PI).abs();

        diff < ACCURACY_BOUNDARY
    }

    /// returns the shortest distance to an other RadianCoordinate along
    /// the surface of the unit-sphere in radians.
    pub fn distance_to(&self, other: &RadianCoordinate) -> f64 {
        // using the haversine formula
        // if the distance between both points is very close to PI, they are
        // antipodal

        let delta_lat = other.lat - self.lat;
        let delta_lon = other.lon - self.lon;

        let a = (delta_lat / 2.0).sin().powi(2)
            + self.lat.cos() * other.lat.cos() * (delta_lon / 2.0).sin().powi(2);
        2.0 * a.sqrt().atan2((1.0_f64 - a).sqrt())
    }

}

impl DegreeCoordinate {
    /// Builds a DegreeCoordinate from latitude and longitude given in
    /// radians.
    pub fn from_radians(lat: f64, lon: f64) -> Self {
        let lat = lat * 90.0 / FRAC_PI_2;
        let lon = normalize_lon(lon) * 180.0 / PI;

        DegreeCoordinate { lat, lon }
    }

    /// returns a SphericalCoordinate parsed from the given string.
    pub fn from_string_tuple(input: &str) -> Result<DegreeCoordinate, CoordinateParsingError> {
        let splits: Vec<&str> = input.split(",").collect();

        if splits.len() != 2 {
            return Err(CoordinateParsingError::NoSemicolon);
        }

        let lat = splits[0];
        let lon = splits[1];

        let lat = match lat.parse::<f64>() {
            Ok(lat) => lat,
            Err(_) => { return Err(CoordinateParsingError::NotAFloat); }
        };

        let lon = match lon.parse::<f64>() {
            Ok(lon) => lon,
            Err(_) => { return Err(CoordinateParsingError::NotAFloat); }
        };

        Ok(DegreeCoordinate{lat, lon})
    }


    /// tries to parse a DegreeCoordinate from a GeoJSON Position
    pub fn from_geojson_position(position: geojson::Position) -> Result<Self, String> {

        if position.len() != 2 {
            return Err("String has more than 2 values".to_string());
        }

        let lat = position[1];
        let lon = position[0];

        Ok(DegreeCoordinate { lat, lon })
    }
}

/// normalizes longitude values given in radians to the range (-PI, PI]
pub fn normalize_lon(lon: f64) -> f64 {
    // restrict values to -/+ TAU
    let mut lon = lon % (TAU);

    if lon <= -PI {
        lon = TAU + lon;
    }

    if lon > PI {
        lon = -TAU + lon;
    }

    lon
}

#[test]
fn test_normalize_lon() {
    assert!((normalize_lon(0.0) - 0.0).abs() < 10e-12);
    assert!((normalize_lon(PI) - PI).abs() < 10e-12);
    assert!((normalize_lon(-PI) - PI).abs() < 10e-12);
    assert!((normalize_lon(TAU) - 0.0).abs() < 10e-12);
    assert!((normalize_lon(PI + FRAC_PI_2) - (-FRAC_PI_2)).abs() < 10e-12);
    assert!((normalize_lon(-PI - FRAC_PI_2) - (FRAC_PI_2)).abs() < 10e-12);
}

#[derive(Clone, Copy, Debug, PartialEq)]
pub enum LongitudeDirection {
    East,
    West,
    None,
}

pub fn east_or_west(from: f64, to: f64) -> LongitudeDirection {
    let delta = normalize_lon(from - to);

    if delta > 0.0 && delta < PI {
        LongitudeDirection::West
    } else if delta < 0.0 && delta > -PI {
        LongitudeDirection::East
    } else {
        LongitudeDirection::None
    }
}

#[derive(Clone, Debug)]
pub enum CoordinateParsingError {
    NoSemicolon,
    NotAFloat,
}

#[test]
fn test_east_or_west() {
    assert_eq!(east_or_west(0.0, 0.0), LongitudeDirection::None);
    assert_eq!(east_or_west(0.0, 1.0), LongitudeDirection::East);
    assert_eq!(east_or_west(0.0, -1.0), LongitudeDirection::West);
    assert_eq!(east_or_west(1.0, -1.0), LongitudeDirection::West);
    assert_eq!(east_or_west(-1.0, 1.0), LongitudeDirection::East);

    assert_eq!(east_or_west(0.5, 1.0), LongitudeDirection::East);
    assert_eq!(east_or_west(-1.0, -0.5), LongitudeDirection::East);

    // wrap around tests
    assert_eq!(east_or_west(3.0, -3.0), LongitudeDirection::East);
    assert_eq!(east_or_west(-3.0, 3.0), LongitudeDirection::West);
}