eorst 1.0.1

#[cfg(test)]
mod test {
    use crate::array_ops::{trimm_array4_owned, fill_nodata_simple};
    use crate::types::{Dimension, RasterDataShape, BlockSize, ImageSize, SamplingMethod};
 
    use crate::data_sources::{DataSource, DataSourceBuilder};
    use crate::metadata::Layer;
    use crate::selection::Stack;
    use crate::rasterdataset::builder::get_overlap;
    use crate::*;
    use std::path::PathBuf;
    use std::collections::BTreeMap;
    use ndarray::Array4;

    macro_rules! td {
        ($rel:expr_2021) => {
            std::path::PathBuf::from(env!("CARGO_MANIFEST_DIR"))
                .join("data")
                .join($rel)
        };
    }

    use num_traits::FromPrimitive;

    use std::ops::{Add, Div};
    #[test]
    fn test_fill_nodata_simple() {
        // Create a 2-band 4x4 array with some nodata values
        let nodata = -999.0;
        let bands = 3;
        let size = 256;

        // Initialize array with sequential values
        let mut data = Array3::<f32>::zeros((bands, size, size));
        for b in 0..bands {
            for r in 0..size {
                for c in 0..size {
                    data[[b, r, c]] = (r * size + c) as f32 + b as f32 * 10000.0;
                }
            }
        }

        // Insert holes in the same place across all bands
        // 1px hole
        data[[0, 50, 50]] = nodata;
        data[[1, 50, 50]] = nodata;
        data[[2, 50, 50]] = nodata;

        // 10px hole (10x1 rectangle)
        for i in 0..10 {
            data[[0, 100 + i, 100]] = nodata;
            data[[1, 100 + i, 100]] = nodata;
            data[[2, 100 + i, 100]] = nodata;
        }

        // 20px hole (20x1 rectangle)
        for i in 0..20 {
            data[[0, 150 + i, 150]] = nodata;
            data[[1, 150 + i, 150]] = nodata;
            data[[2, 150 + i, 150]] = nodata;
        }

        // Run the fill function
        fill_nodata_simple(&mut data, nodata);

        println!("data {data:?}");

        // Call the fill function
        fill_nodata_simple(&mut data, nodata);
        println!("data {data:?}");

        // Assert that no nodata values remain
        for val in data.iter() {
            assert_ne!(*val, nodata, "Found remaining nodata value!");
        }

        // Optionally: check some expected filled values
        // (just as an example for first band, position [0, 2])
        let filled_val = data[[0, 0, 2]];
        assert!(
            filled_val > 0.0 && filled_val < 10.0,
            "Filled value seems wrong"
        );
    }

    #[test]
    fn test_align() {
        let data_source = DataSourceBuilder::from_file(&PathBuf::from(
            "data/cemsre_t55jfm_20200614_sub_abam5.tif",
        ))
        .build();
        let rds_1 = RasterDatasetBuilder::<i16>::from_source(&data_source).build();

        let data_source = DataSourceBuilder::from_file(&PathBuf::from(
            "data/cemsre_t55jfm_20200614_suba_abam5.tif",
        ))
        .build();
        let rds_2 = RasterDatasetBuilder::<i16>::from_source(&data_source)
            .set_template(&rds_1)
            .build();
        assert_eq!(rds_1.metadata.shape, rds_2.metadata.shape);
        assert_eq!(rds_1.metadata.geo_transform, rds_2.metadata.geo_transform);

        let data_source = DataSourceBuilder::from_file(&PathBuf::from(
            "data/cemsre_t55jfm_20200614_suba_abam6.tif",
        ))
        .build();
        let rds_2 = RasterDatasetBuilder::<i16>::from_source(&data_source)
            .set_template(&rds_1)
            .build();
        println!(
            "rds_1 {:?} rds_2{:?}",
            rds_1.metadata.geo_transform, rds_2.metadata.geo_transform
        );
        assert_eq!(rds_1.metadata.shape, rds_2.metadata.shape);
        assert_eq!(rds_1.metadata.geo_transform, rds_2.metadata.geo_transform);
    }

