diff --git a/CHANGES.md b/CHANGES.md index a0a27bc3b..2935d85d3 100644 --- a/CHANGES.md +++ b/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` * * Implement wrapper for `OGR_L_TestCapability` diff --git a/src/alg/mod.rs b/src/alg/mod.rs new file mode 100644 index 000000000..1bba7e150 --- /dev/null +++ b/src/alg/mod.rs @@ -0,0 +1 @@ +pub mod transform; diff --git a/src/alg/transform.rs b/src/alg/transform.rs new file mode 100644 index 000000000..513a9ac2a --- /dev/null +++ b/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, + 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 { + 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, order: Order, reversed: bool) -> Result { + 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::>(); + + 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()); + } + } +} diff --git a/src/lib.rs b/src/lib.rs index 11017979f..359c55e88 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -22,6 +22,7 @@ pub use version::version_info; +pub mod alg; pub mod config; mod dataset; mod driver; diff --git a/src/raster/mod.rs b/src/raster/mod.rs index 8eba798ca..b2d97931d 100644 --- a/src/raster/mod.rs +++ b/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; diff --git a/src/raster/warp.rs b/src/raster/warp.rs index e803fb360..249745139 100644 --- a/src/raster/warp.rs +++ b/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(()) +}