Add introduction comment

.gitignore

@@ -9,6 +9,7 @@ standalone_charm.py
 
 # Other Files
 setting_up_standalone_pyqgis.docx
+symbology-style.db
 
 
 .idea/.name

mtl_buildings_workflow.py

@@ -1,8 +1,21 @@
-from qgis.core import QgsApplication, QgsVectorLayer, QgsProject, QgsProcessingFeedback
+"""
+This script automates the cleaning process of Montreal
+Buildings Geometries.
+Developed by Alireza Adli
+"""
+
+
+from qgis.core import QgsApplication, QgsVectorLayer, QgsProject, \
+  QgsProcessingFeedback
 from qgis.analysis import QgsNativeAlgorithms
 import processing
+import os
+from services_scripts.basic_functions import create_folders, \
+  find_shp_files
 
 
+# This function loads a layer.
+# It's used to count the data records
 def load_layer(path, layer_name):
   the_layer = QgsVectorLayer(path, layer_name, "ogr")
   if not the_layer.isValid():
@@ -23,89 +36,373 @@ qgs.initQgis()
 # Add native algorithms provider
 QgsApplication.processingRegistry().addProvider(QgsNativeAlgorithms())
 
-nrcan_0 = 'C:/Users/a_adli/PycharmProjects/hydroquebec_archetype_gispy/data/nrcan/Autobuilding_QC_VILLE_MONTREAL.shp'
-nrcan, nrcan_name = load_layer(nrcan_0, 'NRCan')
+# Reading the Automatically Extracted Buildings Footprints layer (NRCan)
+# and counting its data records
+nrcan_layer = 'C:/Users/a_adli/PycharmProjects/' \
+              'hydroquebec_archetype_gispy/data/' \
+              'input_data/nrcan/Autobuilding_QC_VILLE_MONTREAL.shp'
+nrcan, nrcan_name = load_layer(nrcan_layer, 'NRCan')
 
 print(f'Loading {nrcan_name}')
 print(f'{nrcan_name} data count: {nrcan.featureCount()}')
 
+# Fixing geometries of the NRCan layer
 print(f'Fixing {nrcan_name} geometries')
+fixed_nrcan_layer = 'C:/Users/a_adli/PycharmProjects/' \
+                    'hydroquebec_archetype_gispy/data/' \
+                    'output_data/fixed_nrcan/fixed_nrcan.shp'
+params_fixing_nrcan = {'INPUT': nrcan_layer,
+                       'METHOD': 0,
+                       'OUTPUT': fixed_nrcan_layer}
+processing.run('native:fixgeometries', params_fixing_nrcan)
 
-fixed_nrcan_0 = 'C:/Users/a_adli/PycharmProjects/hydroquebec_archetype_gispy/data/new_tests/python_fixed_05/py_fixes_05.shp'
-
-params_fixing_nrcan = {'INPUT': nrcan_0, 'METHOD': 0, 'OUTPUT': fixed_nrcan_0}
-fixed_layer_nrcan = processing.run('native:fixgeometries', params_fixing_nrcan)['OUTPUT']
-
-fixed_nrcan, fixed_nrcan_name = load_layer(fixed_layer_nrcan, 'Fixed NRCan')
+fixed_nrcan, fixed_nrcan_name = load_layer(fixed_nrcan_layer, 'Fixed NRCan')
 print(f'{fixed_nrcan_name} data count: {fixed_nrcan.featureCount()}')
 
+# Removing unnecessary parts of the NRCan layer (outside of Montreal) using
+# Administrative boundaries of the Montreal agglomeration data
+clipping_montreal_boundary_layer = 'C:/Users/a_adli/PycharmProjects/' \
+                                   'hydroquebec_archetype_gispy/data/' \
+                                   'input_data/montreal_boundary/' \
+                                   'Montreal_boundary.shp'
+clipped_nrcan_layer = 'C:/Users/a_adli/PycharmProjects/' \
+                      'hydroquebec_archetype_gispy/data/' \
+                      'output_data/clipped_fixed_nrcan/clipped_fixed_nrcan.shp'
 
-clipping_montreal_boundary_layer = 'C:/Users/a_adli/PycharmProjects/hydroquebec_archetype_gispy/data/montreal_boundary|layername=Montreal_boundary'
-clipped_nrcan_layer = 'C:/Users/a_adli/PycharmProjects/hydroquebec_archetype_gispy/data/new_tests/clipped_nrcan_1/clipped_nrcan_1.shp'
-
-params_clipping_nrcan = {'INPUT': fixed_nrcan_0, 'OVERLAY': clipping_montreal_boundary_layer, 'FILTER_EXPRESSION': '', 'FILTER_EXTENT': None, 'OUTPUT': clipped_nrcan_layer}
-clipped_layer_nrcan = processing.run("native:clip", params_clipping_nrcan)['OUTPUT']
-
+params_clipping_nrcan = {'INPUT': fixed_nrcan_layer,
+                         'OVERLAY': clipping_montreal_boundary_layer,
+                         'FILTER_EXPRESSION': '',
+                         'FILTER_EXTENT': None,
+                         'OUTPUT': clipped_nrcan_layer}
+processing.run("native:clip", params_clipping_nrcan)
 print(f'Clipping of {fixed_nrcan_name} is completed.')
 
-clipped_nrcan, clipped_nrcan_name = load_layer(clipped_layer_nrcan, 'Clipped NRCan')
-print(f'{clipped_nrcan_name} data count: {clipped_nrcan.featureCount()}')
+# After each step, data records are counted
+clipped_nrcan, clipped_nrcan_name = \
+  load_layer(clipped_nrcan_layer, 'Clipped NRCan')
+clipped_nrcan_length = clipped_nrcan.featureCount()
+print(f'{clipped_nrcan_name} data count: {clipped_nrcan_length}')
 
 
-params_create_index_nrcan = {'INPUT': clipped_nrcan, 'OUTPUT': 'Output'}
-indexed_layer = processing.run("native:createspatialindex", params_create_index_nrcan)
-
+# Creating spatial index is needed for most of the newly built layers.
+# It adds a file with suffix .qix
+# Without the index, most of the processes cannot be run
+# (smoothly and efficiently)
+params_create_index_nrcan = {'INPUT': clipped_nrcan_layer,
+                             'OUTPUT': 'Output'}
+processing.run("native:createspatialindex", params_create_index_nrcan)
 print(f'Creating spatial index for {clipped_nrcan_name} is completed.')
 
