CUTBLOCK TIMBER TYPE AUTOMATED QUALITY CONTROL AND SPATIAL MISMATCH RESOLUTION
Nov 26, 2024
8 min read
"""Workflow Context:
This script is part of a larger forest inventory process where:
1. Net areas are first delineated from GPS field data
2. Timber type polygons are generated through a separate classification process
3. Cruise plots are established and classified according to type and maturity.
4. This script then validates the geometric relationships and summarizes the attributes of these three layers to ensure data consistency and completeness.This script ensures that timber type polygons exactly match net area geometries and provides summary statistics for cruising and cruise compilation purposes. It ensures that each cutblock is a single polygon or a multipart polygon and that the sum of the timber type areas for each block is equal to the block's net area. It provides cutblock net areas, timber type areas, and number of plots per timber type.The script performs these operations:1. Dissolves Net Areas to ensure one polygon per block2. Updates block IDs in Timber Types based on spatial overlap3. Clips Timber Types layer using dissolved Net Areas4. Updates area fields (areaHa) in both layers5. Validates areas (all areas rounded to 2 decimal places)6. Creates detailed difference analysis layers with zoom capabilities7. Explains and locates area mismatches for each Timber Type polygon8. Analyzes cruise plot distribution, including plot types and maturity classes per timber typeRequired Input Layers: - 'Net Area planned': GPS-derived cutblock boundaries- 'Timber Types': Forest stand polygons with classifications- 'Cruise grid plots': Classified sample plots with PLOT TYPE AND MAT attributesRequired Fields:- 'block' in Net Area planned- 'TT Label' in Timber Types- 'PLOT TYPE AND MAT' in Cruise grid plotsDate: November 12, 2024"""from qgis.core import ( QgsProject, QgsVectorLayer, QgsExpression, QgsField, QgsFeature, QgsExpressionContext, QgsExpressionContextUtils, QgsFeatureRequest, edit)from PyQt5.QtCore import QVariantfrom PyQt5.QtWidgets import QInputDialog, QLineEditimport processingfrom collections import defaultdictimport redef get_layer_name_input(default_name, prompt): """Request layer name input from user with default value""" layer_name, ok = QInputDialog.getText(None, "Input Layer", prompt, QLineEdit.Normal, default_name) if ok and layer_name: return layer_name return default_namedef update_layer_area(layer): """Update area field in hectares for given layer""" area_field_names = ['areaHa', 'areaha'] existing_field = None fields = layer.fields() for field_name in area_field_names: if fields.indexOf(field_name) != -1: existing_field = field_name break if not existing_field: layer.startEditing() field = QgsField('areaHa', QVariant.Double, 'double', 10, 2) layer.dataProvider().addAttributes([field]) layer.updateFields() existing_field = 'areaHa' field_index = layer.fields().indexOf(existing_field) for feature in layer.getFeatures(): area_ha = round(feature.geometry().area() / 10000.0, 2) layer.changeAttributeValue(feature.id(), field_index, area_ha) return existing_fielddef roman_to_int(roman_numeral): """Convert Roman numeral to integer for sorting""" roman_dict = {'I': 1, 'II': 2, 'III': 3, 'IV': 4, 'V': 5, 'VI': 6, 'VII': 7, 'VIII': 8, 'IX': 9, 'X': 10, 'XI': 11, 'XII': 12, 'XIII': 13, 'XIV': 14, 'XV': 15} return roman_dict.get(roman_numeral.upper(), 0)def extract_roman_numeral(text): """Extract Roman numeral from text for sorting""" match = re.search(r'\b([IVXLCDM]+)\b', text.upper()) if match: return match.group(1) return ''def print_area_summaries(timber_layer, net_area_layer, area_field): """Print summary tables of areas and validate against NetArea planned""" print("\n=== AREA SUMMARIES AND VALIDATION ===\n") # Get areas from NetArea planned net_area_blocks = {} for feature in net_area_layer.getFeatures(): block = feature['block'] area = round(feature.geometry().area() / 10000.0, 2) net_area_blocks[block] = area block_areas = defaultdict(float) type_areas = defaultdict(float) block_type_areas = defaultdict(lambda: defaultdict(float)) for feature in timber_layer.getFeatures(): block = feature['block'] timber_type = feature['TT Label'] area = feature[area_field] if block and timber_type: block_areas[block] += area type_areas[timber_type] += area block_type_areas[block][timber_type] += area # Round accumulated areas block_areas = {k: round(v, 2) for k, v in block_areas.items()} type_areas = {k: round(v, 2) for k, v in type_areas.items()} block_type_areas = {k: {tk: round(tv, 2) for tk, tv in v.items()} for k, v in block_type_areas.items()} # Area Validation by Block print("--- Area Validation by Block ---") print(f"{'Block':<10} {'Timber Area (ha)':<20} {'NetArea (ha)':<20} {'Difference (ha)':<20} {'Status':<10}") print("-" * 80) total_timber = 0 total_net = 0 for block in sorted(net_area_blocks.keys()): timber_area = block_areas.get(block, 0) net_area = net_area_blocks[block] diff = round(abs(timber_area - net_area), 2) status = "OK" if diff < 0.1 else "CHECK" print(f"{str(block):<10} {timber_area:<20.2f} {net_area:<20.2f} {diff:<20.2f} {status:<10}") total_timber += timber_area total_net += net_area print("-" * 80) total_timber = round(total_timber, 2) total_net = round(total_net, 2) total_diff = round(abs(total_timber - total_net), 2) print(f"{'Total':<10} {total_timber:<20.2f} {total_net:<20.2f} {total_diff:<20.2f}") # Print Area per Timber Type per Block print("\n--- Area per Timber Type per Block ---") print(f"{'Block':<10} {'Timber Type':<20} {'Area (ha)':<15}") print("-" * 50) for block in sorted(block_type_areas.keys()): # Sort timber types per Roman numerals sorted_timber_types = sorted(block_type_areas[block].keys(), key=lambda x: roman_to_int(extract_roman_numeral(x))) for timber_type in sorted_timber_types: area = block_type_areas[block][timber_type] print(f"{block:<10} {timber_type:<20} {area:<15.2f}")def report_discrepancies(timber_layer, holes_layer, overshoots_layer): """Report and locate area mismatches for each Timber Type polygon""" print("\n--- Detailed Mismatch Report for Each Timber Type Polygon ---") print(f"{'Feature ID':<12} {'TT Label':<20} {'Block':<10} {'Mismatch Type':<30} {'Mismatch Area (ha)':<20}") print("-" * 90) total_mismatch_area = 0 mismatch_count = 0 # Prepare a dictionary to store mismatches per feature mismatches = defaultdict(list) # Process Holes (Missing Timber Types) if holes_layer.featureCount() > 0: for hole in holes_layer.getFeatures(): if hole['area_ha'] > 0.0: request = QgsFeatureRequest().setFilterRect(hole.geometry().boundingBox()) for timber_feature in timber_layer.getFeatures(request): if timber_feature.geometry().intersects(hole.geometry()): mismatches[timber_feature.id()].append(('Missing Coverage (Hole)', round(hole['area_ha'], 2))) # Process Overshoots (Extending Beyond Blocks) if overshoots_layer.featureCount() > 0: for overshoot in overshoots_layer.getFeatures(): if overshoot['area_ha'] > 0.0: request = QgsFeatureRequest().setFilterRect(overshoot.geometry().boundingBox()) for timber_feature in timber_layer.getFeatures(request): if timber_feature.geometry().intersects(overshoot.geometry()): mismatches[timber_feature.id()].append(('Extending Beyond Block (Overshoot)', round(overshoot['area_ha'], 2))) # Iterate through mismatches and print details for feature_id, issues in mismatches.items(): timber_feature = timber_layer.getFeature(feature_id) tt_label = timber_feature['TT Label'] block = timber_feature['block'] for issue_type, issue_area in issues: print(f"{feature_id:<12} {tt_label:<20} {block:<10} {issue_type:<30} {issue_area:<20.2f}") total_mismatch_area += issue_area mismatch_count += 1 if mismatch_count == 0: print("No mismatches found.") else: print("-" * 90) print(f"{'Total':<64} {round(total_mismatch_area,2):<20.2f}") print(f"Number of mismatches: {mismatch_count}")def analyze_cruise_plots(timber_layer, plots_layer): """Analyze cruise plot distribution across timber types""" print("\n=== CRUISE PLOT ANALYSIS ===\n") # Initialize nested counters for each timber type plot_counts = defaultdict(int) plots_outside = [] total_plots = 0 plot_type_by_tt = defaultdict(lambda: defaultdict(int)) # Nested dict for plot type counts by timber type # Spatial join to count plots per timber type joined_result = processing.run("native:joinattributesbylocation", { 'INPUT': plots_layer, 'JOIN': timber_layer, 'PREDICATE': [0], # intersects 'JOIN_FIELDS': ['TT Label', 'areaHa', 'block'], 'METHOD': 0, # one-to-many 'DISCARD_NONMATCHING': False, 'PREFIX': 'tt_', 'OUTPUT': 'memory:joined_plots' }) joined_plots = joined_result['OUTPUT'] # Count plots and identify those outside timber types for plot in joined_plots.getFeatures(): total_plots += 1 tt_label = plot['tt_TT Label'] block = plot['tt_block'] if tt_label: plot_counts[(block, tt_label)] += 1 # Count plot type and mat for this timber type plot_type_mat = plot['PLOT TYPE AND MAT'] if plot_type_mat: plot_type_by_tt[(block, tt_label)][plot_type_mat] += 1 else: plots_outside.append(plot.id()) # Get areas for each timber type per block block_type_areas = defaultdict(float) for feature in timber_layer.getFeatures(): block = feature['block'] tt_label = feature['TT Label'] area = feature['areaHa'] block_type_areas[(block, tt_label)] += area # Print Plot Counts per Timber Type per Block print("--- Number of Plots per Timber Type per Block ---") print(f"{'Block':<10} {'Timber Type':<20} {'Area (ha)':<10} {'Plot Count':<12} {'Plots/ha':<10}") print("-" * 70) # Sort keys based on block and Roman numeral in timber type sorted_keys = sorted(plot_counts.keys(), key=lambda x: (x[0], roman_to_int(extract_roman_numeral(x[1])))) for (block, tt_label) in sorted_keys: count = plot_counts[(block, tt_label)] area = block_type_areas.get((block, tt_label), 0) plots_per_ha = round(count / area, 2) if area > 0 else 0 print(f"{block:<10} {tt_label:<20} {area:<10.2f} {count:<12d} {plots_per_ha:<10.2f}") # Print 'PLOT TYPE AND MAT' Counts per Timber Type per Block print("\n--- 'PLOT TYPE AND MAT' Counts per Timber Type per Block ---") print(f"{'Block':<10} {'Timber Type':<20} {'PLOT TYPE AND MAT':<25} {'Count':<10}") print("-" * 80) for (block, tt_label) in sorted_keys: plot_types = plot_type_by_tt[(block, tt_label)] # Sort plot types alphabetically for plot_type_mat in sorted(plot_types.keys()): count = plot_types[plot_type_mat] print(f"{block:<10} {tt_label:<20} {plot_type_mat:<25} {count:<10d}") print("-" * 80) print(f"Total Plots: {total_plots}") print(f"Plots Outside Timber Types: {len(plots_outside)}") if plots_outside: print("\nPlots Outside Timber Types (Feature IDs):") print(", ".join(map(str, plots_outside)))def dissolve_and_update(): """Main function to process and update forest management layers""" # Get layer names from user input net_area_layer_name = get_layer_name_input("NetArea planned", "Enter cutblock net area layer name:") timber_layer_name = get_layer_name_input("Timber Types", "Enter timber types layer name:") plots_layer_name = get_layer_name_input("Cruise grid plots", "Enter cruise grid plots layer name:") try: # Get project instance and layers project = QgsProject.instance() timber_layers = project.mapLayersByName(timber_layer_name) net_area_layers = project.mapLayersByName(net_area_layer_name) plots_layers = project.mapLayersByName(plots_layer_name) if not timber_layers or not net_area_layers or not plots_layers: print("Error: Could not find one or more required layers") print("Available layers:", [layer.name() for layer in project.mapLayers().values()]) return timber_layer = timber_layers[0] net_area_layer = net_area_layers[0] plots_layer = plots_layers[0] # Step 1: Check and report initial state print("Initial state of NetArea planned:") print(f"{'Block':<6} | {'Number of polygons'}") print("-" * 25) block_counts = {} for feature in net_area_layer.getFeatures(): block = feature['block'] block_counts[block] = block_counts.get(block, 0) + 1 for block, count in sorted(block_counts.items()): print(f"{block:<6} | {count}") # Step 2: Dissolve NetArea planned print("\nDissolving NetArea planned by block...") dissolved_result = processing.run("native:dissolve", { 'INPUT': net_area_layer, 'FIELD': ['block'], 'OUTPUT': 'memory:dissolved' }) dissolved_layer = dissolved_result['OUTPUT'] # Update the original NetArea planned layer net_area_layer.startEditing() net_area_layer.dataProvider().truncate() net_area_layer.dataProvider().addFeatures([f for f in dissolved_layer.getFeatures()]) net_area_layer.commitChanges() # Step 3: Report NetArea planned after dissolve print("\nNetArea planned after dissolve:") print(f"{'Block':<6} | {'Area (ha)'}") print("-" * 25) for feature in net_area_layer.getFeatures(): block = feature['block'] area_ha = round(feature.geometry().area() / 10000.0, 2) print(f"{block:<6} | {area_ha:.2f}") # Step 4: Spatial join with dissolved NetArea planned print("\nUpdating Timber Types with block information...") joined_result = processing.run("native:joinattributesbylocation", { 'INPUT': timber_layer, 'JOIN': net_area_layer, 'PREDICATE': [0], 'JOIN_FIELDS': ['block'], 'METHOD': 1, 'DISCARD_NONMATCHING': False, 'PREFIX': 'join_', 'OUTPUT': 'memory:joined' }) joined_layer = joined_result['OUTPUT'] # Update block field in Timber Types with edit(joined_layer): for feature in joined_layer.getFeatures(): if feature['join_block'] is not None: feature['block'] = feature['join_block'] joined_layer.updateFeature(feature) # Clean up joined fields joined_layer.startEditing() joined_layer.dataProvider().deleteAttributes([joined_layer.fields().indexOf('join_block')]) joined_layer.updateFields() joined_layer.commitChanges() # Step 5: Clip with dissolved NetArea print("Clipping Timber Types...") clipped_result = processing.run("native:clip", { 'INPUT': joined_layer, 'OVERLAY': net_area_layer, 'OUTPUT': 'memory:clipped' }) clipped_layer = clipped_result['OUTPUT'] # Update original layer with clipped results timber_layer.startEditing() timber_layer.dataProvider().truncate() timber_layer.dataProvider().addFeatures([f for f in clipped_layer.getFeatures()]) # Step 6: Update areas print("Updating areas...") area_field = update_layer_area(timber_layer) update_layer_area(net_area_layer) # Commit changes and refresh timber_layer.commitChanges() timber_layer.triggerRepaint() net_area_layer.triggerRepaint() # Step 7: Print area summaries print_area_summaries(timber_layer, net_area_layer, area_field) # Step 8: Difference Analysis print("\n=== SPATIAL DISCREPANCY ANALYSIS ===\n") # Create holes layer print("Analyzing areas where Timber Types don't cover NetArea planned...") holes_result = processing.run("native:difference", { 'INPUT': net_area_layer, 'OVERLAY': timber_layer, 'OUTPUT': 'memory:holes' }) holes_layer = holes_result['OUTPUT'] holes_layer.startEditing() holes_layer.dataProvider().addAttributes([ QgsField('area_ha', QVariant.Double, 'double', 10, 2), QgsField('zoom_scale', QVariant.Double, 'double', 10, 0) # Add zoom scale field ]) holes_layer.updateFields() area_idx = holes_layer.fields().indexOf('area_ha') zoom_idx = holes_layer.fields().indexOf('zoom_scale') for feature in holes_layer.getFeatures(): area_ha = round(feature.geometry().area() / 10000.0, 2) # Set zoom scale based on feature size - smaller features get closer zoom zoom_scale = min(1000, max(100, area_ha * 200)) # Adjust these values as needed holes_layer.changeAttributeValue(feature.id(), area_idx, area_ha) holes_layer.changeAttributeValue(feature.id(), zoom_idx, zoom_scale) holes_layer.commitChanges() # Create overshoots layer print("\nAnalyzing areas where Timber Types extend beyond NetArea planned...") overshoots_result = processing.run("native:difference", { 'INPUT': timber_layer, 'OVERLAY': net_area_layer, 'OUTPUT': 'memory:overshoots' }) overshoots_layer = overshoots_result['OUTPUT'] overshoots_layer.startEditing() overshoots_layer.dataProvider().addAttributes([ QgsField('area_ha', QVariant.Double, 'double', 10, 2), QgsField('zoom_scale', QVariant.Double, 'double', 10, 0) # Add zoom scale field ]) overshoots_layer.updateFields() area_idx_overshoot = overshoots_layer.fields().indexOf('area_ha') zoom_idx_overshoot = overshoots_layer.fields().indexOf('zoom_scale') for feature in overshoots_layer.getFeatures(): area_ha = round(feature.geometry().area() / 10000.0, 2) # Set zoom scale based on feature size - smaller features get closer zoom zoom_scale = min(1000, max(100, area_ha * 200)) # Adjust these values as needed overshoots_layer.changeAttributeValue(feature.id(), area_idx_overshoot, area_ha) overshoots_layer.changeAttributeValue(feature.id(), zoom_idx_overshoot, zoom_scale) overshoots_layer.commitChanges() # Add layers to project project.addMapLayer(holes_layer) holes_layer.setName("Holes_to_fix") project.addMapLayer(overshoots_layer) overshoots_layer.setName("Overshoots_to_fix") # Step 9: Report Discrepancies report_discrepancies(timber_layer, holes_layer, overshoots_layer) # Step 10: Analyze cruise plots analyze_cruise_plots(timber_layer, plots_layer) # Final Report print("\nDiscrepancy layers have been added to the map:") print("- 'Holes_to_fix': Individual areas where Timber Types are missing") print(" * Use block field to identify which block contains the hole") print(" * area_ha field shows the size of each hole") print(" * zoom_scale field can be used to set appropriate zoom level") print("\n- 'Overshoots_to_fix': Individual areas where Timber Types extend beyond blocks") print(" * Contains original Timber Types attributes") print(" * area_ha field shows the size of each overshoot") print(" * zoom_scale field can be used to set appropriate zoom level") print("\nTo use the zoom feature effectively:") print("1. Right-click either layer and go to Properties -> Actions") print("2. Add a new action:") print(" - Type: Python") print(" - Name: Zoom to Feature") print(" - Action: canvas = iface.mapCanvas(); canvas.zoomScale([% \"zoom_scale\" %])") print("\nAll operations completed successfully!") except Exception as e: print(f"Error during operation: {str(e)}") if timber_layer.isEditable(): timber_layer.rollBack() if net_area_layer.isEditable(): net_area_layer.rollBack() print("Changes rolled back due to error") print("Detailed error:", str(e))# Run the functiondissolve_and_update()
Comments