    #[test]
    fn test_iter() {
        let data_source = DataSourceBuilder::from_file(&PathBuf::from(td!(
            "cemsre_t55jfm_20200614_sub_abam5.tif"
        )))
        .build();
        let rds = RasterDatasetBuilder::<i16>::from_source(&data_source).build();
        for iter in rds.iter() {
            let _ = iter.parent.read_block::<i32>(iter.iter_index);
        }
    }
    #[test]
    fn test_stacked_dims() {
        let mut shp_1 = RasterDataShape {
            times: 1,
            layers: 1,
            rows: 1024,
            cols: 1024,
        };
        let shp_2 = RasterDataShape {
            times: 1,
            layers: 1,
            rows: 1024,
            cols: 1024,
        };
        let e = RasterDataShape {
            times: 1,
            layers: 2,
            rows: 1024,
            cols: 1024,
        };
        shp_1.stack(shp_2, Dimension::Layer);
        assert_eq!(shp_1, e);
        let mut shp_1 = RasterDataShape {
            times: 1,
            layers: 1,
            rows: 1024,
            cols: 1024,
        };
        let shp_2 = RasterDataShape {
            times: 1,
            layers: 2,
            rows: 1024,
            cols: 1024,
        };
        let e = RasterDataShape {
            times: 1,
            layers: 3,
            rows: 1024,
            cols: 1024,
        };
        shp_1.stack(shp_2, Dimension::Layer);
        assert_eq!(shp_1, e);
        let mut shp_1 = RasterDataShape {
            times: 1,
            layers: 1,
            rows: 1024,
            cols: 1024,
        };
        let shp_2 = RasterDataShape {
            times: 1,
            layers: 1,
            rows: 1024,
            cols: 1024,
        };
        let e = RasterDataShape {
            times: 2,
            layers: 1,
            rows: 1024,
            cols: 1024,
        };
        shp_1.stack(shp_2, Dimension::Time);
        assert_eq!(shp_1, e);
    }

    #[test]
    fn test_shape_extend() {
        let mut shp_1 = RasterDataShape {
            times: 1,
            layers: 1,
            rows: 1024,
            cols: 1024,
        };
        let shp_2 = RasterDataShape {
            times: 1,
            layers: 1,
            rows: 1024,
            cols: 1024,
        };
        let e = RasterDataShape {
            times: 2,
            layers: 1,
            rows: 1024,
            cols: 1024,
        };
        shp_1.extend(shp_2);
        assert_eq!(shp_1, e);

        let mut shp_1 = RasterDataShape {
            times: 1,
            layers: 2,
            rows: 1024,
            cols: 1024,
        };
        let shp_2 = RasterDataShape {
            times: 1,
            layers: 2,
            rows: 1024,
            cols: 1024,
        };
        let e = RasterDataShape {
            times: 2,
            layers: 2,
            rows: 1024,
            cols: 1024,
        };
        shp_1.extend(shp_2);
        assert_eq!(shp_1, e);
    }
    #[test]
    #[should_panic(expected = "Unable to extend layers")]
    fn test_shape_extend_diff_layer() {
        let mut shp_1 = RasterDataShape {
            times: 1,
            layers: 1,
            rows: 1024,
            cols: 1024,
        };
        let shp_2 = RasterDataShape {
            times: 1,
            layers: 2,
            rows: 1024,
            cols: 1024,
        };
        let _e = RasterDataShape {
            times: 2,
            layers: 2,
            rows: 1024,
            cols: 1024,
        };
        shp_1.extend(shp_2);
    }

    #[test]
    #[should_panic(expected = "Unable to stack layers")]
    fn test_shape_stack_diff_time() {
        let mut shp_1 = RasterDataShape {
            times: 1,
            layers: 1,
            rows: 1024,
            cols: 1024,
        };
        let shp_2 = RasterDataShape {
            times: 2,
            layers: 1,
            rows: 1024,
            cols: 1024,
        };
        shp_1.stack(shp_2, Dimension::Layer);
    }

    #[test]
    #[should_panic(expected = "Unable to stack layers")]
    fn test_shape_stack_diff_time2() {
        let mut shp_1 = RasterDataShape {
            times: 1,
            layers: 1,
            rows: 1024,
            cols: 1024,
        };
        let shp_2 = RasterDataShape {
            times: 2,
            layers: 2,
            rows: 1024,
            cols: 1024,
        };
        shp_1.stack(shp_2, Dimension::Layer);
    }

    #[test]
    fn test_layer_stack_position() {
        let layer_1: Layer = Layer {
            source: PathBuf::from("image_1.img"),
            layer_pos: 0,
            time_pos: 0,
        };
        let mut layer_2: Layer = Layer {
            source: PathBuf::from("image_1.img"),
            layer_pos: 0,
            time_pos: 0,
        };
        let max_time_rds_1 = 3;
        layer_2.stack_position(layer_1, Dimension::Layer, max_time_rds_1);
        assert_eq!(layer_2.layer_pos, 3);
    }

