mmmm

COMBINE.py → src/COMBINE.py

@@ -18,11 +18,10 @@ except ModuleNotFoundError:
 with Flow("COMBINE") as COMBINE:
 
     destination = Parameter("destination")
-    source_rasters = Parameter('sources')
-    testing = Parameter('testing')
+    source_rasters = Parameter("sources")
+    testing = Parameter("testing")
 
     source_rasters = get_source_rasters(source_rasters)
-
     rw_windows = get_rw_windows(source_rasters=source_rasters, sample=testing)
 
     combinations_per_window = get_combinations.map(rw_windows, unmapped(source_rasters))
@@ -43,12 +42,14 @@ with Flow("COMBINE") as COMBINE:
     max_value = get_max_value(combination_dict)
     max_value.set_dependencies(upstream_tasks=[attribute_table_to_file])
 
-    target_raster_written = write_target_raster_chunk.map(rw_windows,
-                                                          unmapped(destination),
-                                                          unmapped(max_value),
-                                                          unmapped(combination_dict),
-                                                          unmapped(source_rasters))
-    target_raster_written.set_dependencies(upstream_tasks=[combination_dict])
+    target_raster_written = write_target_raster_chunk.map(
+        rw_windows,
+        unmapped(destination),
+        unmapped(max_value),
+        unmapped(combination_dict),
+        unmapped(source_rasters),
+    )
+    target_raster_written.set_dependencies(upstream_tasks=[combination_dict, max_value])
 
 
 if __name__ == "__main__":
@@ -61,7 +62,7 @@ if __name__ == "__main__":
         Path to destination *.tif file
     source_raster01: Path to first source raster
     source_raster02: Path to second source raster
-    ...: path to subsequent source  rasters 
+    ...: path to subsequent source  rasters
 
     Returns
     -------
@@ -69,15 +70,21 @@ if __name__ == "__main__":
 
     import argparse
 
-    parser = argparse.ArgumentParser(prog='Open Source COMBINE alternative. By: Hans Roelofsen')
+    parser = argparse.ArgumentParser(
+        prog="Open Source COMBINE alternative. By: Hans Roelofsen"
+    )
     parser.add_argument("destination", type=str, help="path to destination *.tif file")
-    parser.add_argument("sources", type=str, help="path to destination source files", nargs='+')
-    parser.add_argument('--test', action='store_true', help='test run with 20 windows')
+    parser.add_argument(
+        "sources", type=str, help="path to destination source files", nargs="+"
+    )
+    parser.add_argument("--test", action="store_true", help="test run with 20 windows")
     args = parser.parse_args()
 
     print("\n\nStarting COMBINE\n\n")
 
     COMBINE.executor = LocalDaskExecutor()
-    COMBINE.run(parameters=dict(destination=args.destination,
-                                sources=args.sources,
-                                testing=args.test))
+    COMBINE.run(
+        parameters=dict(
+            destination=args.destination, sources=args.sources, testing=args.test
+        )
+    )

attribute_table.py → src/attribute_table.py

@@ -29,7 +29,9 @@ def make_attribute_table(
         index=[v for v in combi_dict.values()],
     )
     tab.to_clipboard(excel=True)
-    tab.to_csv(f"{os.path.splitext(destination)[0]}.csv", index=True, index_label="Value")
+    tab.to_csv(
+        f"{os.path.splitext(destination)[0]}.csv", index=True, index_label="Value"
+    )
 
 
 @task
