use anyhow::Result;
use float_cmp::assert_approx_eq;
use mappers::{
    Ellipsoid,
    projections::{AzimuthalEquidistant, LambertConformalConic, LongitudeLatitude},
};
use mappers_warp::{
    CubicBSpline, MitchellNetravali, RasterBoundsDefinition, Warper, raster_constant_pad,
};
use ndarray::{Array2, Zip};

mod utils;
use utils::*;

#[test]
fn waves() -> Result<()> {
    let src_proj = LongitudeLatitude;
    let tgt_proj = LambertConformalConic::builder()
        .ref_lonlat(80., 24.)
        .standard_parallels(12.472955, 35.1728044444444)
        .ellipsoid(Ellipsoid::WGS84)
        .initialize_projection()?;

    let source_bounds =
        RasterBoundsDefinition::new((60.00, 68.25), (31.75, 40.0), 0.25, 0.25, src_proj)?;
    let target_bounds = RasterBoundsDefinition::new(
        (2_320_000. - 4_000_000., 2_740_000. - 4_000_000.),
        (5_090_000. - 4_000_000., 5_640_000. - 4_000_000.),
        10_000.,
        10_000.,
        tgt_proj,
    )?;

    let warper = Warper::initialize::<CubicBSpline, LongitudeLatitude, LambertConformalConic>(
        &source_bounds,
        &target_bounds,
    )?;

    let source_raster: Array2<f64> = open_nc_data("./tests/data/waves_34.nc")?;
    let ref_raster: Array2<f64> = open_nc_data("./tests/data/waves_ref.nc")?;
    let target_raster = warper.warp_ignore_nodata(&source_raster)?;

    assert_eq!(target_raster.shape(), ref_raster.shape());
    Zip::from(&target_raster)
        .and(&ref_raster)
        .map_collect(|&f, &o| assert_approx_eq!(f64, f, o, epsilon = 1e-6));

    Ok(())
}

#[test]
fn gfs_t2m() -> Result<()> {
    let src_proj = LongitudeLatitude;
    let eu_proj = LambertConformalConic::builder()
        .ref_lonlat(10.0, 52.0)
        .standard_parallels(35.0, 65.0)
        .ellipsoid(Ellipsoid::WGS84)
        .initialize_projection()?;

    let source_domain =
        RasterBoundsDefinition::new((-70.0, 85.0), (17.0, 77.0), 0.25, 0.25, src_proj)?;
    let target_domain = RasterBoundsDefinition::new(
        (-4_120_000., 3_490_000.),
        (-2_750_000., 2_640_000.),
        10_000.,
        10_000.,
        eu_proj,
    )?;

    let warper = Warper::initialize::<CubicBSpline, LongitudeLatitude, LambertConformalConic>(
        &source_domain,
        &target_domain,
    )?;
    let source_raster: Array2<f64> = open_nc_data("./tests/data/gfs_t2m.nc")?;
    let target_raster = warper.warp_ignore_nodata(&source_raster)?;

    target_raster.iter().for_each(|&v| assert!(v.is_finite()));

    assert!(target_raster.max()? <= source_raster.max()?);
    assert!(target_raster.min()? >= source_raster.min()?);

    Ok(())
}

#[test]
fn mitchell() -> Result<()> {
    let src_proj = LongitudeLatitude;
    let eu_proj = LambertConformalConic::builder()
        .ref_lonlat(10.0, 52.0)
        .standard_parallels(35.0, 65.0)
        .ellipsoid(Ellipsoid::WGS84)
        .initialize_projection()?;

    let source_domain =
        RasterBoundsDefinition::new((-70.0, 85.0), (17.0, 77.0), 0.25, 0.25, src_proj)?;
    let target_domain = RasterBoundsDefinition::new(
        (-4_120_000., 3_490_000.),
        (-2_750_000., 2_640_000.),
        10_000.,
        10_000.,
        eu_proj,
    )?;

    let warper = Warper::initialize::<MitchellNetravali, LongitudeLatitude, LambertConformalConic>(
        &source_domain,
        &target_domain,
    )?;
    let source_raster: Array2<f64> = open_nc_data("./tests/data/gfs_t2m.nc")?;
    let target_raster = warper.warp_ignore_nodata(&source_raster)?;

    target_raster.iter().for_each(|&v| assert!(v.is_finite()));

    assert!(target_raster.max()? <= source_raster.max()?);
    assert!(target_raster.min()? >= source_raster.min()?);

    Ok(())
}

#[test]
fn nan_padded_waves() -> Result<()> {
    let src_proj = LongitudeLatitude;
    let tgt_proj = LambertConformalConic::builder()
        .ref_lonlat(80., 24.)
        .standard_parallels(12.472955, 35.1728044444444)
        .ellipsoid(Ellipsoid::WGS84)
        .initialize_projection()?;

    let source_bounds =
        RasterBoundsDefinition::new((59.25, 69.00), (31.00, 40.75), 0.25, 0.25, src_proj)?;
    let target_bounds = RasterBoundsDefinition::new(
        (2_320_000. - 4_000_000., 2_740_000. - 4_000_000.),
        (5_090_000. - 4_000_000., 5_640_000. - 4_000_000.),
        10_000.,
        10_000.,
        tgt_proj,
    )?;

    let warper = Warper::initialize::<CubicBSpline, LongitudeLatitude, LambertConformalConic>(
        &source_bounds,
        &target_bounds,
    )?;

    let source_raster: Array2<f64> = open_nc_data("./tests/data/waves_34.nc")?;
    let source_raster = raster_constant_pad(&source_raster, 3, f64::NAN);
    let ref_raster: Array2<f64> = open_nc_data("./tests/data/waves_ref.nc")?;
    let target_raster = warper.warp_ignore_nodata(&source_raster)?;

    assert_eq!(target_raster.shape(), ref_raster.shape());
    Zip::from(&target_raster)
        .and(&ref_raster)
        .map_collect(|&f, &o| assert_approx_eq!(f64, f, o, epsilon = 1e-6));

    Ok(())
}