-geoindex_layer = 'C:/Users/a_adli/PycharmProjects/hydroquebec_archetype_gispy/data/Geoindex_81670/mamh_usage_predo_2022_s_poly.shp'
+# Reading Shared platform of geospatial data
+# and aerial photographs (GeoIndex) layer
+geoindex_layer = 'C:/Users/a_adli/PycharmProjects/' \
+                 'hydroquebec_archetype_gispy/data/' \
+                 'input_data/Geoindex_81670/mamh_usage_predo_2022_s_poly.shp'
 geoindex, geoindex_name = load_layer(geoindex_layer, 'GeoIndex')
-
 print(f'{geoindex_name} data count: {geoindex.featureCount()}')
 
+params_create_index_geo = {'INPUT': geoindex_layer, 'OUTPUT': 'Output'}
+processing.run("native:createspatialindex", params_create_index_geo)
+
+# Fixing the GeoIndex layer geometries
 print(f'Fixing {geoindex_name} geometries')
+fixed_geoindex = 'C:/Users/a_adli/PycharmProjects/' \
+                 'hydroquebec_archetype_gispy/data/' \
+                 'output_data/fixed_geoindex/fixed_geoindex.shp'
+params_fixing_geoindex = {'INPUT': geoindex,
+                          'METHOD': 0,
+                          'OUTPUT': fixed_geoindex}
+processing.run("native:fixgeometries", params_fixing_geoindex)
 
-fixed_geoindex_0 = 'C:/Users/a_adli/PycharmProjects/hydroquebec_archetype_gispy/data/new_tests/fixed_geoindex_9/fix_geo_9.shp'
+params_create_index_fixed_geo = {'INPUT': fixed_geoindex,
+                                 'OUTPUT': 'Output'}
+processing.run("native:createspatialindex", params_create_index_fixed_geo)
 
-params_fixing_geoindex = {'INPUT': geoindex, 'METHOD': 0, 'OUTPUT': fixed_geoindex_0}
-fixed_layer_geoindex = processing.run("native:fixgeometries", params_fixing_geoindex)['OUTPUT']
+fixed_geoindex_read, fixed_geoindex_name = \
+  load_layer(fixed_geoindex, 'Fixed NRCan')
+print(f'{fixed_geoindex_name} '
+      f'data count: {fixed_geoindex_read.featureCount()}')
 
-fixed_geoindex, fixed_geoindex_name = load_layer(fixed_layer_geoindex, 'Fixed NRCan')
-print(f'{fixed_geoindex_name} data count: {fixed_geoindex.featureCount()}')
-
-
-clipped_geoindex_layer = 'C:/Users/a_adli/PycharmProjects/hydroquebec_archetype_gispy/data/new_tests/clipped_geo_3/clipped_geo_3.shp'
-
-params = {'INPUT': fixed_geoindex, 'OVERLAY': clipping_montreal_boundary_layer, 'FILTER_EXPRESSION': '', 'FILTER_EXTENT': None, 'OUTPUT': clipped_geoindex_layer}
+# Removing unnecessary parts of the GeoIndex layer (outside of Montreal) using
+# Administrative boundaries of the Montreal agglomeration data
+clipped_geoindex_layer = 'C:/Users/a_adli/PycharmProjects/' \
+                         'hydroquebec_archetype_gispy/data/' \
+                         'output_data/clipped_fixed_geoindex/' \
+                         'clipped_fixed_geoindex.shp'
+params = {'INPUT': fixed_geoindex,
+          'OVERLAY': clipping_montreal_boundary_layer,
+          'FILTER_EXPRESSION': '',
+          'FILTER_EXTENT': None,
+          'OUTPUT': clipped_geoindex_layer}
 processing.run("native:clip", params)
-
 print(f'Clipping {fixed_geoindex_name} is completed.')
 
-clipped_geoindex, clipped_geoindex_name = load_layer(clipped_nrcan_layer, 'Clipped GeoIndex')
+clipped_geoindex, clipped_geoindex_name = \
+  load_layer(clipped_geoindex_layer, 'Clipped GeoIndex')
 print(f'{clipped_geoindex_name} data count: {clipped_geoindex.featureCount()}')
 
-params_create_index_geoindex = {'INPUT': clipped_geoindex, 'OUTPUT': 'Output'}
-processing.run("native:createspatialindex", params_create_index_geoindex)
-
+# Creating spatial index fo the GeoIndex layer
+params_create_index_clipped_geo = {'INPUT': clipped_geoindex_layer,
+                                   'OUTPUT': 'Output'}
+processing.run("native:createspatialindex", params_create_index_clipped_geo)
 print(f'Creating spatial index for {clipped_geoindex_name} is completed.')
 
-# Adding the Property Assessment (uniteevaluation) dataset
+# Reading the Property Assessment (AKA uniteevaluationfonciere) dataset
+property_assessment_layer = 'C:/Users/a_adli/PycharmProjects/' \
+                            'hydroquebec_archetype_gispy/data/' \
+                            'input_data/property_assessment/' \
+                            'uniteevaluationfonciere.shp'
 
-uniteevaluationfonciere_layer = 'C:/Users/a_adli/PycharmProjects/hydroquebec_archetype_gispy/data/uniteevaluationfonciere/uniteevaluationfonciere.shp'
+property_assessment_read, property_assessment_layer_name = \
+  load_layer(property_assessment_layer, 'Property Assesment')
+print(f'{property_assessment_layer_name} '
+      f'data count: {property_assessment_read.featureCount()}')
 
