Merge pull request #397 from kartoza/osundwajeff/issue389

Refactor raster_polyline and raster_polygons grid scores
kartoza · Oct 10, 2024 · da32e0c · da32e0c
2 parents f077cdf + 594d7f5
commit da32e0c
Show file tree

Hide file tree

Showing 7 changed files with 208 additions and 89 deletions.
diff --git a/geest/core/polygons_per_grid_cell.py b/geest/core/polygons_per_grid_cell.py
@@ -1,23 +1,34 @@
 import os
 from qgis.PyQt.QtCore import QVariant
 from qgis.core import (
+    QgsGeometry,
     QgsVectorLayer,
     QgsField,
     QgsSpatialIndex,
     QgsProcessingFeedback,
 )
 import processing
-from .create_grids import GridCreator
-from .extents import Extents
+from .utilities import GridAligner
 
 
 class RasterPolygonGridScore:
-    def __init__(self, country_boundary, pixel_size, output_path, crs, input_polygons):
+    def __init__(
+        self,
+        country_boundary,
+        pixel_size,
+        working_dir,
+        crs,
+        input_polygons,
+        output_path,
+    ):
         self.country_boundary = country_boundary
         self.pixel_size = pixel_size
-        self.output_path = output_path
+        self.working_dir = working_dir
         self.crs = crs
         self.input_polygons = input_polygons
+        self.output_path = output_path
+        # Initialize GridAligner with grid size
+        self.grid_aligner = GridAligner(grid_size=100)
 
     def raster_polygon_grid_score(self):
         """
@@ -29,16 +40,52 @@ def raster_polygon_grid_score(self):
         :param input_polygons: Layer of point features to count within each grid cell.
         """
 