#[test]
fn aeqd_to_lcc() -> Result<()> {
    let src_proj = AzimuthalEquidistant::builder()
        .ref_lonlat(19.0926, 52.3469)
        .ellipsoid(Ellipsoid::SPHERE)
        .initialize_projection()?;
    let tgt_proj = LambertConformalConic::builder()
        .ref_lonlat(19.0926, 52.3469)
        .standard_parallels(52.344876525433854, 52.34892347456614)
        .ellipsoid(Ellipsoid::WGS84)
        .initialize_projection()?;

    let source_domain = RasterBoundsDefinition::new(
        (-452500., 452500.),
        (-452500., 452500.),
        1000.,
        1000.,
        src_proj,
    )?;
    let target_domain = RasterBoundsDefinition::new(
        (-449500., 449500.),
        (-449500., 449500.),
        1000.,
        1000.,
        tgt_proj,
    )?;

    let warper = Warper::initialize::<CubicBSpline, AzimuthalEquidistant, LambertConformalConic>(
        &source_domain,
        &target_domain,
    )?;

    let source_raster: Array2<f64> = open_nc_data("./tests/data/aeqd_nan.nc")?;
    let source_raster = raster_constant_pad(&source_raster, 3, f64::NAN);

    let target_raster = warper.warp_ignore_nodata(&source_raster)?;
    let ref_raster: Array2<f64> = open_nc_data("./tests/data/aeqd_ref.nc")?;

    assert_eq!(target_raster.shape(), ref_raster.shape());
    Zip::from(&target_raster)
        .and(&ref_raster)
        .map_collect(|&f, &o| assert_approx_eq!(f64, f, o, epsilon = 1e-6));

    let source_values = source_raster
        .into_iter()
        .filter(|&v| !v.is_nan())
        .collect::<Vec<f64>>();

    let target_values = target_raster
        .into_iter()
        .filter(|&v| !v.is_nan())
        .collect::<Vec<f64>>();

    target_values.iter().for_each(|v| assert!(v.is_finite()));

    let source_max = *source_values.iter().max_by(|a, b| a.total_cmp(b)).unwrap();
    let source_min = *source_values.iter().min_by(|a, b| a.total_cmp(b)).unwrap();
    let target_max = *target_values.iter().max_by(|a, b| a.total_cmp(b)).unwrap();
    let target_min = *target_values.iter().min_by(|a, b| a.total_cmp(b)).unwrap();

    assert!(target_max <= source_max);
    assert_approx_eq!(f64, source_min, target_min, epsilon = 1e-12);

    Ok(())
}

#[test]
fn reverse_waves() -> Result<()> {
    let src_proj = LongitudeLatitude;
    let tgt_proj = LambertConformalConic::builder()
        .ref_lonlat(80., 24.)
        .standard_parallels(12.472955, 35.1728044444444)
        .ellipsoid(Ellipsoid::WGS84)
        .initialize_projection()?;

    let source_bounds =
        RasterBoundsDefinition::new((60.00, 68.25), (31.75, 40.0), 0.25, 0.25, src_proj)?;
    let target_bounds = RasterBoundsDefinition::new(
        (2_320_000. - 4_000_000., 2_740_000. - 4_000_000.),
        (5_090_000. - 4_000_000., 5_640_000. - 4_000_000.),
        10_000.,
        10_000.,
        tgt_proj,
    )?;

    let warper = Warper::initialize::<CubicBSpline, LongitudeLatitude, LambertConformalConic>(
        &source_bounds,
        &target_bounds,
    )?;

    let source_raster: Array2<f64> = open_nc_data("./tests/data/waves_34.nc")?;
    let ref_raster: Array2<f64> = open_nc_data("./tests/data/waves_ref.nc")?;
    let target_raster = warper.warp_ignore_nodata(&source_raster)?;

    assert_eq!(target_raster.shape(), ref_raster.shape());
    Zip::from(&target_raster)
        .and(&ref_raster)
        .map_collect(|&f, &o| assert_approx_eq!(f64, f, o, epsilon = 1e-6));

    let source_bounds_reverse =
        RasterBoundsDefinition::new((63.0, 64.0), (34.0, 36.0), 0.25, 0.25, src_proj)?;

    let reverse_warper =
        Warper::initialize::<CubicBSpline, _, _>(&target_bounds, &source_bounds_reverse)?;

    #[cfg(feature = "multithreading")]
    {
        let parallel_warper = Warper::initialize_parallel::<CubicBSpline, _, _>(
            &target_bounds,
            &source_bounds_reverse,
        )?;
        assert_eq!(reverse_warper, parallel_warper);
    }

    let _ = reverse_warper.warp_reject_nodata(&target_raster)?;

    // for now I don't know a good method for validating the reverse raster so, we just check for errors
    Ok(())
}