-loaded_uniteevaluationfonciere_layer, uniteevaluationfonciere_layer_name = load_layer(uniteevaluationfonciere_layer, 'uniteevaluationfonciere')
-print(f'{uniteevaluationfonciere_layer_name} data count: {loaded_uniteevaluationfonciere_layer.featureCount()}')
+# Creating spatial index for the GeoIndex layer
+params_create_index_property_assessment = {'INPUT': property_assessment_layer,
+                                           'OUTPUT': 'Output'}
+processing.run("native:createspatialindex",
+               params_create_index_property_assessment)
+print(f'Creating spatial index for '
+      f'{property_assessment_layer_name} is completed.')
 
-clipped_nrcan_layer = 'C:/Users/a_adli/PycharmProjects/hydroquebec_archetype_gispy/data/tests/clipped_nrcan_10/clipped_nrcan_10.shp'
-clipped_uniteevaluation_layer = 'C:/Users/a_adli/PycharmProjects/hydroquebec_archetype_gispy/data/tests/clipped_uniteevaluation/clipped_uniteevaluation.shp'
+# For carrying out pairwise clip (clip in QGIS) on
+# property assessment data with Clipped and
+# fixed NRCan data as the overlay, we need to split the latter.
+# This is being done to improve the performance of the clipping task.
+# Otherwise, the task crashes the program.
+# One can decide on the number of this divisions by assigning
+# a number to num_layers. I assume, a number way bigger than
+# 20 can help the performance even more significantly (Will be tested).
 
-params_clip_uniteevaluation = {'INPUT': uniteevaluationfonciere_layer, 'OVERLAY': clipped_nrcan_layer, 'FILTER_EXPRESSION': '', 'FILTER_EXTENT': None, 'OUTPUT': clipped_uniteevaluation_layer}
-processing.run("native:clip", params_clip_uniteevaluation)
+output_layers_folder = 'C:/Users/a_adli/PycharmProjects/' \
+                       'hydroquebec_archetype_gispy/data/' \
+                       'output_data/splitted_clipped_nrcan'
+num_layers = 120
+create_folders(output_layers_folder, num_layers)
 
-print(f'Clipping the {uniteevaluationfonciere_layer_name} is completed.')
+splitting_intervals = clipped_nrcan_length // num_layers
 
-params_create_index_uniteevaluationfonciere = {'INPUT': clipped_uniteevaluation_layer, 'OUTPUT': 'Output'}
-indexed_layer = processing.run("native:createspatialindex", params_create_index_uniteevaluationfonciere)
+for each in range(num_layers):
+  output_layer_path = f'C:/Users/a_adli/PycharmProjects/' \
+                      f'hydroquebec_archetype_gispy/data/' \
+                      f'output_data/splitted_clipped_nrcan/' \
+                      f'layer_{each}/layer_{each}.shp'
+  splitting_params = {'INPUT': clipped_nrcan_layer,
+                      'EXPRESSION': f'$id >= {each * splitting_intervals} '
+                                    f'AND $id < '
+                                    f'{(each + 1) * splitting_intervals}\r\n',
+                      'OUTPUT': output_layer_path}
 
-print(f'Creating spatial index for {uniteevaluationfonciere_layer_name} is completed.')
+  processing.run("native:extractbyexpression", splitting_params)
 