    #[test]
    fn test_reduce() {
        // TODO: Audit that the environment access only happens in single-threaded code.
        unsafe { std::env::set_var("RAYON_NUM_THREADS", "1") };

        fn worker<T>(rdb: &RasterDataBlock<T>, dimension: Dimension) -> Array3<T>
        where
            T: RasterType + FromPrimitive + Add<Output = T> + Div<Output = T>,
        {
            use crate::selection::SumDimension;
            rdb.data.sum_dimension(dimension)
        }
        // Create a ds with 2 time steps
        let source = DataSource {
            source: td!("cemsre_t55jfm_20200614_sub_abam5.tif"),
            bands: [1, 2].to_vec(),
            layer_names: Vec::new(),
        };
        let source_2 = DataSource {
            source: td!("cemsre_t55jfm_20200614_sub_abam5.tif"),
            bands: [1, 2].to_vec(),
            layer_names: Vec::new(),
        };

        let sources = vec![source, source_2];

        let rds = RasterDatasetBuilder::<i16>::from_sources(&sources)
            .block_size(BlockSize {
                rows: 1024,
                cols: 1024,
            })
            .source_composition_dimension(Dimension::Time)
            .band_composition_dimension(Dimension::Layer)
            .build();

        let out_file = PathBuf::from("/tmp/tmp.tif");
        let dimension = Dimension::Layer;
        rds.apply_reduction::<i16>(worker, dimension, 1, &out_file, 32767);
    }

    #[test]
    fn test_overlap() {
        let source = DataSource {
            source: PathBuf::from("data/cemsre_t55jfm_20200614_sub_abam5.tif"),
            bands: [1].to_vec(),
            layer_names: Vec::new(),
        };
        let rds_over_0: RasterDataset<i16> = RasterDatasetBuilder::<i16>::from_source(&source)
            .overlap_size(0)
            .build();

        let source = DataSource {
            source: PathBuf::from("data/cemsre_t55jfm_20200614_sub_abam5.tif"),
            bands: [1].to_vec(),
            layer_names: Vec::new(),
        };
        let rds_over_1: RasterDataset<i16> = RasterDatasetBuilder::<i16>::from_source(&source)
            .overlap_size(1)
            .build();
        println!("rds_over: {rds_over_1}");

        // The idea is that reading with no overlap should give the same result
        // as reading with overlap + trimming the excess
        let n_blocks = 4;
        for block_id in 0..n_blocks {
            //println!("block_id {block_id:?}");
            let array_with_overlap = rds_over_1.read_block::<i32>(block_id);
            //println!("read 1 {array_with_overlap:?}");
            let overlap = rds_over_1.blocks[block_id].overlap;
            //println!("array_with overlap {array_with_overlap:?}");

            let array_with_overlap_trimmed = trimm_array4_owned(&array_with_overlap, &overlap);
            //println!("array_trimmed {array_with_overlap_trimmed:?}");

            // should be the same than with no overlaps!
            let array_with_no_overlap = rds_over_0.read_block::<i32>(block_id);
            //println!("read 2 {array_with_no_overlap:?}");
            assert_eq!(array_with_no_overlap, array_with_overlap_trimmed);
        }
    }
    #[test]
    fn test_get_overlap() {
        let image_size = ImageSize {
            rows: 1580,
            cols: 1514,
        };
        let block_size = BlockSize {
            rows: 1024,
            cols: 1024,
        };
        let overlap_size = 1;
        let block_col_rows = vec![(0, 0), (1, 0), (0, 1), (1, 1)];
        for block_col_row in block_col_rows {
            get_overlap(image_size, block_size, block_col_row, overlap_size);
        }
    }

    #[test]
    fn test_extract_blockwise() {
        let mut e: BTreeMap<i16, Vec<i16>> = BTreeMap::new();
        e.insert(1, vec![1910]);
        e.insert(2, vec![957]);
        e.insert(3, vec![1680]);
        e.insert(4, vec![1847]);
        e.insert(5, vec![1746]);
        e.insert(6, vec![624]);
        e.insert(7, vec![818]);
        e.insert(8, vec![1103]);

        let raster_path = td!("cemsre_t55jfm_20200614_sub_abam5.tif");
        let vector_path = td!("cemsre_t55jfm_20200614_sub_abam5.gpkg");
        let bls: Vec<usize> = vec![16, 32, 64, 128, 256, 512];
        bls.iter().for_each(|b| {
            let data_source = DataSourceBuilder::from_file(&raster_path)
                .bands(vec![1])
                .build();
            let rds = RasterDatasetBuilder::<i16>::from_source(&data_source)
                .block_size(BlockSize { rows: *b, cols: *b })
                .build();
            let data = rds.extract_blockwise(&vector_path, "id", SamplingMethod::Value, None);
            assert_eq!(e, data);
        });
    }

