rstiff 0.2.0

use crate::error::TiffError;

use gdal::raster::Buffer;
use gdal::raster::{GdalDataType, GdalType, RasterizeOptions, rasterize};
use gdal::vector::Geometry;
use gdal::vector::LayerAccess;
use gdal::{Dataset, DriverManager};

// 投影转化
use gdal::spatial_ref::{CoordTransform, SpatialRef};

// 文件大小压缩
use gdal::raster::RasterCreationOptions;

use ndarray::s;
use ndarray::{Array3, Axis};
use std::path::Path;

/// 栅格文件元数据（不加载像素数据）
///
/// 用于快速获取文件信息，计算窗口坐标等
#[derive(Debug, Clone)]
pub struct RasterInfo {
    pub width: usize,
    pub height: usize,
    pub bands: usize,
    pub geo_transform: [f64; 6],
    pub projection: String,
    pub nodata: Option<f64>,
    pub original_type: GdalDataType,
}

impl RasterInfo {
    /// 从文件读取元数据（不加载像素数据）
    pub fn from_file<P: AsRef<Path>>(path: P) -> Result<Self, TiffError> {
        let path_str = path.as_ref().to_string_lossy().to_string();
        let dataset = Dataset::open(&path)?;

        let geo_transform = dataset.geo_transform()?;
        let projection = dataset.projection();
        let (width, height) = dataset.raster_size();
        let bands = dataset.raster_count();

        if bands < 1 {
            return Err(TiffError::BandMissing(path_str, 1));
        }

        let first_band = dataset.rasterband(1)?;
        let nodata = first_band.no_data_value();
        let original_type = first_band.band_type();

        Ok(Self {
            width,
            height,
            bands,
            geo_transform,
            projection,
            nodata,
            original_type,
        })
    }

    /// 获取栅格的地理边界
    ///
    /// # 返回
    /// (min_x, min_y, max_x, max_y)
    pub fn bounds(&self) -> (f64, f64, f64, f64) {
        let gt = self.geo_transform;
        let w = self.width as f64;
        let h = self.height as f64;

        // 四个角点
        let x0 = gt[0];
        let y0 = gt[3];
        let x1 = gt[0] + w * gt[1] + h * gt[2];
        let y1 = gt[3] + w * gt[4] + h * gt[5];

        let min_x = x0.min(x1);
        let max_x = x0.max(x1);
        let min_y = y0.min(y1);
        let max_y = y0.max(y1);

        (min_x, min_y, max_x, max_y)
    }

    /// 获取像素分辨率
    ///
    /// # 返回
    /// (pixel_width, pixel_height) - 注意 height 通常为负值的绝对值
    pub fn res(&self) -> (f64, f64) {
        (self.geo_transform[1].abs(), self.geo_transform[5].abs())
    }

    /// 像素坐标转地理坐标
    ///
    /// # 参数
    /// - `col`: 列号（X 方向）
    /// - `row`: 行号（Y 方向）
    ///
    /// # 返回
    /// (x, y) 地理坐标
    pub fn pixel_to_geo(&self, col: usize, row: usize) -> (f64, f64) {
        let gt = self.geo_transform;
        let x = gt[0] + (col as f64) * gt[1] + (row as f64) * gt[2];
        let y = gt[3] + (col as f64) * gt[4] + (row as f64) * gt[5];
        (x, y)
    }

    /// 地理坐标转像素坐标
    ///
    /// # 参数
    /// - `x`: 地理 X 坐标
    /// - `y`: 地理 Y 坐标
    ///
    /// # 返回
    /// (col, row) 像素坐标
    pub fn geo_to_pixel(&self, x: f64, y: f64) -> (isize, isize) {
        let gt = self.geo_transform;

        // 解方程组：
        // x = gt[0] + col * gt[1] + row * gt[2]
        // y = gt[3] + col * gt[4] + row * gt[5]
        //
        // 对于无旋转情况 (gt[2] = gt[4] = 0):
        // col = (x - gt[0]) / gt[1]
        // row = (y - gt[3]) / gt[5]

        let det = gt[1] * gt[5] - gt[2] * gt[4];
        if det.abs() < 1e-10 {
            // 无旋转的简化情况
            let col = ((x - gt[0]) / gt[1]).floor() as isize;
            let row = ((y - gt[3]) / gt[5]).floor() as isize;
            (col, row)
        } else {
            // 有旋转的一般情况
            let dx = x - gt[0];
            let dy = y - gt[3];
            let col = ((gt[5] * dx - gt[2] * dy) / det).floor() as isize;
            let row = ((-gt[4] * dx + gt[1] * dy) / det).floor() as isize;
            (col, row)
        }
    }

