adding slicing to to raster class

README.rst

@@ -12,6 +12,7 @@ Fauxgeo
         :target: https://readthedocs.org/projects/fauxgeo/?badge=latest
         :alt: Documentation Status
 
+`Documentation at Read the Docs <http://fauxgeo.readthedocs.org/en/latest/>`_
 
 A python library that generates simple OSGeo-supported rasters and vectors.  The primary purpose for this library is to help test geoprocessing functions.
 
@@ -27,6 +28,8 @@ fauxgeo 0.2.0 requires
 * Scipy
 * Matplotlib
 * GDAL
+* PyGeoprocessing
+* Shapely
 
 Installation
 ------------
@@ -50,7 +53,6 @@ The Raster Class
 .. code:: python
 
     import numpy as np
-    from affine import Affine
     import gdal
 
     # set arguments
@@ -98,6 +100,10 @@ The Raster Class
     raster4 = 4.5 - raster3
     raster4 = raster3 ** 4.5
 
+    # get and set slices of a raster
+    print raster4[1:9:2]
+    raster4[1:9:2] = 6
+
     raster4 = raster3.minimum(raster2)
 
     # returns base rasters with same nodata and datatype
@@ -126,7 +132,6 @@ The RasterFactory Class
 
 .. code:: python
 
-    from affine import Affine
     import gdal
 
     # set arguments
@@ -146,7 +151,7 @@ The RasterFactory Class
     test_raster_4 = factory.horizontal_ramp(1, 10)  # interpolated from 1 to 10 across columns
     test_raster_5 = factory.vertical_ramp(1, 10)  # interpolated from 1 to 10 across rows
 
-The Vector Class
+The Vector Class (a work in progress...)
 
 .. code:: python
 
@@ -178,7 +183,9 @@ Planning
 --------
 
 * Add basic visualization functionality
+
 * Raster Operations
+
     * Reclass
     * Overlay - intersection, union, clip
     * Dissolve

fauxgeo/raster.py

@@ -34,15 +34,18 @@ def __init__(self, uri, driver):
         self.dataset = None
 
     @classmethod
-    def from_array(self, array, affine, proj, datatype, nodata_val, driver='GTiff'):
+    def from_array(self, array, affine, proj, datatype, nodata_val, driver='GTiff', filepath=None):
         if len(array.shape) is 2:
             num_bands = 1
         elif len(array.shape) is 3:
             num_bands = len(array)
         else:
             raise ValueError
 
-        dataset_uri = pygeo.geoprocessing.temporary_filename()
+        if filepath:
+            dataset_uri = filepath
+        else:
+            dataset_uri = pygeo.geoprocessing.temporary_filename()
         rows = array.shape[0]
         cols = array.shape[1]
 
@@ -67,7 +70,8 @@ def from_array(self, array, affine, proj, datatype, nodata_val, driver='GTiff'):
         dataset = None
         driver = None
 
-        return Raster(dataset_uri, driver=driver)
+        if not filepath:
+            return Raster(dataset_uri, driver=driver)
 
     @classmethod
     def from_file(self, uri, driver='GTiff'):
@@ -319,16 +323,19 @@ def f(x, y): return np.where((np.not_equal(x, y)), np.fmin(x, y), np.fmin(x, y))
                 return mini
             return self.local_op(raster, min_closure)
 
-    def __getitem__(self):
-        pass  # return numpy slice?  Raster object with sliced numpy array?
+    def __getitem__(self, key):
+        arr = self.get_band(1)
+        return arr[key]
 
-    def __setitem__(self):
-        pass  # set numpy values to raster
+    def __setitem__(self, key, item):
+        arr = self.get_band(1)
+        arr[key] = item
+        self.set_band(1, arr)
 
-    def __getslice__(self):
+    def getslice(self, a, b, c):
         pass
 
-    def __setslice__(self):
+    def setslice(self, a, b, c):
         pass
 
     def __iter__(self):
@@ -589,8 +596,18 @@ def get_cell_area(self):
         a = self.get_affine()
         return abs((a.a * a.e))
 
-    def set_band(self, masked_array):
-        raise NotImplementedError
+    def set_band(self, band_num, masked_array):
+        affine = self.get_affine()
+        proj = self.get_projection()
+        datatype = self.get_datatype(band_num)
+        nodata_val = masked_array.fill_value
+        if np.issubdtype(nodata_val, int):
+            nodata_val = nodata_val.astype(int)
+        else:
+            nodata_val = nodata_val.astype(float)
+        uri = self.uri
+        Raster.from_array(
+            masked_array, affine, proj, datatype, nodata_val, filepath=uri)
 
     def set_bands(self, array):
         raise NotImplementedError

tests/test_raster.py

@@ -263,6 +263,7 @@ def test_resize_pixels(self):
         b = a.resize_pixels(0.5, 'nearest')
         assert(b.get_shape() == (8, 8))
 
+
 class TestRasterReclassMaskedValues(unittest.TestCase):
     def setUp(self):
         self.shape = (4, 4)
@@ -316,6 +317,24 @@ def test_zeros(self):
         assert(ones_raster.get_band(1)[0, 0] == 1)
 
 
+class TestRasterSetBand(unittest.TestCase):
+    def setUp(self):
+        self.shape = (4, 4)
+        self.array = np.ones(self.shape)
+        self.affine = Affine(1, 0, 0, 0, -1, 4)
+        self.proj = 4326
+        self.datatype = gdal.GDT_Float64
+        self.nodata_val = -9999
+        self.factory = RasterFactory(
+            self.proj, self.datatype, self.nodata_val, *self.shape, affine=self.affine)
+
+    def test_set_band(self):
+        raster = self.factory.alternating(-9999, 2.0)
+        new_array = np.ma.ones((4, 4), fill_value=0)
+        raster.set_band(1, new_array)
+        assert(raster[0][0] == 1)
+
+
 class TestRasterSetDatatype(unittest.TestCase):
     def setUp(self):
         self.shape = (4, 4)
@@ -399,7 +418,6 @@ def test_std(self):
 
     def test_sum(self):
         raster = self.factory.alternating(1.0, -9999)
-        print raster
         assert(raster.sum() == 8.0)
 
 
@@ -445,5 +463,27 @@ def test_unique(self):
         assert(raster.unique() == [1, 2, 3, 4])
 
 
+class TestRasterSlice(unittest.TestCase):
+    def setUp(self):
+        self.shape = (10, 10)
+        self.array = np.ones(self.shape)
+        self.affine = Affine(1, 0, 0, 0, -1, 4)
+        self.proj = 4326
+        self.datatype = gdal.GDT_Int16
+        self.nodata_val = -9999
+        self.factory = RasterFactory(
+            self.proj, self.datatype, self.nodata_val, *self.shape, affine=self.affine)
+
+    def test_get_slice(self):
+        raster = self.factory.alternating(0, 1)
+        a = raster[1:3]
+        assert(a[0][0] == 1)
+
+    def test_set_slice(self):
+        raster = self.factory.alternating(0, 1)
+        raster[1:9:2, 1:9:2] = 2
+        assert(raster[1][3] == 2)
+
+
 if __name__ == '__main__':
     unittest.main()