Repository: meteocat/pymica (language: python)

Issue:
Clusters - voronoi
Voronoi and Delaunay operations that were done in the browser with d3 to obtain shp files have to be implemented in Python in clusters.py
Clusters - voronoi
Voronoi and Delaunay operations that were done in the browser with d3 to obtain shp files have to be implemented in Python in clusters.py


Discussion / hints:

Fix patch:
diff --git a/environment_dev.yml b/environment_dev.yml
index f4596ec..b92efb6 100644
--- a/environment_dev.yml
+++ b/environment_dev.yml
@@ -5,7 +5,18 @@ channels:
 dependencies:
   - _libgcc_mutex=0.1=conda_forge
   - _openmp_mutex=4.5=2_gnu
+  - alabaster=0.7.12=pyhd3eb1b0_0
+  - anaconda-anon-usage=0.4.2=py311hfc0e8ea_2
+  - anaconda-client=1.12.1=py311h06a4308_0
+  - anyio=3.7.1=pyhd8ed1ab_0
+  - attrs=23.1.0=pyh71513ae_1
+  - babel=2.11.0=py311h06a4308_0
+  - backports=1.0=pyhd8ed1ab_3
+  - backports.functools_lru_cache=1.6.5=pyhd8ed1ab_0
+  - beautifulsoup4=4.12.2=py311h06a4308_0
   - blosc=1.21.5=h0f2a231_0
+  - boa=0.14.0=pyhd8ed1ab_4
+  - boltons=23.0.0=py311h06a4308_0
   - boost-cpp=1.78.0=h2c5509c_4
   - brotli-python=1.1.0=py311hb755f60_0
   - bzip2=1.0.8=h7f98852_4
@@ -13,13 +24,28 @@ dependencies:
   - ca-certificates=2023.7.22=hbcca054_0
   - cairo=1.16.0=h0c91306_1017
   - certifi=2023.7.22=pyhd8ed1ab_0
+  - cffi=1.15.1=py311h5eee18b_3
   - cfitsio=4.3.0=hbdc6101_0
+  - chardet=4.0.0=py311h06a4308_1003
   - charset-normalizer=3.2.0=pyhd8ed1ab_0
+  - click=8.0.4=py311h06a4308_0
+  - clyent=1.2.2=py311h06a4308_1
   - colorama=0.4.6=pyhd8ed1ab_0
+  - conda=23.7.4=py311h38be061_0
+  - conda-build=3.24.0=py311h38be061_1
+  - conda-index=0.3.0=py311h06a4308_0
+  - conda-package-handling=2.2.0=py311h06a4308_0
+  - conda-package-streaming=0.9.0=py311h06a4308_0
   - coverage=7.3.0=py311h459d7ec_0
+  - cryptography=41.0.3=py311hdda0065_0
   - cython=3.0.2=py311hb755f60_0
+  - dataclasses=0.8=pyhc8e2a94_3
+  - defusedxml=0.7.1=pyhd3eb1b0_0
+  - docutils=0.18.1=py311h06a4308_3
   - exceptiongroup=1.1.3=pyhd8ed1ab_0
   - expat=2.5.0=hcb278e6_1
+  - filelock=3.9.0=py311h06a4308_0
+  - fmt=9.1.0=h924138e_0
   - font-ttf-dejavu-sans-mono=2.37=hab24e00_0
   - font-ttf-inconsolata=3.000=h77eed37_0
   - font-ttf-source-code-pro=2.038=h77eed37_0
@@ -34,13 +60,23 @@ dependencies:
   - geotiff=1.7.1=h22adcc9_11
   - gettext=0.21.1=h27087fc_0
   - giflib=5.2.1=h0b41bf4_3
+  - glob2=0.7=py_0
   - hdf4=4.2.15=h501b40f_6
   - hdf5=1.14.2=nompi_h4f84152_100
   - icu=73.2=h59595ed_0
   - idna=3.4=pyhd8ed1ab_0
+  - imagesize=1.4.1=py311h06a4308_0
+  - importlib_resources=6.1.0=pyhd8ed1ab_0
   - iniconfig=2.0.0=pyhd8ed1ab_0
+  - jinja2=3.1.2=py311h06a4308_0
   - joblib=1.3.2=pyhd8ed1ab_0
   - json-c=0.17=h7ab15ed_0
+  - json5=0.9.14=pyhd8ed1ab_0
+  - jsonpatch=1.32=pyhd3eb1b0_0
+  - jsonpointer=2.1=pyhd3eb1b0_0
+  - jsonschema=4.19.1=pyhd8ed1ab_0
+  - jsonschema-specifications=2023.7.1=pyhd8ed1ab_0
+  - jupyter_core=5.3.0=py311h06a4308_0
   - kealib=1.5.1=hcd42e92_5
   - keyutils=1.6.1=h166bdaf_0
   - krb5=1.21.2=h659d440_0