    /// 将地理边界转换为像素窗口
    ///
    /// # 参数
    /// - `bounds`: (min_x, min_y, max_x, max_y)
    ///
    /// # 返回
    /// (x_off, y_off, width, height)
    pub fn bounds_to_window(
        &self,
        bounds: (f64, f64, f64, f64),
    ) -> Result<(usize, usize, usize, usize), TiffError> {
        let (min_x, min_y, max_x, max_y) = bounds;

        // 转换四个角点
        let (col1, row1) = self.geo_to_pixel(min_x, max_y); // 左上
        let (col2, row2) = self.geo_to_pixel(max_x, min_y); // 右下

        // 整理坐标
        use std::cmp::{max, min};
        let x_start = min(col1, col2);
        let y_start = min(row1, row2);
        let x_end = max(col1, col2);
        let y_end = max(row1, row2);

        // 裁剪到图像范围
        let x_off = x_start.max(0).min(self.width as isize) as usize;
        let y_off = y_start.max(0).min(self.height as isize) as usize;
        let x_end_clamped = x_end.max(0).min(self.width as isize) as usize;
        let y_end_clamped = y_end.max(0).min(self.height as isize) as usize;

        let width = x_end_clamped.saturating_sub(x_off);
        let height = y_end_clamped.saturating_sub(y_off);

        if width == 0 || height == 0 {
            return Err(TiffError::InvalidCropBounds(format!(
                "地理范围与栅格不重叠: bounds=({}, {}, {}, {}), 像素范围: x[{}~{}], y[{}~{}]",
                min_x, min_y, max_x, max_y, x_start, x_end, y_start, y_end
            )));
        }

        Ok((x_off, y_off, width, height))
    }
}

#[derive(Debug, Clone)]
pub struct GeoTiff {
    pub data: Array3<f64>,           //像素矩阵(Bands * Height * Width)
    pub geo_transform: [f64; 6],     //地理变换六参数
    pub projection: String,          //投影信息(WKT)
    pub nodata: Option<f64>,         //无效值标记
    pub original_type: GdalDataType, //原始数据类型
}

impl GeoTiff {
    /// 快速获取文件元数据（不加载像素数据）
    ///
    /// 这是 `RasterInfo::from_file()` 的便捷别名。
    ///
    /// # 示例
    /// ```ignore
    /// let info = GeoTiff::info("large.tif")?;
    /// println!("尺寸: {}x{}", info.width, info.height);
    /// println!("边界: {:?}", info.bounds());
    /// ```
    pub fn info<P: AsRef<Path>>(path: P) -> Result<RasterInfo, TiffError> {
        RasterInfo::from_file(path)
    }

    pub fn read<P: AsRef<Path>>(path: P) -> Result<Self, TiffError> {
        let path_str = path.as_ref().to_string_lossy().to_string(); //用来打印日志

        // 打开文件
        let dataset = Dataset::open(&path)?;

        // 获取关键元数据
        let geo_transform = dataset.geo_transform()?;
        let projection = dataset.projection();

        // 获取波段信息
        let raster_count = dataset.raster_count();
        if raster_count < 1 {
            return Err(TiffError::BandMissing(path_str, 1));
        }
        let (width, height) = dataset.raster_size();

        //以第一景数据参数作为全局数据参数
        let first_band = dataset.rasterband(1)?;
        let mut nodata = first_band.no_data_value();
        let original_type = first_band.band_type();

        // --- 科研补丁：如果元数据没设 NoData，我们通常默认 0.0 是无效背景 ---
        if nodata.is_none() {
            // println!("警告: 影像未定义 NoData，默认将 0.0 视为背景。");
            nodata = Some(f64::NAN);
        }

        let mut flat_data: Vec<f64> = Vec::with_capacity(raster_count * height * width);

        for i in 1..=raster_count {
            let band = dataset.rasterband(i)?;

            let buffer = band.read_as::<f64>(
                (0, 0),          //起点
                (width, height), //读取尺寸
                (width, height), //输出尺寸
                None,            //重采样方式
            )?;

            flat_data.extend_from_slice(buffer.data());
        }

        let data = Array3::from_shape_vec((raster_count, height, width), flat_data)?;

        Ok(Self {
            data,
            geo_transform,
            projection,
            nodata,
            original_type,
        })
    }