    #[test]
    #[cfg(feature = "use_polars")]
    fn extract_blockwise_time() {
        let source = DataSource {
            source: td!("cemsre_t55jfm_20200614_sub_abam5.tif"),
            bands: [1, 3].to_vec(),
            layer_names: Vec::new(),
        };

        let vector_path = td!("cemsre_t55jfm_20200614_sub_abam5.gpkg");
        let overlap_size = 0;
        // two band,
        let rds: RasterDataset<i16> = RasterDatasetBuilder::from_source(&source)
            .overlap_size(overlap_size)
            .build();
        let _extracted_2b =
            rds.extract_blockwise(&vector_path, "id", SamplingMethod::Value, None);
        println!("Extracted: {_extracted_2b:?}");
    }

    #[test]
    fn test_extract_blockwise_padding() {
        let source = DataSource {
            source: td!("cemsre_t55jfm_20200614_sub_abam5.tif"),
            bands: [1, 3].to_vec(),
            layer_names: Vec::new(),
        };

        let vector_path = td!("cemsre_t55jfm_20200614_sub_abam5.gpkg");
        let overlap_size = 0;
        // two band,
        let rds: RasterDataset<i16> = RasterDatasetBuilder::<i16>::from_source(&source)
            .overlap_size(overlap_size)
            .build();
        let _extracted_2b = rds.extract_blockwise(&vector_path, "id", SamplingMethod::Value, None);
        println!("Extracted: {_extracted_2b:?}");

        // Using value
        let source = DataSource {
            source: td!("cemsre_t55jfm_20200614_sub_abam5.tif"),
            bands: [1, 3].to_vec(),
            layer_names: Vec::new(),
        };
        let rds: RasterDataset<i16> = RasterDatasetBuilder::<i16>::from_source(&source)
            .overlap_size(overlap_size)
            .build();
        let extracted_1 = rds.extract_blockwise(&vector_path, "id", SamplingMethod::Value, None);

        let overlap_size = 1;
        let source = DataSource {
            source: td!("cemsre_t55jfm_20200614_sub_abam5.tif"),
            bands: [1, 3].to_vec(),
            layer_names: Vec::new(),
        };
        let rds: RasterDataset<i16> = RasterDatasetBuilder::<i16>::from_source(&source)
            .overlap_size(overlap_size)
            .build();
        let extracted_2 = rds.extract_blockwise(&vector_path, "id", SamplingMethod::Value, None);
        assert_eq!(extracted_1, extracted_2);

        // Using mode
        let overlap_size = 0;
        let source = DataSource {
            source: td!("cemsre_t55jfm_20200614_sub_abam5.tif"),
            bands: [1, 3].to_vec(),
            layer_names: Vec::new(),
        };
        let rds: RasterDataset<i16> = RasterDatasetBuilder::<i16>::from_source(&source)
            .overlap_size(overlap_size)
            .build();
        let extracted_mode = rds.extract_blockwise(&vector_path, "id", SamplingMethod::Value, None);
        assert_eq!(extracted_1, extracted_mode);
    }

    #[test]
    fn test_stack_datasets() {
        let src_fn = td!("cemsre_t55jfm_20200614_sub_abam5.tif");
        let source_1 = DataSourceBuilder::from_file(&src_fn).build();
        let source_2 = source_1.clone();
        let mut rds_1 = RasterDatasetBuilder::<i16>::from_source(&source_1)
            .band_composition_dimension(Dimension::Layer)
            .build();
        let rds_2 = RasterDatasetBuilder::<i16>::from_source(&source_2)
            .band_composition_dimension(Dimension::Layer)
            .build();
        rds_1.stack(&rds_2, Dimension::Layer);

        let block_id = 0;
        let data = rds_1.read_block::<i32>(block_id);
        assert_eq!(data.shape(), [1, 12, 1024, 1024]);
    }

    #[test]
    fn test_extend_datasets() {
        let source = DataSource {
            source: td!("cemsre_t55jfm_20200614_sub_abam5.tif"),
            bands: [1, 2, 3, 4, 5, 6].to_vec(),
            layer_names: Vec::new(),
        };
        let source_2 = source.clone();
        let mut rds_1 = RasterDatasetBuilder::<i16>::from_source(&source)
            .band_composition_dimension(Dimension::Layer)
            .build();
        let rds_2 = RasterDatasetBuilder::<i16>::from_source(&source_2)
            .band_composition_dimension(Dimension::Layer)
            .build();
        rds_1.extend(&rds_2);
        let block_id = 0;
        let data = rds_1.read_block::<i32>(block_id);
        assert_eq!(data.shape(), [2, 6, 1024, 1024]);
    }