-qgs.exitQgis()
+# Creating spatial index for each new layer
+  params_create_index_layers = {'INPUT': output_layer_path,
+                                'OUTPUT': 'Output'}
+  processing.run("native:createspatialindex", params_create_index_layers)
+
+# Putting the remaining features in one last layer
+# (So, overall we're going to have [num_layers + 1] layers).
+remaining_features = num_layers
+os.makedirs(f'C:/Users/a_adli/PycharmProjects/'
+            f'hydroquebec_archetype_gispy/data/'
+            f'output_data/splitted_clipped_nrcan/layer_{remaining_features}')
+output_layer_path = f'C:/Users/a_adli/PycharmProjects/' \
+                    f'hydroquebec_archetype_gispy/data/' \
+                    f'output_data/splitted_clipped_nrcan/' \
+                    f'layer_{remaining_features}/' \
+                    f'layer_{remaining_features}.shp'
+last_splitting_params = {'INPUT': clipped_nrcan_layer,
+                         'EXPRESSION':
+                         f'$id >= {num_layers * splitting_intervals}\r\n',
+                         'OUTPUT': output_layer_path}
+
+processing.run("native:extractbyexpression", last_splitting_params)
+
+# Create spatial index for the last partition (layer)
+params_create_index_last = {'INPUT': output_layer_path,
+                            'OUTPUT': 'Output'}
+processing.run("native:createspatialindex", params_create_index_last)
+
+# Now we clip the property assessment with each of the partitions (layers)
+# from the last process. The process outputs the pairwise clipped version of
+# the Property Assessment data in num_layers (the number is 21
+# in the first run) number of individual layers
+
+# First we make the folders for each clipping process
+clipped_layers_folder = 'C:/Users/a_adli/PycharmProjects/' \
+                        'hydroquebec_archetype_gispy/data/' \
+                        'output_data/clipped_property_assessment_partitions'
+num_clipped_layers = num_layers + 1
+create_folders(clipped_layers_folder, num_clipped_layers)
+
+# Through a for loop, we clip the
+# Property Assessment layer by all the individual overlays
+for each in range(num_clipped_layers):
+  clipping_layer = f'C:/Users/a_adli/PycharmProjects/' \
+                   f'hydroquebec_archetype_gispy/data/' \
+                   f'output_data/splitted_clipped_nrcan/' \
+                   f'layer_{each}/layer_{each}.shp'
+  clipped_layer = f'C:/Users/a_adli/PycharmProjects/' \
+                  f'hydroquebec_archetype_gispy/data/' \
+                  f'output_data/clipped_property_assessment_partitions/' \
+                  f'layer_{each}/layer_{each}.shp'
+
+  pairwise_clipping_params = {'INPUT': property_assessment_layer,
+                              'OVERLAY': clipping_layer,
+                              'FILTER_EXPRESSION': '',
+                              'FILTER_EXTENT': None,
+                              'OUTPUT': clipped_layer}
+  processing.run("native:clip", pairwise_clipping_params)
+
+  # Creating spatial index for each new layer
+  params_create_index_clips = {'INPUT': clipped_layer,
+                               'OUTPUT': 'Output'}
+  processing.run("native:createspatialindex", params_create_index_clips)
+
+print("Pairwise Clipping is completed.")
+
+# Finally, we merge all the clipped layers to make the
+# pairwise clipped Property Assessment layer
+
+# First we extract all the .shp files (clipped layers) from each folder
+pairwise_clipped_layers_path = find_shp_files(clipped_layers_folder)
+
+pairwise_clipped_property_assessment_layer = \
+  'C:/Users/a_adli/PycharmProjects/' \
+  'hydroquebec_archetype_gispy/data/' \
+  'output_data/pairwise_clipped_property_assessment/pairwise_clipped.shp'
+
+merging_params = {'LAYERS': pairwise_clipped_layers_path,
+                  'CRS': None,
+                  'OUTPUT': pairwise_clipped_property_assessment_layer}
+
+processing.run("native:mergevectorlayers", merging_params)
+
+# Create spatial index for the pairwise clipped Property Assessment
+params_create_index_pairwised = \
+  {'INPUT': pairwise_clipped_property_assessment_layer, 'OUTPUT': 'Output'}
+processing.run("native:createspatialindex", params_create_index_pairwised)
+
+pairwise_clipped_property_assessment, \
+ pairwise_clipped_property_assessment_name = \
+ load_layer(pairwise_clipped_property_assessment_layer,
+            'Pairwise Clipped Property Assessment Layer')
+print(f'{pairwise_clipped_property_assessment_name} '
+      f'data count: {pairwise_clipped_property_assessment.featureCount()}')
+
+# -----
+# Delete path and layer fields (This is going to be added after testing the whole program).
+# -----
+
+# In QGIS spatial join is named Join Attributes by Location
+# Spatial join Pairwise Clipped Property Assessment with Clipped NRCan
+property_assessment_nrcan_layer = \
+  'C:/Users/a_adli/PycharmProjects/' \
+  'hydroquebec_archetype_gispy/data/' \
+  'output_data/property_join_nrcan/property_join_nrcan.shp'
+
+spatial_join_nrcan_params = \
+  {'INPUT': pairwise_clipped_property_assessment_layer,
+   'PREDICATE': [0],
+   'JOIN': clipped_nrcan_layer,
+   'JOIN_FIELDS':  [],
+   'METHOD': 0,
+   'DISCARD_NONMATCHING': False,
+   'PREFIX': '',
+   'OUTPUT': property_assessment_nrcan_layer}
+
+feedback = QgsProcessingFeedback()
+processing.run('native:joinattributesbylocation',
+               spatial_join_nrcan_params, feedback=feedback)
+
+params_create_index_joined_nrcan = \
+  {'INPUT': property_assessment_nrcan_layer,
+   'OUTPUT': 'Output'}
+processing.run("native:createspatialindex", params_create_index_joined_nrcan)
+
+property_assessment_nrcan, property_assessment_nrcan_name = \
+  load_layer(property_assessment_nrcan_layer,
+             'Property Assessment Layer Joined with NRCan')
+print(f'{property_assessment_nrcan_name} '
+      f'data count: {property_assessment_nrcan.featureCount()}')
+
+# Spatial join Pairwise Clipped (joined with NRCan) Property Assessment layer with Clipped GeoIndex layer
+property_assessment_nrcan_geo_layer = \
+  'C:/Users/a_adli/PycharmProjects/' \
+  'hydroquebec_archetype_gispy/data/' \
+  'output_data/propertynrcan_join_geo/propertynrcan_join_geo.shp'
+
+
+spatial_join_nrcan_params = {'INPUT': property_assessment_nrcan_layer,
+                             'PREDICATE': [0],
+                             'JOIN': clipped_geoindex_layer,
+                             'JOIN_FIELDS':  [],
+                             'METHOD': 0,
+                             'DISCARD_NONMATCHING': False,
+                             'PREFIX': '',
+                             'OUTPUT': property_assessment_nrcan_geo_layer}
+
+feedback = QgsProcessingFeedback()
+processing.run('native:joinattributesbylocation',
+               spatial_join_nrcan_params, feedback=feedback)
+
+create_index_params = \
+  {'INPUT': property_assessment_nrcan_geo_layer,
+   'OUTPUT': 'Output'}
+processing.run("native:createspatialindex", create_index_params)
+
+property_assessment_nrcan_geo, property_assessment_nrcan_geo_name = \
+  load_layer(property_assessment_nrcan_geo_layer,
+             'Resulting Layer Joined with GeoIndex layer')
+print(f'{property_assessment_nrcan_geo_name} '
+      f'data count: {property_assessment_nrcan_geo.featureCount()}')
+
+# ---
+# There are four steps that will be added later after testing the program:
+# - Aligning GeoIndex layer/features with Property Assessment
+# (We should, firstly, make sure about its benefit)
+# - Count overlapping features
+# - Summarize within
+# - Adding the summarize-within field with a spatial join
+# ---
+
+# In QGIS Delete Identical is named Delete Duplicates
+property_nrcan_geo_deleted_duplicates_layer = \
+  'C:/Users/a_adli/PycharmProjects/' \
+  'hydroquebec_archetype_gispy/data/' \
+  'output_data/deleted_duplicates_property_and_all/delete_duplicates.shp'
+delete_duplicates_params = {'INPUT': property_assessment_nrcan_geo_layer,
+                            'OUTPUT': property_nrcan_geo_deleted_duplicates_layer}
+processing.run("native:deleteduplicategeometries", delete_duplicates_params)
+
+create_index_params = {'INPUT': property_nrcan_geo_deleted_duplicates_layer,
+                       'OUTPUT': 'Output'}
+processing.run("native:createspatialindex", create_index_params)
+
+property_nrcan_geo_deleted_duplicates, \
+ property_nrcan_geo_deleted_duplicates_name = \
+ load_layer(property_nrcan_geo_deleted_duplicates_layer,
+            'Resulting layer with deleted duplicates')
+print(f'{property_nrcan_geo_deleted_duplicates_name} '
+      f'data count: {property_nrcan_geo_deleted_duplicates.featureCount()}')
+
+singled_parts_layer = \
+  'C:/Users/a_adli/PycharmProjects/' \
+  'hydroquebec_archetype_gispy/data/' \
+  'output_data/singled_parts/singled_parts.shp'
+
+singled_parts_params = {'INPUT': property_nrcan_geo_deleted_duplicates_layer,
+                        'OUTPUT': singled_parts_layer}
+processing.run("native:multiparttosingleparts", params)
+
+create_index_params = {'INPUT': singled_parts_layer,
+                       'OUTPUT': 'Output'}
+processing.run("native:createspatialindex", create_index_params)
+
+singled_parts, singled_parts_layer_name = \
+ load_layer(singled_parts_layer,
+            'Property Assessment Layer Joined with NRCan')
+print(f'{singled_parts_layer_name} data count: {singled_parts.featureCount()}')
+
+qgs.exitQgis()

