GCP Transformer and create_and_reproject warp.

Author: michaeljfazio-thales

CHANGES.md

@@ -1,6 +1,8 @@
 # Changes
 
 ## Unreleased
+* Added wrapper for Transformer (currently only implementing GCP Transformer)
+* Added wrapper for GDALCreateAndReprojectImage
 * **Breaking**: Upgrade to `ndarray 0.15`
     * <https://github.com/georust/gdal/pull/175>
 * Implement wrapper for `OGR_L_TestCapability`

src/alg/mod.rs

@@ -0,0 +1 @@
+pub mod transform;

src/alg/transform.rs

@@ -0,0 +1,183 @@
+use core::f64;
+use std::{
+    cell::RefCell,
+    ffi::{CStr, CString},
+};
+
+use libc::c_void;
+
+use crate::{
+    errors::{GdalError, Result},
+    Dataset,
+};
+
+/// Ground Control Point (GCP) used to georeference a dataset.
+#[derive(Debug, Clone)]
+pub struct GCP {
+    id: String,
+    info: Option<String>,
+    pixel_xy: (usize, usize),
+    location_xyz: (f64, f64, f64),
+}
+
+impl GCP {
+    /// Returns a GCP constructed from its GDAL C equivalent.
+    ///
+    /// # Arguments
+    ///
+    /// * `c_gcp` - A valid pointer to a C GDAL GCP representation.
+    pub fn with_c_gcp(c_gcp: *const gdal_sys::GDAL_GCP) -> GCP {
+        unsafe {
+            GCP {
+                id: CStr::from_ptr((*c_gcp).pszId).to_str().unwrap().to_string(),
+                info: CStr::from_ptr((*c_gcp).pszInfo)
+                    .to_str()
+                    .map_or(None, |str| Some(str.to_string())),
+                pixel_xy: ((*c_gcp).dfGCPPixel as usize, (*c_gcp).dfGCPLine as usize),
+                location_xyz: ((*c_gcp).dfGCPX, (*c_gcp).dfGCPY, (*c_gcp).dfGCPZ),
+            }
+        }
+    }
+}
+
+/// Polynomial order.
+#[derive(Copy, Clone)]
+pub enum Order {
+    First = 1,
+    Second = 2,
+    Third = 3,
+}
+/// Transformer used to map between pixel/line/height to longitude/latitude/height.
+pub struct Transformer {
+    c_transformer_ref: RefCell<*mut c_void>,
+}
+
+impl Transformer {
+    /// Construct a ```Transformer``` from a valid C GDAL pointer to a transformer instance.
+    pub(crate) unsafe fn with_c_transformer(c_transformer: *mut c_void) -> Transformer {
+        Transformer {
+            c_transformer_ref: RefCell::new(c_transformer),
+        }
+    }
+
+    /// Extracts the Ground Control Points from the specified dataset.
+    ///
+    /// # Arguments
+    ///
+    /// * `dataset` - The `Dataset` to extract Ground Control Points from.
+    ///
+    /// # Examples
+    ///
+    /// ```no_run
+    /// use std::path::Path;
+    /// use gdal::Dataset;
+    /// use gdal::alg::transform::Transformer;
+    ///
+    /// let path = Path::new("/path/to/dataset");
+    /// let dataset = Dataset::open(&path)?;
+    /// let gcps = Transformer::gcps(&dataset);
+    /// ```
+    pub fn gcps(dataset: &Dataset) -> Vec<GCP> {
+        unsafe {
+            let count = gdal_sys::GDALGetGCPCount(dataset.c_dataset()) as usize;
+            let gcps = gdal_sys::GDALGetGCPs(dataset.c_dataset());
+            let mut result = Vec::with_capacity(count);
+            for i in 0..count {
+                let gcp = GCP::with_c_gcp(gcps.add(i));
+                result.push(gcp);
+            }
+            result
+        }
+    }
+
+    /// Constructs a GCP based polynomial `Transformer`.
+    ///
+    /// # Arguments
+    ///
+    /// * `gcp`    - A vector of Ground Control Points to be used as input.
+    /// * `order`  - The requested polynomial order. Using a third-order polynomial is not recommended due
+    ///              to the potential for numeric instabilities.
+    /// * `reversed` - Set it to `true` to compute the reversed transformation.
+
+    pub fn gcp(gcp: Vec<GCP>, order: Order, reversed: bool) -> Result<Transformer> {
+        let pas_gcp_list = gcp
+            .iter()
+            .map(|p| {
+                let psz_info = match &p.info {
+                    Some(arg) => {
+                        let tmp = CString::new(arg.clone()).unwrap();
+                        tmp.into_raw()
+                    }
+                    None => std::ptr::null_mut(),
+                };
+
+                gdal_sys::GDAL_GCP {
+                    pszId: CString::new(p.id.clone()).unwrap().into_raw(),
+                    pszInfo: psz_info,
+                    dfGCPPixel: p.pixel_xy.0 as f64,
+                    dfGCPLine: p.pixel_xy.1 as f64,
+                    dfGCPX: p.location_xyz.0,
+                    dfGCPY: p.location_xyz.1,
+                    dfGCPZ: p.location_xyz.2,
+                }
+            })
+            .collect::<Vec<_>>();
+
+        let c_transformer = unsafe {
+            gdal_sys::GDALCreateGCPTransformer(
+                gcp.len() as i32,
+                pas_gcp_list.as_ptr(),
+                order as i32,
+                if reversed { 1 } else { 0 },
+            )
+        };
+
+        if c_transformer.is_null() {
+            return Err(GdalError::NullPointer {
+                method_name: "GDALCreateGCPTransformer",
+                msg: "Failed to create GCP Transformer".to_string(),
+            });
+        }
+
+        Ok(unsafe { Transformer::with_c_transformer(c_transformer) })
+    }
+
+    /// Transform a 2D point based on GCP derived polynomial model.
+    ///
+    /// # Arguments
+    ///
+    /// * `x` - The x point (pixel or longitude).
+    /// * `y` - The y point (pixel or longitude).
+    pub fn transform(&self, x: f64, y: f64) -> Option<(f64, f64)> {
+        let mut x: [f64; 1] = [x];
+        let mut y: [f64; 1] = [y];
+        let mut z: [f64; 1] = [0.];
+        let mut r: [i32; 1] = [0];
+        unsafe {
+            gdal_sys::GDALGCPTransform(
+                *self.c_transformer_ref.borrow(),
+                0,
+                1,
+                x.as_mut_ptr(),
+                y.as_mut_ptr(),
+                z.as_mut_ptr(),
+                r.as_mut_ptr(),
+            );
+            if r[0] != 0 {
+                Some((x[0], y[0]))
+            } else {
+                None
+            }
+        }
+    }
+}
+
+/// Cleanup the GCP transformer when it exits scope.
+impl Drop for Transformer {
+    fn drop(&mut self) {
+        unsafe {
+            let c_transformer = self.c_transformer_ref.borrow();
+            gdal_sys::GDALDestroyGCPTransformer(c_transformer.to_owned());
+        }
+    }
+}

