oxigdal-node 0.1.7

//! Algorithm bindings for Node.js
//!
//! This module provides geospatial algorithms including resampling,
//! calculator, terrain analysis, and statistical operations.

use napi::bindgen_prelude::*;
use napi_derive::napi;
use oxigdal_algorithms::raster::compute_zonal_stats as compute_zonal;
use oxigdal_algorithms::raster::{HillshadeParams, hillshade as compute_hillshade};
use oxigdal_algorithms::raster::{
    SlopeAspectConfig, SlopeUnits, aspect_advanced as compute_aspect_advanced,
    slope_advanced as compute_slope_advanced,
};
use oxigdal_algorithms::resampling::{Resampler, ResamplingMethod as CoreResamplingMethod};
use oxigdal_algorithms::vector::{
    AreaMethod as CoreAreaMethod, BufferCapStyle, BufferJoinStyle,
    BufferOptions as CoreBufferOptions, SimplifyMethod as CoreSimplifyMethod, area_polygon,
    buffer_point, buffer_polygon, simplify_linestring,
};
use oxigdal_core::buffer::RasterBuffer;

use crate::buffer::BufferWrapper;
use crate::error::{NodeError, ToNapiResult};
use crate::vector::GeometryWrapper;

/// Resampling methods
#[napi]
pub enum ResamplingMethod {
    /// Nearest neighbor (fast, preserves exact values)
    NearestNeighbor,
    /// Bilinear interpolation (smooth, good for continuous data)
    Bilinear,
    /// Bicubic interpolation (high quality, slower)
    Bicubic,
    /// Lanczos resampling (highest quality, expensive)
    Lanczos,
}

impl From<ResamplingMethod> for CoreResamplingMethod {
    fn from(method: ResamplingMethod) -> Self {
        match method {
            ResamplingMethod::NearestNeighbor => CoreResamplingMethod::Nearest,
            ResamplingMethod::Bilinear => CoreResamplingMethod::Bilinear,
            ResamplingMethod::Bicubic => CoreResamplingMethod::Bicubic,
            ResamplingMethod::Lanczos => CoreResamplingMethod::Lanczos,
        }
    }
}

/// Resamples a raster buffer to a new size
#[napi]
pub fn resample(
    buffer: &BufferWrapper,
    new_width: u32,
    new_height: u32,
    method: ResamplingMethod,
) -> Result<BufferWrapper> {
    let resampler = Resampler::new(method.into());
    let resampled = resampler
        .resample(buffer.inner(), new_width as u64, new_height as u64)
        .to_napi()?;

    Ok(BufferWrapper::from_raster_buffer(resampled))
}

/// Raster calculator - evaluates an expression on raster bands
#[allow(dead_code)]
#[napi]
pub fn calculate(expression: String, bands: Vec<&BufferWrapper>) -> Result<BufferWrapper> {
    if bands.is_empty() {
        return Err(NodeError {
            code: "INVALID_INPUT".to_string(),
            message: "At least one band is required".to_string(),
        }
        .into());
    }

    // Get dimensions from first band
    let width = bands[0].width() as u64;
    let height = bands[0].height() as u64;

    // Verify all bands have same dimensions
    for (i, band) in bands.iter().enumerate() {
        if band.width() as u64 != width || band.height() as u64 != height {
            return Err(NodeError {
                code: "DIMENSION_MISMATCH".to_string(),
                message: format!(
                    "Band {} has different dimensions ({}x{}) than first band ({}x{})",
                    i,
                    band.width(),
                    band.height(),
                    width,
                    height
                ),
            }
            .into());
        }
    }

    // Simple expression evaluation (supports basic operations)
    let result = evaluate_expression(&expression, bands)?;
    Ok(BufferWrapper::from_raster_buffer(result))
}

/// Buffer operations - creates a buffer around geometries
#[napi]
pub fn buffer(geometry: &GeometryWrapper, distance: f64, segments: u32) -> Result<GeometryWrapper> {
    use oxigdal_core::vector::Geometry;

    let options = CoreBufferOptions {
        quadrant_segments: segments as usize,
        cap_style: BufferCapStyle::Round,
        join_style: BufferJoinStyle::Round,
        miter_limit: 5.0,
        simplify_tolerance: 0.0,
    };

    let buffered = match geometry.inner() {
        Geometry::Point(p) => {
            let polygon = buffer_point(p, distance, &options).to_napi()?;
            Geometry::Polygon(polygon)
        }
        Geometry::Polygon(p) => {
            let polygon = buffer_polygon(p, distance, &options).to_napi()?;
            Geometry::Polygon(polygon)
        }
        _ => {
            return Err(NodeError {
                code: "NOT_IMPLEMENTED".to_string(),
                message: "Buffer not implemented for this geometry type".to_string(),
            }
            .into());
        }
    };

    Ok(GeometryWrapper::from_geometry(buffered))
}

