Add RasterStack class

fauxgeo/affine.py

@@ -27,17 +27,14 @@ def __eq__(self, other):
         d = (self.d == other.d)
         e = (self.e == other.e)
         f = (self.f == other.f)
-        if all([a, b, c, d, e, f]):
-            return True
-        else:
-            return False
+        return all([a, b, c, d, e, f])
 
     @classmethod
     def identity(self):
         """Return identify transform."""
         return Affine(1.0, 0.0, 0.0, 0.0, 1.0, 0.0)
 
-    @ classmethod
+    @classmethod
     def from_gdal(self, c, a, b, f, d, e):
         """Convert from a gdal geotransform."""
         return Affine(a, b, c, d, e, f)

fauxgeo/raster.py

@@ -17,6 +17,7 @@
 from shapely.geometry import Polygon
 import shapely
 import pygeoprocessing as pygeo
+import pygeoprocessing.geoprocessing as geoprocess
 
 from fauxgeo.vector import Vector
 from fauxgeo.affine import Affine
@@ -48,7 +49,7 @@ def from_array(array, affine, proj, datatype, nodata_val, resample_method=None,
         if filepath:
             dataset_uri = filepath
         else:
-            dataset_uri = pygeo.geoprocessing.temporary_filename()
+            dataset_uri = geoprocess.temporary_filename()
         rows = array.shape[0]
         cols = array.shape[1]
 
@@ -74,11 +75,12 @@ def from_array(array, affine, proj, datatype, nodata_val, resample_method=None,
         driver = None
 
         if not filepath:
-            return Raster(dataset_uri, resample_method=resample_method, driver=driver)
+            return Raster(
+                dataset_uri, resample_method=resample_method, driver=driver)
 
     @staticmethod
     def from_file(uri, resample_method=None, driver='GTiff'):
-        dataset_uri = pygeo.geoprocessing.temporary_filename()
+        dataset_uri = geoprocess.temporary_filename()
         if not os.path.isabs(uri):
             uri = os.path.join(os.getcwd(), uri)
         # assert existence
@@ -111,7 +113,7 @@ def _delete(self):
         os.remove(self.uri)
 
     def __str__(self):
-        string = '\nRASTER'
+        string = '\nRASTER___'
         string += '\nNumber of Bands: ' + str(self.band_count())
         string += '\nBand 1:\n' + self.get_band(1).__repr__()
         string += self.get_affine().__repr__()
@@ -451,6 +453,13 @@ def get_bands(self):
         self._close_dataset()
         return a
 
+    def get_resample_method(self):
+        if not self.resample_method:
+            raise AttributeError(
+                'Raster object has no assigned resample_method attribute')
+        else:
+            return self.resample_method
+
     def get_nodata(self, band_num):
         nodata_val = None
         self._open_dataset()
@@ -573,7 +582,7 @@ def get_affine(self):
         return Affine.from_gdal(*geotransform)
 
     def get_bounding_box(self):
-        return pygeo.geoprocessing.get_bounding_box(self.uri)
+        return geoprocess.get_bounding_box(self.uri)
 
     def get_aoi(self):
         """May only be suited for non-rotated rasters."""
@@ -623,12 +632,12 @@ def copy(x):
 
         nodata = self.get_nodata(1)
         pixel_op = pixel_op_closure(nodata)
-        dataset_out_uri = pygeo.geoprocessing.temporary_filename()
+        dataset_out_uri = geoprocess.temporary_filename()
         datatype_out = datatype
         nodata_out = nodata
-        pixel_size_out = pygeo.geoprocessing.get_cell_size_from_uri(self.uri)
+        pixel_size_out = geoprocess.get_cell_size_from_uri(self.uri)
 
-        pygeo.geoprocessing.vectorize_datasets(
+        geoprocess.vectorize_datasets(
             dataset_uri_list,
             pixel_op,
             dataset_out_uri,
@@ -659,12 +668,12 @@ def copy(x):
 
         old_nodata = self.get_nodata(1)
         pixel_op = pixel_op_closure(old_nodata, nodata_val)
-        dataset_out_uri = pygeo.geoprocessing.temporary_filename()
-        datatype_out = pygeo.geoprocessing.get_datatype_from_uri(self.uri)
+        dataset_out_uri = geoprocess.temporary_filename()
+        datatype_out = geoprocess.get_datatype_from_uri(self.uri)
         nodata_out = nodata_val
-        pixel_size_out = pygeo.geoprocessing.get_cell_size_from_uri(self.uri)
+        pixel_size_out = geoprocess.get_cell_size_from_uri(self.uri)
 
-        pygeo.geoprocessing.vectorize_datasets(
+        geoprocess.vectorize_datasets(
             dataset_uri_list,
             pixel_op,
             dataset_out_uri,
@@ -695,12 +704,12 @@ def copy(x):
 
         old_nodata = self.get_nodata(1)
         pixel_op = pixel_op_closure(old_nodata, nodata_val)
-        dataset_out_uri = pygeo.geoprocessing.temporary_filename()
+        dataset_out_uri = geoprocess.temporary_filename()
         datatype_out = datatype
         nodata_out = nodata_val
-        pixel_size_out = pygeo.geoprocessing.get_cell_size_from_uri(self.uri)
+        pixel_size_out = geoprocess.get_cell_size_from_uri(self.uri)
 
-        pygeo.geoprocessing.vectorize_datasets(
+        geoprocess.vectorize_datasets(
             dataset_uri_list,
             pixel_op,
             dataset_out_uri,
@@ -718,7 +727,7 @@ def copy(x):
 
     def copy(self, uri=None):
         if not uri:
-            uri = pygeo.geoprocessing.temporary_filename()
+            uri = geoprocess.temporary_filename()
         if not os.path.isabs(uri):
             uri = os.path.join(os.getcwd(), uri)
         shutil.copyfile(self.uri, uri)
@@ -739,13 +748,13 @@ def align(self, raster, resample_method):
 
         def dataset_pixel_op(x, y): return y
         dataset_uri_list = [self.uri, raster.uri]
-        dataset_out_uri = pygeo.geoprocessing.temporary_filename()
-        datatype_out = pygeo.geoprocessing.get_datatype_from_uri(raster.uri)
-        nodata_out = pygeo.geoprocessing.get_nodata_from_uri(raster.uri)
-        pixel_size_out = pygeo.geoprocessing.get_cell_size_from_uri(self.uri)
+        dataset_out_uri = geoprocess.temporary_filename()
+        datatype_out = geoprocess.get_datatype_from_uri(raster.uri)
+        nodata_out = geoprocess.get_nodata_from_uri(raster.uri)
+        pixel_size_out = geoprocess.get_cell_size_from_uri(self.uri)
         bounding_box_mode = "dataset"
 
-        pygeo.geoprocessing.vectorize_datasets(
+        geoprocess.vectorize_datasets(
             dataset_uri_list,
             dataset_pixel_op,
             dataset_out_uri,
@@ -769,13 +778,13 @@ def align_to(self, raster, resample_method):
         def dataset_pixel_op(x, y): return y
 
         dataset_uri_list = [raster.uri, self.uri]
-        dataset_out_uri = pygeo.geoprocessing.temporary_filename()
-        datatype_out = pygeo.geoprocessing.get_datatype_from_uri(self.uri)
-        nodata_out = pygeo.geoprocessing.get_nodata_from_uri(self.uri)
-        pixel_size_out = pygeo.geoprocessing.get_cell_size_from_uri(raster.uri)
+        dataset_out_uri = geoprocess.temporary_filename()
+        datatype_out = geoprocess.get_datatype_from_uri(self.uri)
+        nodata_out = geoprocess.get_nodata_from_uri(self.uri)
+        pixel_size_out = geoprocess.get_cell_size_from_uri(raster.uri)
         bounding_box_mode = "dataset"
 
-        pygeo.geoprocessing.vectorize_datasets(
+        geoprocess.vectorize_datasets(
             dataset_uri_list,
             dataset_pixel_op,
             dataset_out_uri,
@@ -793,13 +802,13 @@ def dataset_pixel_op(x, y): return y
 
     def clip(self, aoi_uri):
         r = None
-        dataset_out_uri = pygeo.geoprocessing.temporary_filename()
+        dataset_out_uri = geoprocess.temporary_filename()
         datatype = self.get_datatype(1)
         nodata = self.get_nodata(1)
         pixel_size = self.get_affine().a
 
         try:
-            pygeo.geoprocessing.clip_dataset_uri(
+            geoprocess.clip_dataset_uri(
                 self.uri,
                 aoi_uri,
                 dataset_out_uri,
@@ -848,21 +857,21 @@ def reproject(self, proj, resample_method, pixel_size=None):
         if pixel_size is None:
             pixel_size = self.get_affine().a
 
-        dataset_out_uri = pygeo.geoprocessing.temporary_filename()
+        dataset_out_uri = geoprocess.temporary_filename()
         srs = osr.SpatialReference()
         srs.ImportFromEPSG(proj)
         wkt = srs.ExportToWkt()
 
-        pygeo.geoprocessing.reproject_dataset_uri(
+        geoprocess.reproject_dataset_uri(
             self.uri, pixel_size, wkt, resample_method, dataset_out_uri)
 
         return Raster.from_tempfile(dataset_out_uri)
 
     def resize_pixels(self, pixel_size, resample_method):
         bounding_box = self.get_bounding_box()
-        output_uri = pygeo.geoprocessing.temporary_filename()
+        output_uri = geoprocess.temporary_filename()
 
-        pygeo.geoprocessing.resize_and_resample_dataset_uri(
+        geoprocess.resize_and_resample_dataset_uri(
             self.uri,
             bounding_box,
             pixel_size,
@@ -901,12 +910,12 @@ def sample_from_raster(self, raster):
 
     def reclass(self, reclass_table, out_nodata=None, out_datatype=None):
         if out_nodata is None:
-            out_nodata = pygeo.geoprocessing.get_nodata_from_uri(self.uri)
+            out_nodata = geoprocess.get_nodata_from_uri(self.uri)
         if out_datatype is None:
-            out_datatype = pygeo.geoprocessing.get_datatype_from_uri(self.uri)
-        dataset_out_uri = pygeo.geoprocessing.temporary_filename()
+            out_datatype = geoprocess.get_datatype_from_uri(self.uri)
+        dataset_out_uri = geoprocess.temporary_filename()
 
-        pygeo.geoprocessing.reclassify_dataset_uri(
+        geoprocess.reclassify_dataset_uri(
             self.uri,
             reclass_table,
             dataset_out_uri,
@@ -951,12 +960,12 @@ def local_op(self, raster, pixel_op_closure, broadcast=False):
 
         nodata = self.get_nodata(1)
         pixel_op = pixel_op_closure(nodata)
-        dataset_out_uri = pygeo.geoprocessing.temporary_filename()
-        datatype_out = pygeo.geoprocessing.get_datatype_from_uri(self.uri)
-        nodata_out = pygeo.geoprocessing.get_nodata_from_uri(self.uri)
-        pixel_size_out = pygeo.geoprocessing.get_cell_size_from_uri(self.uri)
+        dataset_out_uri = geoprocess.temporary_filename()
+        datatype_out = geoprocess.get_datatype_from_uri(self.uri)
+        nodata_out = geoprocess.get_nodata_from_uri(self.uri)
+        pixel_size_out = geoprocess.get_cell_size_from_uri(self.uri)
 
-        pygeo.geoprocessing.vectorize_datasets(
+        geoprocess.vectorize_datasets(
             dataset_uri_list,
             pixel_op,
             dataset_out_uri,

fauxgeo/raster_factory.py

@@ -80,5 +80,6 @@ def create_sample_global_map():
         datatype = 5
         nodata_val = -9999
         proj = 4326
-        factory = RasterFactory(proj, datatype, nodata_val, shape[1], shape[0])
+        factory = RasterFactory(
+            proj, datatype, nodata_val, shape[1], shape[0])
         return factory.uniform(1)

fauxgeo/raster_stack.py

@@ -17,7 +17,9 @@
 from shapely.geometry import Polygon
 import shapely
 import pygeoprocessing as pygeo
+from pygeoprocessing import geoprocessing as geoprocess
 
+from fauxgeo.raster import Raster
 from fauxgeo.vector import Vector
 from fauxgeo.affine import Affine
 
@@ -33,11 +35,66 @@ class RasterStack(object):
     def __init__(self, raster_list):
         self.raster_list = raster_list
 
+    def __del__(self):
+        for i in self.raster_list:
+            del i
+
+    def __str__(self):
+        string =  '\nRasterStack'
+        for i in self.raster_list:
+            string += '\n  ' + i.uri
+        return string
+
+    def get_raster_uri_list(self):
+        """Get Raster URI List."""
+        return [r.uri for r in self.raster_list]
+
     def assert_same_projection(self):
-        pass
+        """Assert Rasters in Same Projection."""
+        return all(self.raster_list[0].get_projection() == r.get_projection() \
+            for r in self.raster_list)
+
+    def assert_same_alignment(self):
+        """Assert Rasters in Same Alignment."""
+        return all(self.raster_list[0].get_affine() == r.get_affine() \
+            for r in self.raster_list)
+
+    def assert_resample_methods_assigned(self):
+        try:
+            for r in self.raster_list:
+                r.get_resample_method()
+            return True
+        except AttributeError:
+            return False
 
     def align(self):
-        pass
+        """Align Rasters."""
+        if not self.assert_resample_methods_assigned():
+            raise ValueError("Resample method not assigned to at least one raster.")
+        dataset_uri_list = self.get_raster_uri_list()
+        dataset_out_uri_list = \
+            [geoprocess.temporary_filename() for _ in self.raster_list]
+        resample_method_list = \
+            [r.get_resample_method() for r in self.raster_list]
+        out_pixel_size = geoprocess.get_cell_size_from_uri(
+            self.raster_list[0].uri)
+        mode = 'dataset'
+        dataset_to_align_index = 0
+
+        raster_uri_list = geoprocess.align_dataset_list(
+            dataset_uri_list,
+            dataset_out_uri_list,
+            resample_method_list,
+            out_pixel_size,
+            mode,
+            dataset_to_align_index,
+            dataset_to_bound_index=0)
+
+        new_raster_list = [Raster.from_file(fp) for fp in dataset_out_uri_list]
+        return RasterStack(new_raster_list)
 
-    def set_standard_nodata_and_dtype(self):
-        pass
+    def set_standard_nodata(self, nodata_val):
+        new_raster_list = []
+        for r in self.raster_list:
+            new_raster_list.append(r.set_nodata(nodata_val))
+        return RasterStack(new_raster_list)