@@ -67,6 +103,10 @@ dependencies:
   - libjpeg-turbo=2.1.5.1=h0b41bf4_0
   - libkml=1.3.0=h37653c0_1015
   - liblapack=3.9.0=18_linux64_openblas
+  - liblief=0.12.3=h6a678d5_0
+  - libllvm14=14.0.6=hcd5def8_4
+  - libmamba=1.4.9=h658169a_0
+  - libmambapy=1.4.9=py311h527f279_0
   - libnetcdf=4.9.2=nompi_h80fb2b6_112
   - libnghttp2=1.52.0=h61bc06f_0
   - libnsl=2.0.0=h7f98852_0
@@ -74,6 +114,7 @@ dependencies:
   - libpng=1.6.39=h753d276_0
   - libpq=15.4=hfc447b1_0
   - librttopo=1.1.0=hb58d41b_14
+  - libsolv=0.7.25=hfc55251_0
   - libspatialite=5.0.1=h15f6e67_28
   - libsqlite=3.43.0=h2797004_0
   - libssh2=1.11.0=h0841786_0
@@ -85,18 +126,30 @@ dependencies:
   - libxml2=2.11.5=h232c23b_1
   - libzip=1.10.1=h2629f0a_2
   - libzlib=1.2.13=hd590300_5
+  - llvmlite=0.40.1=py311ha6695c7_0
   - lz4-c=1.9.4=hcb278e6_0
   - lzo=2.10=h516909a_1000
+  - mamba=1.4.9=py311h3072747_0
+  - markdown-it-py=3.0.0=pyhd8ed1ab_0
+  - markupsafe=2.1.1=py311h5eee18b_0
+  - mdurl=0.1.0=pyhd8ed1ab_0
+  - more-itertools=8.12.0=pyhd3eb1b0_0
+  - nbformat=5.9.2=py311h06a4308_0
   - ncurses=6.4=hcb278e6_0
   - nspr=4.35=h27087fc_0
   - nss=3.92=h1d7d5a4_0
-  - numpy=1.25.2=py311h64a7726_0
+  - numba=0.57.1=py311h96b013e_0
+  - numpy=1.24.4=py311h64a7726_0
   - openjpeg=2.5.0=hfec8fc6_2
-  - openssl=3.1.2=hd590300_0
+  - openssl=3.1.3=hd590300_0
   - packaging=23.1=pyhd8ed1ab_0
+  - patch=2.7.6=h7b6447c_1001
+  - patchelf=0.17.2=h6a678d5_0
   - pcre2=10.40=hc3806b6_0
   - pip=23.2.1=pyhd8ed1ab_0
   - pixman=0.40.0=h36c2ea0_0
+  - pkginfo=1.9.6=py311h06a4308_0
+  - pkgutil-resolve-name=1.3.10=pyhd8ed1ab_1
   - platformdirs=3.10.0=pyhd8ed1ab_0
   - pluggy=1.3.0=pyhd8ed1ab_0
   - pooch=1.7.0=pyha770c72_3
@@ -104,30 +157,71 @@ dependencies:
   - poppler-data=0.4.12=hd8ed1ab_0
   - postgresql=15.4=h8972f4a_0
   - proj=9.2.1=ha643af7_0
+  - prompt-toolkit=3.0.39=pyha770c72_0
+  - prompt_toolkit=3.0.39=hd8ed1ab_0
+  - psutil=5.9.0=py311h5eee18b_0
   - pthread-stubs=0.4=h36c2ea0_1001
+  - py-lief=0.12.3=py311h6a678d5_0
+  - pybind11-abi=4=hd8ed1ab_3
+  - pycosat=0.6.4=py311h5eee18b_0
+  - pycparser=2.21=pyhd3eb1b0_0
+  - pygments=2.15.1=py311h06a4308_1
+  - pyopenssl=23.2.0=py311h06a4308_0
   - pyproj=3.6.0=py311ha169711_1
+  - pyshp=2.3.1=pyhd8ed1ab_0
   - pysocks=1.7.1=pyha2e5f31_6
   - pytest=7.4.1=pyhd8ed1ab_0
   - pytest-cov=4.1.0=pyhd8ed1ab_0
   - python=3.11.5=hab00c5b_0_cpython
+  - python-dateutil=2.8.2=pyhd3eb1b0_0
+  - python-fastjsonschema=2.16.2=py311h06a4308_0
+  - python-libarchive-c=2.9=pyhd3eb1b0_1
   - python_abi=3.11=3_cp311
+  - pytz=2023.3.post1=py311h06a4308_0
+  - pyyaml=6.0=py311h5eee18b_1
   - readline=8.2=h8228510_1
+  - referencing=0.30.2=pyhd8ed1ab_0
+  - reproc=14.2.4.post0=hd590300_0
+  - reproc-cpp=14.2.4.post0=h59595ed_0
   - requests=2.31.0=pyhd8ed1ab_0