/// Calculates the area of a polygon
#[napi]
pub fn area(geometry: &GeometryWrapper, method: String) -> Result<f64> {
    use oxigdal_core::vector::Geometry;

    let area_method = match method.to_lowercase().as_str() {
        "planar" => CoreAreaMethod::Planar,
        "geodetic" => CoreAreaMethod::Geodetic,
        _ => {
            return Err(NodeError {
                code: "INVALID_PARAMETER".to_string(),
                message: format!("Unknown area method: {}", method),
            }
            .into());
        }
    };

    match geometry.inner() {
        Geometry::Polygon(p) => area_polygon(p, area_method).to_napi(),
        _ => Err(NodeError {
            code: "INVALID_GEOMETRY".to_string(),
            message: "Area calculation requires a Polygon geometry".to_string(),
        }
        .into()),
    }
}

/// Simplifies a geometry using the specified method
#[napi]
pub fn simplify(
    geometry: &GeometryWrapper,
    tolerance: f64,
    method: String,
) -> Result<GeometryWrapper> {
    use oxigdal_core::vector::Geometry;

    let simplify_method = match method.to_lowercase().as_str() {
        "douglas-peucker" | "dp" => CoreSimplifyMethod::DouglasPeucker,
        "visvalingam-whyatt" | "vw" => CoreSimplifyMethod::VisvalingamWhyatt,
        _ => {
            return Err(NodeError {
                code: "INVALID_PARAMETER".to_string(),
                message: format!("Unknown simplify method: {}", method),
            }
            .into());
        }
    };

    let simplified = match geometry.inner() {
        Geometry::LineString(ls) => {
            let simple_ls = simplify_linestring(ls, tolerance, simplify_method).to_napi()?;
            Geometry::LineString(simple_ls)
        }
        _ => {
            return Err(NodeError {
                code: "INVALID_GEOMETRY".to_string(),
                message: "Simplify currently only supports LineString geometry".to_string(),
            }
            .into());
        }
    };

    Ok(GeometryWrapper::from_geometry(simplified))
}

/// Validates that a ground-unit pixel size is usable for terrain algorithms.
///
/// A non-positive or non-finite pixel size would silently produce `NaN`/`Inf`
/// slope, aspect, or hillshade values (since it is used as a division/scale
/// factor), so it is rejected explicitly instead of propagating garbage data.
fn validate_pixel_size(pixel_size: f64) -> Result<()> {
    if pixel_size.is_finite() && pixel_size > 0.0 {
        Ok(())
    } else {
        Err(NodeError {
            code: "INVALID_PARAMETER".to_string(),
            message: format!(
                "pixel_size must be a positive, finite number (in the DEM's ground units), got {}",
                pixel_size
            ),
        }
        .into())
    }
}

/// Computes hillshade from a DEM
///
/// `pixel_size` is the ground distance covered by one pixel (e.g. meters or
/// degrees, matching the DEM's coordinate reference system) and must be
/// supplied by the caller since raw pixel buffers carry no georeferencing.
#[napi]
pub fn hillshade(
    dem: &BufferWrapper,
    azimuth: f64,
    altitude: f64,
    z_factor: f64,
    pixel_size: f64,
) -> Result<BufferWrapper> {
    validate_pixel_size(pixel_size)?;

    let params = HillshadeParams {
        azimuth,
        altitude,
        z_factor,
        pixel_size,
        scale: 255.0,
    };

    let result = compute_hillshade(dem.inner(), params).to_napi()?;
    Ok(BufferWrapper::from_raster_buffer(result))
}