    /// 窗口化读取：只读取指定像素范围的数据（节省内存）
    ///
    /// # 参数
    /// - `path`: GeoTiff 文件路径
    /// - `x_off`: 起始列（像素坐标）
    /// - `y_off`: 起始行（像素坐标）
    /// - `width`: 读取宽度（像素）
    /// - `height`: 读取高度（像素）
    ///
    /// # 示例
    /// ```ignore
    /// // 只读取左上角 256x256 的区域
    /// let tile = GeoTiff::read_window("large.tif", 0, 0, 256, 256)?;
    /// ```
    pub fn read_window<P: AsRef<Path>>(
        path: P,
        x_off: usize,
        y_off: usize,
        width: usize,
        height: usize,
    ) -> Result<Self, TiffError> {
        let path_str = path.as_ref().to_string_lossy().to_string();

        // 打开文件
        let dataset = Dataset::open(&path)?;

        // 获取关键元数据
        let geo_transform = dataset.geo_transform()?;
        let projection = dataset.projection();

        // 获取波段信息
        let raster_count = dataset.raster_count();
        if raster_count < 1 {
            return Err(TiffError::BandMissing(path_str.clone(), 1));
        }
        let (raster_width, raster_height) = dataset.raster_size();

        // 检查窗口是否超出范围
        if x_off + width > raster_width || y_off + height > raster_height {
            return Err(TiffError::InvalidCropBounds(format!(
                "窗口超出栅格范围: 窗口({}, {}, {}, {}), 栅格尺寸({}, {})",
                x_off, y_off, width, height, raster_width, raster_height
            )));
        }

        // 以第一景数据参数作为全局数据参数
        let first_band = dataset.rasterband(1)?;
        let mut nodata = first_band.no_data_value();
        let original_type = first_band.band_type();

        if nodata.is_none() {
            nodata = Some(f64::NAN);
        }

        let mut flat_data: Vec<f64> = Vec::with_capacity(raster_count * height * width);

        for i in 1..=raster_count {
            let band = dataset.rasterband(i)?;

            // 关键：指定起始位置和读取尺寸
            let buffer = band.read_as::<f64>(
                (x_off as isize, y_off as isize), // 起点
                (width, height),                  // 读取尺寸
                (width, height),                  // 输出尺寸
                None,                             // 重采样方式
            )?;

            flat_data.extend_from_slice(buffer.data());
        }

        let data = Array3::from_shape_vec((raster_count, height, width), flat_data)?;

        // 更新地理变换参数（左上角坐标偏移）
        let new_geo_transform = [
            geo_transform[0]
                + (x_off as f64) * geo_transform[1]
                + (y_off as f64) * geo_transform[2],
            geo_transform[1],
            geo_transform[2],
            geo_transform[3]
                + (x_off as f64) * geo_transform[4]
                + (y_off as f64) * geo_transform[5],
            geo_transform[4],
            geo_transform[5],
        ];

        Ok(Self {
            data,
            geo_transform: new_geo_transform,
            projection,
            nodata,
            original_type,
        })
    }

    /// 根据地理坐标范围读取数据
    ///
    /// # 参数
    /// - `path`: GeoTiff 文件路径
    /// - `bounds`: 地理坐标范围 (min_x, min_y, max_x, max_y)
    ///
    /// # 示例
    /// ```ignore
    /// // 读取指定经纬度范围的数据
    /// let roi = GeoTiff::read_bounds("large.tif", (116.0, 39.0, 117.0, 40.0))?;
    /// ```
    pub fn read_bounds<P: AsRef<Path>>(
        path: P,
        bounds: (f64, f64, f64, f64), // (min_x, min_y, max_x, max_y)
    ) -> Result<Self, TiffError> {
        // 先获取文件的元数据
        let info = RasterInfo::from_file(&path)?;

        // 计算像素窗口
        let (x_off, y_off, width, height) = info.bounds_to_window(bounds)?;

        // 使用窗口读取
        Self::read_window(path, x_off, y_off, width, height)
    }

