hittekaart 0.2.0

//! GPX data extraction functions.
//!
//! Like the rest of hittekaart, this module is specialized to our use case. That's why this module
//! internally uses [roxmltree](https://github.com/RazrFalcon/roxmltree) for fast XML parsing,
//! instead of the [gpx](https://github.com/georust/gpx) crate that has more features.
//!
//! This module supports reading files that are `gzip` or `brotli` compressed.
//!
//! Note that we throw away all information that we don't care about. Since we need only the
//! coordinates of a track, we simply use a `Vec<Coordinates>` to represent a track.
use std::{
    f64::consts::PI,
    ffi::OsStr,
    fs::{self, File},
    io::{BufReader, Read},
    path::Path,
};

use flate2::bufread::GzDecoder;
use roxmltree::{Document, Node, NodeType};

use super::error::{Error, Result};

/// World coordinates.
#[derive(Debug, Copy, Clone, PartialEq)]
pub struct Coordinates {
    pub longitude: f64,
    pub latitude: f64,
}

impl Coordinates {
    pub fn new(longitude: f64, latitude: f64) -> Self {
        Self {
            longitude,
            latitude,
        }
    }

    /// Calculates the [Web Mercator
    /// projection](https://en.wikipedia.org/wiki/Web_Mercator_projection) of the coordinates.
    ///
    /// Returns the `(x, y)` projection, where both are in the range `[0, 256 * 2^zoom)`.
    ///
    /// Note that coordinates outside the range are silently truncated to the limits.
    pub fn web_mercator(self, zoom: u32) -> (u64, u64) {
        const WIDTH: f64 = super::layer::TILE_WIDTH as f64;
        const HEIGHT: f64 = super::layer::TILE_HEIGHT as f64;

        // We support longitudes from 180° W/E, which is represented as PI in radians.
        // For latitudes, we actually don't go all the way up to 90° N/S, see
        // https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames
        // To quote:
        //   This value will fall in the range (-π, π) for latitudes between 85.0511 °S and
        //   85.0511 °N.
        //   For the curious, the number 85.0511 is the result of arctan(sinh(π)). By using
        //   this bound, the entire map becomes a (very large) square.
        // We compute the exact bound here.
        let max_lon: f64 = PI;
        let max_lat: f64 = PI.sinh().atan();

        let xmax = 2u64.pow(zoom) as f64 * WIDTH;
        let ymax = 2u64.pow(zoom) as f64 * HEIGHT;

        let lambda = self.longitude.to_radians().clamp(-max_lon, max_lon);
        // The coordinate +180° is the same as -180°, so we want to map them to the same point
        let lambda = if lambda == PI { -PI } else { lambda };
        let phi = self.latitude.to_radians().clamp(-max_lat, max_lat);
        let phi = (phi.tan() + (1.0 / phi.cos())).ln();
        let x = xmax / (2.0 * PI) * (lambda + PI);
        let y = ymax * (1.0 - phi / PI) / 2.0;
        (x.floor() as u64, y.floor() as u64)
    }
}

fn is_track_node(node: &Node) -> bool {
    node.node_type() == NodeType::Element && node.tag_name().name() == "trk"
}

fn is_track_segment(node: &Node) -> bool {
    node.node_type() == NodeType::Element && node.tag_name().name() == "trkseg"
}

fn is_track_point(node: &Node) -> bool {
    node.node_type() == NodeType::Element && node.tag_name().name() == "trkpt"
}

/// Extracts a track from the given string.
pub fn extract_from_str(input: &str) -> Result<Vec<Coordinates>> {
    let mut result = Vec::new();
    let document = Document::parse(input)?;
    for node in document.root_element().children().filter(is_track_node) {
        for segment in node.children().filter(is_track_segment) {
            for point in segment.children().filter(is_track_point) {
                let latitude = point
                    .attribute("lat")
                    .ok_or(Error::MissingLatitude)
                    .and_then(|l| l.parse::<f64>().map_err(Error::InvalidLatitude))?;
                let longitude = point
                    .attribute("lon")
                    .ok_or(Error::MissingLongitude)
                    .and_then(|l| l.parse::<f64>().map_err(Error::InvalidLongitude))?;
                result.push(Coordinates {
                    latitude,
                    longitude,
                });
            }
        }
    }
    Ok(result)
}