@@ -40,4 +42,3 @@ def get_max_value(combi_dict: dict) -> int:
     :return:
     """
     return len(combi_dict)
-

combinations.py → src/combinations.py

@@ -11,6 +11,7 @@ except ModuleNotFoundError:
     from src.windows import CombineWindow
     from src.source_rasters import SourceRasters
 
+
 @task
 def get_combinations(
     read_window: CombineWindow,
@@ -25,7 +26,7 @@ def get_combinations(
 
     logger = prefect.context.get("logger")
     logger.info(
-        f"Finding unique combinations in window {read_window.nr:004}"
+        f"Finding unique combinations in window {read_window.nr} out of {read_window.total}"
     )
 
     arrays = [
@@ -34,7 +35,9 @@ def get_combinations(
     ]
 
     unique_combinations = set(list(zip(*arrays)))
-    logger.info(f'lowest value in window {read_window.nr:004}: {np.min(np.array(unique_combinations))}')
+    logger.info(
+        f"lowest value in window {read_window.nr:004}: {np.min(np.array(unique_combinations))}"
+    )
     read_window.set_unique_combinations(uc=unique_combinations)
 
 
@@ -46,5 +49,3 @@ def make_combination_dict(read_windows: list) -> dict:
     return {combi: n for n, combi in enumerate(all_combinations, start=1)}
 
 
-src = r'w:\PROJECTS\Nvk_bdb\b_HDB\b_scenarios\2_regionaal\brondata\Regionaal_Geworteld_Subsector_rel_Y2050\25m\Regionaal_Geworteld_Subsector_rel_Y2050.tif'
-arr = rio.open(src, window=Window(0, 0, 100, 150)).read(1).flatten()
\ No newline at end of file

source_rasters.py → src/source_rasters.py

@@ -7,14 +7,16 @@ from abc import ABC, abstractmethod
 from prefect.engine import signals
 import pandas as pd
 
-class SourceRasters:
 
+class SourceRasters:
     def __init__(self, source_list: list):
 
         self.source_list = [s for s in source_list]
-        self.procedures = [VerifyAsGeospatialRaster(source_list=self.source_list),
-                           VerifyMatchingBounds(source_list=self.source_list)]
-        self.geospatial_profile = getattr(rio.open(self.source_list[0]), 'profile')
+        self.procedures = [
+            VerifyAsGeospatialRaster(source_list=self.source_list),
+            VerifyMatchingBounds(source_list=self.source_list),
+            VerifyMatchingResolution(source_list=self.source_list),
+        ]
 
     def verify_input(self):
         for procedure in self.procedures:
@@ -39,7 +41,7 @@ class VerificationProcedure(ABC):
         pass
 
 
-class VerifyAsGeospatialRaster(ABC):
+class VerifyAsGeospatialRaster(VerificationProcedure):
     def __init__(self, source_list):
         self.rasters = source_list
         self.status = True
@@ -53,10 +55,10 @@ class VerifyAsGeospatialRaster(ABC):
                 logger.info(message)
             else:
                 logger.error(message)
-                self.status=False
+                self.status = False
 
 
-class VerifyMatchingBounds(ABC):
+class VerifyMatchingBounds(VerificationProcedure):
     def __init__(self, source_list):
         self.rasters = source_list
         self.status = True
@@ -72,6 +74,27 @@ class VerifyMatchingBounds(ABC):
                 self.status = False
 
 
+class VerifyMatchingResolution(VerificationProcedure):
+
+    def __init__(self, source_list):
+        self.rasters = source_list
+        self.status = True
+
+    def execute(self):
+        logger = prefect.context.get("logger")
+
+        reference_resolution = getattr(getattr(rio.open(self.rasters[0]), 'transform'), 'a')
+        resolutions = []
+        messages = []
+        for raster in self.rasters:
+            has_resolution, msg = raster_has_resolution(raster, reference_resolution)
+            resolutions.append(has_resolution)
+            messages.append(msg)
+
+        if not all(resolutions):
+            logger.error('Rasters have unequal resolutions')
+            self.status = False
+
 def is_geospatial_raster(file_path: str, msg=False) -> (bool, str):
     """
     Is a file a gdal compatible geospatial raster?
@@ -161,6 +184,33 @@ def rasters_have_identical_bounds(rasters: list) -> (bool, str):
     return out, msg
 
 
+def raster_has_resolution(file_path: str, target_resolution) -> (bool, str):
+    """
+    Check if a raster has target resolution
+    Parameters
+    ----------
+    file_path: target file path
+    target_resolution: expected resolution of the raster
+
+    Returns
+    -------
+
+    """
+
+    if is_geospatial_raster(file_path):
+        raster_resolution = getattr(rio.open(file_path), "res")[0]
+        has_target_resolution = raster_resolution == target_resolution
+        if has_target_resolution:
+            out = True
+            msg = f"{os.path.basename(file_path)} has target resolution."
+        else:
+            out = False
+            msg = f"{os.path.basename(file_path)} resolution {raster_resolution}m does not meet target resolution {target_resolution}."
+        return out, msg
+    else:
+        return False, f"{os.path.basename} is not a geospatial raster"
+
+
 @task
 def get_source_rasters(source_list: list) -> SourceRasters:
 

target_raster.py → src/target_raster.py

@@ -12,31 +12,30 @@ except ModuleNotFoundError:
     from src.source_rasters import SourceRasters
     from src.windows import CombineWindow
 
+
 @task
 def write_target_raster_chunk(
-        rw_window: CombineWindow,
-        destination: str,
-        max_value: int,
-        combi_dict: dict,
-        source_rasters: SourceRasters
+    rw_window: CombineWindow,
+    destination: str,
+    max_value: int,
+    combi_dict: dict,
+    source_rasters: SourceRasters,
 ):