    /// 根据矢量文件范围读取数据（不加载整个栅格）
    ///
    /// # 参数
    /// - `raster_path`: GeoTiff 文件路径
    /// - `vector_path`: 矢量文件路径（KML、Shapefile、GeoJSON 等）
    /// - `apply_mask`: 是否将多边形外的像素设为 NoData
    ///
    /// # 示例
    /// ```ignore
    /// // 只读取 KML 范围内的数据
    /// let roi = GeoTiff::read_by_vector("large.tif", "area.kml", true)?;
    /// ```
    pub fn read_by_vector<P: AsRef<Path>, Q: AsRef<Path>>(
        raster_path: P,
        vector_path: Q,
        apply_mask: bool,
    ) -> Result<Self, TiffError> {
        let vector_path_str = vector_path.as_ref().to_string_lossy().to_string();

        // 1. 读取矢量文件获取范围
        let vector_ds = Dataset::open(&vector_path)?;
        let mut layer = vector_ds.layer(0)?;
        let extent = layer.get_extent()?;

        // 2. 获取栅格元数据
        let info = RasterInfo::from_file(&raster_path)?;

        // 3. 检查投影一致性
        if let Some(vector_srs) = layer.spatial_ref() {
            if let Ok(vector_wkt) = vector_srs.to_wkt() {
                if !info.projection.is_empty()
                    && !vector_wkt.is_empty()
                    && info.projection != vector_wkt
                {
                    eprintln!("警告: 矢量投影与栅格投影可能不一致！建议先进行重投影。");
                }
            }
        }

        // 4. 计算像素窗口
        let bounds = (extent.MinX, extent.MinY, extent.MaxX, extent.MaxY);
        let (x_off, y_off, width, height) = info.bounds_to_window(bounds)?;

        println!(
            "矢量窗口读取: Path={}, x={}, y={}, w={}, h={}",
            vector_path_str, x_off, y_off, width, height
        );

        // 5. 窗口化读取栅格数据
        let mut result = Self::read_window(&raster_path, x_off, y_off, width, height)?;

        // 6. 如果需要掩膜
        if apply_mask {
            let no_data_val = result.nodata.unwrap_or(-9999.0);
            result.nodata = Some(no_data_val);

            // 创建内存数据集用于生成 Mask
            let driver = DriverManager::get_driver_by_name("MEM")?;
            let mut mask_ds = driver.create_with_band_type::<u8, _>("", width, height, 1)?;

            mask_ds.set_geo_transform(&result.geo_transform)?;
            mask_ds.set_projection(&result.projection)?;

            let options = RasterizeOptions {
                all_touched: true,
                merge_algorithm: gdal::raster::MergeAlgorithm::Replace,
                ..Default::default()
            };

            // 收集几何体
            let geometries: Vec<Geometry> = layer
                .features()
                .filter_map(|f| f.geometry().cloned())
                .collect();

            if geometries.is_empty() {
                return Err(TiffError::Gdal(gdal::errors::GdalError::NullPointer {
                    method_name: "rasterize",
                    msg: "矢量文件中没有有效的几何体".into(),
                }));
            }

            // 执行栅格化
            rasterize(&mut mask_ds, &[1], &geometries, &[1.0], Some(options))?;

            // 读取 Mask 数据
            let mask_band = mask_ds.rasterband(1)?;
            let mask_data =
                mask_band.read_as::<u8>((0, 0), (width, height), (width, height), None)?;
            let mask_slice = mask_data.data();

            // 应用掩膜
            for mut band_view in result.data.outer_iter_mut() {
                if let Some(slice) = band_view.as_slice_mut() {
                    for (pixel, &m) in slice.iter_mut().zip(mask_slice.iter()) {
                        if m == 0 {
                            *pixel = no_data_val;
                        }
                    }
                } else {
                    for (pixel, &m) in band_view.iter_mut().zip(mask_slice.iter()) {
                        if m == 0 {
                            *pixel = no_data_val;
                        }
                    }
                }
            }
        }

        Ok(result)
    }

    pub fn write<P: AsRef<Path>>(&self, path: P) -> Result<(), TiffError> {
        let path = path.as_ref();

        // 辅助宏：简化 NoData 计算
        // return (Option<T> [fill_value], Option<f64> [metadata_value])
        // 如果 self.nodata 存在且不是 NaN => 它是有效的 NoData，需要填充且需要标记
        // 如果 self.nodata 是 NaN => 它是无效的 NoData，我们需要填充一个默认值 (如 255/65535)，并且在 metadata 中标记它，以便启用透明
        macro_rules! resolve_nodata {
            ($t:ty, $default:expr) => {{
                if let Some(nd) = self.nodata {
                    if nd.is_finite() {
                        // 有效 nodata: 填充该值，标记该值
                        (Some(nd as $t), Some(nd))
                    } else {
                        // NaN nodata: 填充 max/min 默认值，且标记该值为 NoData (启用透明)
                        (Some($default), Some($default as f64))
                    }
                } else {
                    // None: 不填充，不标记
                    (None, None)
                }
            }};
        }

        match self.original_type {
            GdalDataType::UInt8 => {
                // Byte -> 255 (避免 0 冲突)
                let (fill, meta) = resolve_nodata!(u8, u8::MAX);
                self.write_impl(path, fill, meta, |v| v as u8)
            }
            GdalDataType::UInt16 => {
                // UInt16 -> 65535
                let (fill, meta) = resolve_nodata!(u16, u16::MAX);
                self.write_impl(path, fill, meta, |v| v as u16)
            }
            GdalDataType::Int16 => {
                // Int16 -> -32768 (min)
                let (fill, meta) = resolve_nodata!(i16, i16::MIN);
                self.write_impl(path, fill, meta, |v| v as i16)
            }
            GdalDataType::UInt32 => {
                // UInt32 -> MAX
                let (fill, meta) = resolve_nodata!(u32, u32::MAX);
                self.write_impl(path, fill, meta, |v| v as u32)
            }
            GdalDataType::Int32 => {
                // Int32 -> MIN
                let (fill, meta) = resolve_nodata!(i32, i32::MIN);
                self.write_impl(path, fill, meta, |v| v as i32)
            }
            GdalDataType::Float32 => {
                // 对于浮点数，直接使用 self.nodata 即可 (通常是 NaN 或 -9999.0)
                let (fill, meta) = if let Some(val) = self.nodata {
                    (Some(val as f32), Some(val))
                } else {
                    (None, None)
                };
                self.write_impl(path, fill, meta, |v| v as f32)
            }
            _ => {
                // Default f64
                self.write_impl(path, self.nodata, self.nodata, |v| v)
            }
        }
    }