-        # Create grid
-        self.h_spacing = 100
-        self.v_spacing = 100
-        create_grid = GridCreator(h_spacing=self.h_spacing, v_spacing=self.v_spacing)
-        output_dir = os.path.join("output")
-        merged_output_path = os.path.join(output_dir, "merged_grid.shp")
-        grid_layer = create_grid.create_grids(
-            self.country_boundary, output_dir, self.crs, merged_output_path
+        output_dir = os.path.dirname(self.output_path)
+
+        # Define output directory and ensure it's created
+        os.makedirs(output_dir, exist_ok=True)
+
+        # Load grid layer from the Geopackage
+        geopackage_path = os.path.join(
+            self.working_dir, "study_area", "study_area.gpkg"
+        )
+        if not os.path.exists(geopackage_path):
+            raise ValueError(f"Geopackage not found at {geopackage_path}.")
+
+        grid_layer = QgsVectorLayer(
+            f"{geopackage_path}|layername=study_area_grid", "merged_grid", "ogr"
+        )
+
+        area_layer = QgsVectorLayer(
+            f"{geopackage_path}|layername=study_area_polygons",
+            "study_area_polygons",
+            "ogr",
+        )
+
+        geometries = [feature.geometry() for feature in area_layer.getFeatures()]
+
+        # Combine all geometries into one using unaryUnion
+        area_geometry = QgsGeometry.unaryUnion(geometries)
+
+        # grid_geometry = grid_layer.getGeometry()
+
+        aligned_bbox = self.grid_aligner.align_bbox(
+            area_geometry.boundingBox(), area_layer.extent()
         )
-        grid_layer = QgsVectorLayer(merged_output_path, "merged_grid", "ogr")
+
+        # Extract polylines by location
+        grid_output = processing.run(
+            "native:extractbylocation",
+            {
+                "INPUT": grid_layer,
+                "PREDICATE": [0],
+                "INTERSECT": self.input_polygons,
+                "OUTPUT": "TEMPORARY_OUTPUT",
+            },
+            feedback=QgsProcessingFeedback(),
+        )["OUTPUT"]
+
+        grid_layer = grid_output
 
         # Add score field
         provider = grid_layer.dataProvider()
@@ -114,26 +161,19 @@ def raster_polygon_grid_score(self):
 
         merged_output_vector = os.path.join(output_dir, "merged_grid_vector.shp")
 
-        Merge = processing.run(
+        # Merge the output vector layers
+        merge = processing.run(
             "native:mergevectorlayers",
-            {"LAYERS": [grid_layer], "CRS": None, "OUTPUT": "memory:"},
+            {"LAYERS": [grid_layer], "CRS": self.crs, "OUTPUT": "TEMPORARY_OUTPUT"},
             feedback=QgsProcessingFeedback(),
-        )
-
-        merge = Merge["OUTPUT"]
-
-        extents_processor = Extents(
-            output_dir, self.country_boundary, self.pixel_size, self.crs
-        )
+        )["OUTPUT"]
 
-        # Get the extent of the vector layer
-        country_extent = extents_processor.get_country_extent()
-        xmin, ymin, xmax, ymax = (
-            country_extent.xMinimum(),
-            country_extent.yMinimum(),
-            country_extent.xMaximum(),
-            country_extent.yMaximum(),
-        )
+        xmin, xmax, ymin, ymax = (
+            aligned_bbox.xMinimum(),
+            aligned_bbox.xMaximum(),
+            aligned_bbox.yMinimum(),
+            aligned_bbox.yMaximum(),
+        )  # Extent of the aligned bbox
 
         # Rasterize the clipped grid layer to generate the raster
         rasterize_params = {
@@ -144,13 +184,25 @@ def raster_polygon_grid_score(self):
             "UNITS": 1,
             "WIDTH": self.pixel_size,
             "HEIGHT": self.pixel_size,
-            "EXTENT": f"{xmin},{xmax},{ymin},{ymax}",
-            "NODATA": -9999,
+            "EXTENT": f"{xmin},{ymin},{xmax},{ymax}",
+            "NODATA": None,
             "OPTIONS": "",
             "DATA_TYPE": 5,  # Use Int32 for scores
-            "OUTPUT": self.output_path,
+            "OUTPUT": "TEMPORARY_OUTPUT",
         }
 
-        processing.run(
+        output_file = processing.run(
             "gdal:rasterize", rasterize_params, feedback=QgsProcessingFeedback()
+        )["OUTPUT"]
+
+        processing.run(
+            "gdal:cliprasterbymasklayer",
+            {
+                "INPUT": output_file,
+                "MASK": self.country_boundary,
+                "NODATA": -9999,
+                "CROP_TO_CUTLINE": True,
+                "OUTPUT": self.output_path,
+            },
+            feedback=QgsProcessingFeedback(),
         )
diff --git a/geest/core/polylines_per_grid_cell.py b/geest/core/polylines_per_grid_cell.py
@@ -1,44 +1,90 @@
 import os
 from qgis.PyQt.QtCore import QVariant
 from qgis.core import (
+    QgsGeometry,
     QgsVectorLayer,
     QgsField,
     QgsSpatialIndex,
     QgsProcessingFeedback,
 )
 import processing
-from .create_grids import GridCreator
-from .extents import Extents
+from .utilities import GridAligner
 
 
 class RasterPolylineGridScore:
-    def __init__(self, country_boundary, pixel_size, output_path, crs, input_polylines):
+    def __init__(
+        self,
+        country_boundary,
+        pixel_size,
+        working_dir,
+        crs,
+        input_polylines,
+        output_path,
+    ):
         self.country_boundary = country_boundary
         self.pixel_size = pixel_size
-        self.output_path = output_path
+        self.working_dir = working_dir
         self.crs = crs
         self.input_polylines = input_polylines
+        self.output_path = output_path
+        # Initialize GridAligner with grid size
+        self.grid_aligner = GridAligner(grid_size=100)
 
     def raster_polyline_grid_score(self):
         """
         Generates a raster based on the number of input points within each grid cell.
         :param country_boundary: Layer defining the country boundary to clip the grid.
         :param cellsize: The size of each grid cell.
-        :param output_path: Path to save the output raster.
         :param crs: The CRS in which the grid and raster will be projected.
         :param input_polylines: Layer of point features to count within each grid cell.
         """
 
-        # Create grid
-        self.h_spacing = 100
-        self.v_spacing = 100
-        create_grid = GridCreator(h_spacing=self.h_spacing, v_spacing=self.v_spacing)
-        output_dir = os.path.join("output")
-        merged_output_path = os.path.join(output_dir, "merged_grid.shp")
-        grid_layer = create_grid.create_grids(
-            self.country_boundary, output_dir, self.crs, merged_output_path
+        output_dir = os.path.dirname(self.output_path)
+
+        # Define output directory and ensure it's created
+        os.makedirs(output_dir, exist_ok=True)
+
+        # Load grid layer from the Geopackage
+        geopackage_path = os.path.join(
+            self.working_dir, "study_area", "study_area.gpkg"
+        )
+        if not os.path.exists(geopackage_path):
+            raise ValueError(f"Geopackage not found at {geopackage_path}.")
+
+        grid_layer = QgsVectorLayer(
+            f"{geopackage_path}|layername=study_area_grid", "merged_grid", "ogr"
+        )
+
+        area_layer = QgsVectorLayer(
+            f"{geopackage_path}|layername=study_area_polygons",
+            "study_area_polygons",
+            "ogr",
         )
-        grid_layer = QgsVectorLayer(merged_output_path, "merged_grid", "ogr")
+
+        geometries = [feature.geometry() for feature in area_layer.getFeatures()]
+
+        # Combine all geometries into one using unaryUnion
+        area_geometry = QgsGeometry.unaryUnion(geometries)
+
+        # grid_geometry = grid_layer.getGeometry()
+
+        aligned_bbox = self.grid_aligner.align_bbox(
+            area_geometry.boundingBox(), area_layer.extent()
+        )
+
+        # Extract polylines by location
+        grid_output = processing.run(
+            "native:extractbylocation",
+            {
+                "INPUT": grid_layer,
+                "PREDICATE": [0],
+                "INTERSECT": self.input_polylines,
+                "OUTPUT": "TEMPORARY_OUTPUT",
+            },
+            feedback=QgsProcessingFeedback(),
+        )["OUTPUT"]
+
+        grid_layer = grid_output
 
         # Add score field
         provider = grid_layer.dataProvider()
@@ -48,7 +94,6 @@ def raster_polyline_grid_score(self):
             grid_layer.updateFields()
 
         # Create spatial index for the input points
-        # Reproject the country layer if necessary
         if self.input_polylines.crs() != self.crs:
             self.input_polylines = processing.run(
                 "native:reprojectlayer",
@@ -83,16 +128,10 @@ def raster_polyline_grid_score(self):
             num_polylines = len(unique_intersections)
 
             # Reclassification logic: assign score based on the number of points
-            if num_polylines >= 2:
-                reclass_val = 5
-            elif num_polylines == 1:
-                reclass_val = 3
-            else:
-                reclass_val = 0
-
+            reclass_val = 5 if num_polylines >= 2 else 3 if num_polylines == 1 else 0
             reclass_vals[grid_feat.id()] = reclass_val
 
-        # Step 5: Apply the score values to the grid
+        # Apply the score values to the grid
         grid_layer.startEditing()
         for grid_feat in grid_layer.getFeatures():
             grid_layer.changeAttributeValue(
@@ -102,30 +141,22 @@ def raster_polyline_grid_score(self):
             )
         grid_layer.commitChanges()
 
-        merged_output_vector = os.path.join(output_dir, "merged_grid_vector.shp")
-
-        Merge = processing.run(
+        # Merge the output vector layers
+        merge = processing.run(
             "native:mergevectorlayers",
-            {"LAYERS": [grid_layer], "CRS": None, "OUTPUT": "memory:"},
+            {"LAYERS": [grid_layer], "CRS": self.crs, "OUTPUT": "TEMPORARY_OUTPUT"},
             feedback=QgsProcessingFeedback(),
-        )
-
-        merge = Merge["OUTPUT"]
+        )["OUTPUT"]
 
-        extents_processor = Extents(
-            output_dir, self.country_boundary, self.pixel_size, self.crs
-        )
-
-        # Get the extent of the vector layer
-        country_extent = extents_processor.get_country_extent()
-        xmin, ymin, xmax, ymax = (
-            country_extent.xMinimum(),
-            country_extent.yMinimum(),
-            country_extent.xMaximum(),
-            country_extent.yMaximum(),
-        )
+        xmin, xmax, ymin, ymax = (
+            aligned_bbox.xMinimum(),
+            aligned_bbox.xMaximum(),
+            aligned_bbox.yMinimum(),
+            aligned_bbox.yMaximum(),
+        )  # Extent of the aligned bbox
 
         # Rasterize the clipped grid layer to generate the raster
+        # output_file = os.path.join(output_dir, "rasterized_grid.tif")
         rasterize_params = {
             "INPUT": merge,
             "FIELD": field_name,
@@ -134,13 +165,25 @@ def raster_polyline_grid_score(self):
             "UNITS": 1,
             "WIDTH": self.pixel_size,
             "HEIGHT": self.pixel_size,
-            "EXTENT": f"{xmin},{xmax},{ymin},{ymax}",
-            "NODATA": -9999,
+            "EXTENT": f"{xmin},{ymin},{xmax},{ymax}",
+            "NODATA": None,
             "OPTIONS": "",
             "DATA_TYPE": 5,  # Use Int32 for scores
-            "OUTPUT": self.output_path,
+            "OUTPUT": "TEMPORARY_OUTPUT",
         }
 
-        processing.run(
+        output_file = processing.run(
             "gdal:rasterize", rasterize_params, feedback=QgsProcessingFeedback()
+        )["OUTPUT"]
+
+        processing.run(
+            "gdal:cliprasterbymasklayer",
+            {
+                "INPUT": output_file,
+                "MASK": self.country_boundary,
+                "NODATA": -9999,
+                "CROP_TO_CUTLINE": True,
+                "OUTPUT": self.output_path,
+            },
+            feedback=QgsProcessingFeedback(),
         )
diff --git a/test/test_data/study_area/combined_mask.vrt b/test/test_data/study_area/combined_mask.vrt
@@ -0,0 +1,20 @@
+<VRTDataset rasterXSize="230" rasterYSize="450">
+  <SRS dataAxisToSRSAxisMapping="1,2">PROJCS["WGS 84 / UTM zone 20N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-63],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32620"]]</SRS>
+  <GeoTransform>  7.0720000000000000e+05,  1.0000000000000000e+02,  0.0000000000000000e+00,  1.5609000000000000e+06,  0.0000000000000000e+00, -1.0000000000000000e+02</GeoTransform>
+  <VRTRasterBand dataType="Byte" band="1">
+    <NoDataValue>0</NoDataValue>
+    <ColorInterp>Palette</ColorInterp>
+    <ColorTable>
+      <Entry c1="0" c2="0" c3="0" c4="255" />
+      <Entry c1="255" c2="255" c3="255" c4="255" />
+    </ColorTable>
+    <ComplexSource resampling="nearest">
+      <SourceFilename relativeToVRT="1">saint_lucia_part0.tif</SourceFilename>
+      <SourceBand>1</SourceBand>
+      <SourceProperties RasterXSize="230" RasterYSize="450" DataType="Byte" BlockXSize="230" BlockYSize="282" />
+      <SrcRect xOff="0" yOff="0" xSize="230" ySize="450" />
+      <DstRect xOff="0" yOff="0" xSize="230" ySize="450" />
+      <NODATA>0</NODATA>
+    </ComplexSource>
+  </VRTRasterBand>
+</VRTDataset>
diff --git a/test/test_data/study_area/saint_lucia_part0.tif b/test/test_data/study_area/saint_lucia_part0.tif
diff --git a/test/test_data/study_area/study_area.gpkg b/test/test_data/study_area/study_area.gpkg