-    """
-    Write target raster for a single RW Window
-    """
+
 
     logger = prefect.context.get("logger")
 
     with rio.open(
-            destination,
-            "w",
-            driver="GTiff",
-            width=rw_window.rio_window.width,
-            height=rw_window.rio_window.height,
-            count=1,
-            dtype=rio.dtypes.get_minimum_dtype([0, max_value]),
-            transform=source_rasters.geospatial_profile["transform"],
-            compress="LZW",
-            crs=source_rasters.geospatial_profile["crs"],
+        f"{os.path.splitext(destination)[0]}_{rw_window.nr:05}.tif",
+        "w",
+        driver="GTiff",
+        width=rw_window.rio_window.width,
+        height=rw_window.rio_window.height,
+        count=1,
+        dtype=rio.dtypes.get_minimum_dtype([0, max_value]),
+        transform=rw_window.transform,
+        compress="LZW",
+        crs=source_rasters.geospatial_profile["crs"],
     ) as dest:
         dest.update_tags(
             creation_date=datetime.datetime.now().strftime("%d-%b-%Y_%H:%M:%S"),
@@ -73,6 +72,7 @@ def write_target_raster_chunk(
             indexes=1,
         )
 
+    logger.info(f"Writing chunk {rw_window.nr} out of {rw_window.total} to file.")
 
 
 @task

windows.py → src/windows.py

 import random
 import prefect
 import affine
+import numpy as np
 import rasterio as rio
 from prefect import task
 from rasterio.windows import Window
@@ -16,20 +17,28 @@ class CombineWindow:
     A single read/write window and its associated unique combinations
     """
 
-    def __init__(self, row_start: int, row_stop: int, column_start: int, column_stop: int, nr: int,
-                 top_left_x: int, top_left_y: int, resolution: int):
+    def __init__(
+        self,
+        row_start: int,
+        row_stop: int,
+        column_start: int,
+        column_stop: int,
+        nr: int,
+        totals: int,
+        top_left_x: int,
+        top_left_y: int,
+        resolution: int,
+    ):
         self.rio_window = Window.from_slices(
             (row_start, row_stop), (column_start, column_stop)
         )
         self.unique_combinations = None
         self.nr = nr
-        self.transform = affine.Affine(a=resolution,
-                                      b=0,
-                                      c=top_left_x,
-                                      d=0,
-                                      e=resolution*-1,
-                                      f=top_left_y)
+        self.total = totals
 
+        self.transform = affine.Affine(
+            a=resolution, b=0, c=top_left_x, d=0, e=resolution * -1, f=top_left_y
+        )
 
     def set_unique_combinations(self, uc: set):
         self.unique_combinations = uc
@@ -41,35 +50,49 @@ class CombineWindow:
 @task
 def get_rw_windows(
     source_rasters: SourceRasters,
-    window_rows: int = 2000,
-    window_columns: int = 2000,
+    window_height: int = 2000,
+    window_width: int = 2000,
     sample: bool = False,
 ) -> list:
+    """ """
 
     logger = prefect.context.get("logger")
 
-    row_starts = range(0, (source_rasters.geospatial_profile['height'] - window_rows), window_rows)
-    column_starts = range(0, (source_rasters.geospatial_profile['width'] - window_columns), window_columns)
+    row_starts = range(
+        0, (source_rasters.geospatial_profile["height"] - window_height), window_height
+    )
+    column_starts = range(
+        0, (source_rasters.geospatial_profile["width"] - window_width), window_width
+    )
     windows = []
     n = 0
+    total_windows = np.multiply(len(row_starts), len(column_starts))
     for row_start in row_starts:
         for column_start in column_starts:
             windows.append(
                 CombineWindow(
                     row_start=row_start,
-                    row_stop=row_start + window_rows,
+                    row_stop=row_start + window_height,
                     column_start=column_start,
-                    column_stop=column_start + window_columns,
+                    column_stop=column_start + window_width,
                     nr=n,
-                    top_left_x=(getattr(rio.open(source_rasters.src_list[0]), 'transform') * (column_start, row_start))[0],
-                    top_left_y=(getattr(rio.open(source_rasters.src_list[0]), 'transform') * (column_start, row_start))[1],
-                    resolution=getattr(getattr(rio.open(source_rasters.src_list[0]), 'transform'), 'a'),
+                    totals=total_windows,
+                    top_left_x=(
+                        getattr(rio.open(source_rasters.source_list[0]), "transform")
+                        * (column_start, row_start)
+                    )[0],
+                    top_left_y=(
+                        getattr(rio.open(source_rasters.source_list[0]), "transform")
+                        * (column_start, row_start)
+                    )[1],
+                    resolution=getattr(
+                        getattr(rio.open(source_rasters.source_list[0]), "transform"),
+                        "a",
+                    ),
                 )
             )
             n += 1
-    logger.info(
-                f"Created {n} windows of size {window_rows}X{window_columns}"
-            )
+    logger.info(f"Created {n} windows of size {window_height}X{window_width}")
 
     if sample:
         return random.sample(windows, 20)