    #[test]
    fn test_epsg() {
        fn worker(rdb: &RasterDataBlock<i16>) -> anyhow::Result<Array4<i16>> {
            Ok(rdb.data.to_owned())
        }

        let source = DataSource {
            source: td!("cemsre_t55jfm_20200614_suba_abam5.tif"),
            bands: [1, 2].to_vec(),
            layer_names: Vec::new(),
        };

        let rds = RasterDatasetBuilder::<i16>::from_source(&source)
            .set_epsg(3577)
            .build();
        assert_eq!(
            rds.metadata.shape,
            RasterDataShape {
                times: 1,
                layers: 2,
                rows: 1419,
                cols: 1613
            }
        );

        rds.apply(worker, 1, &PathBuf::from("test.tif")).ok();

        let rds = RasterDatasetBuilder::<i16>::from_source(&source)
            .band_composition_dimension(Dimension::Time)
            .set_epsg(4326)
            .build();
        assert_eq!(
            rds.metadata.shape,
            RasterDataShape {
                times: 2,
                layers: 1,
                rows: 1179,
                cols: 1559
            }
        );
    }

    #[test]
    fn test_layers_from_source() {
        let source = DataSource {
            source: td!("cemsre_t55jfm_20200614_sub_abam5.tif"),
            bands: [1, 2, 3, 4, 5, 6].to_vec(),
            layer_names: Vec::new(),
        };

        let rds = RasterDatasetBuilder::<i16>::from_source(&source)
            .band_composition_dimension(Dimension::Layer)
            .build();
        assert_eq!(
            rds.metadata.shape,
            RasterDataShape {
                times: 1,
                layers: 6,
                rows: 1580,
                cols: 1514
            }
        );

        let source = DataSource {
            source: td!("cemsre_t55jfm_20200614_sub_abam5.tif"),
            bands: [1, 2, 3, 4, 5, 6].to_vec(),
            layer_names: Vec::new(),
        };

        let rds = RasterDatasetBuilder::<i16>::from_source(&source)
            .band_composition_dimension(Dimension::Time)
            .build();
        assert_eq!(
            rds.metadata.shape,
            RasterDataShape {
                times: 6,
                layers: 1,
                rows: 1580,
                cols: 1514
            }
        );
    }

    #[test]
    fn test_zip() {
        let data_source = DataSourceBuilder::from_file(&PathBuf::from(td!(
            "cemsre_t55jfm_20200614_sub_abam5.tif"
        )))
        .build();

        let rds_1 = RasterDatasetBuilder::<i16>::from_source(&data_source).build();
        let rds_2 = RasterDatasetBuilder::<i16>::from_source(&data_source).build();

        for (b_1, b_2) in rds_1.iter().zip(rds_2.iter()) {
            let id = b_1.iter_index;
            let b1_data = b_1.parent.read_block::<i16>(id);
            let b2_data = b_2.parent.read_block::<i16>(id);
            let mask = b2_data.mapv(|v| (v > 1000) as i16);
            let _result = b1_data * mask;
        }
    }

    #[test]
    fn test_raster_data_source_builder() {
        let data_source = DataSourceBuilder::from_file(&PathBuf::from(
            "data/cemsre_t55jfm_20200614_sub_abam5.tif",
        ))
        .build();

        let _ = RasterDatasetBuilder::<i16>::from_source(&data_source).build();
    }

    #[test]
    fn test_data_source_builder() {
        let _data_source = DataSourceBuilder::from_file(&PathBuf::from(
            "data/cemsre_t55jfm_20200614_sub_abam5.tif",
        ))
        .build();
        let _data_source = DataSourceBuilder::from_file(&PathBuf::from(
            "data/cemsre_t55jfm_20200614_sub_abam5.tif",
        ))
        .bands(vec![1, 3])
        .build();
        let _e = [
            "cemsre_t55jfm_20200614_sub_abam5_1",
            "cemsre_t55jfm_20200614_sub_abam5_3",
        ];

        let _data_source = DataSourceBuilder::from_file(&PathBuf::from(
            "data/cemsre_t55jfm_20200614_sub_abam5.tif",
        ))
        .bands(vec![1, 4])
        .build();
    }