services_scripts/adding_area_field.py

@@ -0,0 +1,26 @@
+from qgis.core import QgsProject, QgsExpression, QgsExpressionContext, edit, QgsField, QgsExpressionContextUtils
+
+# Import necessary modules
+from qgis.core import QgsProject, QgsVectorLayer, QgsField
+
+# Define the layer name
+layer_name = "your_layer_name"
+
+# Get the layer by name
+layer = QgsProject.instance().mapLayersByName(layer_name)[0]
+
+# Start editing
+layer.startEditing()
+
+# Add a new field to store the area
+field_name = 'Area'
+field_index = layer.fields().indexFromName(field_name)
+if field_index == -1:
+    layer.dataProvider().addAttributes([QgsField(field_name, QgsField.QVariant.Double)])
+    layer.updateFields()
+
+# Commit changes
+layer.commitChanges()
+
+
+

services_scripts/basic_functions.py

@@ -0,0 +1,52 @@
+import os
+import glob
+import time
+
+
+def create_folders(directory, num_folders):
+  """
+  Create a specified number of folders in the given directory.
+
+  Args:
+  - directory (str): The directory where folders will be created.
+  - num_folders (int): The number of folders to create.
+  """
+  # Check if the directory exists, if not, create it
+  if not os.path.exists(directory):
+    os.makedirs(directory)
+
+  # Create folders
+  for i in range(num_folders):
+    folder_name = f"layer_{i}"
+    folder_path = os.path.join(directory, folder_name)
+    os.makedirs(folder_path)
+    print(f"Created folder: {folder_path}")
+
+
+# In a folder (including several sub folders) this function extract
+# .shp files from each subfolder
+def find_shp_files(root_folder):
+  shp_files = []
+  # Sort folders alphabetically
+  for foldername, _, _ in sorted(os.walk(root_folder)):
+    for filename in sorted(glob.glob(os.path.join(foldername, '*.shp'))):
+      new_file_name = filename.replace('\\', r'/')
+      shp_files.append(new_file_name)
+  return shp_files
+
+
+# Creates a progress bar for a for loop
+def progress_bar(iterable, length=20):
+    total = len(iterable)
+    progress = 0
+    start_time = time.time()
+    for item in iterable:
+        yield item
+        progress += 1
+        percent = progress / total
+        filled_length = int(length * percent)
+        bar = '=' * filled_length + '-' * (length - filled_length)
+        elapsed_time = time.time() - start_time
+        print(f'\r[{bar}] {progress}/{total} ({percent:.0%}) - '
+              f'Elapsed time: {elapsed_time:.2f}s', end='', flush=True)
+    print()

services_scripts/clip_layer.py

@@ -1,4 +1,4 @@
-from qgis.core import QgsApplication, QgsVectorLayer, QgsProcessingFeedback
+from qgis.core import QgsApplication
 from qgis.analysis import QgsNativeAlgorithms
 import processing
 
@@ -12,13 +12,26 @@ qgs.initQgis()
 # Add native algorithms provider
 QgsApplication.processingRegistry().addProvider(QgsNativeAlgorithms())
 
-input_layer = 'C:/Users/a_adli/PycharmProjects/hydroquebec_archetype_gispy/data/tests/python_fixed_05/py_fixes_05.shp'
-clipping_layer = 'C:/Users/a_adli/PycharmProjects/hydroquebec_archetype_gispy/data/montreal_boundary|layername=Montreal_boundary'
-clipped_layer = 'C:/Users/a_adli/PycharmProjects/hydroquebec_archetype_gispy/data/tests/clipped_nrcan_10/clipped_nrcan_10.shp'
+input_layer = \
+  'C:/Users/a_adli/PycharmProjects/' \
+  'hydroquebec_archetype_gispy/data/' \
+  'output_data/python_fixed_05/py_fixes_05.shp'
+clipping_layer = \
+  'C:/Users/a_adli/PycharmProjects/' \
+  'hydroquebec_archetype_gispy/data/' \
+  'input_data/montreal_boundary|layername=Montreal_boundary'
+clipped_layer = \
+  'C:/Users/a_adli/PycharmProjects/' \
+  'hydroquebec_archetype_gispy/data/' \
+  'output_data/clipped_nrcan_10/clipped_nrcan_10.shp'
 
-params = {'INPUT': input_layer, 'OVERLAY': clipping_layer, 'FILTER_EXPRESSION': '', 'FILTER_EXTENT': None, 'OUTPUT': clipped_layer}
+params = {'INPUT': input_layer,
+          'OVERLAY': clipping_layer,
+          'FILTER_EXPRESSION': '',
+          'FILTER_EXTENT': None,
+          'OUTPUT': clipped_layer}
 processing.run("native:clip", params)
 
 print("Clipping is completed.")
 
-qgs.exitQgis()
+qgs.exitQgis()

services_scripts/clip_with_several_overlays.py