    fn write_impl<T: GdalType + Copy + PartialEq>(
        &self,
        path: &Path,
        fill_value: Option<T>,       // 遇到源 NoData/NaN 时，填充到像素中的值
        metadata_value: Option<f64>, // 写入 GDAL Metadata 的值 (如果 None 则不写)
        convert: impl Fn(f64) -> T,
    ) -> Result<(), TiffError> {
        // 实现写入逻辑
        let (bands, height, width) = self.data.dim();

        // 获取GTiff驱动
        let driver = DriverManager::get_driver_by_name("GTiff")?;

        // --- 开启压缩选项 ---
        let mut options = RasterCreationOptions::default();
        // LZW 是无损压缩，PREDICTOR=2 专门优化连续数值（如形变梯度）
        options.set_name_value("COMPRESS", "LZW")?;
        options.set_name_value("PREDICTOR", "2")?;
        options.set_name_value("TILED", "YES")?;
        options.set_name_value("BIGTIFF", "IF_SAFER")?;

        // 创建新文件
        let mut dataset = driver.create_with_band_type::<T, _>(
            path, //路径
            width as usize,
            height as usize,
            bands as usize,
        )?;

        // 3. 设置地理信息
        dataset.set_geo_transform(&self.geo_transform)?;
        dataset.set_projection(&self.projection)?;

        for i in 0..bands {
            let mut band = dataset.rasterband(i + 1)?;
            if let Some(val) = metadata_value {
                band.set_no_data_value(Some(val))?;
            }

            let band_view = self.data.index_axis(Axis(0), i);

            // 数据转换
            // 我们需要处理 NoData 的转换：如果源数据是 NoData，则写入 fill_value
            let src_nodata = self.nodata;
            let map_pixel = |v: &f64| -> T {
                let v = *v;
                let is_nodata = if let Some(snd) = src_nodata {
                    if snd.is_nan() { v.is_nan() } else { v == snd }
                } else {
                    false
                };

                if is_nodata {
                    if let Some(fill) = fill_value {
                        return fill;
                    }
                }
                convert(v)
            };

            let band_vec: Vec<T> = if let Some(slice) = band_view.as_slice() {
                slice.iter().map(map_pixel).collect()
            } else {
                band_view.iter().map(map_pixel).collect()
            };

            let mut buffer = Buffer::new(
                (width, height), // 尺寸
                band_vec,        // 数据 (Vec<T>)
            );

            band.write((0, 0), (width as usize, height as usize), &mut buffer)?;
        }

        dataset.flush_cache()?;

        Ok(())
    }

    // 裁剪基础函数
    pub fn crop(
        &self,
        x_off: usize,
        y_off: usize,
        width: usize,
        height: usize,
    ) -> Result<Self, TiffError> {
        let (_bands, max_h, max_w) = self.data.dim();
        if x_off + width > max_w || y_off + height > max_h {
            let msg = format!(
                "裁剪参数超出范围: x_off={}, y_off={}, width={}, height={}, max_w={}, max_h={}",
                x_off, y_off, width, height, max_w, max_h
            );
            return Err(TiffError::InvalidCropBounds(msg));
        }

        //获取裁剪后的新数据
        let new_data = self
            .data
            .slice(s![.., y_off..y_off + height, x_off..x_off + width])
            .to_owned();

        //更新新的地理变换参数
        let mut new_gt = self.geo_transform;

        let x_shift = x_off as f64;
        let y_shift = y_off as f64;

        new_gt[0] = self.geo_transform[0]
            + (x_shift * self.geo_transform[1])
            + (y_shift * self.geo_transform[2]);
        new_gt[3] = self.geo_transform[3]
            + (x_shift * self.geo_transform[4])
            + (y_shift * self.geo_transform[5]);

        Ok(Self {
            data: new_data,
            geo_transform: new_gt,
            projection: self.projection.clone(), //投影保持不变
            nodata: self.nodata,
            original_type: self.original_type,
        })
    }