    #[test]
    fn test_compositions() {
        let data_source = DataSource {
            source: PathBuf::from("data/cemsre_t55jfm_20200614_sub_abam5.tif"),
            bands: [1, 2, 3].to_vec(),
            layer_names: Vec::new(),
        };
        let data_source_2 = DataSource {
            source: PathBuf::from("data/cemsre_t55jfm_20200614_sub_abam5.tif"),
            bands: [1, 2, 3].to_vec(),
            layer_names: Vec::new(),
        };
        let rds =
            RasterDatasetBuilder::<i16>::from_sources(&vec![data_source_2, data_source]).build();
        assert_eq!(
            rds.metadata.shape,
            RasterDataShape {
                times: 2,
                layers: 3,
                rows: 1580,
                cols: 1514
            }
        );

        let data_source = DataSource {
            source: td!("cemsre_t55jfm_20200614_sub_abam5.tif"),
            bands: [1, 2, 3].to_vec(),
            layer_names: Vec::new(),
        };
        let data_source_2 = DataSource {
            source: td!("cemsre_t55jfm_20200614_sub_abam5.tif"),
            bands: [1, 2, 3].to_vec(),
            layer_names: Vec::new(),
        };
        let rds = RasterDatasetBuilder::<i16>::from_sources(&vec![data_source_2, data_source])
            .band_composition_dimension(Dimension::Layer)
            .source_composition_dimension(Dimension::Layer)
            .build();
        assert_eq!(
            rds.metadata.shape,
            RasterDataShape {
                times: 1,
                layers: 6,
                rows: 1580,
                cols: 1514
            }
        );

        let data_source = DataSource {
            source: td!("cemsre_t55jfm_20200614_sub_abam5.tif"),
            bands: [1, 2, 3].to_vec(),
            layer_names: Vec::new(),
        };
        let data_source_2 = DataSource {
            source: td!("cemsre_t55jfm_20200614_sub_abam5.tif"),
            bands: [1, 2, 3].to_vec(),
            layer_names: Vec::new(),
        };
        let rds = RasterDatasetBuilder::<i16>::from_sources(&vec![data_source_2, data_source])
            .band_composition_dimension(Dimension::Time)
            .source_composition_dimension(Dimension::Time)
            .build();
        assert_eq!(
            rds.metadata.shape,
            RasterDataShape {
                times: 6,
                layers: 1,
                rows: 1580,
                cols: 1514
            }
        );

        let data_source = DataSource {
            source: td!("cemsre_t55jfm_20200614_sub_abam5.tif"),
            bands: [1, 2, 3].to_vec(),
            layer_names: Vec::new(),
        };
        let data_source_2 = DataSource {
            source: td!("cemsre_t55jfm_20200614_sub_abam5.tif"),
            bands: [1, 2, 3].to_vec(),
            layer_names: Vec::new(),
        };
        let rds = RasterDatasetBuilder::<i16>::from_sources(&vec![data_source_2, data_source])
            .band_composition_dimension(Dimension::Time)
            .source_composition_dimension(Dimension::Layer)
            .build();
        assert_eq!(
            rds.metadata.shape,
            RasterDataShape {
                times: 3,
                layers: 2,
                rows: 1580,
                cols: 1514
            }
        );
    }

    // ─── Parallel writer tests ───

    #[test]
    fn test_trimm_array3_asymmetric_trimming() {
        use ndarray::Array3;
        use crate::types::Overlap;

        // Create a 3-band 10x10 array with distinguishable values
        let mut data: Array3<u16> = Array3::zeros((3, 10, 10));
        for b in 0..3 {
            for r in 0..10 {
                for c in 0..10 {
                    data[[b, r, c]] = (r * 10 + c) as u16 + b as u16 * 100;
                }
            }
        }

        // Test 1: zero overlap (no trimming)
        let overlap = Overlap { left: 0, top: 0, right: 0, bottom: 0 };
        let trimmed = crate::array_ops::trimm_array3_asymmetric(&data, &overlap);
        assert_eq!(trimmed.shape(), &[3, 10, 10]);
        assert_eq!(trimmed[[0, 0, 0]], 0);

        // Test 2: symmetric overlap
        let overlap = Overlap { left: 2, top: 2, right: 2, bottom: 2 };
        let trimmed = crate::array_ops::trimm_array3_asymmetric(&data, &overlap);
        assert_eq!(trimmed.shape(), &[3, 6, 6]);
        assert_eq!(trimmed[[0, 0, 0]], data[[0, 2, 2]]);

        // Test 3: asymmetric overlap
        let overlap = Overlap { left: 1, top: 2, right: 3, bottom: 4 };
        let trimmed = crate::array_ops::trimm_array3_asymmetric(&data, &overlap);
        assert_eq!(trimmed.shape(), &[3, 4, 6]);
        assert_eq!(trimmed[[0, 0, 0]], data[[0, 2, 1]]);
    }

    #[test]
    fn test_create_output_geotiff_small() {
        use crate::parallel_writer::create_output_geotiff;
        use crate::types::GeoTransform;
 

        let dir = std::env::temp_dir();
        let path = dir.join(format!("test_output_{}.tif", std::process::id()));
        let _ = std::fs::remove_file(&path);

        let geo_transform = GeoTransform {
            x_ul: 500000.0,
            x_res: 30.0,
            x_rot: 0.0,
            y_ul: 7000000.0,
            y_rot: 0.0,
            y_res: -30.0,
        };

        let result = create_output_geotiff::<u16>(
            &path,
            &geo_transform,
            32755,
            64,
            64,
            3,
            0u16,
        );
        assert!(result.is_ok(), "Failed to create GeoTIFF: {:?}", result.err());
        assert!(path.exists(), "Output file was not created");

        // Verify with GDAL
        let ds = gdal::Dataset::open(&path).expect("Failed to open created GeoTIFF");
        assert_eq!(ds.raster_count(), 3);
        let size = ds.raster_size();
        assert_eq!(size.0, 64);  // cols
        assert_eq!(size.1, 64);  // rows

        let _ = std::fs::remove_file(&path);
    }

