diff --git a/config.json b/config.json index e3e1bf23..5a14ce36 100644 --- a/config.json +++ b/config.json @@ -16,7 +16,7 @@ "author": "Kartoza", "email": "info@kartoza.com", "description": "Gender Enabling Environments Spatial Tool", - "version": "0.5.0", + "version": "0.5.1", "changelog": "", "server": false } diff --git a/geest/__init__.py b/geest/__init__.py index c9cb01d9..f5b46b49 100644 --- a/geest/__init__.py +++ b/geest/__init__.py @@ -87,7 +87,7 @@ def initGui(self): # pylint: disable=missing-function-docstring self.run_action = QAction(icon, "GEEST Settings", self.iface.mainWindow()) self.run_action.triggered.connect(self.run) self.iface.addToolBarIcon(self.run_action) - + self.debug_running = False # Create the dock widget self.dock_widget = GeestDock( parent=self.iface.mainWindow(), @@ -247,6 +247,8 @@ def unload(self): # pylint: disable=missing-function-docstring Unload the plugin from QGIS. Removes all added actions, widgets, and options to ensure a clean unload. """ + + self.kill_debug() # Save geometry before unloading self.save_geometry() @@ -286,11 +288,24 @@ def unload(self): # pylint: disable=missing-function-docstring self.dock_widget.deleteLater() self.dock_widget = None + def kill_debug(self): + # Note that even though this kills the debugpy process I still + # cannot successfully restart the debugger in vscode without restarting QGIS + if self.debug_running: + import psutil + import signal + + """Find the PID of the process listening on the specified port.""" + for conn in psutil.net_connections(kind="tcp"): + if conn.laddr.port == 9000 and conn.status == psutil.CONN_LISTEN: + os.kill(conn.pid, signal.SIGTERM) + def debug(self): """ Enters debug mode. Shows a message to attach a debugger to the process. """ + self.kill_debug() self.display_information_message_box( title="GEEST", message="Close this dialog then open VSCode and start your debug client.", @@ -306,6 +321,8 @@ def debug(self): title="GEEST", message="Visual Studio Code debugger is now attached on port 9000", ) + self.debug_action.setEnabled(False) # prevent user starting it twice + self.debug_running = True def run(self): """ diff --git a/geest/core/algorithms/area_iterator.py b/geest/core/algorithms/area_iterator.py index 020801fc..3b80b23c 100644 --- a/geest/core/algorithms/area_iterator.py +++ b/geest/core/algorithms/area_iterator.py @@ -127,7 +127,12 @@ def __iter__(self) -> Iterator[Tuple[QgsGeometry, QgsGeometry, float]]: return # Iterate over each polygon feature and calculate progress - for index, polygon_feature in enumerate(self.polygon_layer.getFeatures()): + + # Sort polygon features by area in ascending order + sorted_features = sorted( + self.polygon_layer.getFeatures(), key=lambda f: f.geometry().area() + ) + for index, polygon_feature in enumerate(sorted_features): polygon_id: int = polygon_feature.id() # Request the corresponding bbox feature based on the polygon's ID feature_request: QgsFeatureRequest = QgsFeatureRequest().setFilterFid( diff --git a/geest/core/tasks/__init__.py b/geest/core/tasks/__init__.py index b0ba8ce3..d1061167 100644 --- a/geest/core/tasks/__init__.py +++ b/geest/core/tasks/__init__.py @@ -1,2 +1,3 @@ from .study_area import StudyAreaProcessingTask from .ors_checker import OrsCheckerTask +from .grid_from_bbox import GridFromBbox diff --git a/geest/core/tasks/grid_from_bbox.py b/geest/core/tasks/grid_from_bbox.py new file mode 100644 index 00000000..5eea98ef --- /dev/null +++ b/geest/core/tasks/grid_from_bbox.py @@ -0,0 +1,109 @@ +import time +from osgeo import ogr +from qgis.core import QgsTask +from geest.utilities import log_message + + +class GridFromBbox(QgsTask): + """ + A QGIS task to generate grid cells in a bounding box chunk, check intersections, + and store them in memory for later writing. + """ + + def __init__(self, chunk_id, bbox_chunk, geom, cell_size, feedback): + super().__init__(f"CreateGridChunkTask-{chunk_id}", QgsTask.CanCancel) + self.chunk_id = chunk_id + self.bbox_chunk = bbox_chunk # (x_start, x_end, y_start, y_end) + self.geom = geom + self.cell_size = cell_size + self.feedback = feedback + self.run_time = 0.0 + self.features_out = [] # store geometries here + + def run(self): + log_message(f"##################################") + log_message(f"Processing chunk {self.chunk_id}...") + log_message(f"Chunk bbox: {self.bbox_chunk}") + log_message(f"##################################") + start_time = time.time() + x_start, x_end, y_start, y_end = self.bbox_chunk + + # Convert geom to OGR if needed + # If self.geom is an ogr.Geometry, skip this step + # If it's a PyQGIS geometry, convert to WKB or so, e.g.: + # ogr_geom = ogr.CreateGeometryFromWkb(self.geom.asWkb()) + + # We'll assume self.geom is already an ogr.Geometry + skip_intersection_check = False + + chunk_bounds = ogr.CreateGeometryFromWkt( + f"POLYGON(({x_start} {y_start}, {x_end} {y_start}, " + f"{x_end} {y_end}, {x_start} {y_end}, {x_start} {y_start}))" + ) + + # If chunk bounding box is fully outside of the geometry, skip chunk + if not self.geom.Intersects(chunk_bounds): + log_message( + f"Chunk {self.chunk_id} is completely outside geometry, skipping." + ) + return True + + # If the whole chunk is within the geometry, skip intersection check + if self.geom.Contains(chunk_bounds): + log_message( + f"Whole chunk is within the geometry, we will skip check intersection for each feature..." + ) + skip_intersection_check = True + else: + log_message( + f"Whole chunk is NOT within the geometry, we will check intersection for each feature..." + ) + + x = x_start + while x < x_end: + x2 = x + self.cell_size + if x2 <= x: + break + y = y_start + while y < y_end: + y2 = y + self.cell_size + if y2 <= y: + break + # Create cell polygon in memory + ring = ogr.Geometry(ogr.wkbLinearRing) + ring.AddPoint(x, y) + ring.AddPoint(x, y2) + ring.AddPoint(x2, y2) + ring.AddPoint(x2, y) + ring.AddPoint(x, y) + + cell_polygon = ogr.Geometry(ogr.wkbPolygon) + cell_polygon.AddGeometry(ring) + + # Check intersection + if not skip_intersection_check: + if self.geom.Intersects(cell_polygon): + self.features_out.append(cell_polygon) + else: + self.features_out.append(cell_polygon) + + y = y2 + x = x2 + + end_time = time.time() + self.run_time = end_time - start_time + # self.feedback.pushInfo( + log_message( + f"Chunk {self.chunk_id} processed in {end_time - start_time:.2f} s; created {len(self.features_out)} features." + ) + + return True + + def finished(self, result): + # This is called in the main thread after `run` completes + # We do *not* write to the data source here if we want to avoid concurrency issues. + pass + + def cancel(self): + super().cancel() + # clean up if needed diff --git a/geest/core/tasks/study_area.py b/geest/core/tasks/study_area.py index aa457173..3882c601 100644 --- a/geest/core/tasks/study_area.py +++ b/geest/core/tasks/study_area.py @@ -3,38 +3,31 @@ import glob import traceback import datetime +import time +# GDAL / OGR / OSR imports +from osgeo import ogr, osr, gdal from typing import List, Optional from qgis.core import ( QgsTask, - QgsRectangle, QgsFeedback, - QgsFeature, - QgsGeometry, - QgsField, - QgsPointXY, - QgsProcessingContext, - QgsCoordinateTransform, - QgsCoordinateReferenceSystem, - QgsRasterLayer, - QgsSpatialIndex, - QgsWkbTypes, QgsVectorLayer, QgsVectorFileWriter, - QgsFields, - QgsCoordinateTransformContext, - QgsWkbTypes, - Qgis, + QgsApplication, ) -from qgis.PyQt.QtCore import QVariant -import processing # QGIS processing toolbox from geest.utilities import log_message +from .grid_from_bbox import GridFromBbox +from geest.core import setting class StudyAreaProcessingTask(QgsTask): """ A QgsTask subclass for processing study area features. + Processes study-area geometries using GDAL/OGR instead of QGIS API. + It creates bounding boxes, grids, raster masks, a dissolved clip polygon, + and a combined VRT of all masks. + It works through the (multi)part geometries in the input layer, creating bounding boxes and masks. The masks are stored as individual tif files and then a vrt file is created to combine them. The grids are in two forms - the entire bounding box and the individual parts. @@ -42,7 +35,12 @@ class StudyAreaProcessingTask(QgsTask): Any invalid geometries are discarded, and fixed geometries are processed. Args: - QgsTask (_type_): _description_ + layer (QgsVectorLayer): The input vector layer. + field_name (str): The field name in the input layer that holds the study area name. + cell_size_m (float): The cell size for grid spacing in meters. + working_dir (str): The directory path where outputs will be saved. + feedback (QgsFeedback): A feedback object to report progress. + crs (Optional[QgsCoordinateReferenceSystem]): The target CRS. If None, a UTM zone will be computed. Returns: _type_: _description_ @@ -51,920 +49,982 @@ class StudyAreaProcessingTask(QgsTask): def __init__( self, layer: QgsVectorLayer, - field_name: str, - cell_size_m: float, - working_dir: str, - mode: str = "raster", - crs: Optional[QgsCoordinateReferenceSystem] = None, - context: QgsProcessingContext = None, + field_name, + cell_size_m, + working_dir, feedback: QgsFeedback = None, + crs=None, ): """ - Initializes the StudyAreaProcessingTask. - - :param layer: The vector layer containing study area features. - :param field_name: The name of the field containing area names. - :param cell_size: The size of the grid cells in meters. + :param input_vector_path: Path to an OGR-readable vector file (e.g. .gpkg or .shp). + :param field_name: Name of the field that holds the study area name. + :param cell_size_m: Cell size for grid spacing (in meters). :param working_dir: Directory path where outputs will be saved. - :param mode: Processing mode, either 'vector' or 'raster'. Default is raster. - :param crs: Optional CRS for the output CRS. If None, a UTM zone - is calculated based on layer extent or extent of selected features. - :param context: QgsProcessingContext for the task. Default is None. - :param feedback: QgsFeedback object for task cancelling. Default is None. + :param crs_epsg: EPSG code for target CRS. If None, a UTM zone will be computed. """ super().__init__("Study Area Preparation", QgsTask.CanCancel) + self.input_vector_path = self.export_qgs_layer_to_shapefile(layer, working_dir) + self.field_name = field_name + self.cell_size_m = cell_size_m + self.working_dir = working_dir + self.gpkg_path = os.path.join(working_dir, "study_area", "study_area.gpkg") + self.counter = 0 self.feedback = feedback - self.layer: QgsVectorLayer = layer - self.field_name: str = field_name - self.cell_size_m: float = cell_size_m - self.working_dir: str = working_dir - self.mode: str = mode - self.context: QgsProcessingContext = context - self.gpkg_path: str = os.path.join( - self.working_dir, "study_area", "study_area.gpkg" - ) - self.counter: int = 0 - # Remove the GeoPackage if it already exists to start with a clean state + self.metrics = { + "Creating chunks": 0.0, + "Writing chunks": 0.0, + "Complete chunk": 0.0, + "Preparing chunks": 0.0, + } + self.valid_feature_count = 0 + self.current_geom_actual_cell_count = 0 + self.current_geom_cell_count_estimate = 0 + self.total_cells = 0 + self.write_lock = False + # Make sure output directory exists + self.create_study_area_directory(self.working_dir) + + # If GPKG already exists, remove it to start fresh if os.path.exists(self.gpkg_path): try: os.remove(self.gpkg_path) - log_message( - f"Existing GeoPackage removed: {self.gpkg_path}", - tag="Geest", - level=Qgis.Info, - ) + log_message(f"Removed existing GeoPackage: {self.gpkg_path}") except Exception as e: log_message( - f"Error removing existing GeoPackage: {e}", - tag="Geest", - level=Qgis.Critical, + f"Error removing existing GeoPackage: {e}", level="CRITICAL" ) - self.create_study_area_directory(self.working_dir) - - # Determine the correct CRS and layer bounding box based on selected features - selected_features = self.layer.selectedFeatures() - self.layer_bbox = QgsRectangle() - - if selected_features: - # Calculate bounding box based on selected features only - self.layer_bbox: QgsRectangle = QgsRectangle() - for feature in selected_features: - self.layer_bbox.combineExtentWith(feature.geometry().boundingBox()) - else: - # No features selected, use full layer extent - self.layer_bbox: QgsRectangle = self.layer.extent() + # Open the source data using OGR + self.source_ds = ogr.Open(self.input_vector_path, 0) # 0 = read-only + if not self.source_ds: + raise RuntimeError(f"Could not open {self.input_vector_path} with OGR.") + self.source_layer = self.source_ds.GetLayer(0) + if not self.source_layer: + raise RuntimeError("Could not retrieve layer from the data source.") + self.parts_count = self.count_layer_parts() + + # Determine source EPSG (if any) by reading layer's spatial ref + self.src_spatial_ref = self.source_layer.GetSpatialRef() + self.src_epsg = None + if self.src_spatial_ref: + self.src_epsg = self.src_spatial_ref.GetAuthorityCode(None) + + # Compute bounding box from entire layer + # (OGR Envelope: (xmin, xmax, ymin, ymax)) + layer_extent = self.source_layer.GetExtent() + (xmin, xmax, ymin, ymax) = layer_extent + self.layer_bbox = (xmin, xmax, ymin, ymax) - # Determine EPSG code based on provided input or calculated UTM zone if crs is None: - self.epsg_code: int = self.calculate_utm_zone(self.layer_bbox) - self.output_crs: QgsCoordinateReferenceSystem = ( - QgsCoordinateReferenceSystem(f"EPSG:{self.epsg_code}") + # Attempt to pick a suitable UTM zone + self.epsg_code = self.calculate_utm_zone(self.layer_bbox, self.src_epsg) + else: + auth_id = crs.authid() # e.g. "EPSG:4326" + if auth_id.lower().startswith("epsg:"): + epsg_int = int(auth_id.split(":")[1]) + log_message(f"EPSG code is: {epsg_int}") + else: + # Handle case where it's not an EPSG-based CRS + epsg_int = None + raise Exception("CRS is not an EPSG-based ID.") + self.epsg_code = epsg_int + + # Prepare OSR objects for source->target transformation + self.target_spatial_ref = osr.SpatialReference() + self.target_spatial_ref.ImportFromEPSG(self.epsg_code) + + if self.src_spatial_ref: + self.coord_transform = osr.CoordinateTransformation( + self.src_spatial_ref, self.target_spatial_ref ) else: - self.output_crs = crs + self.coord_transform = None + + log_message(f"Using output EPSG:{self.epsg_code}") + # Create aligned bounding box in target CRS space + # We interpret the layer bbox in source CRS (if it has one), transform, and align + self.transformed_layer_bbox = self.transform_and_align_bbox(self.layer_bbox) log_message( - f"Project CRS Set to: {self.output_crs.authid()}", - tag="Geest", - level=Qgis.Info, + f"Transformed layer bbox to target CRS and aligned to grid: {self.transformed_layer_bbox}" ) - # Reproject and align the transformed layer_bbox to a cell_size_m grid and output crs - self.layer_bbox = self.grid_aligned_bbox(self.layer_bbox) + # Tracking table name + self.status_table_name = "study_area_creation_status" - def run(self) -> bool: + def count_layer_parts(self): """ - Runs the task in the background. - """ - try: - result = self.process_study_area() - return result # false if the task was canceled - except Exception as e: - log_message(f"Task failed: {e}", level=Qgis.Critical) - # Print the traceback to the QgsMessageLog - log_message(traceback.format_exc(), level=Qgis.Critical) - # And write it to a text file called error.txt in the working directory - with open(os.path.join(self.working_dir, "error.txt"), "w") as f: - # first the date and time - f.write(f"{datetime.datetime.now()}\n") - # then the traceback - f.write(traceback.format_exc()) - return False + Returns the number of parts in the layer. - def finished(self, result: bool) -> None: + :return: The number of parts in the layer. """ - Called when the task completes successfully. - """ - if result: - log_message( - "Study Area Processing completed successfully.", - tag="Geest", - level=Qgis.Info, - ) - else: - if self.feedback.isCanceled(): - log_message( - "Study Area Processing was canceled by the user.", - tag="Geest", - level=Qgis.Warning, - ) + self.source_layer.ResetReading() + parts_count = 0 + for feature in self.source_layer: + geom = feature.GetGeometryRef() + if not geom: + continue + geometry_type = geom.GetGeometryName() + if geometry_type == "MULTIPOLYGON": + parts_count += geom.GetGeometryCount() else: - log_message( - "Study Area Processing encountered an error.", - tag="Geest", - level=Qgis.Critical, - ) - - def cancel(self) -> None: - """ - Called when the task is canceled. - """ - super().cancel() - self.feedback.cancel() - log_message("Study Area Processing was canceled.", level=Qgis.Warning) + return self.source_layer.GetFeatureCount() + return parts_count - def process_study_area(self) -> None: + def export_qgs_layer_to_shapefile(self, layer, output_dir): """ - Processes each feature in the input layer, creating bounding boxes and grids. - It handles the CRS transformation and calls appropriate processing functions based on geometry type. + Exports a QgsVectorLayer to a Shapefile in output_dir. + Returns the full path to the .shp (main file). """ - # First, write the layer_bbox to the GeoPackage - self.save_to_geopackage( - layer_name="study_area_bbox", - geom=QgsGeometry.fromRect(self.layer_bbox), - area_name="Study Area Bounding Box", + shapefile_path = os.path.join(output_dir, "temp_export.shp") + # Check if the layer has selected features + if layer.selectedFeatureCount() == 0: + has_selection = False + else: + has_selection = True + # Export the selected features of the layer + err_code, err_msg = QgsVectorFileWriter.writeAsVectorFormat( + layer, + shapefile_path, + "UTF-8", # encoding + layer.crs(), # CRS to write + "ESRI Shapefile", # driver name + onlySelected=has_selection, # Export only selected features ) - # Initialize counters for tracking valid and invalid features - invalid_feature_count = 0 - valid_feature_count = 0 - fixed_feature_count = 0 + if err_code != QgsVectorFileWriter.NoError: + raise RuntimeError(f"Failed to export layer to Shapefile: {err_msg}") - # Process individual features - selected_features = self.layer.selectedFeatures() - # count how many features are selected - total_features = len(selected_features) + return shapefile_path - features = selected_features if selected_features else self.layer.getFeatures() + def run(self): + """ + Main entry point (mimics process_study_area from QGIS code). + """ + try: + # 1) Create the bounding box as a single polygon feature + # and save to GeoPackage + self.save_bbox_polygon( + "study_area_bbox", + self.transformed_layer_bbox, + "Study Area Bounding Box", + ) - for feature in features: - if self.feedback.isCanceled(): - return False - geom: QgsGeometry = feature.geometry() - area_name: str = feature[self.field_name] - normalized_name: str = re.sub(r"\s+", "_", area_name.lower()) + # 2) Create the status tracking table + self.create_status_tracking_table() + + # 3) Iterate over features + invalid_feature_count = 0 + self.valid_feature_count = 0 + fixed_feature_count = 0 + self.source_layer.ResetReading() + for feature in self.source_layer: + geom_ref = feature.GetGeometryRef() + if not geom_ref: + continue - # Check for geometry validity - if geom.isEmpty() or not geom.isGeosValid(): - log_message( - f"Feature ID {feature.id()} has an invalid geometry. Attempting to fix.", - tag="Geest", - level=Qgis.Warning, - ) + area_name = feature.GetField(self.field_name) + if not area_name: + area_name = f"area_{feature.GetFID()}" - # Try to fix the geometry - fixed_geom = geom.makeValid() + # Clean up the name + normalized_name = re.sub(r"\s+", "_", area_name.lower()) - # Check if the fixed geometry is valid - if fixed_geom.isEmpty() or not fixed_geom.isGeosValid(): - invalid_feature_count += 1 + # Check validity + if not geom_ref.IsValid(): + # Attempt a fix log_message( - f"Feature ID {feature.id()} could not be fixed. Skipping.", - tag="Geest", - level=Qgis.Critical, + f"Feature {feature.GetFID()} has invalid geometry, attempting to fix." + ) + # OGR >= 3.0 has MakeValid. If unavailable, Buffer(0) can fix many invalid polygons. + try: + geom_ref = geom_ref.MakeValid() + # geom_ref = geom_ref.Buffer(0) + if not geom_ref.IsValid(): + invalid_feature_count += 1 + log_message( + f"Could not fix geometry for feature {feature.GetFID()}. Skipping.", + level="CRITICAL", + ) + continue + else: + fixed_feature_count += 1 + except Exception as e: + invalid_feature_count += 1 + log_message(f"Geometry fix error: {str(e)}", level="CRITICAL") + continue + + self.valid_feature_count += 1 + + # Singlepart vs multipart + # OGR geometry can be geometry collection if multipart + geom_type = ogr.GT_Flatten(geom_ref.GetGeometryType()) + if geom_type == ogr.wkbMultiPolygon: + log_message(f"Processing multipart geometry: {normalized_name}") + self.process_multipart_geometry( + geom_ref, normalized_name, area_name ) - continue - - # Use the fixed geometry if it is valid - geom = fixed_geom - fixed_feature_count += 1 - log_message( - f"Feature ID {feature.id()} geometry fixed successfully.", - tag="Geest", - level=Qgis.Info, - ) - - # Process valid geometry - try: - valid_feature_count += 1 - if geom.isMultipart(): - self.process_multipart_geometry(geom, normalized_name, area_name) else: - self.process_singlepart_geometry(geom, normalized_name, area_name) - except Exception as e: - # Log any unexpected errors during processing - invalid_feature_count += 1 - log_message( - f"Error processing geometry for feature {feature.id()}: {e}", - tag="Geest", - level=Qgis.Critical, - ) - - # Update progress (for example, as percentage) - progress = int((valid_feature_count / self.layer.featureCount()) * 100) - self.setProgress(progress) + log_message(f"Processing singlepart geometry: {normalized_name}") + self.process_singlepart_geometry( + geom_ref, normalized_name, area_name + ) + log_message( + f"Processing complete. Valid: {self.valid_feature_count}, Fixed: {fixed_feature_count}, Invalid: {invalid_feature_count}" + ) + log_message(f"Total cells generated: {self.total_cells}") - # Log the count of valid, fixed, and invalid features processed - log_message( - f"Processing completed. Valid features: {valid_feature_count}, Fixed features: {fixed_feature_count}, Invalid features: {invalid_feature_count}", - tag="Geest", - level=Qgis.Info, - ) - # Create and add the VRT of all generated raster masks if in raster mode - if self.mode == "raster": + # 4) Create a VRT of all generated raster masks self.create_raster_vrt() + + except Exception as e: + log_message(f"Error in run(): {str(e)}") + log_message(traceback.format_exc()) + with open(os.path.join(self.working_dir, "error.txt"), "w") as f: + f.write(f"{datetime.datetime.now()}\n") + f.write(traceback.format_exc()) + return False + return True - def process_singlepart_geometry( - self, geom: QgsGeometry, normalized_name: str, area_name: str - ) -> None: + ########################################################################## + # Table creation logic + ########################################################################## + def create_status_tracking_table(self): """ - Processes a singlepart geometry feature. Creates vector grids or raster masks based on mode. - - :param geom: Geometry of the feature. - :param normalized_name: Name normalized for file storage. - :param area_name: Name of the study area. + Create a table in the GeoPackage to track processing status, + similar to the QGIS version. """ + if not os.path.exists(self.gpkg_path): + # Just create it if no GPKG + driver = ogr.GetDriverByName("GPKG") + ds = driver.CreateDataSource(self.gpkg_path) + if not ds: + raise RuntimeError(f"Could not create GeoPackage {self.gpkg_path}") + ds = None # Close + + # Check if table exists + ds = ogr.Open(self.gpkg_path, 1) # open in update mode + if not ds: + raise RuntimeError(f"Could not open or create {self.gpkg_path} for update.") + try: + # If a layer with status_table_name exists, do nothing + layer = ds.GetLayerByName(self.status_table_name) + if layer: + log_message(f"Table '{self.status_table_name}' already exists.") + return + + # Otherwise, create it + srs = osr.SpatialReference() + srs.ImportFromEPSG(4326) # Arbitrary SRS for table with no geometry + + layer = ds.CreateLayer(self.status_table_name, srs, geom_type=ogr.wkbNone) + layer.CreateField(ogr.FieldDefn("area_name", ogr.OFTString)) + layer.CreateField(ogr.FieldDefn("timestamp_start", ogr.OFTDateTime)) + layer.CreateField(ogr.FieldDefn("timestamp_end", ogr.OFTDateTime)) + layer.CreateField(ogr.FieldDefn("geometry_processed", ogr.OFTInteger)) + layer.CreateField(ogr.FieldDefn("clip_geometry_processed", ogr.OFTInteger)) + layer.CreateField(ogr.FieldDefn("grid_processed", ogr.OFTInteger)) + layer.CreateField(ogr.FieldDefn("mask_processed", ogr.OFTInteger)) + layer.CreateField(ogr.FieldDefn("grid_creation_duration_secs", ogr.OFTReal)) + layer.CreateField( + ogr.FieldDefn("clip_geom_creation_duration_secs", ogr.OFTReal) + ) + layer.CreateField(ogr.FieldDefn("geom_total_duration_secs", ogr.OFTReal)) - # Compute the aligned bounding box based on the transformed geometry - # This will do the CRS transform too - bbox: QgsRectangle = self.grid_aligned_bbox(geom.boundingBox()) - - # Create a feature for the aligned bounding box - # Always save the study area bounding boxes regardless of mode - self.save_to_geopackage( - layer_name="study_area_bboxes", - geom=QgsGeometry.fromRect(bbox), - area_name=normalized_name, - ) - - # Transform the geometry to the output CRS - crs_src: QgsCoordinateReferenceSystem = self.layer.crs() - transform: QgsCoordinateTransform = QgsCoordinateTransform( - crs_src, self.output_crs, self.context.project() - ) - if crs_src != self.output_crs: - geom.transform(transform) - - # Create a feature for the original part - # Always save the study area bounding boxes regardless of mode - self.save_to_geopackage( - layer_name="study_area_polygons", geom=geom, area_name=normalized_name - ) - # Process the geometry based on the selected mode - if self.mode == "vector": - log_message(f"Creating vector grid for {normalized_name}.") - self.create_and_save_grid(geom, bbox) - elif self.mode == "raster": - # Grid must be made first as it is used when creating the clipping vector! - log_message(f"Creating vector grid for {normalized_name}.") - self.create_and_save_grid(geom, bbox) - log_message(f"Creating clip polygon for {normalized_name}.") - self.create_clip_polygon(geom, bbox, normalized_name) - log_message(f"Creating raster mask for {normalized_name}.") - self.create_raster_mask(geom, bbox, normalized_name) - - self.counter += 1 + log_message(f"Table '{self.status_table_name}' created in GeoPackage.") + finally: + ds = None - def process_multipart_geometry( - self, geom: QgsGeometry, normalized_name: str, area_name: str - ) -> None: + def add_row_to_status_tracking_table(self, area_name): """ - Processes each part of a multipart geometry by exploding the parts and delegating - to process_singlepart_geometry, ensuring the part index is included in the output name. - - :param geom: Geometry of the multipart feature. - :param normalized_name: Name normalized for file storage. - :param area_name: Name of the study area. + Adds a new row to the tracking table for area_name. """ - parts: List[QgsGeometry] = geom.asGeometryCollection() - - # Iterate over each part and delegate to process_singlepart_geometry - for part_index, part in enumerate(parts): - # Create a unique name for the part by appending the part index - part_normalized_name = f"{normalized_name}_part{part_index}" - - # Delegate to process_singlepart_geometry - self.process_singlepart_geometry(part, part_normalized_name, area_name) - - def grid_aligned_bbox(self, bbox: QgsRectangle) -> QgsRectangle: + ds = ogr.Open(self.gpkg_path, 1) + if not ds: + raise RuntimeError(f"Could not open {self.gpkg_path} for update.") + layer = ds.GetLayerByName(self.status_table_name) + if not layer: + raise RuntimeError(f"Missing status table layer: {self.status_table_name}") + + feat_defn = layer.GetLayerDefn() + feat = ogr.Feature(feat_defn) + feat.SetField("area_name", area_name) + feat.SetField("timestamp_start", None) + feat.SetField("timestamp_end", None) + feat.SetField("geometry_processed", 0) + feat.SetField("clip_geometry_processed", 0) + feat.SetField("grid_processed", 0) + feat.SetField("mask_processed", 0) + feat.SetField("grid_creation_duration_secs", 0.0) + feat.SetField("clip_geom_creation_duration_secs", 0.0) + feat.SetField("geom_total_duration_secs", 0.0) + layer.CreateFeature(feat) + feat = None + ds = None + + def set_status_tracking_table_value(self, area_name, field_name, value): """ - Transforms and aligns the bounding box to the grid in the output CRS. - The alignment ensures that the bbox aligns with the study area grid, offset by an exact multiple of - the grid size in m. - - :param bbox: The bounding box to be aligned, in the CRS of the input layer. - :return: A new bounding box aligned to the grid, in the output CRS. + Update a field value in the tracking table for the specified area_name. """ - # Transform the bounding box to the output CRS - crs_src: QgsCoordinateReferenceSystem = self.layer.crs() - transform: QgsCoordinateTransform = QgsCoordinateTransform( - crs_src, self.output_crs, self.context.project() + ds = ogr.Open(self.gpkg_path, 1) + if not ds: + raise RuntimeError(f"Could not open {self.gpkg_path} for update.") + layer = ds.GetLayerByName(self.status_table_name) + if not layer: + raise RuntimeError(f"Missing status table layer: {self.status_table_name}") + + layer.SetAttributeFilter(f"area_name = '{area_name}'") + for feature in layer: + feature.SetField(field_name, value) + layer.SetFeature(feature) + layer.ResetReading() + ds = None + log_message( + f"Updated processing status flag for {field_name} for {area_name} to {value}." ) - if crs_src != self.output_crs: - bbox_transformed = transform.transformBoundingBox(bbox) - else: - bbox_transformed = bbox - # Align the bounding box to a grid aligned at 100m intervals, offset by the study area origin - study_area_origin_x = ( - int(self.layer_bbox.xMinimum() // self.cell_size_m) * self.cell_size_m - ) - study_area_origin_y = ( - int(self.layer_bbox.yMinimum() // self.cell_size_m) * self.cell_size_m + ########################################################################## + # Geometry processing + ########################################################################## + def process_singlepart_geometry(self, geom, normalized_name, area_name): + """ + Process a single-part geometry: + 1) Align bounding box + 2) Save bounding box as a feature + 3) Transform geometry + 4) Save geometry + 5) Create vector grid + 6) Create clip polygon + 7) Optionally create raster mask + """ + # Used to track durations + geometry_start_time = time.time() + # Also grab the current timestamp and write it to the tracking table + now_str = datetime.datetime.utcnow().strftime("%Y-%m-%d %H:%M:%S") + self.set_status_tracking_table_value( + normalized_name, "timestamp_start", now_str ) - - # Align bbox to the grid based on the study area origin - x_min = ( - study_area_origin_x - + int( - (bbox_transformed.xMinimum() - study_area_origin_x) // self.cell_size_m - ) - * self.cell_size_m + # Compute aligned bounding box in target CRS + # (We already have a coordinate transformation if the source has a known SRS) + geometry_bbox = geom.GetEnvelope() # (xmin, xmax, ymin, ymax) + aligned_bbox = self.transform_and_align_bbox(geometry_bbox) + + # Save the bounding box for this geometry + self.save_bbox_polygon("study_area_bboxes", aligned_bbox, normalized_name) + + # Add a row to the tracking table + self.add_row_to_status_tracking_table(normalized_name) + + # If needed, transform the geometry to target CRS + # (Only if we have the coordinate transform) + if self.coord_transform: + geom.Transform(self.coord_transform) + + # Save the geometry (in the target CRS) to "study_area_polygons" + self.save_geometry_to_geopackage("study_area_polygons", geom, normalized_name) + self.set_status_tracking_table_value(normalized_name, "geometry_processed", 1) + + # Create the grid + log_message(f"Creating vector grid for {normalized_name}.") + start_time = time.time() + self.create_and_save_grid(normalized_name, geom, aligned_bbox) + self.set_status_tracking_table_value(normalized_name, "grid_processed", 1) + self.set_status_tracking_table_value( + normalized_name, "grid_creation_duration_secs", time.time() - start_time ) - y_min = ( - study_area_origin_y - + int( - (bbox_transformed.yMinimum() - study_area_origin_y) // self.cell_size_m - ) - * self.cell_size_m + # Create clip polygon + log_message(f"Creating clip polygon for {normalized_name}.") + start_time = time.time() + self.create_clip_polygon(geom, aligned_bbox, normalized_name) + self.set_status_tracking_table_value( + normalized_name, "clip_geometry_processed", 1 ) - x_max = ( - study_area_origin_x - + ( - int( - (bbox_transformed.xMaximum() - study_area_origin_x) - // self.cell_size_m - ) - + 1 - ) - * self.cell_size_m + self.set_status_tracking_table_value( + normalized_name, + "clip_geom_creation_duration_secs", + time.time() - start_time, ) - y_max = ( - study_area_origin_y - + ( - int( - (bbox_transformed.yMaximum() - study_area_origin_y) - // self.cell_size_m - ) - + 1 - ) - * self.cell_size_m + # (Optional) create raster mask + log_message(f"Creating raster mask for {normalized_name}.") + self.create_raster_mask(geom, aligned_bbox, normalized_name) + self.set_status_tracking_table_value(normalized_name, "mask_processed", 1) + now_str = datetime.datetime.utcnow().strftime("%Y-%m-%d %H:%M:%S") + self.set_status_tracking_table_value(normalized_name, "timestamp_end", now_str) + self.set_status_tracking_table_value( + normalized_name, + "geom_total_duration_secs", + time.time() - geometry_start_time, ) + self.counter += 1 + progress = int((self.counter / self.parts_count) * 100) + # We use the progress object to notify of progress in the subtask + # And the QgsTask progressChanged signal to track the main task + self.setProgress(progress) + log_message(f"XXXXXXXXXXXX Progress: {progress}% XXXXXXXXXXXXXXXXXXXXXXX") - y_min -= ( - self.cell_size_m - ) # Offset to ensure the grid covers the entire geometry - y_max += ( - self.cell_size_m - ) # Offset to ensure the grid covers the entire geometry - x_min -= ( - self.cell_size_m - ) # Offset to ensure the grid covers the entire geometry - x_max += ( - self.cell_size_m - ) # Offset to ensure the grid covers the entire geometry - - # Return the aligned bbox in the output CRS - return QgsRectangle(x_min, y_min, x_max, y_max) - - def save_to_geopackage( - self, layer_name: str, geom: QgsGeometry, area_name: str - ) -> None: + def process_multipart_geometry(self, geom, normalized_name, area_name): """ - Save features to GeoPackage. Create or append the layer as necessary. - :param layer_name: Name of the layer in the GeoPackage. - :param geom: Feature to append. - :param area_name: Name of the study area. - + Processes each part of a multi-part geometry. """ - self.create_layer_if_not_exists(layer_name) - self.append_to_layer(layer_name, geom, area_name) - - def append_to_layer( - self, layer_name: str, geom: QgsGeometry, area_name: str - ) -> None: + count = geom.GetGeometryCount() + for i in range(count): + part_geom = geom.GetGeometryRef(i) + part_name = f"{normalized_name}_part{i}" + self.process_singlepart_geometry(part_geom, part_name, area_name) + + ########################################################################## + # BBox handling + ########################################################################## + def transform_and_align_bbox(self, bbox): """ - Append feature to an existing layer in the GeoPackage. - - :param layer_name: Name of the layer in the GeoPackage. - :param geom: Feature to append. - :param area_name: Name of the study area. + BBox is (xmin, xmax, ymin, ymax). Transform to target CRS if possible, + then align to cell_size_m grid. + Returns new (xmin, xmax, ymin, ymax) in target CRS. """ - gpkg_layer_path = f"{self.gpkg_path}|layername={layer_name}" - gpkg_layer = QgsVectorLayer(gpkg_layer_path, layer_name, "ogr") - fields = gpkg_layer.fields() - # Create a list of placeholder values based on field count - attributes = [None] * len(fields) - - feature = QgsFeature() - feature.setGeometry(geom) - attributes[fields.indexFromName("area_name")] = area_name - feature.setAttributes(attributes) - - if gpkg_layer.isValid(): - log_message( - f"Appending to existing layer: {layer_name}", - tag="Geest", - level=Qgis.Info, - ) - provider = gpkg_layer.dataProvider() - provider.addFeatures([feature]) - gpkg_layer.updateExtents() - # Do not remove this line, QGIS will crash during the - # study area creation process. - del provider - else: - log_message( - f"Layer '{layer_name}' is not valid for appending.", - tag="Geest", - level=Qgis.Critical, - ) + (xmin, xmax, ymin, ymax) = bbox + + # If we have a coordinate transform, we need to convert min/max + # We'll do a polygon-based approach to ensure correctness + if self.coord_transform: + corner_points = [ + (xmin, ymin), + (xmin, ymax), + (xmax, ymax), + (xmax, ymin), + ] + # Transform each corner + transformed_corners = [] + for x, y in corner_points: + point = ogr.Geometry(ogr.wkbPoint) + point.AddPoint_2D(x, y) + point.Transform(self.coord_transform) + transformed_corners.append((point.GetX(), point.GetY())) + + # Recompute envelope from transformed corners + xs = [pt[0] for pt in transformed_corners] + ys = [pt[1] for pt in transformed_corners] + xmin, xmax = min(xs), max(xs) + ymin, ymax = min(ys), max(ys) + + log_message(f"Transformed bbox pre-alignment: {xmin}, {xmax}, {ymin}, {ymax}") + # Now align to the cell size + cell_size = self.cell_size_m + + def snap_down(value, base): + return (int(value // base)) * base + + def snap_up(value, base): + return (int(value // base) + 1) * base + + # Snap bounding values outward so we always cover the full geometry + x_min_snap = snap_down(xmin, cell_size) - cell_size + y_min_snap = snap_down(ymin, cell_size) - cell_size + x_max_snap = snap_up(xmax, cell_size) + cell_size + y_max_snap = snap_up(ymax, cell_size) + cell_size + log_message( + f"Aligned bbox : {x_min_snap}, {x_max_snap}, {y_min_snap}, {y_max_snap}" + ) + return (x_min_snap, x_max_snap, y_min_snap, y_max_snap) - def create_layer_if_not_exists(self, layer_name: str) -> None: + def save_bbox_polygon(self, layer_name, bbox, area_name): """ - Create a new layer in the GeoPackage if it doesn't already exist. - - It is assumed that all layers have the same structure of + Save a bounding-box polygon to the specified layer (creating it if needed). + """ + # BBox is (xmin, xmax, ymin, ymax) + (xmin, xmax, ymin, ymax) = bbox + ring = ogr.Geometry(ogr.wkbLinearRing) + ring.AddPoint(xmin, ymin) + ring.AddPoint(xmin, ymax) + ring.AddPoint(xmax, ymax) + ring.AddPoint(xmax, ymin) + ring.AddPoint(xmin, ymin) + polygon = ogr.Geometry(ogr.wkbPolygon) + polygon.AddGeometry(ring) + self.save_geometry_to_geopackage(layer_name, polygon, area_name) + + ########################################################################## + # Write geometry to GPKG layers + ########################################################################## + def save_geometry_to_geopackage(self, layer_name, geom, area_name): + """ + Append a single geometry to a (possibly newly created) layer in GPKG. + Each layer has a field 'area_name' (string) and polygon geometry. + """ + self.create_layer_if_not_exists(layer_name) + ds = ogr.Open(self.gpkg_path, 1) + layer = ds.GetLayerByName(layer_name) + if not layer: + raise RuntimeError( + f"Could not open target layer {layer_name} in {self.gpkg_path}" + ) - fid, area_name, geometry (Polygon) + feat_defn = layer.GetLayerDefn() + feature = ogr.Feature(feat_defn) + feature.SetField("area_name", area_name) + feature.SetGeometry(geom) + layer.CreateFeature(feature) + feature = None + ds = None - :param layer_name: Name of the layer to create. + def create_layer_if_not_exists(self, layer_name): """ + Create a GPKG layer if it does not exist. The layer has: + - a string field 'area_name' + - polygon geometry type + - SRS is self.target_spatial_ref + """ + if not os.path.exists(self.gpkg_path): + # Create new GPKG + driver = ogr.GetDriverByName("GPKG") + driver.CreateDataSource(self.gpkg_path) + + ds = ogr.Open(self.gpkg_path, 1) + layer = ds.GetLayerByName(layer_name) + if layer is not None: + ds = None + return # Already exists + + # Create it + layer = ds.CreateLayer( + layer_name, self.target_spatial_ref, geom_type=ogr.wkbPolygon + ) + # area_name field + field_defn = ogr.FieldDefn("area_name", ogr.OFTString) + layer.CreateField(field_defn) + ds = None + + # Helper to update time spent in a named metric block + def track_time(self, metric_name, start_time): + self.metrics[metric_name] += time.time() - start_time + + ########################################################################## + # Create Vector Grid + ########################################################################## + def create_and_save_grid(self, normalized_name, geom, bbox): + """ + Creates a vector grid covering bbox at self.cell_size_m spacing. + Writes those cells that intersect 'geom' to layer 'study_area_grid'. + (In practice, this can be quite large for big extents.) + """ + # ---------------------------- + # Initialize metrics tracking + # ---------------------------- + grid_layer_name = "study_area_grid" + self.create_grid_layer_if_not_exists(grid_layer_name) - fields = QgsFields() - fields.append(QgsField("area_name", QVariant.String)) - geometry_type = QgsWkbTypes.Polygon - - gpkg_layer_path = f"{self.gpkg_path}|layername={layer_name}" - layer = QgsVectorLayer(gpkg_layer_path, layer_name, "ogr") + ds = ogr.Open(self.gpkg_path, 1) # read-write + layer = ds.GetLayerByName(grid_layer_name) + if not layer: + raise RuntimeError(f"Could not open {grid_layer_name} for writing.") - append = True - # Check if the GeoPackage file exists - if not os.path.exists(self.gpkg_path): - append = False - log_message( - f"GeoPackage does not exist. Creating: {self.gpkg_path}", - tag="Geest", - level=Qgis.Info, - ) + xmin, xmax, ymin, ymax = bbox + cell_size = self.cell_size_m - # If the layer doesn't exist, create it - if not layer.isValid(): - log_message( - f"Layer '{layer_name}' does not exist. Creating it.", - tag="Geest", - level=Qgis.Info, - ) - options = QgsVectorFileWriter.SaveVectorOptions() - options.driverName = "GPKG" - options.fileEncoding = "UTF-8" - options.layerName = layer_name - if append: - options.actionOnExistingFile = ( - QgsVectorFileWriter.CreateOrOverwriteLayer - ) + log_message( + f"Creating grid for extents: xmin {xmin}, xmax {xmax}, ymin {ymin}, ymax {ymax}" + ) - # Create a new GeoPackage layer - QgsVectorFileWriter.create( - fileName=self.gpkg_path, - fields=fields, - geometryType=geometry_type, - srs=self.output_crs, - transformContext=QgsCoordinateTransformContext(), - options=options, - ) + # Approx count for logging + x_range_count = int((xmax - xmin) / cell_size) + y_range_count = int((ymax - ymin) / cell_size) + self.current_geom_cell_count_estimate = x_range_count * y_range_count + if self.current_geom_cell_count_estimate == 0: + log_message("No cells to generate.") return + log_message( + f"Estimated total cells to generate: {self.current_geom_cell_count_estimate}" + ) - def create_and_save_grid(self, geom: QgsGeometry, bbox: QgsRectangle) -> None: - """ - Creates a vector grid over the bounding box of a geometry and streams it directly to the GeoPackage. + # OGR geometry intersection can be slow for large grids. + # If this area is huge, consider a more robust approach or indexing. + # For demonstration, we do a naive approach. - :param geom: Geometry to create grid over. - :param bbox: Bounding box for the grid. - """ - grid_layer_name = "study_area_grid" - grid_fields = [QgsField("id", QVariant.Int)] + self.write_lock = False - self.create_layer_if_not_exists(grid_layer_name) - # Access the GeoPackage layer to append features - gpkg_layer_path = f"{self.gpkg_path}|layername={grid_layer_name}" - gpkg_layer = QgsVectorLayer(gpkg_layer_path, grid_layer_name, "ogr") + # worker_tasks = [] + # Get a reference to the global task manager + # task_manager = QgsApplication.taskManager() - # get the highest fid from the vector layer - if gpkg_layer.isValid(): - provider = gpkg_layer.dataProvider() - feature_id = provider.featureCount() - else: - log_message( - f"Failed to access layer '{grid_layer_name}' in the GeoPackage.", - tag="Geest", - level=Qgis.Critical, - ) - return + # Limit concurrency to 8 + # task_manager.setMaxActiveThreadCount(8) + feedback = QgsFeedback() - step = self.cell_size_m # cell_size_mm grid cells - feature_id += 1 - feature_batch = [] - log_message("Checking for meridian wrapping...") - try: - # Handle bounding box for meridian wrapping cases - x_min, x_max = bbox.xMinimum(), bbox.xMaximum() - y_min, y_max = bbox.yMinimum(), bbox.yMaximum() - if x_min > x_max: - # If spanning the meridian, adjust range to handle wraparound - x_ranges = [(x_min, 180), (-180, x_max)] - else: - x_ranges = [(x_min, x_max)] - except Exception as e: + # 1. Chunk the bounding box + # size is squared so 5 will make a 5x5 cell chunk + chunk_size = int(setting(key="chunk_size", default=50)) + bbox_chunks = list( + self.chunk_bbox(xmin, xmax, ymin, ymax, cell_size, chunk_size) + ) + # print out all the chunk bboxes + log_message(f"Chunk count: {len(bbox_chunks)}") + log_message(f"Chunk size: {chunk_size}") + chunk_count = len(bbox_chunks) + for idx, chunk in enumerate(bbox_chunks): + log_message(f"Chunk {idx}: {chunk}") + + self.feedback.setProgress(0) + for idx, chunk in enumerate(bbox_chunks): + start_time = ( + time.time() + ) # used for both create chunk start and total chunk start + task = GridFromBbox(idx, chunk, geom, cell_size, feedback) + self.track_time("Creating chunks", start_time) + # Not running in thread for now, see note below + task.run() + + self.write_chunk(layer, task, normalized_name) + # We use the progress object to notify of progress in the subtask + # And the QgsTask progressChanged signal to track the main task log_message( - f"Error handling meridian wrapping: {str(e)}", level=Qgis.Critical + f"XXXXXX Chunks Progress: {int((idx + 1 / chunk_count) * 100)}% XXXXXX" ) - import traceback - - log_message(traceback.format_exc(), level=Qgis.Critical) - raise - log_message("Meridian wrapping handled.") - - step = self.cell_size_m # Grid cell size - total_cells = ((x_max - x_min) // step) * ((y_max - y_min) // step) - cell_count = 0 - features_per_batch = 10000 - progress_log_interval = 10000 - log_message("Preparing spatial index...") - - # Spatial index for efficient intersection checks - spatial_index = QgsSpatialIndex() - feature = QgsFeature() - feature.setGeometry(geom) - spatial_index.addFeature(feature) - log_message("Spatial index prepared.") - + self.feedback.setProgress(int((idx + 1 / chunk_count) * 100)) + + # This is blocking, but we're in a thread + # Crashes QGIS, needs to be refactored to use QgsTask subtasks + # task.taskCompleted.connect(write_grids) + # worker_tasks.append(task) + # log_message(f"Task {idx} created for chunk {chunk}") + # log_message(f"{len(worker_tasks)} tasks queued.") + # task_manager.addTask(task) + self.track_time("Complete chunk", start_time) + ds = None + # ---------------------------- + # Print out metrics summary + # ---------------------------- + log_message("=== Metrics Summary ===") + for k, v in self.metrics.items(): + log_message(f"{k}: {v:.4f} seconds") + self.total_cells += self.current_geom_actual_cell_count + log_message(f"Grid creation completed for area {normalized_name}.") + + def write_chunk(self, layer, task, normalized_name): + start_time = time.time() + # Write locking is intended for a future version where we might have multiple threads + # currently I am just using the grid_from_bbox task to generate the geometries in a + # single thread by calling its run method directly. + # The reason for that is that QgsTask subtasks are only called after the parent's + # run method is completed, so I can't use them to write the geometries to the layer + # If write_lock is true, wait for the lock to be released + while self.write_lock: + log_message("Waiting for write lock...") + time.sleep(0.1) + log_message("Write lock released.") + log_message(f"Writing {len(task.features_out)} features to layer.") + self.track_time("Preparing chunks", task.run_time) + self.write_lock = True + feat_defn = layer.GetLayerDefn() try: - log_message("Creating grid features...") - # Iterate through defined ranges for x to handle wraparounds - for x_range in x_ranges: - for x in range(int(x_range[0]), int(x_range[1]), step): - for y in range(int(y_min), int(y_max), step): - if self.feedback.isCanceled(): - return False - - # Create cell rectangle - rect = QgsRectangle(x, y, x + step, y + step) - grid_geom = QgsGeometry.fromRect(rect) - # log_message(f"Creating rect: {rect}") - # Intersect check - if spatial_index.intersects( - grid_geom.boundingBox() - ) and grid_geom.intersects(geom): - # Create the grid feature - feature = QgsFeature() - feature.setGeometry(grid_geom) - feature.setAttributes([feature_id]) - feature_batch.append(feature) - feature_id += 1 - # log_message(f"Feature ID {feature_id} created.") - - # Commit in batches - if len(feature_batch) >= features_per_batch: - provider.addFeatures(feature_batch) - # gpkg_layer.updateExtents() - feature_batch.clear() - log_message( - f"Committed {feature_id} features to {grid_layer_name}." - ) - gpkg_layer.commitChanges() - # Reset provider to free memory - # del provider - # gpkg_layer = QgsVectorLayer( - # gpkg_layer_path, grid_layer_name, "ogr" - # ) - # provider = gpkg_layer.dataProvider() - - # Log progress periodically - cell_count += 1 - if cell_count % progress_log_interval == 0: - progress = (cell_count / total_cells) * 100 - log_message( - f"Processed {cell_count}/{total_cells} grid cells ({progress:.2f}% complete)." - ) - self.setProgress(int(progress)) - - # Commit any remaining features - if feature_batch: - provider.addFeatures(feature_batch) - gpkg_layer.commitChanges() - gpkg_layer.updateExtents() - provider.forceReload() - del gpkg_layer - log_message( - f"Final commit: {feature_id} features written to {grid_layer_name}." - ) - + for geometry in task.features_out: + feature = ogr.Feature(feat_defn) + feature.SetField("grid_id", self.current_geom_actual_cell_count) + feature.SetField("area_name", normalized_name) + feature.SetGeometry(geometry) + layer.CreateFeature(feature) + feature = None + + self.current_geom_actual_cell_count += 1 + if self.current_geom_actual_cell_count % 10000 == 0: + try: + log_message( + f" Cell count: {self.current_geom_actual_cell_count}" + ) + log_message( + f" Total cells: {self.current_geom_cell_count_estimate}" + ) + percent = ( + self.current_geom_actual_cell_count + / float(self.current_geom_cell_count_estimate) + ) * 100.0 + log_message(f" Percent complete: {percent:.2f}%") + log_message( + f"Grid creation for part {normalized_name}: {self.current_geom_actual_cell_count}/{self.current_geom_cell_count_estimate} ({percent:.2f}%)" + ) + except ZeroDivisionError: + pass + # commit changes + layer.SyncToDisk() + self.track_time("Writing chunks", start_time) + self.write_lock = False except Exception as e: - log_message(f"Error during grid creation: {str(e)}", level=Qgis.Critical) - import traceback - - log_message(traceback.format_exc(), level=Qgis.Critical) - return + log_message(f"write_grids: {str(e)}") + log_message(f"write_grids: {traceback.format_exc()}") + self.write_lock = False - log_message(f"Grid creation completed. Total features written: {feature_id}.") - - def create_clip_polygon( - self, geom: QgsGeometry, aligned_box: QgsRectangle, normalized_name: str - ) -> str: + def create_grid_layer_if_not_exists(self, layer_name): """ - Creates a polygon like the study area geometry, but with each - area extended to completely cover the grid cells that intersect with the geometry. - - We do this by combining the area polygon and all the cells that occur along - its edge. This is necessary because gdalwarp clip only includes cells that - have 50% coverage or more and it has no -at option like gdalrasterize has. - - The cells and the original area polygon are dissolved into a single polygon - representing the outer boundary of all the edge cells and the original area polygon. - - :param geom: Geometry to be processed. - :param aligned_box: Aligned bounding box for the geometry. - :param normalized_name: Name for the output area - - :return: None. + Create a grid layer with 'grid_id' as integer field + and a polygon geometry if it does not exist. """ - # Create a memory layer to hold the geometry - temp_layer = QgsVectorLayer( - f"Polygon?crs={self.output_crs.authid()}", "temp_mask_layer", "memory" - ) - temp_layer_data_provider = temp_layer.dataProvider() - # get the geometry as a linestring - linestring = geom.coerceToType(QgsWkbTypes.LineString)[0] - # The linestring and the original polygon geom should have the same bbox - if linestring.boundingBox() != geom.boundingBox(): - raise Exception("Bounding boxes of linestring and polygon do not match.") - - # we are going to select all grid cells that intersect the linestring - # so that we can ensure gdal rasterize does not miss any land areas - gpkg_layer_path = f"{self.gpkg_path}|layername=study_area_grid" - gpkg_layer = QgsVectorLayer(gpkg_layer_path, "study_area_grid", "ogr") - - # Create a spatial index for efficient spatial querying - spatial_index = QgsSpatialIndex(gpkg_layer.getFeatures()) - - # Get feature IDs of candidates that may intersect with the multiline geometry - candidate_ids = spatial_index.intersects(geom.boundingBox()) - - if len(candidate_ids) == 0: - raise Exception( - f"No candidate cells on boundary of the geometry: {self.working_dir}" + if not os.path.exists(self.gpkg_path): + driver = ogr.GetDriverByName("GPKG") + driver.CreateDataSource(self.gpkg_path) + + ds = ogr.Open(self.gpkg_path, 1) + layer = ds.GetLayerByName(layer_name) + if layer is None: + layer = ds.CreateLayer( + layer_name, self.target_spatial_ref, geom_type=ogr.wkbPolygon ) + field_defn = ogr.FieldDefn("grid_id", ogr.OFTInteger) + layer.CreateField(field_defn) + field_defn = ogr.FieldDefn("area_name", ogr.OFTString) + layer.CreateField(field_defn) + ds = None + + ########################################################################## + # Create Clip Polygon + ########################################################################## + def create_clip_polygon(self, geom, aligned_box, normalized_name): + """ + Creates a polygon that includes the original geometry plus all grid cells + that intersect the boundary of the geometry. Then dissolves them into one polygon. + """ + # 1) We load the grid from GPKG + grid_ds = ogr.Open(self.gpkg_path, 0) + grid_layer = grid_ds.GetLayerByName("study_area_grid") + if not grid_layer: + raise RuntimeError("Missing study_area_grid layer.") + + # We'll do a bounding box filter for performance + (xmin, ymin, xmax, ymax) = aligned_box + grid_layer.SetSpatialFilterRect(xmin, ymin, xmax, ymax) + + # 2) We'll gather all grid cells that intersect *the boundary* of geom + # In OGR, we can do: + boundary = geom.GetBoundary() # line geometry for polygon boundary + + union_geom = ogr.Geometry(ogr.wkbPolygon) + union_geom.Destroy() # We'll handle it differently—see below. + + # For union accumulation, start with a null geometry + dissolved_geom = None + + # For clarity, transform boundary to the same SRS if needed (already is). + # We'll just do an Intersects check with each cell. + + grid_layer.ResetReading() + count = 0 + for f in grid_layer: + cell_geom = f.GetGeometryRef() + if not cell_geom: + continue + if boundary.Intersects(cell_geom): + # We'll union + if dissolved_geom is None: + dissolved_geom = cell_geom.Clone() + else: + dissolved_geom = dissolved_geom.Union(cell_geom) + count += 1 + grid_layer.ResetReading() + + # Also union the original geom itself + if dissolved_geom is None: + # No boundary cells found, fallback + dissolved_geom = geom.Clone() + else: + dissolved_geom = dissolved_geom.Union(geom) - # Filter candidates by precise geometry intersection - intersecting_ids = [] - for feature_id in candidate_ids: - feature = gpkg_layer.getFeature(feature_id) - if feature.geometry().intersects(linestring): - intersecting_ids.append(feature_id) - - if len(intersecting_ids) == 0: - raise Exception("No cells on boundary of the geometry") - # Select intersecting features in the layer - gpkg_layer.selectByIds(intersecting_ids) - - log_message( - f"Selected {len(intersecting_ids)} features that intersect with the multiline geometry." - ) - - # Define a field to store the area name - temp_layer_data_provider.addAttributes( - [QgsField(self.field_name, QVariant.String)] - ) - temp_layer.updateFields() - - # Add the geometry to the memory layer - temp_feature = QgsFeature() - temp_feature.setGeometry(geom) - temp_feature.setAttributes( - [normalized_name] - ) # Setting an arbitrary value for the mask - - # Add the main geometry for this part of the country - temp_layer_data_provider.addFeature(temp_feature) - - # Now all the grid cells get added that intersect with the geometry border - # since gdal rasterize only includes cells that have 50% coverage or more - # by the looks of things. - selected_features = gpkg_layer.selectedFeatures() - new_features = [] - - for feature in selected_features: - # Create a new feature for emp_layer - new_feature = QgsFeature() - new_feature.setGeometry(feature.geometry()) - new_feature.setAttributes([normalized_name]) - new_features.append(new_feature) - - # Add the features to temp_layer - temp_layer_data_provider.addFeatures(new_features) - - # commit all changes - temp_layer.updateExtents() - temp_layer.commitChanges() - # check how many features we have - feature_count = temp_layer.featureCount() - log_message( - f"Added {feature_count} features to the temp layer for mask creation." - ) - - output = processing.run( - "native:dissolve", - { - "INPUT": temp_layer, - "FIELD": [], - "SEPARATE_DISJOINT": False, - "OUTPUT": "TEMPORARY_OUTPUT", - }, - )["OUTPUT"] - # get the first feature from the output layer - feature = next(output.getFeatures()) - # get the geometry - clip_geom = feature.geometry() - - self.save_to_geopackage( - layer_name="study_area_clip_polygons", - geom=clip_geom, - area_name=normalized_name, + # dissolved_geom is now the final clip polygon + self.save_geometry_to_geopackage( + "study_area_clip_polygons", dissolved_geom, normalized_name ) log_message(f"Created clip polygon: {normalized_name}") - return - def create_raster_mask( - self, geom: QgsGeometry, aligned_box: QgsRectangle, mask_name: str - ) -> str: + ########################################################################## + # Split the bbox into chunks for parallel processing + ########################################################################## + def chunk_bbox(self, xmin, xmax, ymin, ymax, cell_size, chunk_size=1000): """ - Creates a 1-bit raster mask for a single geometry. - - :param geom: Geometry to be rasterized. - :param aligned_box: Aligned bounding box for the geometry. - :param mask_name: Name for the output raster file. + Generator that yields bounding box chunks. Each chunk is a tuple: + (x_start, x_end, y_start, y_end). - :return: The path to the created raster mask. + `chunk_size` indicates how many cells in the X-direction + (and optionally also Y-direction) you want per chunk. """ - mask_filepath = os.path.join(self.working_dir, "study_area", f"{mask_name}.tif") - # Create a memory layer to hold the geometry - temp_layer = QgsVectorLayer( - f"Polygon?crs={self.output_crs.authid()}", "temp_mask_layer", "memory" - ) - temp_layer_data_provider = temp_layer.dataProvider() + x_range_count = int((xmax - xmin) / cell_size) + y_range_count = int((ymax - ymin) / cell_size) - # Define a field to store the mask value - temp_layer_data_provider.addAttributes( - [QgsField(self.field_name, QVariant.String)] - ) - temp_layer.updateFields() - - # Add the geometry to the memory layer - temp_feature = QgsFeature() - temp_feature.setGeometry(geom) - temp_feature.setAttributes(["1"]) # Setting an arbitrary value for the mask - - # Add the main geometry for this part of the country - temp_layer_data_provider.addFeature(temp_feature) - - # commit all changes - temp_layer.updateExtents() - temp_layer.commitChanges() - - # Ensure resolution parameters are properly formatted as float values - x_res = self.cell_size_m # 100m pixel size in X direction - y_res = self.cell_size_m # 100m pixel size in Y direction - - # Define rasterization parameters for the temporary layer - params = { - "INPUT": temp_layer, - "FIELD": None, - "BURN": 1, - "USE_Z": False, - "UNITS": 1, - "WIDTH": x_res, - "HEIGHT": y_res, - "EXTENT": f"{aligned_box.xMinimum()},{aligned_box.xMaximum()}," - f"{aligned_box.yMinimum()},{aligned_box.yMaximum()}", # Extent of the aligned bbox - "NODATA": 0, - "OPTIONS": "", - "DATA_TYPE": 0, # byte - "INIT": None, - "INVERT": False, - "EXTRA": "-co NBITS=1 -at", # -at is for all touched cells - "OUTPUT": mask_filepath, - } - processing.run("gdal:rasterize", params) - log_message(f"Created raster mask: {mask_filepath}") - return mask_filepath + x_blocks = range(0, x_range_count, chunk_size) + y_blocks = range(0, y_range_count, chunk_size) + for x_block_start in x_blocks: + log_message(f"Processing chunk {x_block_start} of {x_range_count}") + x_block_end = min(x_block_start + chunk_size, x_range_count) - def calculate_utm_zone(self, bbox: QgsRectangle) -> int: - """ - Calculates the UTM zone based on the centroid of the bounding box. - The calculation is based on transforming the centroid to WGS84 (EPSG:4326). + # Convert from cell index to real coords + x_start_coord = xmin + x_block_start * cell_size + x_end_coord = xmin + x_block_end * cell_size - :param bbox: The bounding box of the layer. - :return: The EPSG code for the UTM zone. - """ - # Get the source CRS (from the input layer) - crs_src: QgsCoordinateReferenceSystem = self.layer.crs() - crs_wgs84: QgsCoordinateReferenceSystem = QgsCoordinateReferenceSystem( - "EPSG:4326" - ) + for y_block_start in y_blocks: + log_message(f"Processing chunk {y_block_start} of {y_range_count}") + y_block_end = min(y_block_start + chunk_size, y_range_count) - # Create the transform object - transform: QgsCoordinateTransform = QgsCoordinateTransform( - crs_src, crs_wgs84, self.context.project() - ) + # Convert from cell index to real coords + y_start_coord = ymin + y_block_start * cell_size + y_end_coord = ymin + y_block_end * cell_size - # Transform the centroid to WGS84 - centroid: QgsPointXY = bbox.center() - centroid_wgs84 = transform.transform(centroid) + log_message( + f"Created Chunk bbox: {x_start_coord}, {x_end_coord}, {ymin}, {ymax}" + ) + yield (x_start_coord, x_end_coord, y_start_coord, y_end_coord) - lon = centroid_wgs84.x() - lat = centroid_wgs84.y() + ########################################################################## + # Create Raster Mask + ########################################################################## + def create_raster_mask(self, geom, aligned_box, mask_name): + """ + Creates a 1-bit raster mask for a single geometry using gdal.Rasterize. + """ + mask_filepath = os.path.join(self.working_dir, "study_area", f"{mask_name}.tif") - # Calculate the UTM zone based on the longitude - utm_zone = int((lon + 180) / 6) + 1 + driver_mem = ogr.GetDriverByName("Memory") + mem_ds = driver_mem.CreateDataSource("temp") + mem_lyr = mem_ds.CreateLayer( + "temp_mask_layer", self.target_spatial_ref, geom_type=ogr.wkbPolygon + ) - # EPSG code for the UTM zone - if lat >= 0: - epsg_code = 32600 + utm_zone # Northern Hemisphere - else: - epsg_code = 32700 + utm_zone # Southern Hemisphere + # Create a field to burn + field_def = ogr.FieldDefn("burnval", ogr.OFTInteger) + mem_lyr.CreateField(field_def) + + # Create feature + feat_defn = mem_lyr.GetLayerDefn() + feat = ogr.Feature(feat_defn) + feat.SetField("burnval", 1) + feat.SetGeometry(geom.Clone()) + mem_lyr.CreateFeature(feat) + feat = None + + # Now call gdal.Rasterize + x_res = self.cell_size_m + y_res = self.cell_size_m + + (xmin, xmax, ymin, ymax) = aligned_box + + # For pixel width/height, we can compute: + # width in coordinate space: (xmax - xmin) + width = int((xmax - xmin) / x_res) + height = int((ymax - ymin) / y_res) + if width < 1 or height < 1: + log_message("Extent is too small for raster creation. Skipping mask.") + return - return epsg_code + # Create the raster + # NB: gdal.GetDriverByName('GTiff').Create() expects col, row order + target_ds = gdal.GetDriverByName("GTiff").Create( + mask_filepath, + width, + height, + 1, # 1 band + gdal.GDT_Byte, + options=["NBITS=1"], # 1-bit + ) + if not target_ds: + raise RuntimeError(f"Could not create raster {mask_filepath}") + + # Set geotransform (origin x, pixel width, rotation, origin y, rotation, pixel height) + # Note y_res is negative if north-up. We'll use negative for correct alignment in typical north-up data + geotransform = (xmin, x_res, 0.0, ymax, 0.0, -y_res) + target_ds.SetGeoTransform(geotransform) + target_ds.SetProjection(self.target_spatial_ref.ExportToWkt()) + + # Rasterize + err = gdal.RasterizeLayer( + target_ds, + [1], # bands + mem_lyr, # layer + burn_values=[1], # burn value for the feature + options=["ALL_TOUCHED=TRUE"], + ) + if err != 0: + log_message(f"Error in RasterizeLayer: {err}", level="CRITICAL") - def create_study_area_directory(self, working_dir: str) -> None: - """ - Creates the 'study_area' directory if it doesn't already exist. + target_ds.FlushCache() + target_ds = None + mem_ds = None - :param working_dir: Path to the working directory. - """ - study_area_dir = os.path.join(working_dir, "study_area") - if not os.path.exists(study_area_dir): - try: - os.makedirs(study_area_dir) - log_message( - f"Created study area directory: {study_area_dir}", - tag="Geest", - level=Qgis.Info, - ) - except Exception as e: - log_message(f"Error creating directory: {e}", level=Qgis.Critical) + log_message(f"Created raster mask: {mask_filepath}") + return mask_filepath - def create_raster_vrt(self, output_vrt_name: str = "combined_mask.vrt") -> None: + ########################################################################## + # Create VRT + ########################################################################## + def create_raster_vrt(self, output_vrt_name="combined_mask.vrt"): """ - Creates a VRT file from all generated raster masks and adds it to the QGIS map. - - :param output_vrt_name: The name of the VRT file to create. + Creates a VRT file from all .tif masks in the 'study_area' dir using + gdal.BuildVRT (Python API) or gdalbuildvrt approach. """ - log_message( - f"Creating VRT of masks '{output_vrt_name}' layer to the map.", - tag="Geest", - level=Qgis.Info, - ) - # Directory containing raster masks raster_dir = os.path.join(self.working_dir, "study_area") raster_files = glob.glob(os.path.join(raster_dir, "*.tif")) if not raster_files: - log_message( - "No raster masks found to combine into VRT.", - tag="Geest", - level=Qgis.Warning, - ) + log_message("No raster masks found to build VRT.") return vrt_filepath = os.path.join(raster_dir, output_vrt_name) + log_message(f"Building VRT: {vrt_filepath}") - # Define the VRT parameters - params = { - "INPUT": raster_files, - "RESOLUTION": 0, # Use highest resolution among input files - "SEPARATE": False, # Combine all input rasters as a single band - "OUTPUT": vrt_filepath, - "PROJ_DIFFERENCE": False, - "ADD_ALPHA": False, - "ASSIGN_CRS": None, - "RESAMPLING": 0, - "SRC_NODATA": "0", - "EXTRA": "", - } + # Use gdal.BuildVRT + vrt = gdal.BuildVRT(vrt_filepath, raster_files, separate=False) + if vrt is None: + log_message("Failed to create VRT.", level="CRITICAL") + return + vrt.FlushCache() + vrt = None - # Run the gdal:buildvrt processing algorithm to create the VRT - processing.run("gdal:buildvirtualraster", params) log_message(f"Created VRT: {vrt_filepath}") - # Add the VRT to the QGIS map - vrt_layer = QgsRasterLayer(vrt_filepath, "Combined Mask VRT") + ########################################################################## + # CRS / UTM calculation + ########################################################################## + def calculate_utm_zone(self, bbox, source_epsg=None): + """ + Determine a UTM zone from the centroid of (xmin, xmax, ymin, ymax), + reprojected into WGS84 if possible. Return EPSG code. + """ + (xmin, xmax, ymin, ymax) = bbox + cx = 0.5 * (xmin + xmax) + cy = 0.5 * (ymin + ymax) - if vrt_layer.isValid(): - # Not thread safe, use signal instead - # QgsProject.instance().addMapLayer(vrt_layer) - log_message("Added VRT layer to the map.") + # If there's no source SRS, we'll assume it's already lat/lon + if not source_epsg: + # fallback if no known EPSG + log_message( + "Source has no EPSG, defaulting to a naive assumption of WGS84 bounding box." + ) + lon, lat = cx, cy + else: + # We have a known EPSG, so transform centroid to WGS84 + src_ref = osr.SpatialReference() + src_ref.ImportFromEPSG(int(source_epsg)) + wgs84_ref = osr.SpatialReference() + wgs84_ref.ImportFromEPSG(4326) + ct = osr.CoordinateTransformation(src_ref, wgs84_ref) + point = ogr.Geometry(ogr.wkbPoint) + point.AddPoint(cx, cy) + point.Transform(ct) + lon = point.GetX() + lat = point.GetY() + + # Standard formula for UTM zone + utm_zone = int((lon + 180) / 6) + 1 + # We guess north or south + if lat >= 0: + return 32600 + utm_zone # Northern Hemisphere else: - log_message("Failed to add VRT layer to the map.", level=Qgis.Critical) + return 32700 + utm_zone # Southern Hemisphere + + ########################################################################## + # Directory creation + ########################################################################## + def create_study_area_directory(self, working_dir): + """ + Create 'study_area' subdir if not exist + """ + study_area_dir = os.path.join(working_dir, "study_area") + if not os.path.exists(study_area_dir): + os.makedirs(study_area_dir) + log_message(f"Created directory {study_area_dir}") + + ############################################################################## diff --git a/geest/core/workflows/single_point_buffer_workflow.py b/geest/core/workflows/single_point_buffer_workflow.py index 23aa64b5..0c6a4668 100644 --- a/geest/core/workflows/single_point_buffer_workflow.py +++ b/geest/core/workflows/single_point_buffer_workflow.py @@ -60,9 +60,13 @@ def __init__( log_message("single_buffer_point_layer not valid", level=Qgis.Critical) log_message(f"Layer Source: {layer_source}", level=Qgis.Critical) return False - + default_buffer_distance = int( + self.attributes.get("default_single_buffer_distance", "5000") + ) self.buffer_distance = int( - self.attributes.get("single_buffer_point_layer_distance", "5000") + self.attributes.get( + "single_buffer_point_layer_distance", default_buffer_distance + ) ) def _process_features_for_area( diff --git a/geest/core/workflows/workflow_base.py b/geest/core/workflows/workflow_base.py index 623bdd0e..a0d6b0d1 100644 --- a/geest/core/workflows/workflow_base.py +++ b/geest/core/workflows/workflow_base.py @@ -27,7 +27,7 @@ combine_rasters_to_vrt, ) from geest.core.constants import GDAL_OUTPUT_DATA_TYPE -from geest.utilities import log_message +from geest.utilities import log_message, log_layer_count class WorkflowBase(QObject): @@ -56,6 +56,8 @@ def __init__( :working_directory: Folder containing study_area.gpkg and where the outputs will be placed. If not set will be taken from QSettings. """ super().__init__() + log_layer_count() # For performance tuning, write the number of open layers to a log file + # we will log the layer count again at then end of the workflow self.item = item # ⭐️ This is a reference - whatever you change in this item will directly update the tree self.cell_size_m = cell_size_m self.feedback = feedback @@ -252,6 +254,7 @@ def execute(self) -> bool: return False area_iterator = AreaIterator(self.gpkg_path) + log_layer_count() # For performance tuning, write the number of open layers to a log file try: for index, (current_area, clip_area, current_bbox, progress) in enumerate( area_iterator @@ -326,6 +329,7 @@ def execute(self) -> bool: ) self.attributes["execution_end_time"] = datetime.datetime.now().isoformat() self.attributes["error_file"] = None + log_layer_count() # For performance tuning, write the number of open layers to a log file return True except Exception as e: @@ -354,6 +358,7 @@ def execute(self) -> bool: f.write(traceback.format_exc()) self.attributes["error_file"] = error_path self.attributes["error"] = f"Failed to process {self.workflow_name}: {e}" + log_layer_count() # For performance tuning, write the number of open layers to a log file return False def _create_workflow_directory(self) -> str: diff --git a/geest/gui/dialogs/analysis_aggregation_dialog.py b/geest/gui/dialogs/analysis_aggregation_dialog.py index 9de7e4ea..16564799 100644 --- a/geest/gui/dialogs/analysis_aggregation_dialog.py +++ b/geest/gui/dialogs/analysis_aggregation_dialog.py @@ -19,7 +19,7 @@ QLineEdit, ) from qgis.PyQt.QtGui import QPixmap, QDesktopServices -from qgis.PyQt.QtCore import Qt, QUrl, QSettings +from qgis.PyQt.QtCore import Qt, QUrl, QSettings, QByteArray from qgis.core import Qgis, QgsMapLayerProxyModel, QgsProject from geest.utilities import ( resources_path, @@ -27,6 +27,7 @@ setting, is_qgis_dark_theme_active, get_ui_class, + log_window_geometry, ) from qgis.gui import QgsMapLayerComboBox from geest.gui.widgets import CustomBannerLabel @@ -191,12 +192,17 @@ def restore_dialog_geometry(self): Restore the dialog geometry from QSettings. """ settings = QSettings() - geometry = settings.value("AnalysisAggregationDialog/geometry") + # Restore geometry (fall back to empty QByteArray if setting not found) + geometry = settings.value("AnalysisAggregationDialog/geometry", QByteArray()) + log_window_geometry(geometry) if geometry: log_message("Restoring dialog geometry") try: - self.resize(geometry.width(), geometry.height()) - except AttributeError: + self.restoreGeometry(geometry) + except Exception as e: + log_message( + "Restoring geometry failed", tag="Geest", level=Qgis.Warning + ) pass else: log_message("No saved geometry found, resizing dialog") @@ -216,6 +222,12 @@ def save_geometry(self): settings = QSettings() settings.setValue("AnalysisAggregationDialog/geometry", self.geometry()) + def closeEvent(self, event): + # Save geometry before closing + settings = QSettings() + settings.setValue("AnalysisAggregationDialog/geometry", self.saveGeometry()) + super().closeEvent(event) + def setup_table(self): """ Set up the QTableWidget to display the analysis dimensions. diff --git a/geest/gui/geest_settings.py b/geest/gui/geest_settings.py index 23712382..df81e9d7 100644 --- a/geest/gui/geest_settings.py +++ b/geest/gui/geest_settings.py @@ -59,9 +59,12 @@ def __init__(self, parent=None): self.ors_key_line_edit.setText(ors_key) else: self.ors_key_line_edit.setPlaceholderText("Enter your ORS API key here") - ors_request_size = int(setting(key="ors_request_size", default=100)) + ors_request_size = int(setting(key="ors_request_size", default=5)) self.ors_request_size.setValue(ors_request_size) + chunk_size = int(setting(key="chunk_size", default=50)) + self.chunk_size.setValue(chunk_size) + def apply(self): """Process the animation sequence. @@ -84,6 +87,7 @@ def apply(self): set_setting(key="ors_key", value=self.ors_key_line_edit.text()) set_setting(key="ors_request_size", value=self.ors_request_size.value()) + set_setting(key="chunk_size", value=self.chunk_size.value()) class GeestOptionsFactory(QgsOptionsWidgetFactory): diff --git a/geest/gui/panels/create_project_panel.py b/geest/gui/panels/create_project_panel.py index ffe164be..cc972e16 100644 --- a/geest/gui/panels/create_project_panel.py +++ b/geest/gui/panels/create_project_panel.py @@ -1,6 +1,7 @@ import os import json import shutil +import traceback from PyQt5.QtWidgets import ( QWidget, QFileDialog, @@ -13,8 +14,8 @@ QgsProject, Qgis, QgsProject, - QgsProcessingContext, QgsFeedback, + QgsLayerTreeGroup, ) from qgis.PyQt.QtCore import QSettings, pyqtSignal @@ -85,6 +86,7 @@ def initUI(self): self.load_boundary_button.clicked.connect(self.load_boundary) self.progress_bar.setVisible(False) + self.child_progress_bar.setVisible(False) def on_previous_button_clicked(self): self.switch_to_previous_tab.emit() @@ -162,6 +164,7 @@ def create_project(self): shutil.copy(default_model_path, model_path) except Exception as e: QMessageBox.critical(self, "Error", f"Failed to copy model.json: {e}") + self.enable_widgets() return # open the model.json to set the analysis cell size, then close it again with open(model_path, "r") as f: @@ -172,9 +175,9 @@ def create_project(self): # Create the processor instance and process the features debug_env = int(os.getenv("GEEST_DEBUG", 0)) - context = QgsProcessingContext() - context.setProject(QgsProject.instance()) - feedback = QgsFeedback() + feedback = ( + QgsFeedback() + ) # Used to cancel tasks and measure subtask progress try: processor = StudyAreaProcessingTask( @@ -183,29 +186,65 @@ def create_project(self): cell_size_m=self.cell_size_spinbox.value(), crs=crs, working_dir=self.working_dir, - context=context, feedback=feedback, ) # Hook up the QTask feedback signal to the progress bar - self.progress_updated(0) + # Measure overall task progress from the task object itself processor.progressChanged.connect(self.progress_updated) processor.taskCompleted.connect(self.on_task_completed) - + # Measure subtask progress from the feedback object + feedback.progressChanged.connect(self.subtask_progress_updated) + self.disable_widgets() if debug_env: processor.process_study_area() else: self.queue_manager.add_task(processor) self.queue_manager.start_processing() except Exception as e: - QMessageBox.critical(self, "Error", f"Error processing study area: {e}") + trace = traceback.format_exc() + QMessageBox.critical( + self, "Error", f"Error processing study area: {e}\n{trace}" + ) + self.enable_widgets() return self.settings.setValue("last_working_directory", self.working_dir) self.set_working_directory.emit(self.working_dir) + self.enable_widgets() + + def disable_widgets(self): + """Disable all widgets in the panel.""" + for widget in self.findChildren(QWidget): + widget.setEnabled(False) - def progress_updated(self, progress): + def enable_widgets(self): + """Enable all widgets in the panel.""" + for widget in self.findChildren(QWidget): + widget.setEnabled(True) + + # Slot that listens for changes in the study_area task object which is used to measure overall task progress + def progress_updated(self, progress: float): """Slot to be called when the task progress is updated.""" + log_message(f"\n\n\n\n\n\n\Progress: {progress}\n\n\n\n\n\n\n\n") self.progress_bar.setVisible(True) + self.progress_bar.setEnabled(True) self.progress_bar.setValue(int(progress)) + # This is a sneaky hack to show the exact progress in the label + # since QProgressBar only takes ints. See Qt docs for more info. + # Use the 'setFormat' method to display the exact float: + float_value_as_string = f"Total progress: {progress}%" + self.progress_bar.setFormat(float_value_as_string) + self.add_bboxes_to_map() # will just refresh them if already there + + # Slot that listens for changes in the progress object which is used to measure subtask progress + def subtask_progress_updated(self, progress: float): + self.child_progress_bar.setVisible(True) + self.child_progress_bar.setEnabled(True) + self.child_progress_bar.setValue(int(progress)) + # This is a sneaky hack to show the exact progress in the label + # since QProgressBar only takes ints. See Qt docs for more info. + # Use the 'setFormat' method to display the exact float: + float_value_as_string = f"Current geometry progress: {progress}%" + self.child_progress_bar.setFormat(float_value_as_string) def on_task_completed(self): """Slot to be called when the task completes successfully.""" @@ -215,6 +254,8 @@ def on_task_completed(self): level=Qgis.Info, ) self.progress_bar.setVisible(False) + self.child_progress_bar.setVisible(False) + self.enable_widgets() self.switch_to_next_tab.emit() def update_recent_projects(self, directory): @@ -241,13 +282,13 @@ def resizeEvent(self, event): def set_font_size(self): # Scale the font size to fit the text in the available space - log_message(f"Description Label Width: {self.description.rect().width()}") + # log_message(f"Description Label Width: {self.description.rect().width()}") # scale the font size linearly from 16 pt to 8 ps as the width of the panel decreases font_size = int( linear_interpolation(self.description.rect().width(), 12, 16, 400, 600) ) - log_message(f"Description Label Font Size: {font_size}") + # log_message(f"Description Label Font Size: {font_size}") self.description.setFont(QFont("Arial", font_size)) self.description2.setFont(QFont("Arial", font_size)) self.description3.setFont(QFont("Arial", font_size)) @@ -256,3 +297,69 @@ def set_font_size(self): self.cell_size_spinbox.setFont(QFont("Arial", font_size)) self.layer_combo.setFont(QFont("Arial", font_size)) self.field_combo.setFont(QFont("Arial", font_size)) + + def add_bboxes_to_map(self): + """Add the study area layers to the map. + + If it is already there we will just refrehs it. + + This provides the user with visual feedback as each geometry gets processed. + + This method is a cut and paste of the same method in tree_panel.py + but only adds two layers (well one layer and one table) to the map. + + """ + gpkg_path = os.path.join(self.working_dir, "study_area", "study_area.gpkg") + project = QgsProject.instance() + + # Check if 'Geest' group exists, otherwise create it + root = project.layerTreeRoot() + geest_group = root.findGroup("Geest Study Area") + if geest_group is None: + geest_group = root.insertGroup( + 0, "Geest Study Area" + ) # Insert at the top of the layers panel + + layers = [ + "study_area_bboxes", + "study_area_creation_status", + ] + for layer_name in layers: + gpkg_layer_path = f"{gpkg_path}|layername={layer_name}" + layer = QgsVectorLayer(gpkg_layer_path, layer_name, "ogr") + + if not layer.isValid(): + log_message( + f"Failed to add '{layer_name}' layer to the map.", + tag="Geest", + level=Qgis.Critical, + ) + continue + + source_qml = resources_path("resources", "qml", f"{layer_name}.qml") + result = layer.loadNamedStyle(source_qml) + if result[0]: # loadNamedStyle returns (success, error_message) + print(f"Successfully applied QML style to layer '{layer_name}'") + else: + print(f"Failed to apply QML style: {result[1]}") + + # Check if a layer with the same data source exists in the correct group + existing_layer = None + for child in geest_group.children(): + if isinstance(child, QgsLayerTreeGroup): + continue + if child.layer().source() == gpkg_layer_path: + existing_layer = child.layer() + break + + # If the layer exists, refresh it instead of removing and re-adding + if existing_layer is not None: + log_message(f"Refreshing existing layer: {existing_layer.name()}") + existing_layer.reload() + else: + # Add the new layer to the appropriate subgroup + QgsProject.instance().addMapLayer(layer, False) + layer_tree_layer = geest_group.addLayer(layer) + log_message( + f"Added layer: {layer.name()} to group: {geest_group.name()}" + ) diff --git a/geest/gui/panels/credits_panel.py b/geest/gui/panels/credits_panel.py index 5b188027..d0a04278 100644 --- a/geest/gui/panels/credits_panel.py +++ b/geest/gui/panels/credits_panel.py @@ -63,11 +63,11 @@ def resizeEvent(self, event): def set_font_size(self): # Scale the font size to fit the text in the available space - log_message(f"Label Width: {self.description.rect().width()}") + # log_message(f"Label Width: {self.description.rect().width()}") # scale the font size linearly from 16 pt to 8 ps as the width of the panel decreases font_size = int( linear_interpolation(self.description.rect().width(), 12, 16, 400, 600) ) - log_message(f"Label Font Size: {font_size}") + # log_message(f"Label Font Size: {font_size}") self.description.setFont(QFont("Arial", font_size)) self.description.repaint() diff --git a/geest/gui/panels/intro_panel.py b/geest/gui/panels/intro_panel.py index fa03bf01..842d265e 100644 --- a/geest/gui/panels/intro_panel.py +++ b/geest/gui/panels/intro_panel.py @@ -50,11 +50,10 @@ def resizeEvent(self, event): def set_font_size(self): # Scale the font size to fit the text in the available space - log_message(f"Intro Label Width: {self.intro_label.rect().width()}") + # log_message(f"Intro Label Width: {self.intro_label.rect().width()}") # scale the font size linearly from 16 pt to 8 ps as the width of the panel decreases font_size = int( linear_interpolation(self.intro_label.rect().width(), 12, 16, 400, 600) ) - - log_message(f"Intro Label Font Size: {font_size}") + # log_message(f"Intro Label Font Size: {font_size}") self.intro_label.setFont(QFont("Arial", font_size)) diff --git a/geest/gui/panels/open_project_panel.py b/geest/gui/panels/open_project_panel.py index d530a730..bf5ddc5e 100644 --- a/geest/gui/panels/open_project_panel.py +++ b/geest/gui/panels/open_project_panel.py @@ -165,10 +165,10 @@ def resizeEvent(self, event): def set_font_size(self): # Scale the font size to fit the text in the available space - log_message(f"Label Width: {self.label.rect().width()}") + # log_message(f"Label Width: {self.label.rect().width()}") # scale the font size linearly from 16 pt to 8 ps as the width of the panel decreases font_size = int( linear_interpolation(self.label.rect().width(), 12, 16, 400, 600) ) - log_message(f"Label Font Size: {font_size}") + # log_message(f"Label Font Size: {font_size}") self.label.setFont(QFont("Arial", font_size)) diff --git a/geest/gui/panels/ors_panel.py b/geest/gui/panels/ors_panel.py index 3a20a709..96c0eec3 100644 --- a/geest/gui/panels/ors_panel.py +++ b/geest/gui/panels/ors_panel.py @@ -96,10 +96,10 @@ def resizeEvent(self, event): def set_font_size(self): # Scale the font size to fit the text in the available space - log_message(f"Label Width: {self.description.rect().width()}") + # log_message(f"Label Width: {self.description.rect().width()}") # scale the font size linearly from 16 pt to 8 ps as the width of the panel decreases font_size = int( linear_interpolation(self.description.rect().width(), 12, 16, 400, 600) ) - log_message(f"Label Font Size: {font_size}") + # log_message(f"Label Font Size: {font_size}") self.description.setFont(QFont("Arial", font_size)) diff --git a/geest/gui/panels/setup_panel.py b/geest/gui/panels/setup_panel.py index 108c14bc..42f6b5bd 100644 --- a/geest/gui/panels/setup_panel.py +++ b/geest/gui/panels/setup_panel.py @@ -60,10 +60,10 @@ def resizeEvent(self, event): def set_font_size(self): # Scale the font size to fit the text in the available space - log_message(f"Label Width: {self.description.rect().width()}") + # log_message(f"Label Width: {self.description.rect().width()}") # scale the font size linearly from 16 pt to 8 ps as the width of the panel decreases font_size = int( linear_interpolation(self.description.rect().width(), 12, 16, 400, 600) ) - log_message(f"Label Font Size: {font_size}") + # log_message(f"Label Font Size: {font_size}") self.description.setFont(QFont("Arial", font_size)) diff --git a/geest/gui/panels/tree_panel.py b/geest/gui/panels/tree_panel.py index 71efccd7..2e1cc8d2 100644 --- a/geest/gui/panels/tree_panel.py +++ b/geest/gui/panels/tree_panel.py @@ -488,7 +488,6 @@ def update_action_text(): menu.addAction(clear_item_action) menu.addAction(clear_results_action) menu.addAction(run_item_action) - menu.addAction(add_to_map_action) add_wee_score = QAction("Add WEE Score to Map") add_wee_score.triggered.connect( lambda: self.add_to_map( @@ -929,6 +928,7 @@ def add_study_area_to_map(self): "study_area_grid", "study_area_bboxes", "study_area_bbox", + "study_area_creation_status", ] for layer_name in layers: gpkg_layer_path = f"{gpkg_path}|layername={layer_name}" @@ -1574,3 +1574,13 @@ def expand_all_nodes(self, index=None): for row in range(row_count): child_index = self.treeView.model().index(row, 0, index) self.expand_all_nodes(child_index) + + def disable_widgets(self): + """Disable all widgets in the panel.""" + for widget in self.findChildren(QWidget): + widget.setEnabled(False) + + def enable_widgets(self): + """Enable all widgets in the panel.""" + for widget in self.findChildren(QWidget): + widget.setEnabled(True) diff --git a/geest/gui/widgets/configuration_widgets/safety_raster_configuration_widget.py b/geest/gui/widgets/configuration_widgets/safety_raster_configuration_widget.py index 417f1c57..4aa3fb5b 100644 --- a/geest/gui/widgets/configuration_widgets/safety_raster_configuration_widget.py +++ b/geest/gui/widgets/configuration_widgets/safety_raster_configuration_widget.py @@ -3,6 +3,7 @@ ) from .base_configuration_widget import BaseConfigurationWidget from geest.utilities import log_message +from qgis.core import Qgis class SafetyRasterConfigurationWidget(BaseConfigurationWidget): diff --git a/geest/gui/widgets/custom_banner_label.py b/geest/gui/widgets/custom_banner_label.py index dc5b8025..f906c8c8 100644 --- a/geest/gui/widgets/custom_banner_label.py +++ b/geest/gui/widgets/custom_banner_label.py @@ -27,10 +27,10 @@ def paintEvent(self, event): # Scale the font size to fit the text in the available space font_size = 16 threshold = 430 - log_message(f"Banner Label Width: {self.rect().width()}") + # log_message(f"Banner Label Width: {self.rect().width()}") if self.rect().width() < threshold: font_size = int(14 * (self.rect().width() / threshold)) - log_message(f"Font Size: {font_size}") + # log_message(f"Font Size: {font_size}") painter.setFont(QFont("Arial", font_size)) painter.drawText(text_rect, Qt.AlignCenter | Qt.AlignBottom, self.text) diff --git a/geest/resources/qml/study_area_polygons.qml b/geest/resources/qml/study_area_polygons.qml index 38802d9e..5f38d2eb 100644 --- a/geest/resources/qml/study_area_polygons.qml +++ b/geest/resources/qml/study_area_polygons.qml @@ -1,194 +1,50 @@ - - - 1 - 1 - 1 - 0 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + - + - + @@ -198,351 +54,135 @@ + + + + + + + + + + + + + + + + + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 0 0 - 1 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 0 - - - 0 - generatedlayout - - - - - - - - - - - - - - - "area_name" - 2 diff --git a/geest/ui/create_project_panel_base.ui b/geest/ui/create_project_panel_base.ui index c591b223..bcf64b1e 100644 --- a/geest/ui/create_project_panel_base.ui +++ b/geest/ui/create_project_panel_base.ui @@ -13,7 +13,7 @@ Form - + @@ -251,8 +251,8 @@ - - + + @@ -282,8 +282,8 @@ - - + + 0 @@ -295,7 +295,7 @@ - + @@ -325,6 +325,19 @@ + + + + + 0 + 0 + + + + 0 + + + diff --git a/geest/ui/geest_settings_base.ui b/geest/ui/geest_settings_base.ui index 93e8de6d..fea76bb2 100644 --- a/geest/ui/geest_settings_base.ui +++ b/geest/ui/geest_settings_base.ui @@ -13,7 +13,7 @@ Geest Settings - + @@ -27,7 +27,7 @@ - + Open Route Service Options @@ -66,7 +66,36 @@ - + + + + Study Area Preparation Options + + + + + + Chunk size (10 would process in chunks of 10x10 cells) + + + + + + + 10000 + + + 10 + + + 100 + + + + + + + Qt::Vertical @@ -79,7 +108,7 @@ - + Advanced Options diff --git a/geest/ui/project_setup_progress_panel_base.ui b/geest/ui/project_setup_progress_panel_base.ui new file mode 100644 index 00000000..57eb5ec8 --- /dev/null +++ b/geest/ui/project_setup_progress_panel_base.ui @@ -0,0 +1,200 @@ + + + SetupPanelBase + + + + 0 + 0 + 604 + 938 + + + + Form + + + + + + Qt::Horizontal + + + + + + + + 16 + + + + Your project is being created. This could take some time if your analysis covers a large area. For each polygon in your boundary layer, an analysis grid and a set of other geometries needed for the analysis will be created. + + + Qt::MarkdownText + + + Qt::AlignJustify|Qt::AlignTop + + + true + + + + + + + + + + ../resources/geest-banner.png + + + true + + + + + + + + + + 0 + 0 + + + + + 80 + 40 + + + + + 80 + 40 + + + + + 16 + + + + ◀️ + + + + + + + + 0 + 0 + + + + 0 + + + + + + + + 0 + 0 + + + + + 80 + 40 + + + + + 80 + 40 + + + + + 16 + + + + ▶️ + + + + + + + + 0 + 0 + + + + 0 + + + + + + + + + + 300 + 300 + + + + + 300 + 300 + + + + TextLabel + + + + + + + <html><head/><body><p align="center"><span style=" font-size:16pt; font-weight:600;">GEEST Project Setup Progress</span></p></body></html> + + + Qt::RichText + + + Qt::AlignCenter + + + true + + + + + + + Qt::Vertical + + + + 20 + 40 + + + + + + + + + diff --git a/geest/utilities.py b/geest/utilities.py index 7a2b8c94..207ea125 100644 --- a/geest/utilities.py +++ b/geest/utilities.py @@ -21,13 +21,152 @@ import os import logging import inspect -from qgis.PyQt.QtCore import QUrl, QSettings +from datetime import datetime +import tempfile +import re +import platform +import subprocess + +from qgis.PyQt.QtCore import QUrl, QSettings, QRect from qgis.PyQt import uic -from qgis.core import QgsMessageLog, Qgis, QgsProject +from qgis.core import QgsMessageLog, Qgis, QgsProject, QgsLayerTreeGroup from qgis.PyQt.QtWidgets import QApplication +from qgis.core import QgsProject, Qgis from geest.core import setting +def log_window_geometry(geometry): + """ + Creates an ASCII-art diagram of the dialog's dimensions based on the + given geometry (a QRect) and logs it with log_message in QGIS. + + Example output: + + +-------------------- 500 px -------------------+ + | | + | 300 px + | | + +-----------------------------------------------+ + + """ + try: + if type(geometry) == QRect: + rect = geometry + else: + rect = geometry.rect() + except AttributeError: + log_message("Could not get geometry from dialog", level=Qgis.Warning) + log_message(type(geometry), level=Qgis.Warning) + return + + w = rect.width() + h = rect.height() + char_width = 20 - len(str(w)) + top_line = f"\n+{'-'*char_width} {w} px {'-'*20}+" + middle_line = f"|{' ' * 47}{h} px" + bottom_line = f"+{'-'*47}+\n" + + diagram = ( + f"{top_line}\n" + f"| |\n" + f"{middle_line}\n" + f"| |\n" + f"{bottom_line}" + ) + + log_message(diagram) + + +def get_free_memory_mb(): + """ + Attempt to return the free system memory in MB (approx). + Uses only modules from the Python standard library. + """ + system = platform.system() + + # --- Windows --- + if system == "Windows": + try: + import ctypes.wintypes + + class MEMORYSTATUSEX(ctypes.Structure): + _fields_ = [ + ("dwLength", ctypes.wintypes.DWORD), + ("dwMemoryLoad", ctypes.wintypes.DWORD), + ("ullTotalPhys", ctypes.c_ulonglong), + ("ullAvailPhys", ctypes.c_ulonglong), + ("ullTotalPageFile", ctypes.c_ulonglong), + ("ullAvailPageFile", ctypes.c_ulonglong), + ("ullTotalVirtual", ctypes.c_ulonglong), + ("ullAvailVirtual", ctypes.c_ulonglong), + ("ullAvailExtendedVirtual", ctypes.c_ulonglong), + ] + + memoryStatus = MEMORYSTATUSEX() + memoryStatus.dwLength = ctypes.sizeof(MEMORYSTATUSEX) + ctypes.windll.kernel32.GlobalMemoryStatusEx(ctypes.byref(memoryStatus)) + return memoryStatus.ullAvailPhys / (1024 * 1024) + except Exception: + pass + + # --- Linux --- + elif system == "Linux": + # /proc/meminfo is a common place to get memory info on Linux + try: + with open("/proc/meminfo") as f: + meminfo = f.read() + match = re.search(r"^MemAvailable:\s+(\d+)\skB", meminfo, re.MULTILINE) + if match: + return float(match.group(1)) / 1024.0 + except Exception: + pass + + # --- macOS (Darwin) --- + elif system == "Darwin": + # One approach is to parse the output of the 'vm_stat' command + try: + vm_stat = subprocess.check_output(["vm_stat"]).decode("utf-8") + page_size = 4096 # Usually 4096 bytes + free_pages = 0 + # Look for "Pages free: " + match = re.search(r"Pages free:\s+(\d+).", vm_stat) + if match: + free_pages = int(match.group(1)) + return free_pages * page_size / (1024.0 * 1024.0) + except Exception: + pass + + # If none of the above worked or on an unsupported OS, return 0.0 + return 0.0 + + +def log_layer_count(): + """ + Append the number of layers in the project and a timestamp to a text file, + along with free system memory (approximate), using only standard library dependencies. + """ + # Count QGIS layers + layer_count = len(QgsProject.instance().mapLayers()) + + # Gather system free memory (MB) + free_memory_mb = get_free_memory_mb() + + # Create a timestamp + timestamp = datetime.now().strftime("%Y-%m-%d %H:%M:%S") + + # Compose the log entry text + log_entry = f"{timestamp} - Layer count: {layer_count} - Free memory: {free_memory_mb:.2f} MB\n" + + # Send to QGIS log (optional) + log_message(log_entry, level=Qgis.Info, tag="LayerCount") + + # Also write to a log file in the system temp directory + tmp_dir = tempfile.gettempdir() + log_file_path = os.path.join(tmp_dir, "layer_count_log.txt") + with open(log_file_path, "a") as log_file: + log_file.write(log_entry) + + def resources_path(*args): """Get the path to our resources folder. diff --git a/shell.nix b/shell.nix index d2276e87..6443e98e 100644 --- a/shell.nix +++ b/shell.nix @@ -20,6 +20,7 @@ in pkgs.mkShell rec { python3Packages.jsonschema python3Packages.pandas python3Packages.odfpy + python3Packages.psutil # This executes some shell code to initialize a venv in $venvDir before # dropping into the shell diff --git a/start_qgis.sh b/start_qgis.sh index 9fa8d82c..82ba8907 100755 --- a/start_qgis.sh +++ b/start_qgis.sh @@ -10,5 +10,5 @@ esac # Running on local used to skip tests that will not work in a local dev env nix-shell -p \ - 'qgis.override { extraPythonPackages = (ps: [ ps.pyqtwebengine ps.jsonschema ps.debugpy ps.future ]);}' \ + 'qgis.override { extraPythonPackages = (ps: [ ps.pyqtwebengine ps.jsonschema ps.debugpy ps.future ps.psutil ]);}' \ --command "GEEST_DEBUG=${DEBUG_MODE} RUNNING_ON_LOCAL=1 qgis --profile GEEST2" diff --git a/test/test_data/admin/fake_admin_epsg8859.gpkg b/test/test_data/admin/fake_admin_epsg8859.gpkg new file mode 100644 index 00000000..e476da38 Binary files /dev/null and b/test/test_data/admin/fake_admin_epsg8859.gpkg differ diff --git a/test/test_grid_from_bbox.py b/test/test_grid_from_bbox.py new file mode 100644 index 00000000..8e6e9a6d --- /dev/null +++ b/test/test_grid_from_bbox.py @@ -0,0 +1,140 @@ +#!/usr/bin/env python +""" +Test suite for grid_from_bbox.py. + +Version Changed: 2025-01-24 +""" + +from geest.core.tasks import GridFromBbox +from qgis.core import QgsRectangle, QgsCoordinateReferenceSystem + + +def example_bbox(): + """Returns a standard bounding box.""" + return QgsRectangle(0, 0, 100, 100) + + +def example_crs(): + """Returns the Coordinate Reference System.""" + return QgsCoordinateReferenceSystem("EPSG:4326") + + +def test_create_grid_from_bbox_standard(): + """Test grid generation with standard parameters.""" + bbox = example_bbox() + crs = example_crs() + + for cell_width, cell_height in [(10, 10), (5, 5), (20, 20)]: + grid = GridFromBbox( + bbox=bbox, cell_width=cell_width, cell_height=cell_height, crs=crs + ) + + assert grid is not None, "Grid creation failed to return a valid result." + assert grid.isValid(), "Generated grid is not valid." + assert grid.featureCount() > 0, "Grid should contain at least one feature." + + +def test_create_grid_from_bbox_edge_cases(): + """Test edge cases for bounding boxes.""" + crs = example_crs() + + # Zero-area bbox + zero_area_bbox = QgsRectangle(10, 10, 10, 10) + grid = GridFromBbox(bbox=zero_area_bbox, cell_width=5, cell_height=5, crs=crs) + assert grid.featureCount() == 0, "Grid for zero-area bbox should have no features." + + # Inverted bbox + inverted_bbox = QgsRectangle(100, 100, 0, 0) + try: + GridFromBbox(bbox=inverted_bbox, cell_width=10, cell_height=10, crs=crs) + assert False, "Expected ValueError for inverted bbox, but none was raised." + except ValueError as e: + assert "Invalid bounding box dimensions" in str( + e + ), f"Unexpected error message: {e}" + + +def test_invalid_cell_dimensions(): + """Test invalid cell dimensions.""" + bbox = example_bbox() + crs = example_crs() + + # Negative cell dimensions + try: + GridFromBbox(bbox=bbox, cell_width=-10, cell_height=10, crs=crs) + assert ( + False + ), "Expected ValueError for negative cell width, but none was raised." + except ValueError as e: + assert "Cell dimensions must be positive" in str( + e + ), f"Unexpected error message: {e}" + + try: + GridFromBbox(bbox=bbox, cell_width=10, cell_height=-10, crs=crs) + assert ( + False + ), "Expected ValueError for negative cell height, but none was raised." + except ValueError as e: + assert "Cell dimensions must be positive" in str( + e + ), f"Unexpected error message: {e}" + + # Zero cell dimensions + try: + GridFromBbox(bbox=bbox, cell_width=0, cell_height=10, crs=crs) + assert False, "Expected ValueError for zero cell width, but none was raised." + except ValueError as e: + assert "Cell dimensions must be positive" in str( + e + ), f"Unexpected error message: {e}" + + +def test_large_bbox(): + """Test handling of a large bounding box.""" + crs = example_crs() + large_bbox = QgsRectangle(-180, -90, 180, 90) + grid = GridFromBbox(bbox=large_bbox, cell_width=1, cell_height=1, crs=crs) + + assert ( + grid.featureCount() > 10000 + ), "Large bounding box should generate a substantial number of features." + + +def test_crs_mismatch(): + """Test behavior when CRS is mismatched or invalid.""" + bbox = QgsRectangle(0, 0, 100, 100) + + invalid_crs = QgsCoordinateReferenceSystem() + try: + GridFromBbox(bbox=bbox, cell_width=10, cell_height=10, crs=invalid_crs) + assert False, "Expected ValueError for invalid CRS, but none was raised." + except ValueError as e: + assert "Invalid CRS provided" in str(e), f"Unexpected error message: {e}" + + mismatched_crs = QgsCoordinateReferenceSystem("EPSG:3857") + result = GridFromBbox(bbox=bbox, cell_width=10, cell_height=10, crs=mismatched_crs) + + assert result.isValid(), "Grid generation should handle CRS mismatches gracefully." + + +def run_tests(): + """Runs all tests sequentially.""" + test_create_grid_from_bbox_standard() + print("test_create_grid_from_bbox_standard passed.") + + test_create_grid_from_bbox_edge_cases() + print("test_create_grid_from_bbox_edge_cases passed.") + + test_invalid_cell_dimensions() + print("test_invalid_cell_dimensions passed.") + + test_large_bbox() + print("test_large_bbox passed.") + + test_crs_mismatch() + print("test_crs_mismatch passed.") + + +if __name__ == "__main__": + run_tests() diff --git a/test/test_study_area_processing_task.py b/test/test_study_area_processing_task.py index 1d1fabb1..8885fd7e 100644 --- a/test/test_study_area_processing_task.py +++ b/test/test_study_area_processing_task.py @@ -1,204 +1,207 @@ +#!/usr/bin/env python +""" +Test suite for study_area.py + +versionadded: 2025-01-24 +""" + import os import unittest -from qgis.core import ( - QgsVectorLayer, - QgsProcessingContext, - QgsFeedback, -) -from geest.core.tasks import ( - StudyAreaProcessingTask, -) # Adjust the import path as necessary +import shutil + +from osgeo import ogr, gdal +from qgis.core import QgsVectorLayer, QgsFeedback +from geest.core.tasks import StudyAreaProcessingTask from utilities_for_testing import prepare_fixtures -class TestStudyAreaProcessingTask(unittest.TestCase): - """Test suite for the StudyAreaProcessingTask class.""" +class TestStudyAreaProcessor(unittest.TestCase): + """ + Comprehensive test suite for StudyAreaProcessor class. + """ @classmethod def setUpClass(cls): - """Set up shared resources for the test suite.""" - + """ + Setup for the entire test suite. + """ + cls.cell_size_m = 1000 cls.test_data_directory = prepare_fixtures() - cls.working_directory = os.path.join(cls.test_data_directory, "output") - - # Create the output directory if it doesn't exist - if not os.path.exists(cls.working_directory): - os.makedirs(cls.working_directory) - - cls.layer_path = os.path.join( - cls.test_data_directory, "admin", "fake_admin0.gpkg" + cls.input_admin_path = os.path.join( + cls.test_data_directory, "admin", "Admin0.shp" ) + cls.layer = QgsVectorLayer(cls.input_admin_path, "Admin0", "ogr") cls.field_name = "Name" - cls.cell_size_m = 100 - cls.context = QgsProcessingContext() - cls.feedback = QgsFeedback() - - # Ensure the working directory exists - if not os.path.exists(cls.working_directory): - os.makedirs(cls.working_directory) + # Define working directories + cls.working_dir = os.path.join(cls.test_data_directory, "output") + cls.gpkg_path = os.path.join(cls.working_dir, "study_area", "study_area.gpkg") - # Load the test layer - cls.layer = QgsVectorLayer(cls.layer_path, "Test Layer", "ogr") - if not cls.layer.isValid(): - raise RuntimeError(f"Failed to load test layer from {cls.layer_path}") + @classmethod + def tearDownClass(cls): + """ + Cleanup after all tests. + """ + pass def setUp(self): - """Set up test-specific resources.""" - # Clean up the output directory before each test - for filename in os.listdir(self.working_directory): - file_path = os.path.join(self.working_directory, filename) - if os.path.isfile(file_path): - os.remove(file_path) - - def test_initialization(self): - """Test initialization of the task.""" - task = StudyAreaProcessingTask( - layer=self.layer, - field_name=self.field_name, - cell_size_m=self.cell_size_m, - working_dir=self.working_directory, - context=self.context, - feedback=self.feedback, - ) - - self.assertEqual(task.layer, self.layer) - self.assertEqual(task.field_name, self.field_name) - self.assertEqual(task.cell_size_m, self.cell_size_m) - self.assertEqual(task.working_dir, self.working_directory) - self.assertTrue( - os.path.exists(os.path.join(self.working_directory, "study_area")) - ) - - def test_calculate_utm_zone(self): - """Test UTM zone calculation.""" - task = StudyAreaProcessingTask( - layer=self.layer, - field_name=self.field_name, - cell_size_m=self.cell_size_m, - working_dir=self.working_directory, - context=self.context, - feedback=self.feedback, - ) - - bbox = self.layer.extent() - utm_zone = task.calculate_utm_zone(bbox) - - # Validate the calculated UTM zone (adjust based on test data location) - self.assertTrue(utm_zone >= 32600 or utm_zone >= 32700) - - def test_process_study_area(self): - """Test processing of study area features.""" - task = StudyAreaProcessingTask( - layer=self.layer, - field_name=self.field_name, - cell_size_m=self.cell_size_m, - working_dir=self.working_directory, - context=self.context, - feedback=self.feedback, - ) - - result = task.process_study_area() - self.assertTrue(result) - - # Validate output GeoPackage - gpkg_path = os.path.join( - self.working_directory, "study_area", "study_area.gpkg" - ) - self.assertTrue( - os.path.exists(gpkg_path), - msg=f"GeoPackage not created in {self.working_directory}", - ) + """ + Set up the environment for the test, loading the test data layers. + """ + # Create the output directory if it doesn't exist + if not os.path.exists(self.working_dir): + # print(f"Creating working directory {self.working_dir}") + os.makedirs(self.working_dir) - def test_process_singlepart_geometry(self): - """Test processing of singlepart geometry.""" - task = StudyAreaProcessingTask( + # Define paths to test layers + self.processor = StudyAreaProcessingTask( layer=self.layer, field_name=self.field_name, + working_dir=self.working_dir, cell_size_m=self.cell_size_m, - working_dir=self.working_directory, - context=self.context, - feedback=self.feedback, - ) + crs=None, + ) + + def tearDown(self): + # Recursively delete everything in the working_dir + shutil.rmtree(self.working_dir) + + def test_create_layer_if_not_exists(self): + """ + Test creating a layer if it does not exist. + """ + layer_name = "test_layer" + self.processor.create_layer_if_not_exists(layer_name) + + ds = ogr.Open(self.gpkg_path) + layer = ds.GetLayerByName(layer_name) + self.assertIsNotNone(layer, "Layer was not created.") + ds = None + + @unittest.skip("Skipping test for now") + def test_create_and_save_grid(self): + """ + Test grid creation and saving. + """ + bbox = (-10, 10, -10, 10) + geom = ogr.Geometry(ogr.wkbPolygon) + ring = ogr.Geometry(ogr.wkbLinearRing) + ring.AddPoint(-10, -10) + ring.AddPoint(10, -10) + ring.AddPoint(10, 10) + ring.AddPoint(-10, 10) + ring.AddPoint(-10, -10) + geom.AddGeometry(ring) + + self.processor.create_and_save_grid("test_grid", geom, bbox) + + ds = ogr.Open(self.gpkg_path) + layer = ds.GetLayerByName("study_area_grid") + self.assertIsNotNone(layer, "Grid layer was not created.") + self.assertGreater(layer.GetFeatureCount(), 0, "Grid layer has no features.") + ds = None + + @unittest.skip("Skipping test for now") + def test_create_raster_mask(self): + """ + Test raster mask creation. + """ + bbox = (-10, 10, -10, 10) + geom = ogr.Geometry(ogr.wkbPolygon) + ring = ogr.Geometry(ogr.wkbLinearRing) + ring.AddPoint(-10, -10) + ring.AddPoint(10, -10) + ring.AddPoint(10, 10) + ring.AddPoint(-10, 10) + ring.AddPoint(-10, -10) + geom.AddGeometry(ring) + + mask_path = self.processor.create_raster_mask(geom, bbox, "test_mask") + + self.assertTrue(os.path.exists(mask_path), "Raster mask was not created.") - feature = next(self.layer.getFeatures()) - task.process_singlepart_geometry(feature.geometry(), "test_area", "Test Area") + def test_calculate_utm_zone(self): + """ + Test UTM zone calculation. + """ + bbox = (-10, 10, -10, 10) + utm_code = self.processor.calculate_utm_zone(bbox) + self.assertEqual(utm_code, 32631, "UTM zone calculation is incorrect.") + + def test_create_study_area_directory(self): + """ + Test study area directory creation. + """ + self.processor.create_study_area_directory(self.working_dir) - # Validate GeoPackage outputs - gpkg_path = os.path.join( - self.working_directory, "study_area", "study_area.gpkg" - ) - self.assertTrue( - os.path.exists(gpkg_path), - msg=f"GeoPackage not created in {self.working_directory}", - ) - # Validate mask is a valid file - mask_path = os.path.join( - self.working_directory, "study_area", "saint_lucia_part0.tif" - ) self.assertTrue( - os.path.exists(mask_path), - msg=f"mask saint_lucia_part0.tif not created in {mask_path}", - ) - - def test_grid_aligned_bbox(self): - """Test grid alignment of bounding boxes.""" - task = StudyAreaProcessingTask( - layer=self.layer, - field_name=self.field_name, - cell_size_m=self.cell_size_m, - working_dir=self.working_directory, - context=self.context, - feedback=self.feedback, - ) - - bbox = self.layer.extent() - aligned_bbox = task.grid_aligned_bbox(bbox) - - # Validate grid alignment - self.assertAlmostEqual( - (aligned_bbox.xMaximum() - aligned_bbox.xMinimum()) % self.cell_size_m, 0 - ) - self.assertAlmostEqual( - (aligned_bbox.yMaximum() - aligned_bbox.yMinimum()) % self.cell_size_m, 0 - ) - + os.path.exists(os.path.join(self.working_dir, "study_area")), + "Study area directory was not created.", + ) + + def test_track_time(self): + """ + Test the track_time helper method. + """ + self.processor.metrics = {"test_metric": 0} + start_time = 0 + self.processor.track_time("test_metric", start_time) + self.assertIn("test_metric", self.processor.metrics, "Metric was not tracked.") + + @unittest.skip("Skipping test for now") + def test_write_chunk(self): + """ + Test the write_chunk method for layer writing. + """ + # Create dummy layer and task + self.processor.create_layer_if_not_exists("chunk_layer") + ds = ogr.Open(self.gpkg_path, 1) + layer = ds.GetLayerByName("chunk_layer") + self.assertIsNotNone(layer, "Chunk layer creation failed.") + task = type("DummyTask", (), {"features_out": [ogr.Geometry(ogr.wkbPolygon)]})() + task.features_out[0].AddGeometry(ogr.Geometry(ogr.wkbLinearRing)) + self.processor.write_chunk(layer, task, "test_chunk") + self.assertGreater(layer.GetFeatureCount(), 0, "Chunk writing failed.") + + @unittest.skip("Skipping test for now") def test_create_raster_vrt(self): - """Test creation of a VRT from raster masks.""" - task = StudyAreaProcessingTask( - layer=self.layer, - field_name=self.field_name, - cell_size_m=self.cell_size_m, - working_dir=self.working_directory, - context=self.context, - feedback=self.feedback, - ) - - # Generate raster masks - task.process_study_area() - - # Create VRT - task.create_raster_vrt() - - # Validate VRT file - vrt_path = os.path.join( - self.working_directory, "study_area", "combined_mask.vrt" - ) + """ + Test VRT creation. + """ + vrt_name = "test_vrt.vrt" + self.processor.create_raster_vrt(output_vrt_name=vrt_name) self.assertTrue( - os.path.exists(vrt_path), - msg=f"VRT file not created in {self.working_directory}", + os.path.exists(os.path.join(self.working_dir, "study_area", vrt_name)), + "VRT file was not created.", + ) + + @unittest.skip("Skipping test for now") + def test_create_clip_polygon(self): + """ + Test clip polygon creation. + """ + bbox = (-10, 10, -10, 10) + aligned_box = (-12, 12, -12, 12) + + # Create a test geometry + geom = ogr.Geometry(ogr.wkbPolygon) + ring = ogr.Geometry(ogr.wkbLinearRing) + ring.AddPoint(-10, -10) + ring.AddPoint(10, -10) + ring.AddPoint(10, 10) + ring.AddPoint(-10, 10) + ring.AddPoint(-10, -10) + geom.AddGeometry(ring) + + # Call create_clip_polygon + self.processor.create_clip_polygon(geom, aligned_box, "test_clip") + + ds = ogr.Open(self.gpkg_path) + layer = ds.GetLayerByName("study_area_clip_polygons") + self.assertIsNotNone(layer, "Clip polygon layer was not created.") + self.assertGreater( + layer.GetFeatureCount(), 0, "Clip polygon layer has no features." ) - @classmethod - def tearDownClass(cls): - """Clean up shared resources.""" - cleanup = False - if os.path.exists(cls.working_directory) and cleanup: - for root, dirs, files in os.walk(cls.working_directory, topdown=False): - for name in files: - os.remove(os.path.join(root, name)) - for name in dirs: - os.rmdir(os.path.join(root, name)) - if __name__ == "__main__": unittest.main()