/// Computes slope from a DEM
///
/// `pixel_size` is the ground distance covered by one pixel (e.g. meters or
/// degrees, matching the DEM's coordinate reference system). When
/// `as_percent` is `true`, the result is percent rise instead of degrees.
#[napi]
pub fn slope(
    dem: &BufferWrapper,
    pixel_size: f64,
    z_factor: f64,
    as_percent: bool,
) -> Result<BufferWrapper> {
    validate_pixel_size(pixel_size)?;

    let config = SlopeAspectConfig {
        slope_units: if as_percent {
            SlopeUnits::Percent
        } else {
            SlopeUnits::Degrees
        },
        z_factor,
        ..Default::default()
    };

    let result = compute_slope_advanced(dem.inner(), pixel_size, &config).to_napi()?;
    Ok(BufferWrapper::from_raster_buffer(result))
}

/// Computes aspect from a DEM
///
/// `pixel_size` is the ground distance covered by one pixel (e.g. meters or
/// degrees, matching the DEM's coordinate reference system). Aspect (the
/// compass direction a slope faces) does not depend on the vertical
/// exaggeration factor since it is a ratio of the x/y gradients.
#[napi]
pub fn aspect(dem: &BufferWrapper, pixel_size: f64) -> Result<BufferWrapper> {
    validate_pixel_size(pixel_size)?;

    let result = compute_aspect_advanced(dem.inner(), pixel_size, &SlopeAspectConfig::default())
        .to_napi()?;
    Ok(BufferWrapper::from_raster_buffer(result))
}

/// Computes zonal statistics
#[napi]
pub fn zonal_stats(raster: &BufferWrapper, zones: &BufferWrapper) -> Result<Vec<ZonalStatistics>> {
    let core_stats = compute_zonal(raster.inner(), zones.inner()).to_napi()?;

    let results = core_stats
        .into_iter()
        .map(|s| ZonalStatistics {
            zone_id: s.zone_id,
            count: s.count as u32,
            min: s.min,
            max: s.max,
            mean: s.mean,
            stddev: s.std_dev,
            sum: s.sum,
        })
        .collect();

    Ok(results)
}

/// Zonal statistics result
#[napi(object)]
pub struct ZonalStatistics {
    pub zone_id: i32,
    pub count: u32,
    pub min: f64,
    pub max: f64,
    pub mean: f64,
    pub stddev: f64,
    pub sum: f64,
}

// Helper functions

impl GeometryWrapper {
    #[allow(dead_code)]
    pub(crate) fn from_geometry(geom: oxigdal_core::vector::Geometry) -> Self {
        Self { inner: geom }
    }
}

/// Simple expression evaluator for raster calculator
#[allow(dead_code)]
fn evaluate_expression(expr: &str, bands: Vec<&BufferWrapper>) -> Result<RasterBuffer> {
    // For now, use the RasterCalculator from oxigdal-algorithms
    // Simple single-band pass-through
    if expr.len() == 1 && expr.chars().next().is_some_and(|c| c.is_ascii_uppercase()) {
        let band_idx = (expr.as_bytes()[0] - b'A') as usize;
        if band_idx >= bands.len() {
            return Err(NodeError {
                code: "INVALID_EXPRESSION".to_string(),
                message: format!("Band {} not found", expr),
            }
            .into());
        }
        Ok(bands[band_idx].inner().clone())
    } else {
        Err(NodeError {
            code: "NOT_IMPLEMENTED".to_string(),
            message: "Complex expressions not yet supported".to_string(),
        }
        .into())
    }
}

#[cfg(test)]
#[allow(clippy::unwrap_used, clippy::expect_used, clippy::panic)]
mod tests {
    use super::*;
    use oxigdal_core::types::RasterDataType;

    /// Builds a 5x5 float32 DEM ramping linearly along x: elevation = x * 10.0.
    fn ramp_dem() -> BufferWrapper {
        let mut buf = RasterBuffer::zeros(5, 5, RasterDataType::Float32);
        for y in 0..5 {
            for x in 0..5 {
                buf.set_pixel(x, y, (x as f64) * 10.0)
                    .expect("set_pixel should succeed for in-bounds coordinates");
            }
        }
        BufferWrapper::from_raster_buffer(buf)
    }

    /// Asserts a `Result<BufferWrapper>` is an error and returns its message.
    ///
    /// `BufferWrapper` does not implement `Debug` (it wraps a large pixel
    /// buffer), so `Result::expect_err`/`unwrap_err` cannot be used directly.
    fn expect_error_message(result: Result<BufferWrapper>) -> String {
        match result {
            Ok(_) => panic!("expected an error, got Ok(..)"),
            Err(e) => e.to_string(),
        }
    }