    // 根据矢量数据进行裁剪
    pub fn crop_by_vector<P: AsRef<Path>>(
        &self,
        vector_path: P,
        apply_mask: bool,
    ) -> Result<Self, TiffError> {
        let path_str = vector_path.as_ref().to_string_lossy().to_string();

        // 1. 读取矢量
        let dataset = Dataset::open(&vector_path).map_err(TiffError::Gdal)?;
        let mut layer = dataset.layer(0).map_err(TiffError::Gdal)?;

        // --- 关键步骤：检查投影是否一致 ---
        // 这一步非常重要，否则裁剪框会飞到十万八千里外
        let raster_srs = self.projection.clone();
        if let Some(vector_srs) = layer.spatial_ref() {
            if let Ok(vector_wkt) = vector_srs.to_wkt() {
                // 简单的字符串比对，虽然不完美但能拦截大部分错误
                // 实际工程中可能需要用 osr::SpatialReference::is_same
                if !raster_srs.is_empty() && !vector_wkt.is_empty() && raster_srs != vector_wkt {
                    eprintln!("警告: 矢量投影与栅格投影可能不一致！建议先进行重投影。");
                    // 可以在这里返回 Err，或者继续尝试
                }
            }
        }

        // 2. 获取范围
        let extent = layer.get_extent().map_err(TiffError::Gdal)?;

        // 3. 计算裁剪窗口
        let (x_off, y_off, width, height) = self.compute_crop_window(extent, self.geo_transform)?;

        println!(
            "矢量裁剪: Path={}, x={}, y={}, w={}, h={}",
            path_str, x_off, y_off, width, height
        );

        // 4. 执行基础裁剪 (这一步得到的 cropped_tif 是 Array3)
        let mut cropped_tif = self.crop(x_off, y_off, width, height)?;

        // 5. 如果需要掩膜 (Masking)
        if apply_mask {
            // 确定 NoData 值，默认取 -9999.0 或 NAN
            let no_data_val = cropped_tif.nodata.unwrap_or(-9999.0);
            cropped_tif.nodata = Some(no_data_val);

            // 创建内存数据集用于生成 Mask (单波段 u8)
            let driver = DriverManager::get_driver_by_name("MEM").map_err(TiffError::Gdal)?;
            let mut mask_ds = driver
                .create_with_band_type::<u8, _>("", width, height, 1)
                .map_err(TiffError::Gdal)?;

            // 必须设置地理参考，否则栅格化位置不对
            mask_ds
                .set_geo_transform(&cropped_tif.geo_transform)
                .map_err(TiffError::Gdal)?;
            mask_ds
                .set_projection(&cropped_tif.projection)
                .map_err(TiffError::Gdal)?;

            let options = RasterizeOptions {
                all_touched: true, // true: 只要碰到像素就涂色; false: 中心点在多边形内才涂色
                merge_algorithm: gdal::raster::MergeAlgorithm::Replace,
                ..Default::default()
            };

            // 收集几何体
            let geometries: Vec<Geometry> = layer
                .features()
                .filter_map(|f| f.geometry().cloned())
                .collect();

            if geometries.is_empty() {
                return Err(TiffError::Gdal(gdal::errors::GdalError::NullPointer {
                    method_name: "rasterize",
                    msg: "矢量文件中没有有效的几何体".into(),
                }));
            }

            // 执行栅格化：多边形内部为 1，外部默认为 0
            rasterize(&mut mask_ds, &[1], &geometries, &[1.0], Some(options))?;

            // 读取 Mask 数据
            let mask_band = mask_ds.rasterband(1).map_err(TiffError::Gdal)?;
            let mask_data = mask_band
                .read_as::<u8>((0, 0), (width, height), (width, height), None)
                .map_err(TiffError::Gdal)?;
            let mask_slice = mask_data.data();

            // --- 核心修复：针对三维数组应用掩膜 ---
            // cropped_tif.data 是 Array3 [Bands, Height, Width]
            // 我们必须遍历每一个波段，分别应用同一个 2D Mask

            for mut band_view in cropped_tif.data.outer_iter_mut() {
                // band_view 现在是 2D 视图 (Height, Width)
                if let Some(slice) = band_view.as_slice_mut() {
                    // 快速路径 (内存连续)
                    // slice 长度 = H*W, mask_slice 长度 = H*W，这下 zip 对齐了
                    for (pixel, &m) in slice.iter_mut().zip(mask_slice.iter()) {
                        if m == 0 {
                            *pixel = no_data_val;
                        }
                    }
                } else {
                    // 慢速路径 (内存不连续)
                    for (pixel, &m) in band_view.iter_mut().zip(mask_slice.iter()) {
                        if m == 0 {
                            *pixel = no_data_val;
                        }
                    }
                }
            }
        }

        Ok(cropped_tif)
    }

    fn compute_crop_window(
        &self,
        extent: gdal::vector::Envelope,
        geo_transform: [f64; 6],
    ) -> Result<(usize, usize, usize, usize), TiffError> {
        let px1 = ((extent.MinX - geo_transform[0]) / geo_transform[1] + 0.5) as isize;
        let px2 = ((extent.MaxX - geo_transform[0]) / geo_transform[1] + 0.5) as isize;

        let py1 = ((extent.MinY - geo_transform[3]) / geo_transform[5] + 0.5) as isize;
        let py2 = ((extent.MaxY - geo_transform[3]) / geo_transform[5] + 0.5) as isize;

        // 整理坐标
        use std::cmp::{max, min};
        let x_start = min(px1, px2);
        let y_start = min(py1, py2);
        let x_end = max(px1, px2);
        let y_end = max(py1, py2);

        let (_, h, w) = self.data.dim();

        // 裁剪图像
        let x_off = x_start.max(0).min(w as isize) as usize;
        let y_off = y_start.max(0).min(h as isize) as usize;

        let x_end_clamped = x_end.max(0).min(w as isize) as usize;
        let y_end_clamped = y_end.max(0).min(h as isize) as usize;
        let width = x_end_clamped.saturating_sub(x_off); //无符号整型不能为负数，该函数更完备
        let height = y_end_clamped.saturating_sub(y_off);

        if width == 0 || height == 0 {
            return Err(TiffError::InvalidCropBounds(format!(
                "矢量范围与图像不重叠或无效。矢量像素范围: x[{}~{}], y[{}~{}]",
                x_start, x_end, y_start, y_end
            )));
        }
        Ok((x_off, y_off, width, height))
    }

