Merge pull request #735 from kartoza/timlinux/issue704

Geenrate aggregate for polygon work opportunities layer

geest/core/algorithms/__init__.py

@@ -1,4 +1,13 @@
 from .area_iterator import AreaIterator
 from .population_processor import PopulationRasterProcessingTask
-from .wee_score_processor import WEEByPopulationScoreProcessingTask
+from .wee_by_population_score_processor import WEEByPopulationScoreProcessingTask
 from .subnational_aggregation_processor import SubnationalAggregationProcessingTask
+from .opportunities_mask_processor import OpportunitiesMaskProcessor
+from .utilities import (
+    assign_crs_to_raster_layer,
+    assign_crs_to_vector_layer,
+    subset_vector_layer,
+    geometry_to_memory_layer,
+    check_and_reproject_layer,
+    combine_rasters_to_vrt,
+)

geest/core/algorithms/population_processor.py

@@ -1,7 +1,7 @@
 import os
 import traceback
 import shutil
-from typing import Optional, Tuple
+from typing import Optional
 import subprocess
 import platform
 
@@ -10,7 +10,6 @@
     QgsTask,
     QgsProcessingContext,
     QgsFeedback,
-    QgsCoordinateReferenceSystem,
     QgsRasterLayer,
     QgsRasterDataProvider,
     QgsVectorLayer,
@@ -34,7 +33,6 @@ class PopulationRasterProcessingTask(QgsTask):
         study_area_gpkg_path (str): Path to the GeoPackage containing study area masks.
         output_dir (str): Directory to save the output rasters.
         cell_size_m (float): Cell size for the output rasters.
-        crs (Optional[QgsCoordinateReferenceSystem]): CRS for the output rasters. Will use the CRS of the clip layer if not provided.
         context (Optional[QgsProcessingContext]): QGIS processing context.
         feedback (Optional[QgsFeedback]): QGIS feedback object.
         force_clear (bool): Flag to force clearing of all outputs before processing.