+  - requests-toolbelt=1.0.0=py311h06a4308_0
+  - rich=13.6.0=pyhd8ed1ab_0
+  - ripgrep=13.0.0=h2f28480_2
+  - rpds-py=0.10.3=py311h46250e7_1
+  - ruamel.yaml=0.17.21=py311h5eee18b_0
   - scikit-learn=1.3.0=py311hc009520_0
   - scipy=1.11.2=py311h64a7726_0
   - setuptools=68.1.2=pyhd8ed1ab_0
+  - shapely=2.0.1=py311he06c224_3
+  - six=1.16.0=pyhd3eb1b0_1
   - snappy=1.1.10=h9fff704_0
+  - sniffio=1.3.0=pyhd8ed1ab_0
+  - snowballstemmer=2.2.0=pyhd3eb1b0_0
+  - soupsieve=2.5=py311h06a4308_0
+  - sphinx=5.0.2=py311h06a4308_0
+  - sphinxcontrib-applehelp=1.0.2=pyhd3eb1b0_0
+  - sphinxcontrib-devhelp=1.0.2=pyhd3eb1b0_0
+  - sphinxcontrib-htmlhelp=2.0.0=pyhd3eb1b0_0
+  - sphinxcontrib-jsmath=1.0.1=pyhd3eb1b0_0
+  - sphinxcontrib-qthelp=1.0.3=pyhd3eb1b0_0
+  - sphinxcontrib-serializinghtml=1.1.5=pyhd3eb1b0_0
   - sqlite=3.43.0=h2c6b66d_0
   - threadpoolctl=3.2.0=pyha21a80b_0
   - tiledb=2.16.3=hf0b6e87_1
   - tk=8.6.12=h27826a3_0
   - toml=0.10.2=pyhd8ed1ab_0
   - tomli=2.0.1=pyhd8ed1ab_0
+  - toolz=0.12.0=pyhd8ed1ab_0
+  - tqdm=4.65.0=py311h92b7b1e_0
+  - traitlets=5.7.1=py311h06a4308_0
+  - truststore=0.8.0=py311h06a4308_0
   - typing-extensions=4.7.1=hd8ed1ab_0
   - typing_extensions=4.7.1=pyha770c72_0
   - tzcode=2023c=h0b41bf4_0
   - tzdata=2023c=h71feb2d_0
   - urllib3=2.0.4=pyhd8ed1ab_0
+  - watchgod=0.8.2=pyhd8ed1ab_0
+  - wcwidth=0.2.8=pyhd8ed1ab_0
   - wheel=0.41.2=pyhd8ed1ab_0
   - xerces-c=3.2.4=hac6953d_3
   - xorg-kbproto=1.0.7=h7f98852_1002
@@ -142,6 +236,15 @@ dependencies:
   - xorg-xextproto=7.3.0=h0b41bf4_1003
   - xorg-xproto=7.0.31=h7f98852_1007
   - xz=5.2.6=h166bdaf_0
+  - yaml=0.2.5=h7b6447c_0
+  - yaml-cpp=0.7.0=h59595ed_3
+  - zipp=3.17.0=pyhd8ed1ab_0
   - zlib=1.2.13=hd590300_5
+  - zstandard=0.19.0=py311h5eee18b_0
   - zstd=1.5.5=hfc55251_0
+  - pip:
+    - accessible-pygments==0.0.4
+    - pydata-sphinx-theme==0.14.1
+    - sphinx-book-theme==1.0.1
+    - sphinxcontrib-autoprogram==0.1.8
 prefix: /home/ecm/anaconda3/envs/pymica-dev
diff --git a/pymica/utils/clusters.py b/pymica/utils/clusters.py
index ce7cfcc..1865059 100644
--- a/pymica/utils/clusters.py
+++ b/pymica/utils/clusters.py
@@ -1,147 +1,284 @@
 """Tool for creating clusters with the K-Means method.
 """
 import json
+from glob import glob
+from os import remove
 from os.path import dirname
 
 import numpy as np
+import pyproj
+import shapefile
 from osgeo import gdal, ogr, osr
 from scipy.ndimage import gaussian_filter
+from scipy.spatial import Voronoi
+from shapely.geometry import MultiPolygon, shape
 from sklearn.cluster import KMeans
 
-from pymica.utils.geotools import get_utm_epsg_from_lonlat, reproject_point
+from pymica.utils.geotools import reproject_point
 
 
-def create_clusters(locations: list | str, n_clusters: int, geojson_file: str) -> None:
-    """Group locations in a specified number of clusters using K-Means algorithm.
+def create_clusters(
+    locations: list | str,
+    n_clusters: int,
+    shp_file: str,
+    bounding_box: list,
+    epsg_code: int,
+) -> None:
+    """Group locations into a specified number of clusters and save them as an ESRI
+    Shapefile.
+
+    The Shapefile will have the extent of the specified `bounding_box`, and it's
+    crucial that the coordinates of the bounding box are in the same projection as the
+    locations.
 
     Args:
