diff --git a/Cargo.toml b/Cargo.toml
index f932a997..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,6 +22,7 @@ geo-types = "0.3"
 # gdal-sys = { path = "gdal-sys", version = "0.2"}
 gdal-sys = "0.2"
 num-traits = "0.2"
+ndarray = {version = "0.12.1", optional = true }
 
 [workspace]
 members = ["gdal-sys"]
diff --git a/src/errors.rs b/src/errors.rs
index 77dfbcd2..078bc7cf 100644
--- a/src/errors.rs
+++ b/src/errors.rs
@@ -11,13 +11,16 @@ pub struct Error {
     inner: Context<ErrorKind>,
 }
 
-#[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),
+    #[cfg(feature = "ndarray")]
+    #[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 +112,10 @@ impl From<std::str::Utf8Error> for Error {
         Error { inner: Context::new(ErrorKind::StrUtf8Error(err)) }
     }
 }
+
+#[cfg(feature = "ndarray")]
+impl From<ndarray::ShapeError> 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..1a76c2c2 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -28,6 +28,9 @@ extern crate geo_types;
 extern crate libc;
 extern crate num_traits;
 
+#[cfg(feature = "ndarray")]
+extern crate ndarray;
+
 pub use version::version_info;
 
 pub mod config;
diff --git a/src/raster/dataset.rs b/src/raster/dataset.rs
index a424fda2..c67d3aa5 100644
--- a/src/raster/dataset.rs
+++ b/src/raster/dataset.rs
@@ -10,6 +10,9 @@ 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::*;
 
 pub type GeoTransform = [c_double; 6];
@@ -32,7 +35,6 @@ impl Drop for Dataset {
     }
 }
 
-
 impl Dataset {
     pub fn open(path: &Path) -> Result<Dataset> {
         _register_drivers();
@@ -69,6 +71,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);
@@ -91,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))?;
@@ -102,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<GeoTransform> {
-        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()
             )
         };
 
@@ -115,7 +148,7 @@ impl Dataset {
         if rv != CPLErr::CE_None {
             Err(_last_cpl_err(rv))?;
         }
-        Ok(tr)
+        Ok(transformation)
     }
 
     pub fn create_copy(
@@ -192,6 +225,24 @@ impl Dataset {
         self.rasterband(band_index)?.read_as(window, window_size, size)
     }
 
+    #[cfg(feature = "ndarray")]
+    /// Read a 'Array2<T>' 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<T: Copy + GdalType>(
+        &self,
+        band_index: isize,
+        window: (isize, isize),
+        window_size: (usize, usize),
+        array_size: (usize, usize),
+    ) -> Result<Array2<T>>
+    {
+        self.rasterband(band_index)?.read_as_array(window, window_size, array_size)
+    }
+
     /// Write a 'Buffer<T>' into a 'Dataset'.
     /// # Arguments
     /// * band_index - the band_index
diff --git a/src/raster/rasterband.rs b/src/raster/rasterband.rs
index 703429dc..092bb4a6 100644
--- a/src/raster/rasterband.rs
+++ b/src/raster/rasterband.rs
@@ -6,6 +6,9 @@ use metadata::Metadata;
 use gdal_sys::{self, CPLErr, GDALDataType, GDALMajorObjectH, GDALRasterBandH, GDALRWFlag};
 use utils::_last_cpl_err;
 
+#[cfg(feature = "ndarray")]
+use ndarray::{Array2};
+
 use errors::*;
 
 pub struct RasterBand<'a> {
@@ -22,6 +25,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<T>' from a 'Dataset'. T implements 'GdalType'
     /// # Arguments
     /// * band_index - the band_index
@@ -65,6 +85,52 @@ impl <'a> RasterBand<'a> {
         Ok(Buffer{size, data})
     }
 
+    #[cfg(feature = "ndarray")]
+    /// Read a 'Array2<T>' 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<T: Copy + GdalType>(
+        &self,
+        window: (isize, isize),
+        window_size: (usize, usize),
+        array_size: (usize, usize),
+    ) -> Result<Array2<T>>
+    {
+        let pixels = (array_size.0 * array_size.1) as usize;
+        let mut data: Vec<T> = 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<T>'.
     /// # Arguments
     /// * band_index - the band_index
@@ -80,6 +146,41 @@ impl <'a> RasterBand<'a> {
         )
     }
 
+    #[cfg(feature = "ndarray")]
+    /// Read a 'Array2<T>' 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<T: Copy + GdalType>(
+        &self,
+        block_index: (usize, usize)
+    ) -> Result<Array2<T>>
+    {
+        let size = self.block_size();
+        let pixels = (size.0 * size.1) as usize;
+        let mut data: Vec<T> = 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<T>' into a 'Dataset'.
     /// # Arguments
     /// * band_index - the band_index
diff --git a/src/raster/tests.rs b/src/raster/tests.rs
index 32cb6c7d..a502c091 100644
--- a/src/raster/tests.rs
+++ b/src/raster/tests.rs
@@ -1,8 +1,11 @@
 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;
+
 
 macro_rules! fixture {
     ($name:expr) => (
@@ -34,6 +37,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 +201,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 +213,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 +229,30 @@ fn test_read_raster_as() {
     assert_eq!(dataset.band_type(1).unwrap(), GDALDataType::GDT_Byte);
 }
 
+#[test]
+#[cfg(feature = "ndarray")]
+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::<u8>(
+        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 +261,56 @@ fn test_read_full_raster_as() {
     assert_eq!(rv.size.1, 50);
 }
 
+#[test]
+#[cfg(feature = "ndarray")]
+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::<u8>(block_index);
+    assert!(result.is_ok());
+}
+
+#[test]
+#[cfg(feature = "ndarray")]
+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::<u8>(block).unwrap();
+    let dimension = (1, 100);
+    assert_eq!(array.dim(), dimension);
+}
+
+#[test]
+#[cfg(feature = "ndarray")]
+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::<u8>(block).unwrap();
+    let dimension = (1, 100);
+    assert_eq!(array.dim(), dimension);
+}
+
+#[test]
+#[cfg(feature = "ndarray")]
+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::<u8>(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();