@@ -0,0 +1,51 @@
+from qgis.core import QgsApplication
+from qgis.analysis import QgsNativeAlgorithms
+import processing
+from services_scripts.basic_functions import create_folders
+
+clipped_layers_folder = \
+  'C:/Users/a_adli/PycharmProjects/' \
+  'hydroquebec_archetype_gispy/' \
+  'data/aEndeavor/clip_all'
+num_clipped_layers = 21
+create_folders(clipped_layers_folder, num_clipped_layers)
+
+# Set the path to QGIS installation
+QgsApplication.setPrefixPath("C:/Program Files/QGIS 3.34.1/apps/qgis", True)
+
+# Initialize QGIS application
+qgs = QgsApplication([], False)
+qgs.initQgis()
+
+# Add native algorithms provider
+QgsApplication.processingRegistry().addProvider(QgsNativeAlgorithms())
+
+input_layer = \
+  'C:/Users/a_adli/PycharmProjects/' \
+  'hydroquebec_archetype_gispy/data/' \
+  'input_data/uniteevaluationfonciere/uniteevaluationfonciere.shp'
+
+for each in range(num_clipped_layers):
+  clipping_layer = \
+    f'C:/Users/a_adli/PycharmProjects/' \
+    f'hydroquebec_archetype_gispy/data/' \
+    f'output_data/divided_all/layer_{each}/layer_{each}.shp'
+  clipped_layer = \
+      f'C:/Users/a_adli/PycharmProjects/' \
+      f'hydroquebec_archetype_gispy/data/' \
+      f'output_data/clip_all/layer_{each}/layer_{each}.shp'
+
+  params = {'INPUT': input_layer,
+            'OVERLAY': clipping_layer,
+            'FILTER_EXPRESSION': '',
+            'FILTER_EXTENT': None,
+            'OUTPUT': clipped_layer}
+  processing.run("native:clip", params)
+
+  # Creating spatian index for each new layer
+  create_spatial_indext_params = {'INPUT': clipped_layer, 'OUTPUT': 'Output'}
+  processing.run("native:createspatialindex", create_spatial_indext_params)
+
+print("Clipping is completed.")
+
+qgs.exitQgis()

services_scripts/count_records.py

@@ -1,4 +1,4 @@
-from qgis.core import QgsApplication, QgsVectorLayer, QgsProcessingFeedback
+from qgis.core import QgsApplication
 from qgis.analysis import QgsNativeAlgorithms
 from load_layer import load_layer
 
@@ -13,9 +13,12 @@ qgs.initQgis()
 # Add native algorithms provider
 QgsApplication.processingRegistry().addProvider(QgsNativeAlgorithms())
 
-layer = 'C:/Users/a_adli/PycharmProjects/hydroquebec_archetype_gispy/data/tests/clipped_nrcan_10/clipped_nrcan_10.shp'
+layer = \
+  'C:/Users/a_adli/PycharmProjects/' \
+  'hydroquebec_archetype_gispy/data/' \
+  'output_data/clipped_nrcan_10/clipped_nrcan_10.shp'
 layer, layer_name = load_layer(layer, 'NRCan')
 
 print(f'{layer_name} data count: {layer.featureCount()}')
 
-qgs.exitQgis()
+qgs.exitQgis()

services_scripts/create_spatial_index.py

@@ -1,4 +1,4 @@
-from qgis.core import QgsApplication, QgsVectorLayer, QgsProcessingFeedback
+from qgis.core import QgsApplication
 from qgis.analysis import QgsNativeAlgorithms
 import processing
 
@@ -11,12 +11,16 @@ qgs.initQgis()
 
 # Add native algorithms provider
 QgsApplication.processingRegistry().addProvider(QgsNativeAlgorithms())
-layer = 'C:/Users/a_adli/PycharmProjects/hydroquebec_archetype_gispy/data/tests/fixed_geoindex_3/fix_geo.shp'
+layer = \
+  'C:/Users/a_adli/PycharmProjects/' \
+  'hydroquebec_archetype_gispy/data/' \
+  'output_data/merged_all_delete_columns/merge_divisions.shp'
 
-params = {'INPUT': layer, 'OUTPUT': 'Output'}
+params = {'INPUT': layer,
+          'OUTPUT': 'Output'}
 processing.run("native:createspatialindex", params)
 
 
 print("Creating spatial index is completed.")
 
-qgs.exitQgis()
+qgs.exitQgis()

services_scripts/delete_duplicates.py

@@ -0,0 +1,28 @@
+from qgis.core import QgsApplication
+from qgis.analysis import QgsNativeAlgorithms
+import processing
+
+# In QGIS Delete Identical is named Delete Duplicates
+# Set the path to QGIS installation
+QgsApplication.setPrefixPath("C:/Program Files/QGIS 3.34.1/apps/qgis", True)
+
+# Initialize QGIS application
+qgs = QgsApplication([], False)
+qgs.initQgis()
+
+# Add native algorithms provider
+QgsApplication.processingRegistry().addProvider(QgsNativeAlgorithms())
+input_layer = \
+  'C:/Users/a_adli/PycharmProjects/' \
+  'hydroquebec_archetype_gispy/data/' \
+  'output_data/qgis_spatial_pariwiseunitNrcan_geo|layername=spatial_with_geo'
+output_layer = \
+  'C:/Users/a_adli/PycharmProjects/' \
+  'hydroquebec_archetype_gispy/data/' \
+  'output_data/empty/delete_duplicates_02.shp'
+
+params = {'INPUT': input_layer,
+          'OUTPUT': output_layer}
+processing.run("native:deleteduplicategeometries", params)
+
+qgs.exitQgis()

services_scripts/delete_field.py

@@ -0,0 +1,44 @@
+from qgis.core import QgsApplication, QgsVectorLayer, QgsProject, edit
+from qgis.analysis import QgsNativeAlgorithms
+
+
+def load_layer(path, layer_name):
+  the_layer = QgsVectorLayer(path, layer_name, "ogr")
+  if not the_layer.isValid():
+    print(f'{layer_name} failed to load!')
+  else:
+    QgsProject.instance().addMapLayer(the_layer)
+  return the_layer, layer_name
+
+
+# Set the path to QGIS installation
+QgsApplication.setPrefixPath("C:/Program Files/QGIS 3.34.1/apps/qgis", True)
+
+# Initialize QGIS application
+qgs = QgsApplication([], False)
+qgs.initQgis()
+
+# Add native algorithms provider
+QgsApplication.processingRegistry().addProvider(QgsNativeAlgorithms())
+
+layer_path = \
+  'C:/Users/a_adli/PycharmProjects/' \
+  'hydroquebec_archetype_gispy/data/' \
+  'output_data/merged_all_delete_columns/merge_divisions.shp'
+output_layer = load_layer(layer_path, 'Merged Pairwise Clipped')[0]
+
+column_name = 'path'  # and 'layer'
+
+# Start editing
+with edit(output_layer):
+  # Get the index of the column to delete
+  idx = output_layer.fields().indexFromName(column_name)
+
+  # Delete the field
+  output_layer.deleteAttribute(idx)
+
+  # Update layer fields
+  output_layer.updateFields()
+
+
+qgs.exitQgis()