/// Compression format of the data.
#[derive(Debug, Copy, Clone, PartialEq, Eq, Hash)]
pub enum Compression {
    /// Indicates that no compression is applied, and the file is plain GPX.
    None,
    /// Indicates that the file is gzip compressed.
    Gzip,
    /// Indicates that the file is brotli compressed.
    Brotli,
}

impl Compression {
    /// Suggests a [`Compression`] from the given path name.
    ///
    /// This will suggest [`Compression::Brotli`] for files ending in `.br`, [`Compression::Gzip`]
    /// for files ending with `.gz` or `.gzip`, and [`Compression::None`] for files ending with
    /// `.gpx`.
    ///
    /// If the file does not end with any of the aforementioned extensions, `None` is returned
    /// instead.
    pub fn suggest_from_path<P: AsRef<Path>>(path: P) -> Option<Compression> {
        let Some(ext) = path.as_ref().extension() else {
            return None;
        };
        if OsStr::new("br") == ext {
            Some(Compression::Brotli)
        } else if [OsStr::new("gz"), OsStr::new("gzip")].contains(&ext) {
            Some(Compression::Gzip)
        } else if OsStr::new("gpx") == ext {
            Some(Compression::None)
        } else {
            None
        }
    }
}

/// Extracts the relevant GPX data from the given file.
///
/// Note that the content must be valid UTF-8, as that is what our parser expects.
pub fn extract_from_file<P: AsRef<Path>>(
    path: P,
    compression: Compression,
) -> Result<Vec<Coordinates>> {
    let content = match compression {
        Compression::None => fs::read_to_string(path).map_err(|e| Error::Io("reading file", e))?,
        Compression::Gzip => {
            let mut result = String::new();
            GzDecoder::new(BufReader::new(
                File::open(path).map_err(|e| Error::Io("opening file", e))?,
            ))
            .read_to_string(&mut result)
            .map_err(|e| Error::Io("reading file", e))?;
            result
        }
        Compression::Brotli => {
            let mut result = Vec::new();
            brotli::BrotliDecompress(
                &mut BufReader::new(File::open(path).map_err(|e| Error::Io("opening file", e))?),
                &mut result,
            )
            .map_err(|e| Error::Io("decompressing", e))?;
            String::from_utf8(result)?
        }
    };
    extract_from_str(&content)
}

#[cfg(test)]
mod test {
    use super::*;
    use rstest::*;

    #[rstest]
    #[case((0.0, 0.0), 0, (128, 128))]
    #[case((-180.0, 0.0), 0, (0, 128))]
    #[case((180.0, 0.0), 0, (0, 128))]
    #[case((179.99, 0.0), 0, (255, 128))]
    #[case((0.0, 90.0), 0, (128, 0))]
    #[case((0.0, -90.0), 0, (128, 255))]
    #[case((0.0, 0.0), 4, (2048, 2048))]
    #[case((-180.0, 0.0), 4, (0, 2048))]
    #[case((180.0, 0.0), 4, (0, 2048))]
    #[case((179.99, 0.0), 4, (4095, 2048))]
    #[case((0.0, 90.0), 4, (2048, 0))]
    #[case((0.0, -90.0), 4, (2048, 4095))]
    #[case((-90.0, 0.0), 0, (64, 128))]
    #[case((90.0, 0.0), 0, (192, 128))]
    // Example taken from
    // https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Example:_Convert_a_GPS_coordinate_to_a_pixel_position_in_a_Web_Mercator_tile
    #[case((139.7006793, 35.6590699), 18, (232798 * 256 + 238, 103246 * 256 + 105))]
    fn web_mercator(#[case] input: (f64, f64), #[case] zoom: u32, #[case] expected: (u64, u64)) {
        assert_eq!(
            Coordinates::new(input.0, input.1).web_mercator(zoom),
            expected
        );
    }
}