Rasterize polygon functionality

Fixes #390

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(),
  )

test/test_polygons_per_grid_cell.py

@@ -23,11 +23,11 @@ def test_raster_polygon_grid_score(self):
 
  # Load the input data (polygons and country boundary layers)
  self.polygon_layer = QgsVectorLayer(
- os.path.join(self.test_data_dir, "polygons/blocks.shp"),
+ os.path.join(self.test_data_dir, "polygons", "blocks.shp"),
  "test_polygons",
  "ogr",
  )
- self.country_boundary = os.path.join(self.test_data_dir, "admin/Admin0.shp")
+ self.country_boundary = os.path.join(self.test_data_dir, "admin", "Admin0.shp")
 
  self.assertTrue(self.polygon_layer.isValid(), "The polygon layer is not valid.")
 
@@ -45,9 +45,10 @@ def test_raster_polygon_grid_score(self):
  rasterizer = RasterPolygonGridScore(
  country_boundary=self.country_boundary,
  pixel_size=self.pixel_size,
- output_path=self.output_path,
+ working_dir=self.test_data_dir,
  crs=self.crs,
  input_polygons=self.polygon_layer,
+ output_path=self.output_path,
  )
 
  # Run the raster_polygon_grid_score method