-        locations (list | str): List of dictionaries including location data or the
-                                file path with data as json format. Data for each
-                                location is a dictionary with 'id', 'alt', 'lon' and
-                                'lat' as keys.
-        n_clusters (int): Number of clusters to create.
-        geojson_file (str): File path where clusters are saved.
+        locations (list or str): List of dictionaries containing location data or the
+            file path with data in JSON format. Each dictionary should include keys
+            'id', 'alt', 'lon', and 'lat' for each location.
+        n_clusters (int): The number of clusters to group the locations into.
+        shp_file (str): The path to the output ESRI Shapefile.
+        bounding_box (list): A list representing the extent of the Shapefile in the
+            format [min_x, min_y, max_x, max_y]. The coordinates in this bounding box
+            must be in the same projection as the locations.
+        epsg_code (int): The EPSG code specifying the coordinate system for the
+            provided locations.
+
+    Returns:
+        None
     """
     if isinstance(locations, str):
         f_p = open(locations)
         locations = json.load(f_p)
 
-    positions = np.zeros([len(locations), 2])
-    utm_proj4 = get_utm_epsg_from_lonlat(locations[0]["lon"], locations[0]["lat"])
-
+    x_coords, y_coords = [], []
     positions_list = []
-    for i, loc_value in enumerate(locations):
-        positions[i][0] = loc_value["lon"]
-        positions[i][1] = loc_value["lat"]
-
-        utm_coordinates = reproject_point(
-            (loc_value["lon"], loc_value["lat"]), "EPSG:4326", utm_proj4
+    for loc_value in locations:
+        x_coord, y_coord = reproject_point(
+            (loc_value["lon"], loc_value["lat"]), 4326, epsg_code
         )
-        positions_list.append(utm_coordinates)
+        x_coords.append(x_coord)
+        y_coords.append(y_coord)
+
+        positions_list.append((x_coord, y_coord))
+
     positions_list = np.array(positions_list)
 
     kmeans = KMeans(n_clusters=n_clusters, random_state=0).fit(positions_list)
+    cluster_labels = kmeans.fit_predict(positions_list)
+
+    positions_list = np.append(
+        positions_list,
+        np.array(
+            [
+                [0.9 * min(x_coords), 0.9 * min(y_coords)],
+                [0.9 * min(x_coords), 1.1 * max(y_coords)],
+                [1.1 * max(x_coords), 0.9 * min(y_coords)],
+                [1.1 * max(x_coords), 1.1 * max(y_coords)],
+            ]
+        ),
+        axis=0,
+    )
 
-    out_geojson = {"type": "FeatureCollection", "features": []}
-    for i, labels_value in enumerate(kmeans.labels_):
-        locations[i]["cluster"] = labels_value
-        out_geojson["features"].append(
-            {
-                "type": "Feature",
-                "geometry": {
-                    "type": "Point",
-                    "coordinates": [locations[i]["lon"], locations[i]["lat"]],
-                },
-                "properties": {
-                    "cluster": int(labels_value),
-                    "id": locations[i]["id"],
-                    "alt": locations[i]["alt"],
-                },
-            }
-        )
+    # Calculate Voronoi diagram
+    vor = Voronoi(positions_list)
 
-    try:
-        with open(geojson_file, "w") as f_p:
-            json.dump(out_geojson, f_p)
-    except FileNotFoundError:
-        raise FileNotFoundError(dirname(geojson_file) + " directory does not exist.")
+    # Create a dictionary to store Voronoi cells for each cluster
+    cluster_cells = {cluster_id: [] for cluster_id in range(n_clusters)}
 
+    # Assign Voronoi cells to clusters
+    for station_id, region_id in enumerate(vor.point_region):
+        if region_id >= 0:
+            try:
+                cluster_id = cluster_labels[station_id]
+            except Exception:
+                cluster_id = 0
+            vertices = vor.regions[region_id]
 
-def create_reprojected_geometries(geojson_file: str, epsg: int) -> ogr.DataSource:
-    """Reprojection of an ogr file to specified projection.
-    Taken from: https://pcjericks.github.io/py-gdalogr-cookbook/projection.html
+            # Filter out invalid vertices (negative indices)
+            vertices = [v for v in vertices if v >= 0]
 
-    Args:
-        geojson_file (str): Path to a GeoJSON file with geometries.
-        epsg (int): EPSG code of the new projection
+            if vertices:
+                polygon = [vor.vertices[v] for v in vertices]
 
-    Raises:
-        ValueError: If EPSG code is wrong.
-        IOError: If `geojson_file` does not exist.
+                # Ensure the polygon is closed by adding the first point to the end
+                if np.array_equal(polygon[0], polygon[-1]):
+                    polygon.append(polygon[0])
 