    pub fn reproject(&self, target_epsg: i32) -> Result<GeoTiff, TiffError> {
        // 获取源和目标投影
        let mut src_srs = SpatialRef::from_wkt(&self.projection)?;
        let mut target_srs = SpatialRef::from_epsg(target_epsg as u32)?;

        // --- 核心修复：强制使用传统的 [经度, 纬度] 顺序 ---
        // TraditionalGisOrder 保证了坐标总是 [东/经, 北/纬]，这也是 GeoTransform 的存储方式
        src_srs
            .set_axis_mapping_strategy(gdal::spatial_ref::AxisMappingStrategy::TraditionalGisOrder);
        target_srs
            .set_axis_mapping_strategy(gdal::spatial_ref::AxisMappingStrategy::TraditionalGisOrder);

        // 尝试从源投影中提取 EPSG 代码字符串
        if let Ok(src_epsg) = src_srs.auth_code() {
            if src_epsg == target_epsg {
                // 进一步确认 Authority 是 EPSG
                if let Some(auth_name) = src_srs.auth_name() {
                    if auth_name == "EPSG" {
                        println!("源投影与目标投影一致 (EPSG:{}), 跳过转换。", target_epsg);
                        return Ok(self.clone());
                    }
                }
            }
        }

        // 计算坐标转化
        let transform = CoordTransform::new(&src_srs, &target_srs)?;

        // 4. 计算重投影后的新边界 (调用我们之前讨论的辅助函数)
        let (new_width, new_height, new_gt) = self.compute_reprojected_bounds(&transform)?;
        // --- 5. 创建源数据集 (内存 MEM) ---
        let (bands, h, w) = self.data.dim();
        let driver = DriverManager::get_driver_by_name("MEM")?;

        let mut src_ds = driver.create_with_band_type::<f64, _>("", w, h, bands)?;
        src_ds.set_geo_transform(&self.geo_transform)?;
        src_ds.set_projection(&self.projection)?;

        for b in 0..bands {
            let mut band = src_ds.rasterband(b + 1)?;
            if let Some(nd) = self.nodata {
                band.set_no_data_value(Some(nd))?;
            }

            // 提取单波段数据并写入
            let data_vec: Vec<f64> = self
                .data
                .index_axis(ndarray::Axis(0), b)
                .iter()
                .cloned()
                .collect();
            let mut buffer = Buffer::new((w, h), data_vec);
            band.write((0, 0), (w, h), &mut buffer)?;
        }

        // --- 6. 创建目标数据集 (内存 MEM) ---
        let mut dst_ds =
            driver.create_with_band_type::<f64, _>("", new_width, new_height, bands)?;
        dst_ds.set_geo_transform(&new_gt)?;
        let target_wkt = target_srs.to_wkt()?;
        dst_ds.set_projection(&target_wkt)?;

        // 初始化目标盘的 NoData 背景
        for b in 1..=bands {
            let mut band = dst_ds.rasterband(b)?;
            if let Some(nd) = self.nodata {
                // 1. 告诉 GDAL 这个波段的无效值是多少
                // 注意：在某些 gdal 版本中是 set_no_data_value(nd)，有些是 Some(nd)
                // 如果报错，请尝试将 Some(nd) 改为 nd
                band.set_no_data_value(Some(nd))?;

                // 2. 预填充背景
                // 第一个参数是实部 (nd)，第二个参数是虚部 (0.0)
                // 即使你的数据不是复数，也必须提供第二个参数
                band.fill(nd, Some(0.0))?;
            } else {
                // 科学研究建议：如果没有定义 nodata，重投影默认会产生 0 黑边
                // 建议给 InSAR 数据默认一个 nan 或 -9999.0
                band.fill(0.0, Some(0.0))?;
            }
        }

        // --- 7. 执行重投影 (Warp) ---
        // 注意：这里建议显式使用 gdal::raster::reproject
        gdal::raster::reproject(&src_ds, &dst_ds)?;

        // --- 8. 读取结果并封装回 Array3 ---
        let mut new_flat_data = Vec::with_capacity(bands * new_height * new_width);
        for b in 1..=bands {
            let band = dst_ds.rasterband(b)?;
            let buffer = band.read_as::<f64>(
                (0, 0),
                (new_width, new_height),
                (new_width, new_height),
                None,
            )?;
            new_flat_data.extend_from_slice(buffer.data());
        }

        let new_data = Array3::from_shape_vec((bands, new_height, new_width), new_flat_data)
            .map_err(|e| TiffError::InvalidCropBounds(e.to_string()))?;

        Ok(GeoTiff {
            data: new_data,
            geo_transform: new_gt,
            projection: target_wkt,
            nodata: self.nodata,
            original_type: self.original_type,
        })
    }

