From dbc4b7e3875ac6bd4059d1627eae4101b97ca92f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9=20Rocha?= Date: Sun, 10 Mar 2019 16:53:37 +0100 Subject: [PATCH 1/4] Add from RasterBand to Ndarray. --- Cargo.toml | 1 + src/errors.rs | 10 ++++- src/lib.rs | 3 +- src/raster/dataset.rs | 29 +++++++++++- src/raster/rasterband.rs | 97 ++++++++++++++++++++++++++++++++++++++++ src/raster/tests.rs | 81 ++++++++++++++++++++++++++++++++- 6 files changed, 216 insertions(+), 5 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index f932a997..67d46d40 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -21,6 +21,7 @@ geo-types = "0.3" # gdal-sys = { path = "gdal-sys", version = "0.2"} gdal-sys = "0.2" num-traits = "0.2" +ndarray = "0.12.1" [workspace] members = ["gdal-sys"] diff --git a/src/errors.rs b/src/errors.rs index 77dfbcd2..3c0bc4a0 100644 --- a/src/errors.rs +++ b/src/errors.rs @@ -11,13 +11,15 @@ pub struct Error { inner: Context, } -#[derive(Clone, Eq, PartialEq, Debug, Fail)] +#[derive(Clone, PartialEq, Debug, Fail)] pub enum ErrorKind { #[fail(display = "FfiNulError")] FfiNulError(#[cause] std::ffi::NulError), #[fail(display = "StrUtf8Error")] StrUtf8Error(#[cause] std::str::Utf8Error), + #[fail(display = "NdarrayShapeError")] + NdarrayShapeError(#[cause] ndarray::ShapeError), #[fail(display = "CPL error class: '{:?}', error number: '{}', error msg: '{}'", class, number, msg)] CplError { class: CPLErr::Type, @@ -109,3 +111,9 @@ impl From for Error { Error { inner: Context::new(ErrorKind::StrUtf8Error(err)) } } } + +impl From for Error { + fn from(err: ndarray::ShapeError) -> Error { + Error { inner: Context::new(ErrorKind::NdarrayShapeError(err)) } + } +} \ No newline at end of file diff --git a/src/lib.rs b/src/lib.rs index 016976ea..13b51d91 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -21,12 +21,11 @@ #![crate_type = "lib"] extern crate failure; -#[macro_use] -extern crate failure_derive; extern crate gdal_sys; extern crate geo_types; extern crate libc; extern crate num_traits; +extern crate ndarray; pub use version::version_info; diff --git a/src/raster/dataset.rs b/src/raster/dataset.rs index a424fda2..68be1582 100644 --- a/src/raster/dataset.rs +++ b/src/raster/dataset.rs @@ -9,6 +9,7 @@ use raster::types::GdalType; use gdal_major_object::MajorObject; use metadata::Metadata; use gdal_sys::{self, CPLErr, GDALAccess, GDALDatasetH, GDALDataType, GDALMajorObjectH}; +use ndarray::Array2; use errors::*; @@ -32,7 +33,6 @@ impl Drop for Dataset { } } - impl Dataset { pub fn open(path: &Path) -> Result { _register_drivers(); @@ -69,6 +69,16 @@ impl Dataset { (size_x, size_y) } + /// Get block size from a 'Dataset'. + /// # Arguments + /// * band_index - the band_index + /* + pub fn size_block(&self, band_index: isize) -> (usize, usize) { + let band = self.rasterband(band_index)?; + band.size_block() + } + */ + pub fn driver(&self) -> Driver { unsafe { let c_driver = gdal_sys::GDALGetDatasetDriver(self.c_dataset); @@ -192,6 +202,23 @@ impl Dataset { self.rasterband(band_index)?.read_as(window, window_size, size) } + /// Read a 'Array2' from a 'Dataset'. T implements 'GdalType'. + /// # Arguments + /// * band_index - the band_index + /// * 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' + pub fn read_as_array( + &self, + band_index: isize, + window: (isize, isize), + window_size: (usize, usize), + array_size: (usize, usize), + ) -> Result> + { + self.rasterband(band_index)?.read_as_array(window, window_size, array_size) + } + /// Write a 'Buffer' into a 'Dataset'. /// # Arguments /// * band_index - the band_index diff --git a/src/raster/rasterband.rs b/src/raster/rasterband.rs index 703429dc..0553d917 100644 --- a/src/raster/rasterband.rs +++ b/src/raster/rasterband.rs @@ -4,6 +4,7 @@ use raster::types::GdalType; use gdal_major_object::MajorObject; use metadata::Metadata; use gdal_sys::{self, CPLErr, GDALDataType, GDALMajorObjectH, GDALRasterBandH, GDALRWFlag}; +use ndarray::{Array2}; use utils::_last_cpl_err; use errors::*; @@ -22,6 +23,23 @@ impl <'a> RasterBand<'a> { RasterBand { c_rasterband, owning_dataset } } + /// Get block size from a 'Dataset'. + /// # Arguments + /// * band_index - the band_index + pub fn block_size(&self) -> (usize, usize) { + let mut size_x = 0; + let mut size_y = 0; + + unsafe { + gdal_sys::GDALGetBlockSize( + self.c_rasterband, + &mut size_x, + &mut size_y + ) + }; + (size_x as usize, size_y as usize) + } + /// Read a 'Buffer' from a 'Dataset'. T implements 'GdalType' /// # Arguments /// * band_index - the band_index @@ -65,6 +83,51 @@ impl <'a> RasterBand<'a> { Ok(Buffer{size, data}) } + /// Read a 'Array2' from a 'Dataset'. T implements 'GdalType'. + /// # Arguments + /// * 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' + /// # Docs + /// The Matrix shape is (rows, cols) and raster shape is (cols in x-axis, rows in y-axis). + pub fn read_as_array( + &self, + window: (isize, isize), + window_size: (usize, usize), + array_size: (usize, usize), + ) -> Result> + { + let pixels = (array_size.0 * array_size.1) as usize; + let mut data: Vec = Vec::with_capacity(pixels); + + let values = unsafe { + gdal_sys::GDALRasterIO( + self.c_rasterband, + GDALRWFlag::GF_Read, + window.0 as c_int, + window.1 as c_int, + window_size.0 as c_int, + window_size.1 as c_int, + data.as_mut_ptr() as GDALRasterBandH, + array_size.0 as c_int, + array_size.1 as c_int, + T::gdal_type(), + 0, + 0 + ) + }; + if values != CPLErr::CE_None { + Err(_last_cpl_err(values))?; + } + + unsafe { + data.set_len(pixels); + }; + + // Matrix shape is (rows, cols) and raster shape is (cols in x-axis, rows in y-axis) + Array2::from_shape_vec((array_size.1, array_size.0) , data).map_err(Into::into) + } + /// Read a full 'Dataset' as 'Buffer'. /// # Arguments /// * band_index - the band_index @@ -80,6 +143,40 @@ impl <'a> RasterBand<'a> { ) } + /// Read a 'Array2' from a 'Dataset' block. T implements 'GdalType' + /// # Arguments + /// * block_index - the block index + /// # Docs + /// The Matrix shape is (rows, cols) and raster shape is (cols in x-axis, rows in y-axis). + pub fn read_block( + &self, + block_index: (usize, usize) + ) -> Result> + { + let size = self.block_size(); + let pixels = (size.0 * size.1) as usize; + let mut data: Vec = Vec::with_capacity(pixels); + + //let no_data: + let rv = unsafe { + gdal_sys::GDALReadBlock( + self.c_rasterband, + block_index.0 as c_int, + block_index.1 as c_int, + data.as_mut_ptr() as GDALRasterBandH + ) + }; + if rv != CPLErr::CE_None { + Err(_last_cpl_err(rv))?; + } + + unsafe { + data.set_len(pixels); + }; + + Array2::from_shape_vec((size.1, size.0), data).map_err(Into::into) + } + // Write a 'Buffer' into a 'Dataset'. /// # Arguments /// * band_index - the band_index diff --git a/src/raster/tests.rs b/src/raster/tests.rs index 32cb6c7d..d8bef946 100644 --- a/src/raster/tests.rs +++ b/src/raster/tests.rs @@ -2,6 +2,7 @@ use std::path::Path; use super::{ByteBuffer, Driver, Dataset}; use gdal_sys::GDALDataType; use metadata::Metadata; +use ndarray::arr2; macro_rules! fixture { @@ -34,6 +35,15 @@ fn test_get_raster_size() { assert_eq!(size_y, 50); } +#[test] +fn test_get_raster_block_size() { + let band_index = 1; + let dataset = Dataset::open(fixture!("tinymarble.png")).unwrap(); + let rasterband = dataset.rasterband(band_index).unwrap(); + let (size_x, size_y) = rasterband.block_size(); + assert_eq!(size_x, 100); + assert_eq!(size_y, 1); +} #[test] fn test_get_raster_count() { @@ -189,7 +199,6 @@ fn test_geo_transform() { assert_eq!(dataset.geo_transform().unwrap(), transform); } - #[test] fn test_get_driver_by_name() { let missing_driver = Driver::get("wtf"); @@ -202,6 +211,7 @@ fn test_get_driver_by_name() { assert_eq!(driver.long_name(), "GeoTIFF"); } + #[test] fn test_read_raster_as() { let dataset = Dataset::open(fixture!("tinymarble.png")).unwrap(); @@ -217,6 +227,29 @@ fn test_read_raster_as() { assert_eq!(dataset.band_type(1).unwrap(), GDALDataType::GDT_Byte); } +#[test] +fn test_read_raster_as_array() { + let band_index = 1; + let (left, top) = (19, 5); + let (window_size_x, window_size_y) = (3, 4); + let (array_size_x, array_size_y) = (3, 4); + let dataset = Dataset::open(fixture!("tinymarble.png")).unwrap(); + let values = dataset.read_as_array::( + band_index, + (left, top), + (window_size_x, window_size_y), + (array_size_x, array_size_y) + ).unwrap(); + + let data = arr2(&[[226, 225, 157], + [215, 222, 225], + [213, 231, 229], + [171, 189, 192]]); + + assert_eq!(values, data); + assert_eq!(dataset.band_type(band_index).unwrap(), GDALDataType::GDT_Byte); +} + #[test] fn test_read_full_raster_as() { let dataset = Dataset::open(fixture!("tinymarble.png")).unwrap(); @@ -225,6 +258,52 @@ fn test_read_full_raster_as() { assert_eq!(rv.size.1, 50); } +#[test] +fn test_read_block_as_array() { + let band_index = 1; + let block_index = (0, 0); + let dataset = Dataset::open(fixture!("tinymarble.png")).unwrap(); + let rasterband = dataset.rasterband(band_index).unwrap(); + let result = rasterband.read_block::(block_index); + assert!(result.is_ok()); +} + +#[test] +fn test_read_block_dimension() { + let band_index = 1; + let block = (0, 0); + let dataset = Dataset::open(fixture!("tinymarble.png")).unwrap(); + let rasterband = dataset.rasterband(band_index).unwrap(); + let array = rasterband.read_block::(block).unwrap(); + let dimension = (1, 100); + assert_eq!(array.dim(), dimension); +} + +#[test] +fn test_read_block_last_dimension() { + let band_index = 1; + let block = (0, 49); + let dataset = Dataset::open(fixture!("tinymarble.png")).unwrap(); + let rasterband = dataset.rasterband(band_index).unwrap(); + let array = rasterband.read_block::(block).unwrap(); + let dimension = (1, 100); + assert_eq!(array.dim(), dimension); +} + +#[test] +fn test_read_block_data() { + let band_index = 1; + let block = (0, 0); + let dataset = Dataset::open(fixture!("tinymarble.png")).unwrap(); + let rasterband = dataset.rasterband(band_index).unwrap(); + let array = rasterband.read_block::(block).unwrap(); + assert_eq!(array[[0, 0]], 0); + assert_eq!(array[[0, 1]], 9); + assert_eq!(array[[0, 98]], 24); + assert_eq!(array[[0, 99]], 51); +} + + #[test] fn test_get_band_type() { let driver = Driver::get("MEM").unwrap(); From 06985f857bf40990614d5f302f6d945c00b749a7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9=20Rocha?= Date: Sun, 10 Mar 2019 17:00:06 +0100 Subject: [PATCH 2/4] Reinsert failure_derive. --- src/lib.rs | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/lib.rs b/src/lib.rs index 13b51d91..c1b9cf6d 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -21,6 +21,8 @@ #![crate_type = "lib"] extern crate failure; +#[macro_use] +extern crate failure_derive; extern crate gdal_sys; extern crate geo_types; extern crate libc; From c68ab7a40be010d0de0d2ad85a1f73da15cb2f85 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9=20Rocha?= Date: Thu, 14 Mar 2019 22:06:02 +0100 Subject: [PATCH 3/4] Ndarray placed inside a feature. --- Cargo.toml | 3 ++- src/errors.rs | 2 ++ src/lib.rs | 2 ++ src/raster/dataset.rs | 36 ++++++++++++++++++++++++++++++------ src/raster/rasterband.rs | 6 +++++- src/raster/tests.rs | 33 ++++++++++++++++++++++++++++++++- 6 files changed, 73 insertions(+), 9 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index 67d46d40..4761344c 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -12,6 +12,7 @@ documentation = "https://georust.github.io/gdal/" [features] bindgen = ["gdal-sys/bindgen"] +array = ["ndarray"] [dependencies] failure = "0.1" @@ -21,7 +22,7 @@ geo-types = "0.3" # gdal-sys = { path = "gdal-sys", version = "0.2"} gdal-sys = "0.2" num-traits = "0.2" -ndarray = "0.12.1" +ndarray = {version = "0.12.1", optional = true } [workspace] members = ["gdal-sys"] diff --git a/src/errors.rs b/src/errors.rs index 3c0bc4a0..078bc7cf 100644 --- a/src/errors.rs +++ b/src/errors.rs @@ -18,6 +18,7 @@ pub enum ErrorKind { FfiNulError(#[cause] std::ffi::NulError), #[fail(display = "StrUtf8Error")] StrUtf8Error(#[cause] std::str::Utf8Error), + #[cfg(feature = "ndarray")] #[fail(display = "NdarrayShapeError")] NdarrayShapeError(#[cause] ndarray::ShapeError), #[fail(display = "CPL error class: '{:?}', error number: '{}', error msg: '{}'", class, number, msg)] @@ -112,6 +113,7 @@ impl From for Error { } } +#[cfg(feature = "ndarray")] impl From for Error { fn from(err: ndarray::ShapeError) -> Error { Error { inner: Context::new(ErrorKind::NdarrayShapeError(err)) } diff --git a/src/lib.rs b/src/lib.rs index c1b9cf6d..1a76c2c2 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -27,6 +27,8 @@ extern crate gdal_sys; extern crate geo_types; extern crate libc; extern crate num_traits; + +#[cfg(feature = "ndarray")] extern crate ndarray; pub use version::version_info; diff --git a/src/raster/dataset.rs b/src/raster/dataset.rs index 68be1582..c67d3aa5 100644 --- a/src/raster/dataset.rs +++ b/src/raster/dataset.rs @@ -9,6 +9,8 @@ use raster::types::GdalType; use gdal_major_object::MajorObject; use metadata::Metadata; use gdal_sys::{self, CPLErr, GDALAccess, GDALDatasetH, GDALDataType, GDALMajorObjectH}; + +#[cfg(feature = "ndarray")] use ndarray::Array2; use errors::*; @@ -101,10 +103,23 @@ impl Dataset { Ok(()) } - pub fn set_geo_transform(&self, tr: &GeoTransform) -> Result<()> { - assert_eq!(tr.len(), 6); + /// Affine transformation called geotransformation. + /// + /// This is like a linear transformation preserves points, straight lines and planes. + /// Also, sets of parallel lines remain parallel after an affine transformation. + /// # Arguments + /// * transformation - coeficients of transformations + /// + /// x-coordinate of the top-left corner pixel (x-offset) + /// width of a pixel (x-resolution) + /// row rotation (typically zero) + /// y-coordinate of the top-left corner pixel + /// column rotation (typically zero) + /// height of a pixel (y-resolution, typically negative) + pub fn set_geo_transform(&self, transformation: &GeoTransform) -> Result<()> { + assert_eq!(transformation.len(), 6); let rv = unsafe { - gdal_sys::GDALSetGeoTransform(self.c_dataset, tr.as_ptr() as *mut f64) + gdal_sys::GDALSetGeoTransform(self.c_dataset, transformation.as_ptr() as *mut f64) }; if rv != CPLErr::CE_None { Err(_last_cpl_err(rv))?; @@ -112,12 +127,20 @@ impl Dataset { Ok(()) } + /// Get affine transformation coefficients. + /// + /// x-coordinate of the top-left corner pixel (x-offset) + /// width of a pixel (x-resolution) + /// row rotation (typically zero) + /// y-coordinate of the top-left corner pixel + /// column rotation (typically zero) + /// height of a pixel (y-resolution, typically negative) pub fn geo_transform(&self) -> Result { - let mut tr = GeoTransform::default(); + let mut transformation = GeoTransform::default(); let rv = unsafe { gdal_sys::GDALGetGeoTransform( self.c_dataset, - tr.as_mut_ptr() + transformation.as_mut_ptr() ) }; @@ -125,7 +148,7 @@ impl Dataset { if rv != CPLErr::CE_None { Err(_last_cpl_err(rv))?; } - Ok(tr) + Ok(transformation) } pub fn create_copy( @@ -202,6 +225,7 @@ impl Dataset { self.rasterband(band_index)?.read_as(window, window_size, size) } + #[cfg(feature = "ndarray")] /// Read a 'Array2' from a 'Dataset'. T implements 'GdalType'. /// # Arguments /// * band_index - the band_index diff --git a/src/raster/rasterband.rs b/src/raster/rasterband.rs index 0553d917..092bb4a6 100644 --- a/src/raster/rasterband.rs +++ b/src/raster/rasterband.rs @@ -4,9 +4,11 @@ use raster::types::GdalType; use gdal_major_object::MajorObject; use metadata::Metadata; use gdal_sys::{self, CPLErr, GDALDataType, GDALMajorObjectH, GDALRasterBandH, GDALRWFlag}; -use ndarray::{Array2}; use utils::_last_cpl_err; +#[cfg(feature = "ndarray")] +use ndarray::{Array2}; + use errors::*; pub struct RasterBand<'a> { @@ -83,6 +85,7 @@ impl <'a> RasterBand<'a> { Ok(Buffer{size, data}) } + #[cfg(feature = "ndarray")] /// Read a 'Array2' from a 'Dataset'. T implements 'GdalType'. /// # Arguments /// * window - the window position from top left @@ -143,6 +146,7 @@ impl <'a> RasterBand<'a> { ) } + #[cfg(feature = "ndarray")] /// Read a 'Array2' from a 'Dataset' block. T implements 'GdalType' /// # Arguments /// * block_index - the block index diff --git a/src/raster/tests.rs b/src/raster/tests.rs index d8bef946..576971d5 100644 --- a/src/raster/tests.rs +++ b/src/raster/tests.rs @@ -1,7 +1,9 @@ use std::path::Path; -use super::{ByteBuffer, Driver, Dataset}; +use super::{ByteBuffer, Driver, Dataset, warp}; use gdal_sys::GDALDataType; use metadata::Metadata; + +#[cfg(feature = "ndarray")] use ndarray::arr2; @@ -228,6 +230,7 @@ fn test_read_raster_as() { } #[test] +#[cfg(feature = "ndarray")] fn test_read_raster_as_array() { let band_index = 1; let (left, top) = (19, 5); @@ -259,6 +262,7 @@ fn test_read_full_raster_as() { } #[test] +#[cfg(feature = "ndarray")] fn test_read_block_as_array() { let band_index = 1; let block_index = (0, 0); @@ -269,6 +273,7 @@ fn test_read_block_as_array() { } #[test] +#[cfg(feature = "ndarray")] fn test_read_block_dimension() { let band_index = 1; let block = (0, 0); @@ -280,6 +285,7 @@ fn test_read_block_dimension() { } #[test] +#[cfg(feature = "ndarray")] fn test_read_block_last_dimension() { let band_index = 1; let block = (0, 49); @@ -291,6 +297,7 @@ fn test_read_block_last_dimension() { } #[test] +#[cfg(feature = "ndarray")] fn test_read_block_data() { let band_index = 1; let block = (0, 0); @@ -358,3 +365,27 @@ fn test_get_offset() { let offset = rasterband.offset(); assert_eq!(offset, Some(0.0)); } + +#[test] +fn test_reproject() { + let driver = Driver::get("MEM").unwrap(); + let size_x = 100; + let size_y = 50; + let bands = 3; + let filename = ""; // Don't need concrete filename because is a file in memory + let destiny = driver.create(filename, size_x, size_y, bands).unwrap(); + let source = Dataset::open(fixture!("tinybluemarble.tif")).unwrap(); + + println!("{:?}", source.geo_transform().expect("Geotransform vector.")); + + // c = x-coordinate of the upper-left corner of the upper-left pixel + // a = width of a pixel + // b = row rotation (typically zero) + // f = y-coordinate of the of the upper-left corner of the upper-left pixel + // d = column rotation (typically zero) + // e = height of a pixel (typically negative) + let transform = [0.0, 10.0, 0.0, 0.0, 0.0, -10.0]; + assert!(destiny.set_geo_transform(&transform).is_ok()); + assert!(warp::reproject(&source, &destiny).is_ok()); + println!("{:?}", destiny.geo_transform().expect("Geotransform vector.")); +} From 52468129ca8c1dfc43d0727acccd4aa6aa61ab1d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9=20Rocha?= Date: Thu, 14 Mar 2019 22:17:45 +0100 Subject: [PATCH 4/4] Without reproject test. --- src/raster/tests.rs | 24 ------------------------ 1 file changed, 24 deletions(-) diff --git a/src/raster/tests.rs b/src/raster/tests.rs index 576971d5..a502c091 100644 --- a/src/raster/tests.rs +++ b/src/raster/tests.rs @@ -365,27 +365,3 @@ fn test_get_offset() { let offset = rasterband.offset(); assert_eq!(offset, Some(0.0)); } - -#[test] -fn test_reproject() { - let driver = Driver::get("MEM").unwrap(); - let size_x = 100; - let size_y = 50; - let bands = 3; - let filename = ""; // Don't need concrete filename because is a file in memory - let destiny = driver.create(filename, size_x, size_y, bands).unwrap(); - let source = Dataset::open(fixture!("tinybluemarble.tif")).unwrap(); - - println!("{:?}", source.geo_transform().expect("Geotransform vector.")); - - // c = x-coordinate of the upper-left corner of the upper-left pixel - // a = width of a pixel - // b = row rotation (typically zero) - // f = y-coordinate of the of the upper-left corner of the upper-left pixel - // d = column rotation (typically zero) - // e = height of a pixel (typically negative) - let transform = [0.0, 10.0, 0.0, 0.0, 0.0, -10.0]; - assert!(destiny.set_geo_transform(&transform).is_ok()); - assert!(warp::reproject(&source, &destiny).is_ok()); - println!("{:?}", destiny.geo_transform().expect("Geotransform vector.")); -}