From 594d7f57bff32b6efe8cab7b848bd16605fe0d83 Mon Sep 17 00:00:00 2001 From: Jeff Osundwa Date: Thu, 10 Oct 2024 01:48:58 +0300 Subject: [PATCH] Rasterize polygon functionality Fixes #390 --- geest/core/polygons_per_grid_cell.py | 120 +++++++++++++++++++-------- test/test_polygons_per_grid_cell.py | 7 +- 2 files changed, 90 insertions(+), 37 deletions(-) diff --git a/geest/core/polygons_per_grid_cell.py b/geest/core/polygons_per_grid_cell.py index 7467a12..c9c3d18 100644 --- 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/test/test_polygons_per_grid_cell.py b/test/test_polygons_per_grid_cell.py index a11444f..c85e455 100644 --- a/test/test_polygons_per_grid_cell.py +++ b/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