diff --git a/services_scripts/split_layer_by_expression.py b/services_scripts/split_layer_by_expression.py
new file mode 100644
index 0000000..7a3291e
--- /dev/null
+++ b/services_scripts/split_layer_by_expression.py
@@ -0,0 +1,80 @@
+from qgis.core import QgsApplication, QgsVectorLayer, QgsProject, QgsProcessingFeedback
+from qgis.analysis import QgsNativeAlgorithms
+import processing
+import os
+
+
+# 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
+
+
+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}")
+
+
+# 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/aEndeavor/divided_all'
+num_layers = 3
+create_folders(output_layers_folder, num_layers)
+
+layer_path = 'C:/Users/a_adli/PycharmProjects/hydroquebec_archetype_gispy/data/aEndeavor/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/hydroquebec_archetype_gispy/data/aEndeavor/divided_all/layer_{each}/layer_{each}.shp'
+  params = {'INPUT': layer_path,
+            'EXPRESSION': f'$id >= {each * intervals} AND $id < {(each + 1) * intervals}\r\n',
+            'OUTPUT': output_layer_path}
+
+  processing.run("native:extractbyexpression", params)
+
+
+remaining_features = num_layers
+os.makedirs(f'C:/Users/a_adli/PycharmProjects/hydroquebec_archetype_gispy/data/aEndeavor/divided_all/layer_{remaining_features}')
+output_layer_path = f'C:/Users/a_adli/PycharmProjects/hydroquebec_archetype_gispy/data/aEndeavor/divided_all/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()