    /// 辅助函数：计算重投影后的图像边界和分辨率
    fn compute_reprojected_bounds(
        &self,
        transform: &CoordTransform,
    ) -> Result<(usize, usize, [f64; 6]), TiffError> {
        let (_, h, w) = self.data.dim();
        let gt = self.geo_transform;

        // 定义图像的四个角点 (像素坐标)
        let corners_pixel = vec![
            (0.0, 0.0),           // 左上
            (w as f64, 0.0),      // 右上
            (w as f64, h as f64), // 右下
            (0.0, h as f64),      // 左下
        ];

        // 1. 将像素坐标 -> 源地理坐标 -> 目标地理坐标
        let mut xs = Vec::new();
        let mut ys = Vec::new();

        for (px, py) in corners_pixel {
            // Pixel -> Geo (Source)
            let (x, y) = self.pixel_to_projected(px, py, gt, transform)?;
            xs.push(x);
            ys.push(y);
        }

        // 2. 获取目标坐标系下的 Bounding Box
        let min_x = xs.iter().cloned().fold(f64::INFINITY, f64::min);
        let max_x = xs.iter().cloned().fold(f64::NEG_INFINITY, f64::max);
        let min_y = ys.iter().cloned().fold(f64::INFINITY, f64::min);
        let max_y = ys.iter().cloned().fold(f64::NEG_INFINITY, f64::max);

        // 3. 估算新的分辨率 (Pixel Size)
        // 简单方法：利用中心点附近的映射比例来估算
        // 取中心点 (cx, cy) 和 (cx+1, cy+1)
        let cx = w as f64 / 2.0;
        let cy = h as f64 / 2.0;

        let p0 = self.pixel_to_projected(cx, cy, gt, transform)?;
        let p1 = self.pixel_to_projected(cx + 1.0, cy, gt, transform)?; // 右移一像素
        let p2 = self.pixel_to_projected(cx, cy + 1.0, gt, transform)?; // 下移一像素

        // 计算目标坐标系下的距离
        let res_x = ((p1.0 - p0.0).powi(2) + (p1.1 - p0.1).powi(2)).sqrt();
        let res_y = ((p2.0 - p0.0).powi(2) + (p2.1 - p0.1).powi(2)).sqrt();

        // 4. 计算新的宽高
        let width_geo = max_x - min_x;
        let height_geo = max_y - min_y;

        let new_w = (width_geo / res_x).ceil() as usize;
        let new_h = (height_geo / res_y).ceil() as usize;

        // 5. 构建新的 GeoTransform
        // [top_left_x, res_x, 0, top_left_y, 0, -res_y] (假设无旋转，且为北向)
        // 注意：max_y 通常是左上角的 Y (在北半球/大多数投影中)
        let new_gt = [min_x, res_x, 0.0, max_y, 0.0, -res_y];

        Ok((new_w, new_h, new_gt))
    }

    // 小助手：将源像素坐标转换为目标投影坐标
    fn pixel_to_projected(
        &self,
        px: f64,
        py: f64,
        gt: [f64; 6],
        transform: &CoordTransform,
    ) -> Result<(f64, f64), TiffError> {
        // 1. 像素坐标 -> 源地理坐标 (GeoTransform)
        let geo_x = gt[0] + px * gt[1] + py * gt[2];
        let geo_y = gt[3] + px * gt[4] + py * gt[5];

        // 2. 源地理坐标 -> 目标地理坐标 (CoordTransform)
        // 准备输入数据（X, Y, Z）
        let mut x = [geo_x];
        let mut y = [geo_y];
        let mut z = [0.0];

        // 使用 transform_coords，它会直接修改数组中的值
        transform.transform_coords(&mut x, &mut y, &mut z)?;

        Ok((x[0], y[0]))
    }
}

#[test]
fn test_reproject() -> Result<(), TiffError> {
    let p = std::path::Path::new("./data/Hawaiin.tif");
    let tif = GeoTiff::read(p)?;
    println!("tif_nodata: {:?}", tif.nodata);

    // --- 测试重投影：从 WGS84 (4326) 转到 UTM 4N (32604) ---
    // 对于 InSAR 博士生来说，将地理坐标转为投影坐标是计算形变量（米）的基础
    let target_epsg = 32604;
    println!("正在开始重投影至 EPSG:{}...", target_epsg);

    let tif_projected = tif.reproject(target_epsg as i32)?;

    println!("重投影完成！新尺寸: {:?}", tif_projected.data.dim());
    println!("新投影 WKT: {}", &tif_projected.projection[..100]); // 打印前100个字符查看

    let path = std::path::Path::new("./data/Hawaiin_UTM4N.tif");
    tif_projected.write(path)?;
    println!("结果已保存至: {}", path.display());

    Ok(())
}