services_scripts/fix_geometries.py

@@ -1,4 +1,4 @@
-from qgis.core import QgsApplication, QgsVectorLayer, QgsProcessingFeedback
+from qgis.core import QgsApplication
 from qgis.analysis import QgsNativeAlgorithms
 import processing
 
@@ -11,8 +11,14 @@ qgs.initQgis()
 
 # Add native algorithms provider
 QgsApplication.processingRegistry().addProvider(QgsNativeAlgorithms())
-layer = 'C:/Users/a_adli/PycharmProjects/hydroquebec_archetype_gispy/data/Geoindex_81670/mamh_usage_predo_2022_s_poly.shp'
-fixed_layer = 'C:/Users/a_adli/PycharmProjects/hydroquebec_archetype_gispy/data/tests/fixed_geoindex_7/py_fixes_07.shp'
+layer = \
+  'C:/Users/a_adli/PycharmProjects/' \
+  'hydroquebec_archetype_gispy/data/' \
+  'input_data/Geoindex_81670/mamh_usage_predo_2022_s_poly.shp'
+fixed_layer = \
+  'C:/Users/a_adli/PycharmProjects/' \
+  'hydroquebec_archetype_gispy/data/' \
+  'output_data/fixed_geoindex_7/py_fixes_07.shp'
 
 params = {'INPUT': layer, 'METHOD': 0, 'OUTPUT': fixed_layer}
 fix_layer = processing.run("native:fixgeometries", params)['OUTPUT']

services_scripts/intersection.py

@@ -0,0 +1,42 @@
+# This services is not used in the main workflow
+from qgis.core import QgsApplication
+from qgis.analysis import QgsNativeAlgorithms
+import processing
+
+# Set the path to QGIS installation
+QgsApplication.setPrefixPath('C:/Program Files/QGIS 3.34.1/apps/qgis', True)
+
+# Initialize QGIS application
+qgs = QgsApplication([], False)
+qgs.initQgis()
+
+
+# Add native algorithms provider
+QgsApplication.processingRegistry().addProvider(QgsNativeAlgorithms())
+
+layer_1_path = \
+  'C:/Users/a_adli/PycharmProjects/' \
+  'hydroquebec_archetype_gispy/data/' \
+  'output_data/clipped_nrcan_1/clipped_nrcan_1.shp'
+layer_2_path = \
+  'C:/Users/a_adli/PycharmProjects/' \
+  'hydroquebec_archetype_gispy/data/' \
+  'input_data/property_assessment/uniteevaluationfonciere.shp'
+output_layer_path = \
+  'C:/Users/a_adli/PycharmProjects/' \
+  'hydroquebec_archetype_gispy/data/' \
+  'output_data/intersection_nrcan_property/intersection_nrcan_property_02.shp'
+
+
+params = {'INPUT': layer_1_path,
+          'OVERLAY': layer_2_path,
+          'INPUT_FIELDS': [],
+          'OVERLAY_FIELDS': [],
+          'OVERLAY_FIELDS_PREFIX': '',
+          'OUTPUT': output_layer_path,
+          'GRID_SIZE': None}
+
+processing.run("native:intersection", params)
+
+# Exit QGIS application
+qgs.exitQgis()

services_scripts/join_attributes_by_location.py

@@ -0,0 +1,45 @@
+# Join attributes by location in QGIS is the same as Spatial Join in ArcGIS
+
+from qgis.core import QgsApplication, QgsProcessingFeedback
+from qgis.analysis import QgsNativeAlgorithms
+import processing
+
+# Set the path to QGIS installation
+QgsApplication.setPrefixPath('C:/Program Files/QGIS 3.34.1/apps/qgis', True)
+
+# Initialize QGIS application
+qgs = QgsApplication([], False)
+qgs.initQgis()
+
+
+# Add native algorithms provider
+QgsApplication.processingRegistry().addProvider(QgsNativeAlgorithms())
+
+layer_1_path = \
+  'C:/Users/a_adli/PycharmProjects/' \
+  'hydroquebec_archetype_gispy/data/' \
+  'output_data/clipped_nrcan_1/clipped_nrcan_1.shp'
+layer_2_path = \
+  'C:/Users/a_adli/PycharmProjects/' \
+  'hydroquebec_archetype_gispy/data/' \
+  'output_data/intersection_nrcan_property/intersection_nrcan_property_02.shp'
+output_layer_path = \
+  'C:/Users/a_adli/PycharmProjects/' \
+  'hydroquebec_archetype_gispy/data/' \
+  'output_data/join_nrcan_intersected_nrcan/join_nrcan_intersected_nrcan.shp'
+
+
+params = {'INPUT': layer_1_path,
+          'PREDICATE': [0],
+          'JOIN': layer_2_path,
+          'JOIN_FIELDS':  [],
+          'METHOD': 0,
+          'DISCARD_NONMATCHING': False,
+          'PREFIX': '',
+          'OUTPUT': output_layer_path}
+
+feedback = QgsProcessingFeedback()
+processing.run('native:joinattributesbylocation', params, feedback=feedback)
+
+# Exit QGIS application
+qgs.exitQgis()

services_scripts/load_layer.py

@@ -1,7 +1,7 @@
 from qgis.core import *
 
 # Supply the path to the qgis install location
-QgsApplication.setPrefixPath(prefixPath='C:/Program Files/QGIS 3.34.1/apps/qgis', useDefaultPaths=True)
+QgsApplication.setPrefixPath('C:/Program Files/QGIS 3.34.1/apps/qgis', True)
 
 qgs = QgsApplication([], False)
 