    #[test]
    fn slope_rejects_non_positive_pixel_size() {
        let dem = ramp_dem();

        let message = expect_error_message(slope(&dem, 0.0, 1.0, false));
        assert!(message.contains("INVALID_PARAMETER"));

        let message = expect_error_message(slope(&dem, -5.0, 1.0, false));
        assert!(message.contains("INVALID_PARAMETER"));

        let message = expect_error_message(slope(&dem, f64::NAN, 1.0, false));
        assert!(message.contains("INVALID_PARAMETER"));
    }

    #[test]
    fn hillshade_rejects_non_positive_pixel_size() {
        let dem = ramp_dem();
        let message = expect_error_message(hillshade(&dem, 315.0, 45.0, 1.0, 0.0));
        assert!(message.contains("INVALID_PARAMETER"));
    }

    #[test]
    fn aspect_rejects_non_positive_pixel_size() {
        let dem = ramp_dem();
        let message = expect_error_message(aspect(&dem, f64::INFINITY));
        assert!(message.contains("INVALID_PARAMETER"));
    }

    #[test]
    fn slope_scales_inversely_with_pixel_size() {
        let dem = ramp_dem();

        let slope_1m = slope(&dem, 1.0, 1.0, false).expect("slope at 1.0 pixel_size");
        let slope_10m = slope(&dem, 10.0, 1.0, false).expect("slope at 10.0 pixel_size");

        // At an interior pixel, the same elevation ramp spread over a larger
        // pixel_size means a gentler rise/run gradient, hence a smaller slope
        // angle. This is exactly the bug: previously pixel_size was hardcoded
        // to 1.0 regardless of what was passed in, so both calls would have
        // produced identical output.
        let center_1m = slope_1m.get_pixel(2, 2).expect("center pixel readable");
        let center_10m = slope_10m.get_pixel(2, 2).expect("center pixel readable");

        assert!(
            center_1m > center_10m,
            "slope with pixel_size=1.0 ({center_1m}) should exceed slope with pixel_size=10.0 ({center_10m})"
        );
    }

    #[test]
    fn slope_as_percent_flag_changes_units() {
        let dem = ramp_dem();

        let degrees = slope(&dem, 1.0, 1.0, false).expect("slope in degrees");
        let percent = slope(&dem, 1.0, 1.0, true).expect("slope in percent");

        let degrees_value = degrees.get_pixel(2, 2).expect("center pixel readable");
        let percent_value = percent.get_pixel(2, 2).expect("center pixel readable");

        // Previously `as_percent` was discarded entirely, so both calls would
        // yield identical (degree-based) output regardless of the flag.
        assert!(
            (degrees_value - percent_value).abs() > 1e-6,
            "as_percent=true must change the output units: degrees={degrees_value}, percent={percent_value}"
        );

        // tan(degrees) * 100 == percent, cross-checking against the known
        // conversion formula rather than just asserting inequality. The
        // buffers are float32, so allow for that reduced precision rather
        // than requiring f64-level exactness.
        let expected_percent = degrees_value.to_radians().tan() * 100.0;
        assert!(
            (expected_percent - percent_value).abs() < 1e-2,
            "percent value {percent_value} should equal tan(degrees) * 100 = {expected_percent}"
        );
    }

    #[test]
    fn hillshade_pixel_size_affects_output() {
        let dem = ramp_dem();

        let shade_1m = hillshade(&dem, 315.0, 45.0, 1.0, 1.0).expect("hillshade at 1.0");
        let shade_10m = hillshade(&dem, 315.0, 45.0, 1.0, 10.0).expect("hillshade at 10.0");

        let value_1m = shade_1m.get_pixel(2, 2).expect("center pixel readable");
        let value_10m = shade_10m.get_pixel(2, 2).expect("center pixel readable");

        assert!(
            (value_1m - value_10m).abs() > 1e-6,
            "hillshade should differ when pixel_size changes: {value_1m} vs {value_10m}"
        );
    }

    #[test]
    fn aspect_pixel_size_accepted_for_valid_input() {
        let dem = ramp_dem();
        let result = aspect(&dem, 30.0);
        assert!(result.is_ok(), "aspect should accept a valid pixel_size");
    }
}