    #[test]
    fn test_write_block_roundtrip() {
        use crate::parallel_writer::{self, ParallelGeoTiffWriter};
        use crate::types::{GeoTransform, Offset, ReadWindow, Size};
        use ndarray::Array3;

        let dir = std::env::temp_dir();
        let path = dir.join(format!("test_write_{}.tif", std::process::id()));
        let _ = std::fs::remove_file(&path);

        let geo_transform = GeoTransform {
            x_ul: 500000.0,
            x_res: 30.0,
            x_rot: 0.0,
            y_ul: 7000000.0,
            y_rot: 0.0,
            y_res: -30.0,
        };

        // Create output file
        parallel_writer::create_output_geotiff::<u16>(
            &path,
            &geo_transform,
            32755,
            128,    // total cols
            128,    // total rows
            2,      // 2 bands
            0u16,
        ).expect("Failed to create output");

        let writer = ParallelGeoTiffWriter::new(
            path.clone(),
            geo_transform,
            32755,
            128,
            128,
            2,
        );

        // Write a block at offset (32, 32) with size (64, 64)
        let mut block_data: Array3<u16> = Array3::zeros((2, 64, 64));
        // Fill with distinguishable values
        for b in 0..2 {
            for r in 0..64 {
                for c in 0..64 {
                    block_data[[b, r, c]] = (r * 64 + c) as u16 + b as u16 * 10000;
                }
            }
        }

        let window = ReadWindow {
            offset: Offset { cols: 32, rows: 32 },
            size: Size { cols: 64, rows: 64 },
        };

        parallel_writer::write_block(
            &writer,
            block_data.view(),
            window,
        ).expect("Failed to write block");

        // Drop writer to flush/close
        drop(writer);

        // Verify: read back the data and check it matches
        let ds = gdal::Dataset::open(&path).expect("Failed to open for verification");
        assert_eq!(ds.raster_count(), 2);

        for band_idx in 1..=2 {
            let band = ds.rasterband(band_idx).expect("Failed to get band");
            let data = band.read_as::<u16>(
                (32, 32),    // offset
                (64, 64),    // window size
                (64, 64),    // output size
                None,
            ).expect("Failed to read band").to_array().expect("Failed to convert to array");

            for r in 0..64 {
                for c in 0..64 {
                    let expected = (r * 64 + c) as u16 + (band_idx - 1) as u16 * 10000;
                    assert_eq!(data[[r, c]], expected,
                        "Mismatch at band={}, row={}, col={}", band_idx, r, c);
                }
            }
        }

        let _ = std::fs::remove_file(&path);
    }

    // ─── COG Validation Tests ───

    #[test]
    fn test_build_overviews() {
        use crate::parallel_writer::{self, ParallelGeoTiffWriter};
        use crate::types::GeoTransform;

        let dir = std::env::temp_dir();
        let path = dir.join(format!("test_overviews_{}.tif", std::process::id()));
        let _ = std::fs::remove_file(&path);

        let geo_transform = GeoTransform {
            x_ul: 500000.0,
            x_res: 30.0,
            x_rot: 0.0,
            y_ul: 7000000.0,
            y_rot: 0.0,
            y_res: -30.0,
        };

        // Create output file
        parallel_writer::create_output_geotiff::<u16>(
            &path,
            &geo_transform,
            32755,
            128,
            128,
            1,
            0u16,
        ).expect("Failed to create output");

        let writer = ParallelGeoTiffWriter::new(
            path.clone(),
            geo_transform,
            32755,
            128,
            128,
            1,
        );

        // Build overviews
        writer.build_overviews("CUBIC", &[2, 4, 8]).expect("Failed to build overviews");
        drop(writer);

        // Verify overviews exist
        let ds = gdal::Dataset::open(&path).expect("Failed to open for verification");
        let band = ds.rasterband(1).expect("Failed to get band");
        let overview_count = band.overview_count().expect("Failed to get overview count");
        assert!(
            overview_count > 0,
            "Expected overviews to exist, but got count={}",
            overview_count
        );

        let _ = std::fs::remove_file(&path);
    }