-    Returns:
-        ogr.DataSource: Reprojected geometries.
-    """
-    out_proj = osr.SpatialReference()
-    out_proj.ImportFromEPSG(epsg)
-    if len(out_proj.ExportToPrettyWkt()) <= 1:
-        raise ValueError("Wrong EPSG code: {}".format(epsg))
+                cluster_cells[cluster_id].append(polygon)
 
-    in_proj = osr.SpatialReference()
-    in_proj.ImportFromEPSG(4326)
+    # Create a shapefile
+    auxiliary_shapefile = dirname(shp_file) + "/auxiliary_voronoi.shp"
+    sf = shapefile.Writer(auxiliary_shapefile, shapeType=shapefile.POLYGON)
 
-    transf = osr.CoordinateTransformation(in_proj, out_proj)
+    # Define the fields in the shapefile
+    sf.field("ClusterID", "N")
 
-    in_ds = ogr.Open(geojson_file)
-    if in_ds is None:
-        raise IOError("File {} doesn't exist".format(geojson_file))
-    in_layer = in_ds.GetLayer()
+    # Add polygons to the shapefile
+    for cluster_id, polygons in cluster_cells.items():
+        for polygon in polygons:
+            sf.poly([[(x, y) for x, y in polygon]])
+            sf.record(cluster_id)
 
-    mem_driver = ogr.GetDriverByName("MEMORY")
-    proj_ds = mem_driver.CreateDataSource("memData")
-    proj_layer = proj_ds.CreateLayer(
-        "clusters", out_proj, geom_type=ogr.wkbMultiPolygon
-    )
+    # Save the shapefile
+    sf.close()
+
+    with open(auxiliary_shapefile.replace(".shp", ".prj"), "w") as prj_file:
+        prj_file.write(pyproj.CRS.from_epsg(epsg_code).to_wkt(version="WKT1_ESRI"))
+
+    # Re-open the shapefile to merge same cluster ID polygons
+    sf = shapefile.Reader(auxiliary_shapefile)
+
+    # Extract the shapes and records (attributes)
+    shapes = sf.shapes()
+    records = sf.records()
+
+    # Create a dictionary to store merged polygons
+    merged_polygons = {}
+
+    # Iterate through the shapes and records
+    for shp, record in zip(shapes, records):
+        attribute_value = record[sf.fields[1].index("ClusterID")]
 
-    in_layer_def = in_layer.GetLayerDefn()
-    for i in range(0, in_layer_def.GetFieldCount()):
-        field_def = in_layer_def.GetFieldDefn(i)
-        proj_layer.CreateField(field_def)
-
-    out_layer_def = proj_layer.GetLayerDefn()
-
-    feature = in_layer.GetNextFeature()
-    while feature:
-        geom = feature.GetGeometryRef()
-        geom.Transform(transf)
-        proj_feat = ogr.Feature(out_layer_def)
-        proj_feat.SetGeometry(geom)
-        for i in range(0, out_layer_def.GetFieldCount()):
-            proj_feat.SetField(
-                out_layer_def.GetFieldDefn(i).GetNameRef(), feature.GetField(i)
+        # Convert the shape to a Shapely Polygon
+        polygon = shape(shp)
+
+        if attribute_value not in merged_polygons:
+            merged_polygons[attribute_value] = polygon
+        else:
+            merged_polygons[attribute_value] = merged_polygons[attribute_value].union(
+                polygon
             )
-        proj_layer.CreateFeature(proj_feat)
-        proj_feat = None
-        feature = in_layer.GetNextFeature()
 
-    return proj_ds
+    auxiliary_merged_shapefile = dirname(shp_file) + "/auxiliary_merged.shp"
+    sf_merged = shapefile.Writer(
+        auxiliary_merged_shapefile, shapeType=shapefile.POLYGON
+    )
+
+    # Copy the fields from the original shapefile to the merged shapefile
+    for field in sf.fields[1:]:
+        sf_merged.field(*field)
+
+    # Add merged polygons and their attributes to the merged shapefile
+    for attribute_value, merged_polygon in merged_polygons.items():
+        if isinstance(merged_polygon, MultiPolygon):
+            for (
+                geom
+            ) in merged_polygon.geoms:  # Iterate through the geometries in MultiPolygon
+                sf_merged.poly([list(geom.exterior.coords)])
+                sf_merged.record(attribute_value)
+        else:
+            sf_merged.poly([list(merged_polygon.exterior.coords)])
+            sf_merged.record(attribute_value)
+
+    with open(auxiliary_merged_shapefile.replace(".shp", ".prj"), "w") as prj_file:
+        prj_file.write(pyproj.CRS.from_epsg(epsg_code).to_wkt(version="WKT1_ESRI"))
+
+    # Save the merged shapefile
+    sf_merged.close()
+
+    clip_clusters_extent(auxiliary_merged_shapefile, shp_file, bounding_box, epsg_code)
+
+    with open(shp_file.replace(".shp", ".prj"), "w") as prj_file:
+        prj_file.write(pyproj.CRS.from_epsg(epsg_code).to_wkt(version="WKT1_ESRI"))
+
+    # Remove auxiliary files
+    for cluster_file in glob(auxiliary_merged_shapefile[:-3] + "*"):
+        remove(cluster_file)
+
+    for cluster_file in glob(auxiliary_shapefile[:-3] + "*"):
+        remove(cluster_file)
+
+
+def clip_clusters_extent(
+    input_shapefile_path: str,
+    output_shapefile_path: str,
+    bounding_box: list,
+    epsg_code: int,
+) -> None:
+    """Clip a shapefile to a defined extent. Bounding box coordinates must be in the
+    same projection as the input shapefile. The clipped shapfile will have the same
+    input projection but with a smaller extent.
+
+    Args:
+        input_shapefile_path (str): Path to a shapefile to clip.
+        output_shapefile_path (str): Path where clipped shapefile is saved.
+        bounding_box (list): Coordinates extent following (xmin, ymin, xmax, ymax).
+        epsg_code (int): EPSG code.
+    """
+    # Specify the extent as a bounding box (xmin, ymin, xmax, ymax)
+    xmin, ymin, xmax, ymax = bounding_box
+
+    # Open the input Shapefile
+    input_ds = ogr.Open(input_shapefile_path)
+
+    # Get the layer from the dataset
+    layer = input_ds.GetLayer()
+
+    # Create a geometry representing the clipping extent
+    clip_extent_geometry = ogr.Geometry(ogr.wkbPolygon)
+    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)
+    clip_extent_geometry.AddGeometry(ring)
+
+    # Set the spatial reference for the clipping extent
+    clip_extent_srs = osr.SpatialReference()
+    clip_extent_srs.ImportFromEPSG(epsg_code)
+    clip_extent_geometry.AssignSpatialReference(clip_extent_srs)
+
+    # Create an output Shapefile with the same structure as the input
+    output_driver = ogr.GetDriverByName("ESRI Shapefile")
+    output_ds = output_driver.CreateDataSource(output_shapefile_path)
+
+    # Create a layer in the output dataset
+    output_layer = output_ds.CreateLayer("clipped", geom_type=ogr.wkbMultiPolygon)
+
+    # Create a feature definition based on the original layer
+    output_layer.CreateFields(layer.schema)
+
+    # Perform the clipping operation
+    layer.SetSpatialFilter(clip_extent_geometry)
+    for feature in layer:
+        clipped_geometry = feature.GetGeometryRef().Intersection(clip_extent_geometry)
+        if clipped_geometry:
+            new_feature = ogr.Feature(output_layer.GetLayerDefn())
+            new_feature.SetGeometry(clipped_geometry)
+            for field_index in range(feature.GetFieldCount()):
+                field_name = feature.GetFieldDefnRef(field_index).GetName()
+                field_value = feature.GetField(field_name)
+                new_feature.SetField(field_name, field_value)
+            output_layer.CreateFeature(new_feature)
+            new_feature = None
+
+    # Clean up and close datasets
+    input_ds = None
+    output_ds = None
 
 
 def rasterize_clusters(
-    ds_in: ogr.DataSource, raster_config: dict, sigma: float = 15
+    shapefile_path: str, raster_config: dict, sigma: float = 15
 ) -> None:
-    """Rasterize clusters from the GeoJSON file generated by
+    """
+    Rasterize clusters from a GeoJSON file and save the result as a raster image.
 
     Args:
-        ds_in (ogr.DataSource): Input layer.
-        raster_config (dict): The output raster characteristics. Must include
-            'out_file' (path to save the rasterized clusters), 'size'
-            (raster size (x, y)) and 'geotransform' ([ul_x, x_res, x_rot,
-            ul_y, y_rot, y_res]) keys.
+        shapefile_path (str): The path to the GeoJSON file containing cluster
+            information.
+        raster_config (dict): A dictionary with the following required keys:
+            - 'out_file' (str): Path to save the rasterized clusters.
+            - 'size' (tuple): Raster size (x, y) in pixels.
+            - 'geotransform' (list): GeoTransform information [ul_x, x_res, x_rot,
+                                                               ul_y, y_rot, y_res].
         sigma (float, optional): Sigma parameter for a Gaussian filter. Defaults to 15.
 
+    Returns:
+        None
+
     Raises:
-        ValueError: If `raster_config` does not have all required keys.
+        ValueError: If `raster_config` does not have all the required keys.
     """
     if raster_config.keys() < {"out_file", "size", "geotransform"}:
         raise KeyError(
@@ -149,6 +286,8 @@ def rasterize_clusters(
             "parameter."
         )
 
+    ds_in = ogr.Open(shapefile_path)
+
     layer = ds_in.GetLayer()
     num_layers = layer.GetFeatureCount()
     proj = layer.GetSpatialRef()
@@ -165,7 +304,7 @@ def rasterize_clusters(
     ds_out.SetProjection(proj.ExportToWkt())
 
     for i in range(num_layers):
-        layer.SetAttributeFilter("cluster={}".format(i))
+        layer.SetAttributeFilter("ClusterID={}".format(float(i)))
         gdal.RasterizeLayer(ds_out, [i + 1], layer, burn_values=[1])
 
     data = ds_out.ReadAsArray().astype(np.float32)


Test patch:
diff --git a/pymica_tests/clusters_test.py b/pymica_tests/clusters_test.py
index db653db..7214f45 100644
--- a/pymica_tests/clusters_test.py
+++ b/pymica_tests/clusters_test.py
@@ -1,17 +1,17 @@
 """Tests creation of clusters files (shapefiles and rasters)
 """
+import json
 import unittest
+from glob import glob
 from os import remove
 from os.path import exists
 from tempfile import gettempdir
 
 from osgeo import gdal, ogr
-import json
 
 from pymica.utils.clusters import (
-    create_reprojected_geometries,
-    rasterize_clusters,
     create_clusters,
+    rasterize_clusters,
 )
 
 
@@ -21,23 +21,8 @@ class TestClusters(unittest.TestCase):
     with open("pymica_tests/data/sample_station_metadata.json", "r") as f_p:
         sample_data = json.load(f_p)
 
-    def test_create_reprojected_geometries(self):
-        """Test creation of reprojected geometries"""
-        d_s = create_reprojected_geometries("pymica_tests/data/clusters.json", 25831)
-
-        self.assertIsInstance(d_s, ogr.DataSource)
-        layer = d_s.GetLayer()
-
-        self.assertEqual(layer.GetFeatureCount(), 3)
-        self.assertTrue(layer.GetSpatialRef() is not None)
-
-        with self.assertRaises(IOError) as err:
-            create_reprojected_geometries("/bad_file.json", 25831)
-        self.assertEqual("File /bad_file.json doesn't exist", str(err.exception))
-
-        with self.assertRaises(ValueError) as err:
-            create_reprojected_geometries("pymica_tests/data/clusters.json", -999)
-        self.assertEqual("Wrong EPSG code: -999", str(err.exception))
+    clusters_test_file = "pymica_tests/data/test_clusters_3.shp"
+    bounding_box = [260000, 4488100, 530000, 4750000]
 
     def test_rasterize_clusters(self):
         """Test rasterization of clusters"""
@@ -46,9 +31,8 @@ def test_rasterize_clusters(self):
         if exists(out_file):
             remove(out_file)
 
-        layer = create_reprojected_geometries("pymica_tests/data/clusters.json", 25831)
         rasterize_clusters(
-            layer,
+            self.clusters_test_file,
             {"size": [1000, 1000], "geotransform": geotransform, "out_file": out_file},
         )
         d_s = gdal.Open(out_file)
@@ -60,7 +44,7 @@ def test_rasterize_clusters(self):
         # Errors
         with self.assertRaises(KeyError) as err:
             rasterize_clusters(
-                layer,
+                self.clusters_test_file,
                 {
                     "size": [1000, 1000],
                     "geotransform": geotransform,
@@ -74,7 +58,7 @@ def test_rasterize_clusters(self):
 
         with self.assertRaises(KeyError) as err:
             rasterize_clusters(
-                layer,
+                self.clusters_test_file,
                 {"out_file": gettempdir() + "/test.tiff", "geotransform": geotransform},
             )
         self.assertEqual(
@@ -85,7 +69,8 @@ def test_rasterize_clusters(self):
 
         with self.assertRaises(KeyError) as err:
             rasterize_clusters(
-                layer, {"out_file": gettempdir() + "/test.tiff", "size": [1000, 1000]}
+                self.clusters_test_file,
+                {"out_file": gettempdir() + "/test.tiff", "size": [1000, 1000]},
             )
         self.assertEqual(
             "`out_file`, `size`, `geotransform` must be in the `raster_config` "
@@ -95,71 +80,49 @@ def test_rasterize_clusters(self):
 
     def test_create_clusters(self):
         """Test creation of clusters with Kmeans"""
-        create_clusters(self.sample_data, 6, "pymica_tests/data/test_clusters_6.json")
-
-        self.assertTrue(exists("pymica_tests/data/test_clusters_6.json"))
+        clusters_file_path = "pymica_tests/data/test_clusters_file_3.shp"
 
-        with open("pymica_tests/data/test_clusters_6.json", "r") as f_p:
-            clusters = json.load(f_p)
-
-        self.assertEqual(
-            clusters["features"][0]["geometry"]["coordinates"], [2.18091, 41.39004]
-        )
-        self.assertEqual(
-            clusters["features"][0]["properties"],
-            {"cluster": 1, "id": "AN", "alt": 7.5},
+        create_clusters(
+            self.sample_data, 3, clusters_file_path, self.bounding_box, 25831
         )
 
-        self.assertEqual(
-            clusters["features"][-1]["geometry"]["coordinates"], [1.89716, 42.32211]
-        )
-        self.assertEqual(
-            clusters["features"][-1]["properties"],
-            {"cluster": 2, "id": "ZD", "alt": 2478},
-        )
+        self.assertTrue(exists(clusters_file_path))
+
+        clusters = ogr.Open(clusters_file_path)
+        layer = clusters.GetLayer()
 
-        if exists("pymica_tests/data/test_clusters_6.json"):
-            remove("pymica_tests/data/test_clusters_6.json")
+        self.assertEqual(layer.GetFeatureCount(), 3)
+        self.assertEqual(layer.GetFeature(0).ClusterID, 0.0)
+        self.assertEqual(layer.GetFeature(2).ClusterID, 2.0)
 
     def test_create_clusters_from_file(self):
         """Test creation of clusters with Kmeans from file"""
+        clusters_file_path = "pymica_tests/data/test_clusters_file_6.shp"
+
         create_clusters(
             "pymica_tests/data/sample_station_metadata.json",
             6,
-            "pymica_tests/data/test_clusters_6.json",
+            clusters_file_path,
+            self.bounding_box,
+            25831,
         )
 
-        self.assertTrue(exists("pymica_tests/data/test_clusters_6.json"))
+        self.assertTrue(exists(clusters_file_path))
 
-        with open("pymica_tests/data/test_clusters_6.json", "r") as f_p:
-            clusters = json.load(f_p)
-
-        self.assertEqual(
-            clusters["features"][0]["geometry"]["coordinates"], [2.18091, 41.39004]
-        )
-        self.assertEqual(
-            clusters["features"][0]["properties"],
-            {"cluster": 1, "id": "AN", "alt": 7.5},
-        )
+        clusters = ogr.Open(clusters_file_path)
+        layer = clusters.GetLayer()
 
-        self.assertEqual(
-            clusters["features"][-1]["geometry"]["coordinates"], [1.89716, 42.32211]
-        )
-        self.assertEqual(
-            clusters["features"][-1]["properties"],
-            {"cluster": 2, "id": "ZD", "alt": 2478},
-        )
+        self.assertEqual(layer.GetFeatureCount(), 6)
+        self.assertEqual(layer.GetFeature(0).ClusterID, 0.0)
+        self.assertEqual(layer.GetFeature(2).ClusterID, 2.0)
 
-        if exists("pymica_tests/data/test_clusters_6.json"):
-            remove("pymica_tests/data/test_clusters_6.json")
+    def tearDown(self) -> None:
+        cluster_files = glob("pymica_tests/data/test_clusters_file_3.*")
+        for cluster_file in cluster_files:
+            remove(cluster_file)
 
-    def test_create_clusters_wrong_dir(self):
-        """Test error raise from creation of clusters"""
-        with self.assertRaises(FileNotFoundError) as err:
-            create_clusters(
-                "pymica_tests/data/sample_station_metadata.json",
-                6,
-                "not/found/test_clusters_6.json",
-            )
+        cluster_files = glob("pymica_tests/data/test_clusters_file_6.*")
+        for cluster_file in cluster_files:
+            remove(cluster_file)
 
-        self.assertEqual(err.exception.args[0], "not/found directory does not exist.")
+        return super().tearDown()
diff --git a/pymica_tests/data/test_clusters_3.dbf b/pymica_tests/data/test_clusters_3.dbf
new file mode 100644
index 0000000..2ab6cef
Binary files /dev/null and b/pymica_tests/data/test_clusters_3.dbf differ
diff --git a/pymica_tests/data/test_clusters_3.prj b/pymica_tests/data/test_clusters_3.prj
new file mode 100644
index 0000000..382d168
--- /dev/null
+++ b/pymica_tests/data/test_clusters_3.prj
@@ -0,0 +1,1 @@
+PROJCS["ETRS_1989_UTM_Zone_31N",GEOGCS["GCS_ETRS_1989",DATUM["D_ETRS_1989",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",3.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]
\ No newline at end of file
diff --git a/pymica_tests/data/test_clusters_3.shp b/pymica_tests/data/test_clusters_3.shp
new file mode 100644
index 0000000..6095aa8
Binary files /dev/null and b/pymica_tests/data/test_clusters_3.shp differ
diff --git a/pymica_tests/data/test_clusters_3.shx b/pymica_tests/data/test_clusters_3.shx
new file mode 100644
index 0000000..bb55d91
Binary files /dev/null and b/pymica_tests/data/test_clusters_3.shx differ