src/lib.rs

@@ -22,6 +22,7 @@
 
 pub use version::version_info;
 
+pub mod alg;
 pub mod config;
 mod dataset;
 mod driver;

src/raster/mod.rs

@@ -6,7 +6,7 @@ mod warp;
 
 pub use rasterband::{Buffer, ByteBuffer, ColorInterpretation, RasterBand};
 pub use types::{GDALDataType, GdalType};
-pub use warp::reproject;
+pub use warp::{create_and_reproject, reproject};
 
 #[cfg(test)]
 mod tests;

src/raster/warp.rs

@@ -1,7 +1,11 @@
-use crate::dataset::Dataset;
 use crate::utils::_last_cpl_err;
+use crate::{dataset::Dataset, spatial_ref::SpatialRef, Driver};
 use gdal_sys::{self, CPLErr, GDALResampleAlg};
-use std::ptr::{null, null_mut};
+use std::{
+    ffi::CString,
+    path::Path,
+    ptr::{null, null_mut},
+};
 
 use crate::errors::*;
 
@@ -25,3 +29,36 @@ pub fn reproject(src: &Dataset, dst: &Dataset) -> Result<()> {
     }
     Ok(())
 }
+
+pub fn create_and_reproject(
+    src: &Dataset,
+    dst_path: &Path,
+    dst_srs: &SpatialRef,
+    dst_driver: &Driver,
+) -> Result<()> {
+    let psz_dst_filename = dst_path
+        .to_str()
+        .expect("Destination path must be supplied.");
+    let psz_dst_wkt = dst_srs.to_wkt().expect("Failed to obtain WKT for SRS.");
+
+    let rv = unsafe {
+        gdal_sys::GDALCreateAndReprojectImage(
+            src.c_dataset(),
+            null(),
+            CString::new(psz_dst_filename)?.as_ptr(),
+            CString::new(psz_dst_wkt)?.as_ptr(),
+            dst_driver.c_driver(),
+            null_mut(),
+            GDALResampleAlg::GRA_Bilinear,
+            0.0,
+            0.0,
+            None,
+            null_mut(),
+            null_mut(),
+        )
+    };
+    if rv != CPLErr::CE_None {
+        return Err(_last_cpl_err(rv));
+    }
+    Ok(())
+}