    #[test]
    fn test_translate_to_cog() {
        use crate::gdal_utils::translate_to_cog;
        use crate::parallel_writer::{self, ParallelGeoTiffWriter};
        use crate::types::{GeoTransform, Offset, ReadWindow, Size};
        use ndarray::Array3;

        let dir = std::env::temp_dir();
        let intermediate = dir.join(format!("test_cog_inter_{}.tif", std::process::id()));
        let cog_output = dir.join(format!("test_cog_{}.tif", std::process::id()));
        let _ = std::fs::remove_file(&intermediate);
        let _ = std::fs::remove_file(&cog_output);

        let geo_transform = GeoTransform {
            x_ul: 500000.0,
            x_res: 30.0,
            x_rot: 0.0,
            y_ul: 7000000.0,
            y_rot: 0.0,
            y_res: -30.0,
        };

        // Use 1024x1024 to ensure COG driver auto-generates overviews
        let img_size = 1024usize;
        let block_size = 256usize;

        // Create intermediate file and write some data
        parallel_writer::create_output_geotiff::<u16>(
            &intermediate,
            &geo_transform,
            32755,
            img_size,
            img_size,
            1,
            0u16,
        ).expect("Failed to create intermediate");

        let writer = ParallelGeoTiffWriter::new(
            intermediate.clone(),
            geo_transform,
            32755,
            img_size,
            img_size,
            1,
        );

        // Write a block of data to make the file non-empty
        let mut block_data: Array3<u16> = Array3::zeros((1, block_size, block_size));
        for r in 0..block_size {
            for c in 0..block_size {
                block_data[[0, r, c]] = (r * block_size + c) as u16;
            }
        }
        let window = ReadWindow {
            offset: Offset { cols: block_size as isize, rows: block_size as isize },
            size: Size { cols: block_size as isize, rows: block_size as isize },
        };
        parallel_writer::write_block(&writer, block_data.view(), window)
            .expect("Failed to write block");
        drop(writer);

        // Translate to COG via in-process create_copy (AUTO overviews generated by COG driver)
        translate_to_cog(
            &intermediate, &cog_output,
            "DEFLATE", "NEAREST",
        ).expect("Failed to translate to COG");

        // Verify COG has overviews (auto-generated by COG driver)
        let ds = gdal::Dataset::open(&cog_output).expect("Failed to open COG");
        let band = ds.rasterband(1).expect("Failed to get band");
        let overview_count = band.overview_count().expect("Failed to get overview count");
        assert!(
            overview_count > 0,
            "COG should have auto-generated overviews, but got count={}",
            overview_count
        );

        // Verify COG has correct size
        let size = ds.raster_size();
        assert_eq!(size.0, img_size); // cols
        assert_eq!(size.1, img_size); // rows

        // Verify data integrity — read back the written block
        let band = ds.rasterband(1).expect("Failed to get band");
        let data = band.read_as::<u16>(
            (block_size as isize, block_size as isize),
            (block_size, block_size),
            (block_size, block_size),
            None,
        ).expect("Failed to read band").to_array().expect("Failed to convert to array");
        assert_eq!(data[[0, 0]], 0u16);
        assert_eq!(data[[1, 1]], (block_size + 1) as u16);

        // Verify overviews contain valid data (not all zeros)
        let ov = band.overview(0).expect("Failed to get overview");
        let ov_size_x = ov.x_size();
        let ov_size_y = ov.y_size();
        // Read from center of overview where data should exist
        let read_offset = ((ov_size_x as isize) / 2 - 8, (ov_size_y as isize) / 2 - 8);
        let ov_data = ov.read_as::<u16>(
            read_offset,
            (16, 16),
            (16, 16),
            None,
        ).expect("Failed to read overview").to_array().expect("Failed to convert overview to array");
        assert!(
            ov_data.iter().any(|v| *v != 0),
            "Overview should contain non-zero data"
        );

        let _ = std::fs::remove_file(&intermediate);
        let _ = std::fs::remove_file(&cog_output);
    }

    #[test]
    fn test_output_config_defaults() {
        use crate::types::{OutputConfig, OutputFormat};

        let config = OutputConfig::default();
        assert_eq!(config.format, OutputFormat::GeoTiff);
        assert_eq!(config.compression, "LZW");
        assert_eq!(config.overview_resampling, "CUBIC");
        assert_eq!(config.overview_levels, vec![2, 4, 8, 16, 32]);
    }

    #[test]
    fn test_output_format_variants() {
        use crate::types::OutputFormat;

        // Test all variants can be created
        let _geotiff = OutputFormat::GeoTiff;
        let _geotiff_overviews = OutputFormat::GeoTiffOverviews;
        let _cog = OutputFormat::COG;

        // Test default is GeoTiff
        assert_eq!(OutputFormat::default(), OutputFormat::GeoTiff);

        // Test equality
        assert_eq!(OutputFormat::GeoTiff, OutputFormat::GeoTiff);
        assert_ne!(OutputFormat::GeoTiff, OutputFormat::COG);
    }
}