surtgis_core/
resample.rs

1//! Raster resampling: align a raster to a reference grid.
2//!
3//! Resamples a source raster to match the exact grid (origin, cell size,
4//! dimensions) of a reference raster. Essential for stacking rasters from
5//! different sources (e.g., Copernicus DEM at 30m + Sentinel-2 at 10m).
6
7use ndarray::Array2;
8
9use crate::error::Result;
10use crate::raster::Raster;
11
12/// Interpolation method for resampling.
13#[derive(Debug, Clone, Copy, PartialEq, Eq)]
14pub enum ResampleMethod {
15    /// Nearest neighbor — fast, preserves discrete values (classification, SCL)
16    NearestNeighbor,
17    /// Bilinear interpolation — smooth, best for continuous data (elevation, indices)
18    Bilinear,
19}
20
21/// Resample a raster to match the grid of a reference raster.
22///
23/// The output has the same dimensions, transform, and CRS as `reference`.
24/// Each output pixel is computed by mapping its geographic coordinate back
25/// to the source raster and interpolating.
26///
27/// # Arguments
28/// * `source` - Raster to resample
29/// * `reference` - Reference raster defining the target grid
30/// * `method` - Interpolation method
31///
32/// # Example
33/// ```ignore
34/// // Resample Sentinel-2 band (10m) to match DEM grid (30m)
35/// let aligned = resample_to_grid(&s2_band, &dem, ResampleMethod::Bilinear)?;
36/// ```
37pub fn resample_to_grid(
38    source: &Raster<f64>,
39    reference: &Raster<f64>,
40    method: ResampleMethod,
41) -> Result<Raster<f64>> {
42    let (out_rows, out_cols) = reference.shape();
43    let ref_gt = reference.transform();
44    let src_gt = source.transform();
45    let (src_rows, src_cols) = source.shape();
46    let src_data = source.data();
47    let nodata = source.nodata();
48
49    let mut output = Array2::<f64>::from_elem((out_rows, out_cols), f64::NAN);
50
51    for out_row in 0..out_rows {
52        for out_col in 0..out_cols {
53            // Geographic coordinate of output pixel center
54            let geo_x = ref_gt.origin_x + (out_col as f64 + 0.5) * ref_gt.pixel_width;
55            let geo_y = ref_gt.origin_y + (out_row as f64 + 0.5) * ref_gt.pixel_height;
56
57            // Map to source pixel coordinates (continuous)
58            let src_col_f = (geo_x - src_gt.origin_x) / src_gt.pixel_width - 0.5;
59            let src_row_f = (geo_y - src_gt.origin_y) / src_gt.pixel_height - 0.5;
60
61            // Check if source pixel is within bounds
62            if src_col_f < -0.5
63                || src_row_f < -0.5
64                || src_col_f >= src_cols as f64 - 0.5
65                || src_row_f >= src_rows as f64 - 0.5
66            {
67                // Log first few out-of-bounds pixels for debugging
68                if out_row == 0 && out_col < 5 {
69                    eprintln!(
70                        "    [resample] out({},{}) geo=({:.1},{:.1}) src=({:.1},{:.1}) OOB (src_bounds=0..{},0..{})",
71                        out_row, out_col, geo_x, geo_y, src_col_f, src_row_f, src_cols, src_rows
72                    );
73                }
74                continue; // NaN (outside source extent)
75            }
76
77            // Log first few valid pixels
78            if out_row == 0 && (out_col < 3 || out_col == 42 || out_col == 84) {
79                eprintln!(
80                    "    [resample] out({},{}) geo=({:.1},{:.1}) src=({:.1},{:.1}) VALID",
81                    out_row, out_col, geo_x, geo_y, src_col_f, src_row_f
82                );
83            }
84
85            output[[out_row, out_col]] = match method {
86                ResampleMethod::NearestNeighbor => {
87                    let r = src_row_f.round().max(0.0) as usize;
88                    let c = src_col_f.round().max(0.0) as usize;
89                    let r = r.min(src_rows - 1);
90                    let c = c.min(src_cols - 1);
91                    let v = src_data[[r, c]];
92                    if is_nodata(v, nodata) { f64::NAN } else { v }
93                }
94                ResampleMethod::Bilinear => {
95                    let r0 = src_row_f.floor().max(0.0) as usize;
96                    let c0 = src_col_f.floor().max(0.0) as usize;
97                    let r1 = (r0 + 1).min(src_rows - 1);
98                    let c1 = (c0 + 1).min(src_cols - 1);
99
100                    let dr = src_row_f - r0 as f64;
101                    let dc = src_col_f - c0 as f64;
102                    let dr = dr.clamp(0.0, 1.0);
103                    let dc = dc.clamp(0.0, 1.0);
104
105                    let v00 = src_data[[r0, c0]];
106                    let v01 = src_data[[r0, c1]];
107                    let v10 = src_data[[r1, c0]];
108                    let v11 = src_data[[r1, c1]];
109
110                    // Weighted average of valid neighbors (NaN-tolerant bilinear).
111                    // This prevents NaN borders at internal COG tile edges from
112                    // producing stripe artifacts during resampling.
113                    let neighbors = [
114                        (v00, (1.0 - dc) * (1.0 - dr)),
115                        (v01, dc * (1.0 - dr)),
116                        (v10, (1.0 - dc) * dr),
117                        (v11, dc * dr),
118                    ];
119                    let mut wsum = 0.0;
120                    let mut wtotal = 0.0;
121                    for &(v, w) in &neighbors {
122                        if !is_nodata(v, nodata) && v.is_finite() {
123                            wsum += v * w;
124                            wtotal += w;
125                        }
126                    }
127                    if wtotal > 0.0 {
128                        wsum / wtotal
129                    } else {
130                        f64::NAN
131                    }
132                }
133            };
134        }
135    }
136
137    let mut result = Raster::from_array(output);
138    result.set_transform(ref_gt.clone());
139    result.set_nodata(Some(f64::NAN));
140    if let Some(crs) = reference.crs() {
141        result.set_crs(Some(crs.clone()));
142    }
143    Ok(result)
144}
145
146fn is_nodata(v: f64, nodata: Option<f64>) -> bool {
147    if !v.is_finite() {
148        return true;
149    }
150    if let Some(nd) = nodata {
151        if nd.is_finite() && (v - nd).abs() < 1e-10 {
152            return true;
153        }
154    }
155    false
156}
157
158#[cfg(test)]
159mod tests {
160    use super::*;
161    use crate::raster::GeoTransform;
162
163    fn make_raster(rows: usize, cols: usize, value: f64, gt: GeoTransform) -> Raster<f64> {
164        let arr = Array2::from_elem((rows, cols), value);
165        let mut r = Raster::from_array(arr);
166        r.set_transform(gt);
167        r.set_nodata(Some(f64::NAN));
168        r
169    }
170
171    #[test]
172    fn test_resample_same_grid() {
173        // Same grid → output equals input
174        let gt = GeoTransform::new(0.0, 100.0, 10.0, -10.0);
175        let src = make_raster(10, 10, 42.0, gt);
176        let reference = make_raster(10, 10, 0.0, gt);
177
178        let result = resample_to_grid(&src, &reference, ResampleMethod::NearestNeighbor).unwrap();
179        assert!((result.data()[[5, 5]] - 42.0).abs() < 1e-10);
180    }
181
182    #[test]
183    fn test_resample_downsample() {
184        // Source 10m → reference 30m (3x downsample)
185        let src_gt = GeoTransform::new(0.0, 300.0, 10.0, -10.0);
186        let ref_gt = GeoTransform::new(0.0, 300.0, 30.0, -30.0);
187
188        // Source: 30x30 at 10m, all value 100.0
189        let src = make_raster(30, 30, 100.0, src_gt);
190        // Reference: 10x10 at 30m
191        let reference = make_raster(10, 10, 0.0, ref_gt);
192
193        let result = resample_to_grid(&src, &reference, ResampleMethod::Bilinear).unwrap();
194        assert_eq!(result.shape(), (10, 10));
195        // All values should be 100.0 (uniform source)
196        assert!((result.data()[[5, 5]] - 100.0).abs() < 1e-10);
197    }
198
199    #[test]
200    fn test_resample_with_offset() {
201        // Source and reference have different origins (the real-world case)
202        let src_gt = GeoTransform::new(100.0, 500.0, 10.0, -10.0);
203        let ref_gt = GeoTransform::new(115.0, 485.0, 30.0, -30.0); // offset by 15m,15m
204
205        // Source: gradient where value = col index
206        let mut src_data = Array2::<f64>::zeros((40, 40));
207        for r in 0..40 {
208            for c in 0..40 {
209                src_data[[r, c]] = c as f64;
210            }
211        }
212        let mut src = Raster::from_array(src_data);
213        src.set_transform(src_gt);
214        src.set_nodata(Some(f64::NAN));
215
216        let reference = make_raster(10, 10, 0.0, ref_gt);
217
218        let result = resample_to_grid(&src, &reference, ResampleMethod::Bilinear).unwrap();
219        assert_eq!(result.shape(), (10, 10));
220        // First output pixel center is at (115+15, 485-15) = (130, 470)
221        // In source: col = (130 - 100) / 10 - 0.5 = 2.5 → bilinear between col 2 and 3
222        let v = result.data()[[0, 0]];
223        assert!(v.is_finite());
224        assert!((v - 2.5).abs() < 0.5, "Expected ~2.5, got {}", v);
225    }
226
227    #[test]
228    fn test_resample_outside_bounds() {
229        // Reference extends beyond source → NaN for out-of-bounds pixels
230        let src_gt = GeoTransform::new(100.0, 200.0, 10.0, -10.0);
231        let ref_gt = GeoTransform::new(0.0, 300.0, 10.0, -10.0); // extends far beyond source
232
233        let src = make_raster(10, 10, 42.0, src_gt);
234        let reference = make_raster(30, 30, 0.0, ref_gt);
235
236        let result = resample_to_grid(&src, &reference, ResampleMethod::NearestNeighbor).unwrap();
237        // Pixel at (0,0) maps to geo (5, 295), source origin is (100, 200) → outside
238        assert!(result.data()[[0, 0]].is_nan());
239        // Pixel near source area should have value
240        // Source covers x=[100,200], y=[100,200]. ref pixel at col=10,row=10 → geo(105,195)
241        let v = result.data()[[10, 10]];
242        assert!(v.is_finite());
243    }
244}
surtgis_core/resample.rs

surtgis_core/
resample.rs
