Select resampling algo

CHANGES.md

@@ -1,6 +1,14 @@
 # Changes
 
 ## Unreleased
+
+* **Breaking**: Add support to select a resampling algorithm when reading a raster
+    * <https://github.com/georust/gdal/pull/141>
+
+    Now, it is necessary to provide a `Option<ResampleAlg>` when reading a raster.
+    If `None`, it uses `ResampleAlg::NearestNeighbour` which was the
+    default behavior.
+
 * Implement wrapper for `OGR_L_TestCapability`
     * <https://github.com/georust/gdal/pull/160>
 * **Breaking**: Use `DatasetOptions` to pass as `Dataset::open_ex` parameters and
@@ -76,8 +84,11 @@
   `SpatialRef::axis_mapping_strategy` instead.
 * Add support for reading and setting rasterband colour interpretations
     * <https://github.com/georust/gdal/pull/144>
+<<<<<<< HEAD
 * Fixed memory leak in `Geometry::from_wkt`
     * <https://github.com/georust/gdal/pull/172>
+=======
+>>>>>>> Added entry to CHANGES
 
 ## 0.7.1
 * fix docs.rs build for gdal-sys

examples/rasterband.rs

@@ -13,7 +13,7 @@ fn main() {
     println!("rasterband type: {:?}", rasterband.band_type());
     println!("rasterband scale: {:?}", rasterband.scale());
     println!("rasterband offset: {:?}", rasterband.offset());
-    if let Ok(rv) = rasterband.read_as::<u8>((20, 30), (2, 3), (2, 3)) {
+    if let Ok(rv) = rasterband.read_as::<u8>((20, 30), (2, 3), (2, 3), None) {
         println!("{:?}", rv.data);
     }
 }

src/raster/rasterband.rs

@@ -3,7 +3,10 @@ use crate::gdal_major_object::MajorObject;
 use crate::metadata::Metadata;
 use crate::raster::{GDALDataType, GdalType};
 use crate::utils::{_last_cpl_err, _last_null_pointer_err, _string};
-use gdal_sys::{self, CPLErr, GDALColorInterp, GDALMajorObjectH, GDALRWFlag, GDALRasterBandH};
+use gdal_sys::{
+    self, CPLErr, GDALColorInterp, GDALMajorObjectH, GDALRWFlag, GDALRasterBandH,
+    GDALRasterIOExtraArg,
+};
 use libc::c_int;
 use std::ffi::CString;
 
@@ -12,6 +15,100 @@ use ndarray::Array2;
 
 use crate::errors::*;
 
+/// Resampling algorithms, map GDAL defines
+#[derive(Debug, Copy, Clone)]
+pub enum ResampleAlg {
+    /// Nearest neighbour
+    NearestNeighbour,
+    /// Bilinear (2x2 kernel)
+    Bilinear,
+    /// Cubic Convolution Approximation (4x4 kernel)
+    Cubic,
+    /// Cubic B-Spline Approximation (4x4 kernel)
+    CubicSpline,
+    /// Lanczos windowed sinc interpolation (6x6 kernel)
+    Lanczos,
+    /// Average
+    Average,
+    /// Mode (selects the value which appears most often of all the sampled points)
+    Mode,
+    /// Gauss blurring
+    Gauss,
+}
+
+fn map_resample_alg(alg: &ResampleAlg) -> u32 {
+    match alg {
+        ResampleAlg::NearestNeighbour => 0,
+        ResampleAlg::Bilinear => 1,
+        ResampleAlg::Cubic => 2,
+        ResampleAlg::CubicSpline => 3,
+        ResampleAlg::Lanczos => 4,
+        ResampleAlg::Average => 5,
+        ResampleAlg::Mode => 6,
+        ResampleAlg::Gauss => 7,
+    }
+}
+
+/// Extra options used to read a raster.
+///
+/// For documentation, see `gdal_sys::GDALRasterIOExtraArg`.
+#[derive(Debug)]
+pub struct RasterIOExtraArg {
+    pub n_version: usize,
+    pub e_resample_alg: ResampleAlg,
+    pub pfn_progress: gdal_sys::GDALProgressFunc,
+    p_progress_data: *mut libc::c_void,
+    pub b_floating_point_window_validity: usize,
+    pub df_x_off: f64,
+    pub df_y_off: f64,
+    pub df_x_size: f64,
+    pub df_y_size: f64,
+}
+
+impl Default for RasterIOExtraArg {
+    fn default() -> Self {
+        Self {
+            n_version: 1,
+            pfn_progress: None,
+            p_progress_data: std::ptr::null_mut(),
+            e_resample_alg: ResampleAlg::NearestNeighbour,
+            b_floating_point_window_validity: 0,
+            df_x_off: 0.0,
+            df_y_off: 0.0,
+            df_x_size: 0.0,
+            df_y_size: 0.0,
+        }
+    }
+}
+
+impl Into<GDALRasterIOExtraArg> for RasterIOExtraArg {
+    fn into(self) -> GDALRasterIOExtraArg {
+        let RasterIOExtraArg {
+            n_version,
+            e_resample_alg,
+            pfn_progress,
+            p_progress_data,
+            b_floating_point_window_validity,
+            df_x_off,
+            df_y_off,
+            df_x_size,
+            df_y_size,
+        } = self;
+
+        GDALRasterIOExtraArg {
+            nVersion: n_version as c_int,
+            eResampleAlg: map_resample_alg(&e_resample_alg),
+            pfnProgress: pfn_progress,
+            pProgressData: p_progress_data,
+            bFloatingPointWindowValidity: b_floating_point_window_validity as c_int,
+            dfXOff: df_x_off,
+            dfYOff: df_y_off,
+            dfXSize: df_x_size,
+            dfYSize: df_y_size,
+        }
+    }
+}
+
 /// Represents a single band of a dataset.
 ///
 /// This object carries the lifetime of the dataset that
@@ -73,19 +170,30 @@ impl<'a> RasterBand<'a> {
     /// * window_size - the window size (GDAL will interpolate data if window_size != buffer_size)
     /// * size - the desired size to read
     /// * buffer - a slice to hold the data (length must equal product of size parameter)
+    /// * e_resample_alg - the resample algorithm used for the interpolation
     pub fn read_into_slice<T: Copy + GdalType>(
         &self,
         window: (isize, isize),
         window_size: (usize, usize),
         size: (usize, usize),
         buffer: &mut [T],
+        e_resample_alg: Option<ResampleAlg>,
     ) -> Result<()> {
         let pixels = (size.0 * size.1) as usize;
         assert_eq!(buffer.len(), pixels);
 
-        //let no_data:
+        let resample_alg = e_resample_alg.unwrap_or(ResampleAlg::NearestNeighbour);
+
+        let mut options: GDALRasterIOExtraArg = RasterIOExtraArg {
+            e_resample_alg: resample_alg,
+            ..Default::default()
+        }
+        .into();
+
+        let options_ptr: *mut GDALRasterIOExtraArg = &mut options;
+
         let rv = unsafe {
-            gdal_sys::GDALRasterIO(
+            gdal_sys::GDALRasterIOEx(
                 self.c_rasterband,
                 GDALRWFlag::GF_Read,
                 window.0 as c_int,
@@ -98,6 +206,7 @@ impl<'a> RasterBand<'a> {
                 T::gdal_type(),
                 0,
                 0,
+                options_ptr,
             )
         };
         if rv != CPLErr::CE_None {
@@ -113,11 +222,13 @@ impl<'a> RasterBand<'a> {
     /// * window - the window position from top left
     /// * window_size - the window size (GDAL will interpolate data if window_size != buffer_size)
     /// * buffer_size - the desired size of the 'Buffer'
+    /// * e_resample_alg - the resample algorithm used for the interpolation
     pub fn read_as<T: Copy + GdalType>(
         &self,
         window: (isize, isize),
         window_size: (usize, usize),
         size: (usize, usize),
+        e_resample_alg: Option<ResampleAlg>,
     ) -> Result<Buffer<T>> {
         let pixels = (size.0 * size.1) as usize;
 
@@ -131,7 +242,7 @@ impl<'a> RasterBand<'a> {
         unsafe {
             data.set_len(pixels);
         };
-        self.read_into_slice(window, window_size, size, &mut data)?;
+        self.read_into_slice(window, window_size, size, &mut data, e_resample_alg)?;
 
         Ok(Buffer { size, data })
     }
@@ -143,15 +254,17 @@ impl<'a> RasterBand<'a> {
     /// * window - the window position from top left
     /// * window_size - the window size (GDAL will interpolate data if window_size != array_size)
     /// * array_size - the desired size of the 'Array'
+    /// * e_resample_alg - the resample algorithm used for the interpolation
     /// # Docs
     /// The Matrix shape is (rows, cols) and raster shape is (cols in x-axis, rows in y-axis).
     pub fn read_as_array<T: Copy + GdalType>(
         &self,
         window: (isize, isize),
         window_size: (usize, usize),
         array_size: (usize, usize),
+        e_resample_alg: Option<ResampleAlg>,
     ) -> Result<Array2<T>> {
-        let data = self.read_as::<T>(window, window_size, array_size)?;
+        let data = self.read_as::<T>(window, window_size, array_size, e_resample_alg)?;
 
         // Matrix shape is (rows, cols) and raster shape is (cols in x-axis, rows in y-axis)
         Ok(Array2::from_shape_vec(
@@ -169,6 +282,7 @@ impl<'a> RasterBand<'a> {
             (0, 0),
             (size.0 as usize, size.1 as usize),
             (size.0 as usize, size.1 as usize),
+            None,
         )
     }
 

src/raster/tests.rs

@@ -1,5 +1,6 @@
 use crate::dataset::Dataset;
 use crate::metadata::Metadata;
+use crate::raster::rasterband::ResampleAlg;
 use crate::raster::{ByteBuffer, ColorInterpretation};
 use crate::Driver;
 use gdal_sys::GDALDataType;
@@ -73,17 +74,50 @@ fn test_get_projection() {
 fn test_read_raster() {
     let dataset = Dataset::open(fixture!("tinymarble.png")).unwrap();
     let rb = dataset.rasterband(1).unwrap();
-    let rv = rb.read_as::<u8>((20, 30), (2, 3), (2, 3)).unwrap();
+    let rv = rb.read_as::<u8>((20, 30), (2, 3), (2, 3), None).unwrap();
     assert_eq!(rv.size.0, 2);
     assert_eq!(rv.size.1, 3);
     assert_eq!(rv.data, vec!(7, 7, 7, 10, 8, 12));
 
     let mut buf = rv;
-    rb.read_into_slice((20, 30), (2, 3), (2, 3), &mut buf.data)
+    rb.read_into_slice((20, 30), (2, 3), (2, 3), &mut buf.data, None)
         .unwrap();
     assert_eq!(buf.data, vec!(7, 7, 7, 10, 8, 12));
 }
 
+#[test]
+fn test_read_raster_with_default_resample() {
+    let dataset = Dataset::open(fixture!("tinymarble.png")).unwrap();
+    let rb = dataset.rasterband(1).unwrap();
+    let rv = rb.read_as::<u8>((20, 30), (4, 4), (2, 2), None).unwrap();
+    assert_eq!(rv.size.0, 2);
+    assert_eq!(rv.size.1, 2);
+    assert_eq!(rv.data, vec!(10, 4, 6, 11));
+
+    let mut buf = rv;
+    rb.read_into_slice((20, 30), (4, 4), (2, 2), &mut buf.data, None)
+        .unwrap();
+    assert_eq!(buf.data, vec!(10, 4, 6, 11));
+}
+
+#[test]
+fn test_read_raster_with_average_resample() {
+    let dataset = Dataset::open(fixture!("tinymarble.png")).unwrap();
+    let rb = dataset.rasterband(1).unwrap();
+    let resample_alg = ResampleAlg::Average;
+    let rv = rb
+        .read_as::<u8>((20, 30), (4, 4), (2, 2), Some(resample_alg))
+        .unwrap();
+    assert_eq!(rv.size.0, 2);
+    assert_eq!(rv.size.1, 2);
+    assert_eq!(rv.data, vec!(8, 6, 8, 12));
+
+    let mut buf = rv;
+    rb.read_into_slice((20, 30), (4, 4), (2, 2), &mut buf.data, Some(resample_alg))
+        .unwrap();
+    assert_eq!(buf.data, vec!(8, 6, 8, 12));
+}
+
 #[test]
 fn test_write_raster() {
     let driver = Driver::get("MEM").unwrap();
@@ -103,11 +137,11 @@ fn test_write_raster() {
     assert!(res.is_ok());
 
     // read a pixel from the left side
-    let left = rb.read_as::<u8>((5, 5), (1, 1), (1, 1)).unwrap();
+    let left = rb.read_as::<u8>((5, 5), (1, 1), (1, 1), None).unwrap();
     assert_eq!(left.data[0], 50u8);
 
     // read a pixel from the right side
-    let right = rb.read_as::<u8>((15, 5), (1, 1), (1, 1)).unwrap();
+    let right = rb.read_as::<u8>((15, 5), (1, 1), (1, 1), None).unwrap();
     assert_eq!(right.data[0], 20u8);
 }
 
@@ -240,7 +274,7 @@ fn test_get_driver_by_name() {
 fn test_read_raster_as() {
     let dataset = Dataset::open(fixture!("tinymarble.png")).unwrap();
     let rb = dataset.rasterband(1).unwrap();
-    let rv = rb.read_as::<u8>((20, 30), (2, 3), (2, 3)).unwrap();
+    let rv = rb.read_as::<u8>((20, 30), (2, 3), (2, 3), None).unwrap();
     assert_eq!(rv.data, vec!(7, 7, 7, 10, 8, 12));
     assert_eq!(rv.size.0, 2);
     assert_eq!(rv.size.1, 3);
@@ -261,6 +295,7 @@ fn test_read_raster_as_array() {
             (left, top),
             (window_size_x, window_size_y),
             (array_size_x, array_size_y),
+            None,
         )
         .unwrap();