@@ -15,4 +15,4 @@ def load_layer(path, layer_name):
         print(f'{layer_name} failed to load!')
     else:
         QgsProject.instance().addMapLayer(the_layer)
-    return the_layer, layer_name
+    return the_layer, layer_name

services_scripts/merge_layers.py

@@ -0,0 +1,32 @@
+from qgis.core import QgsApplication
+from qgis.analysis import QgsNativeAlgorithms
+import processing
+from services_scripts.basic_functions import find_shp_files
+
+
+root_folder = \
+  'C:/Users/a_adli/PycharmProjects/' \
+  'hydroquebec_archetype_gispy/data/' \
+  'output_data/clip_all'
+layers_path = find_shp_files(root_folder)
+
+# Set the path to QGIS installation
+QgsApplication.setPrefixPath("C:/Program Files/QGIS 3.34.1/apps/qgis", True)
+
+# Initialize QGIS application
+qgs = QgsApplication([], False)
+qgs.initQgis()
+
+# Add native algorithms provider
+QgsApplication.processingRegistry().addProvider(QgsNativeAlgorithms())
+
+merged_layer = \
+  'C:/Users/a_adli/PycharmProjects/' \
+  'hydroquebec_archetype_gispy/data/' \
+  'output_data/merged_all/merge_divisions.shp'
+
+params = {'LAYERS': layers_path,
+          'CRS': None,
+          'OUTPUT': merged_layer}
+
+processing.run("native:mergevectorlayers", params)

services_scripts/multipart_to_singleparts.py

@@ -0,0 +1,29 @@
+from qgis.core import QgsApplication
+from qgis.analysis import QgsNativeAlgorithms
+import processing
+
+# Set the path to QGIS installation
+QgsApplication.setPrefixPath("C:/Program Files/QGIS 3.34.1/apps/qgis", True)
+
+# Initialize QGIS application
+qgs = QgsApplication([], False)
+qgs.initQgis()
+
+# Add native algorithms provider
+QgsApplication.processingRegistry().addProvider(QgsNativeAlgorithms())
+
+
+input_layer = \
+  'C:/Users/a_adli/PycharmProjects/' \
+  'hydroquebec_archetype_gispy/data/' \
+  'output_data/delete_dups_4_partitioned|layername=delete_dups__partitioned'
+singled_parts_layer = \
+  'C:/Users/a_adli/PycharmProjects/' \
+  'hydroquebec_archetype_gispy/data/' \
+  'output_data/multi_to_single/multi_to_single2.shp'
+
+params = {'INPUT': input_layer,
+          'OUTPUT': singled_parts_layer}
+processing.run("native:multiparttosingleparts", params)
+
+qgs.exitQgis()

services_scripts/split_layer_by_expression.py

@@ -0,0 +1,81 @@
+from qgis.core import QgsApplication, QgsVectorLayer, QgsProject
+from qgis.analysis import QgsNativeAlgorithms
+import processing
+import os
+from basic_functions import create_folders
+
+
+# This function loads a layer. It's used to count the data records
+def load_layer(path, layer_name):
+  the_layer = QgsVectorLayer(path, layer_name, "ogr")
+  if not the_layer.isValid():
+    print(f'{layer_name} failed to load!')
+  else:
+    QgsProject.instance().addMapLayer(the_layer)
+  return the_layer, layer_name
+
+
+# Set the path to QGIS installation
+QgsApplication.setPrefixPath('C:/Program Files/QGIS 3.34.1/apps/qgis', True)
+
+# Initialize QGIS application
+qgs = QgsApplication([], False)
+qgs.initQgis()
+
+
+# Add native algorithms provider
+QgsApplication.processingRegistry().addProvider(QgsNativeAlgorithms())
+
+output_layers_folder = \
+  'C:/Users/a_adli/PycharmProjects/' \
+  'hydroquebec_archetype_gispy/data/' \
+  'output_data/divided_all'
+num_layers = 20
+create_folders(output_layers_folder, num_layers)
+
+layer_path = \
+  'C:/Users/a_adli/PycharmProjects/' \
+  'hydroquebec_archetype_gispy/data/' \
+  'output_data/02_clipped_nrcan_mtl_boundary/clipped_nrcan.shp'
+
+fixed_nrcan, fixed_nrcan_name = load_layer(layer_path, 'Clipped NRCan')
+layer_length = fixed_nrcan.featureCount()
+print(f'{fixed_nrcan_name} data count: {layer_length}')
+
+intervals = layer_length // num_layers
+
+for each in range(num_layers):
+  output_layer_path = \
+    f'C:/Users/a_adli/PycharmProjects/' \
+    f'hydroquebec_archetype_gispy/data/' \
+    f'output_data/divided_all/layer_{each}/layer_{each}.shp'
+  params = {'INPUT': layer_path,
+            'EXPRESSION': f'$id >= {each * intervals} '
+                          f'AND $id < {(each + 1) * intervals}\r\n',
+            'OUTPUT': output_layer_path}
+
+  processing.run("native:extractbyexpression", params)
+
+# Creating spatian index for each new layer
+  create_spatial_indext_params = {'INPUT': output_layer_path,
+                                  'OUTPUT': 'Output'}
+  processing.run("native:createspatialindex", create_spatial_indext_params)
+
+
+remaining_features = num_layers
+os.makedirs(f'C:/Users/a_adli/PycharmProjects/'
+            f'hydroquebec_archetype_gispy/data/'
+            f'output_data/divided_all/layer_{remaining_features}')
+output_layer_path = \
+  f'C:/Users/a_adli/PycharmProjects/' \
+  f'hydroquebec_archetype_gispy/data/' \
+  f'output_data/divided_all/' \
+  f'layer_{remaining_features}/layer_{remaining_features}.shp'
+params = {'INPUT': layer_path,
+          'EXPRESSION': f'$id >= {num_layers * intervals}\r\n',
+          'OUTPUT': output_layer_path}
+
+processing.run("native:extractbyexpression", params)
+
+# Exit QGIS application
+qgs.exitQgis()