@@ -46,7 +44,6 @@ def __init__(
         study_area_gpkg_path: str,
         working_directory: str,
         cell_size_m: float,
-        target_crs: Optional[QgsCoordinateReferenceSystem] = None,
         context: Optional[QgsProcessingContext] = None,
         feedback: Optional[QgsFeedback] = None,
         force_clear: bool = False,
@@ -61,15 +58,14 @@ def __init__(
                 os.remove(os.path.join(self.output_dir, file))
         self.cell_size_m = cell_size_m
         os.makedirs(self.output_dir, exist_ok=True)
-        self.target_crs = target_crs
-        if not self.target_crs:
-            layer: QgsVectorLayer = QgsVectorLayer(
-                f"{self.study_area_gpkg_path}|layername=study_area_clip_polygons",
-                "study_area_clip_polygons",
-                "ogr",
-            )
-            self.target_crs = layer.crs()
-            del layer
+
+        layer: QgsVectorLayer = QgsVectorLayer(
+            f"{self.study_area_gpkg_path}|layername=study_area_clip_polygons",
+            "study_area_clip_polygons",
+            "ogr",
+        )
+        self.target_crs = layer.crs()
+        del layer
         self.context = context
         self.feedback = feedback
         self.global_min = float("inf")
@@ -132,13 +128,10 @@ def clip_population_rasters(self) -> None:
             clip_layer.commitChanges()
 
             layer_name = f"{index}.tif"
-            output_path = os.path.join(self.output_dir, f"clipped_{layer_name}")
-            log_message(f"Processing mask {output_path}")
-
-            if not self.force_clear and os.path.exists(output_path):
-                log_message(f"Reusing existing clipped raster: {output_path}")
-                self.clipped_rasters.append(output_path)
-                continue
+            phase1_output = os.path.join(
+                self.output_dir, f"clipped_phase1_{layer_name}"
+            )
+            log_message(f"Processing mask {phase1_output}")
 
             # Clip the population raster using the mask
             params = {
@@ -151,21 +144,84 @@ def clip_population_rasters(self) -> None:
                 "KEEP_RESOLUTION": True,
                 "OPTIONS": "",
                 "DATA_TYPE": 5,  # Float32
-                "OUTPUT": output_path,
+                "OUTPUT": phase1_output,
             }
 
+            if not self.force_clear and os.path.exists(phase1_output):
+                log_message(f"Reusing existing clip phase 1 raster: {phase1_output}")
+            else:
+                result = processing.run("gdal:cliprasterbymasklayer", params)
+                if not result["OUTPUT"]:
+                    log_message(
+                        f"Failed to do phase1 clip raster for mask: {layer_name}"
+                    )
+                    continue
+
+            del clip_layer
+
+            clipped_layer = QgsRasterLayer(
+                phase1_output, f"Phase1 Clipped {layer_name}"
+            )
+            if not clipped_layer.isValid():
+                log_message(f"Invalid clipped raster layer for phase1: {layer_name}")
+                continue
+            del clipped_layer
+
+            log_message("Expanding clip layer to area bbox now ....")
+            # Now we need to expand the raster to the area_bbox so that it alighns
+            # with the clipped products produced by workflows
+            phase2_output = os.path.join(
+                self.output_dir, f"clipped_phase2_{layer_name}"
+            )
+
+            if not self.force_clear and os.path.exists(phase2_output):
+                log_message(f"Reusing existing phase2 clipped raster: {phase2_output}")
+                self.clipped_rasters.append(phase2_output)
+                continue
+
+            clip_layer = QgsVectorLayer("Polygon", "clip", "memory")
+            clip_layer.setCrs(self.target_crs)
+            clip_layer.startEditing()
+            feature = QgsFeature()
+            feature.setGeometry(current_bbox)
+            clip_layer.addFeature(feature)
+            clip_layer.commitChanges()
+
+            params = {
+                "INPUT": phase1_output,
+                "MASK": clip_layer,
+                "SOURCE_CRS": None,
+                "TARGET_CRS": self.target_crs,
+                "TARGET_EXTENT": clip_area.boundingBox(),
+                "NODATA": None,
+                "ALPHA_BAND": False,
+                "CROP_TO_CUTLINE": False,
+                "KEEP_RESOLUTION": False,
+                "SET_RESOLUTION": False,
+                "X_RESOLUTION": self.cell_size_m,
+                "Y_RESOLUTION": self.cell_size_m,
+                "MULTITHREADING": False,
+                "OPTIONS": "",
+                "DATA_TYPE": 0,
+                "EXTRA": "",
+                "OUTPUT": phase2_output,
+            }
             result = processing.run("gdal:cliprasterbymasklayer", params)
+            del clip_layer
 
             if not result["OUTPUT"]:
-                log_message(f"Failed to clip raster for mask: {layer_name}")
+                log_message(f"Failed to do phase2 clip raster for mask: {layer_name}")
                 continue
 
-            clipped_layer = QgsRasterLayer(output_path, f"Clipped {layer_name}")
+            clipped_layer = QgsRasterLayer(
+                phase2_output, f"Phase2 Clipped {layer_name}"
+            )
             if not clipped_layer.isValid():
-                log_message(f"Invalid clipped raster layer for mask: {layer_name}")
+                log_message(f"Invalid clipped raster layer for phase2: {layer_name}")
                 continue
+            del clipped_layer
 
-            self.clipped_rasters.append(output_path)
+            self.clipped_rasters.append(phase2_output)
 
             log_message(f"Processed mask {layer_name}")
 
@@ -217,7 +273,7 @@ def resample_population_rasters(self) -> None:
                 return
 
             layer_name = f"{index}.tif"
-            input_path = os.path.join(self.output_dir, f"clipped_{layer_name}")
+            input_path = os.path.join(self.output_dir, f"clipped_phase2_{layer_name}")
             output_path = os.path.join(self.output_dir, f"resampled_{layer_name}")
 
             log_message(f"Resampling {output_path} from {input_path}")

geest/core/algorithms/subnational_aggregation_processor.py

@@ -63,7 +63,7 @@ class SubnationalAggregationProcessingTask(QgsTask):
     | ![#0000FF](#) `#0000FF` | Highly enabling, medium population        |
     | ![#0000FF](#) `#0000FF` | Highly enabling, high population          |
 
-    See the wee_score_processor.py module for more details on how this is computed.
+    See the wee_by_population_score_processor.py module for more details on how this is computed.
 
     📒 The majority score for each subnational area is calculated by counting the
        number of pixels in each WEE Score and WEE x Population Score class as per the
@@ -112,15 +112,17 @@ def __init__(
 
         # These folders should already exist from the aggregation analysis and population raster processing
         self.population_folder = os.path.join(working_directory, "population")
-        self.wee_folder = os.path.join(working_directory, "wee_score")
+        self.wee_by_population_folder = os.path.join(
+            working_directory, "wee_by_population_score"
+        )
 
         if not os.path.exists(self.population_folder):
             raise Exception(
                 f"Population folder not found:\n{self.population_folder}\nPlease run population raster processing first."
             )
-        if not os.path.exists(self.wee_folder):
+        if not os.path.exists(self.wee_by_population_folder):
             raise Exception(
-                f"WEE folder not found.\n{self.wee_folder}\nPlease run WEE raster processing first."
+                f"WEE folder not found.\n{self.wee_by_population_folder}\nPlease run WEE raster processing first."
             )
 
         self.force_clear = force_clear
@@ -175,7 +177,7 @@ def aggregate(self) -> None:
         params = {
             "INPUT": output,
             "INPUT_RASTER": os.path.join(
-                self.wee_folder, "wee_by_population_score.vrt"
+                self.wee_by_population_folder, "wee_by_population_score.vrt"
             ),
             "RASTER_BAND": 1,
             "COLUMN_PREFIX": "_",