1#!/usr/bin/env python3
  2# fmt: off
  3"""👋🌎 🌲🔥
  4# Raster clustering
  5## Usage
  6### Overview
  71. Choose your raster files
  82. Configure nodata, scaling strategies and weights in the `config.toml` file
  93. Choose "distance threshold" (or "number of clusters") for the [Agglomerative](https://scikit-learn.org/stable/modules/generated/sklearn.cluster.AgglomerativeClustering.html) clustering algorithm. Recommended:
 10   - Start with a distance threshold of 10.0 and decrease for more or increase for less clusters
 11   - After calibrating the distance threshold; 
 12   - [Sieve](https://gdal.org/en/latest/programs/gdal_sieve.html) small clusters (merge them to the biggest neighbor) with the `--sieve integer_pixels_size` option 
 13
 14### Execution
 15```bash
 16# get command line help
 17python -m fire2a.agglomerative_clustering -h
 18python -m fire2a.agglomerative_clustering --help
 19
 20# activate your qgis dev environment
 21source ~/pyqgisdev/bin/activate 
 22# execute 
 23(qgis) $ python -m fire2a.agglomerative_clustering -d 10.0
 24
 25# windows💩 users should use QGIS's python
 26C:\\PROGRA~1\\QGIS33~1.3\\bin\\python-qgis.bat -m fire2a.agglomerative_clustering -d 10.0
 27```
 28[More info on: How to windows 💩 using qgis's python](https://github.com/fire2a/fire2a-lib/tree/main/qgis-launchers)
 29
 30### Preparation
 31#### 1. Choose your raster files
 32- Any [GDAL compatible](https://gdal.org/en/latest/drivers/raster/index.html) raster will be read
 33- Place them all in the same directory where the script will be executed
 34- "Quote them" if they have any non alphanumerical chars [a-zA-Z0-9]
 35
 36#### 2. Preprocessing configuration
 37See the `config.toml` file for example of the configuration of the preprocessing steps. The file is structured as follows:
 38
 39```toml
 40["filename.tif"]
 41no_data_strategy = "most_frequent"
 42scaling_strategy = "onehot"
 43fill_value = 0
 44weight = 1
 45```
 46
 471. __scaling_strategy__
 48   - can be "standard", "robust", "onehot"
 49   - default is "robust"
 50   - [Standard](https://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.StandardScaler.html): (x-mean)/stddev
 51   - [Robust](https://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.RobustScaler.html): same but droping the tails of the distribution
 52   - [OneHot](https://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.OneHotEncoder.html): __for CATEGORICAL DATA__
 53
 542. __no_data_strategy__
 55   - can be "mean", "median", "most_frequent", "constant"
 56   - default is "mean"
 57   - categorical data should use "most_frequent" or "constant"
 58   - "constant" will use the value in __fill_value__ (see below)
 59   - [SimpleImputer](https://scikit-learn.org/stable/modules/generated/sklearn.impute.SimpleImputer.html)
 60
 613. __fill_value__
 62   - used when __no_data_strategy__ is "constant"
 63   - default is 0
 64   - [SimpleImputer](https://scikit-learn.org/stable/modules/generated/sklearn.impute.SimpleImputer.html)
 65
 664. __weight__
 67   - default is 1
 68   - used to give more importance to some features than others
 69   - This is done after the nodata imputation and scaling steps, before clustering
 70
 71
 72#### 3. Clustering configuration
 73
 741. __Agglomerative__ clustering algorithm is used. The following parameters are muttually exclusive:
 75- `-n` or `--n_clusters`: The number of clusters to form as well as the number of centroids to generate.
 76- `-d` or `--distance_threshold`: The linkage distance threshold above which, clusters will not be merged. When scaling start with 10.0 and downward (0.0 is compute the whole algorithm).
 77- More [parameters](https://scikit-learn.org/stable/modules/generated/sklearn.cluster.AgglomerativeClustering.html) for clustering can be passed directly into the pipelie method as keyword arguments
 78
 792. __Sieve filter__ is applied to remove small clusters. The sieve filter is applied using the [GDAL sieve library](https://gdal.org/en/latest/programs/gdal_sieve.html#gdal-sieve)
 80- `--sieve`: Use GDAL sieve filter to merge small clusters (number of pixels) into the biggest neighbor
 81
 82#### 4. Post-processing
 83Outputs can be:
 84- A raster file with the cluster labels and a polygon file with the cluster polygons
 85- A polygon file with the cluster polygons, with attribute being the number of pixels in each cluster
 86- A plot of the input data distributions, the rescaled data distributions, and the cluster size history and histogram (crashes QGIS in windows)
 87
 88Or use the `--script` option to return the label_map and the pipeline object for further processing in another python script:
 89    ```python
 90    from fire2a.agglomerative_clustering import main
 91    label_map, pipe1, pipe2 = main(["-d", "10.0", "-s"])
 92    ```
 93"""
 94# fmt: on
 95# from IPython.terminal.embed import InteractiveShellEmbed
 96# InteractiveShellEmbed()()
 97import logging
 98import sys
 99from pathlib import Path
100
101import numpy as np
102from osgeo import gdal, ogr, osr
103from sklearn.base import BaseEstimator, TransformerMixin
104from sklearn.cluster import AgglomerativeClustering
105from sklearn.compose import ColumnTransformer
106from sklearn.impute import SimpleImputer
107from sklearn.neighbors import radius_neighbors_graph
108from sklearn.pipeline import Pipeline
109from sklearn.preprocessing import OneHotEncoder, RobustScaler, StandardScaler
110
111from fire2a.utils import fprint, read_toml
112
113try:
114    GDT = gdal.GDT_Int64
115except:
116    GDT = gdal.GDT_Int32
117try:
118    OFT = ogr.OFTInteger64
119except:
120    OFT = ogr.OFTInteger
121
122logger = logging.getLogger(__name__)
123
124
125def check_shapes(data_list):
126    """Check if all data arrays have the same shape and are 2D.
127    Returns the shape of the data arrays if they are all equal.
128    """
129    from functools import reduce
130
131    def equal_or_error(x, y):
132        """Check if x and y are equal, returns x if equal else raises a ValueError."""
133        if x == y:
134            return x
135        else:
136            raise ValueError("All data arrays must have the same shape")
137
138    shape = reduce(equal_or_error, (data.shape for data in data_list))
139    if len(shape) != 2:
140        raise ValueError("All data arrays must be 2D")
141    height, width = shape
142    return height, width
143
144
145def get_map_neighbors(height, width, num_neighbors=8):
146    """Get the neighbors of each cell in a 2D grid.
147    n_jobs=-1 uses all available cores.
148    """
149
150    grid_points = np.indices((height, width)).reshape(2, -1).T
151
152    nb4 = radius_neighbors_graph(grid_points, radius=1, metric="manhattan", include_self=False, n_jobs=-1)
153    nb8 = radius_neighbors_graph(grid_points, radius=2 ** (1 / 2), metric="euclidean", include_self=False, n_jobs=-1)
154
155    # assert nb4.shape[0] == width * height
156    # assert nb8.shape[1] == width * height
157    # for n in range(width * height):
158    #     _, neighbors = np.nonzero(nb4[n])
159    #     assert 2<= len(neighbors) <= 4, f"{n=} {neighbors=}"
160    #     assert 3<= len(neighbors) <= 8, f"{n=} {neighbors=}"
161    return nb4, nb8
162
163
164class NoDataImputer(BaseEstimator, TransformerMixin):
165    """A custom Imputer that treats a specified nodata_value as np.nan and supports different strategies per column"""
166
167    def __init__(self, no_data_values, strategies, constants):
168        self.no_data_values = no_data_values
169        self.strategies = strategies
170        self.constants = constants
171        self.imputers = []
172        for no_data_value, strategy, constant in zip(no_data_values, strategies, constants):
173            if no_data_value:
174                self.imputers += [SimpleImputer(strategy=strategy, missing_values=no_data_value, fill_value=constant)]
175            else:
176                self.imputers += [SimpleImputer(strategy=strategy, fill_value=constant)]
177
178    def fit(self, X, y=None):
179        for i, imputer in enumerate(self.imputers):
180            imputer.fit(X[:, [i]], y)
181        return self
182
183    def transform(self, X):
184        for i, imputer in enumerate(self.imputers):
185            X[:, [i]] = imputer.transform(X[:, [i]])
186        self.output_data = X
187        return X
188
189
190class RescaleAllToCommonRange(BaseEstimator, TransformerMixin):
191    """A custom transformer that rescales all features to a common range [0, 1]"""
192
193    def __init__(self, weight_map):
194        self.weight_map = weight_map
195
196    def fit(self, X, y=None):
197        # Determine the combined range of all scaled features
198        self.min_val = [x.min() for x in X.T]
199        self.max_val = [x.max() for x in X.T]
200        return self
201
202    def transform(self, X):
203        # Rescale all features to match the common range
204        for i, (x, mi, ma) in enumerate(zip(X.T, self.min_val, self.max_val)):
205            if ma - mi == 0:
206                X.T[i] = x * self.weight_map[i]
207            else:
208                X.T[i] = (x - mi) / (ma - mi) * self.weight_map[i]
209        return X
210
211
212class CustomAgglomerativeClustering(BaseEstimator, TransformerMixin):
213    def __init__(self, height, width, neighbors=4, **kwargs):
214        self.height = height
215        self.width = width
216        self.neighbors = neighbors
217
218        self.grid_points = np.indices((height, width)).reshape(2, -1).T
219        if neighbors == 4:
220            connectivity = radius_neighbors_graph(
221                self.grid_points, radius=1, metric="manhattan", include_self=False, n_jobs=-1
222            )
223        elif neighbors == 8:
224            connectivity = radius_neighbors_graph(
225                self.grid_points, radius=2 ** (1 / 2), metric="euclidean", include_self=False, n_jobs=-1
226            )
227
228        self.connectivity = connectivity
229        self.kwargs = kwargs
230        self.model = AgglomerativeClustering(connectivity=self.connectivity, **self.kwargs)
231
232    def fit(self, X, y=None):
233        logger.debug("not sure why, but this method is never called alas needed")
234        self.model.fit(X)
235        return self
236
237    def fit_predict(self, X, y=None):
238        self.input_data = X
239        return self.model.fit_predict(X)
240
241
242def pipelie(observations, info_list, height, width, **kwargs):
243    """A scipy pipeline to achieve Agglomerative Clustering with connectivity on 2d matrix
244    Steps are:
245    1. Impute missing values
246    2. Scale the features
247    3. Rescale all features to a common range
248    4. Cluster the data using Agglomerative Clustering with connectivity
249    5. Reshape the labels back to the original spatial map shape
250    6. Return the labels and the pipeline object
251
252    Args:
253        observations (np.ndarray): The input data to cluster (n_samples, n_features) shaped
254        info_list (list): A list of dictionaries containing information about each feature
255        height (int): The height of the spatial map
256        width (int): The width of the spatial map
257        kwargs: Additional keyword arguments for AgglomerativeClustering, at least one of n_clusters or distance_threshold
258
259    Returns:
260        np.ndarray: The labels of the clusters, reshaped to the original 2d spatial map shape
261        Pipeline: The pipeline object containing all the steps of the pipeline
262    """
263    # kwargs = {"n_clusters": args.n_clusters, "distance_threshold": args.distance_threshold}
264
265    # imputer strategies
266    no_data_values = [info["NoDataValue"] for info in info_list]
267    no_data_strategies = [info["no_data_strategy"] for info in info_list]
268    fill_values = [info["fill_value"] for info in info_list]
269    weights = [info["weight"] for info in info_list]
270    # scaling_strategies = [info["scaling_strategy"] for info in info_list]
271
272    # scaling strategies
273    index_map = {}
274    for strategy in ["robust", "standard", "onehot"]:
275        index_map[strategy] = [i for i, info in enumerate(info_list) if info["scaling_strategy"] == strategy]
276    # index_map
277    # !cat config.toml
278
279    # Create transformers for each type
280    robust_transformer = Pipeline(steps=[("robust_step", RobustScaler())])
281    standard_transformer = Pipeline(steps=[("standard_step", StandardScaler())])
282    onehot_transformer = Pipeline(steps=[("onehot_step", OneHotEncoder(sparse_output=False))])
283    # OneHotEncoder._n_features_outs):
284
285    # Combine transformers using ColumnTransformer
286    feature_scaler = ColumnTransformer(
287        transformers=[
288            ("robust", robust_transformer, index_map["robust"]),
289            ("standard", standard_transformer, index_map["standard"]),
290            ("onehot", onehot_transformer, index_map["onehot"]),
291        ]
292    )
293
294    # # Create a temporary directory for caching calculations
295    # # FOR ACCESING STEPS LATER ON VERY LARGE DATASETS
296    # import tempfile
297    # import joblib
298    # temp_dir = tempfile.mkdtemp()
299    # memory = joblib.Memory(location=temp_dir, verbose=0)
300
301    # Create and apply the pipeline
302    # part 1 until feature scaling
303    pipe1 = Pipeline(
304        # n_features_in_ : int
305        # feature_names_in_ : ndarray of shape (`n_features_in_`,)
306        steps=[
307            ("no_data_imputer", NoDataImputer(no_data_values, no_data_strategies, fill_values)),
308            ("feature_scaling", feature_scaler),
309        ],
310        # memory=memory,
311        verbose=True,
312    )
313    # map weights to new columns (onehot feature scaler creates one column per category)
314    obs1 = pipe1.fit_transform(observations)
315    cat_names = pipe1.named_steps["feature_scaling"]["onehot"].get_feature_names_out()
316    split_names = [name.split("_")[0] for name in cat_names]
317    cat_count = np.unique(split_names, return_counts=True)[1]
318    onehot_map = {}
319    for i, key in enumerate(index_map["onehot"]):
320        onehot_map[key] = cat_count[i]
321    # onehot_map = {key: cat_count[i] for i, key in enumerate(index_map["onehot"])}
322    weight_map = []
323    for name, idxs in index_map.items():
324        for idx in idxs:
325            if name == "onehot":
326                weight_map += [weights[idx]] * onehot_map[idx]
327                continue
328            weight_map += [weights[idx]]
329    # part 2 use weight_map and cluster
330    pipe2 = Pipeline(
331        steps=[
332            ("common_rescaling", RescaleAllToCommonRange(weight_map)),
333            ("agglomerative_clustering", CustomAgglomerativeClustering(height, width, neighbors=4, **kwargs)),
334        ],
335        # memory=memory,
336        verbose=True,
337    )
338
339    # apply pipeLIE
340    labels = pipe2.fit_predict(obs1)
341
342    # Reshape the labels back to the original spatial map shape
343    labels_reshaped = labels.reshape(height, width)
344    return labels_reshaped, pipe1, pipe2
345
346
347def write(
348    label_map,
349    width,
350    height,
351    output_raster="",
352    output_poly="output.shp",
353    authid="EPSG:3857",
354    geotransform=(0, 1, 0, 0, 0, 1),
355    nodata=None,
356    feedback=None,
357):
358
359    from fire2a.processing_utils import get_output_raster_format, get_vector_driver_from_filename
360
361    # setup drivers for raster and polygon output formats
362    if output_raster == "":
363        raster_driver = "MEM"
364    else:
365        try:
366            raster_driver = get_output_raster_format(output_raster, feedback=feedback)
367        except Exception:
368            raster_driver = "GTiff"
369    try:
370        poly_driver = get_vector_driver_from_filename(output_poly)
371    except Exception:
372        poly_driver = "ESRI Shapefile"
373
374    # create raster output
375    src_ds = gdal.GetDriverByName(raster_driver).Create(output_raster, width, height, 1, GDT)
376    src_ds.SetGeoTransform(geotransform)  # != 0 ?
377    src_ds.SetProjection(authid)  # != 0 ?
378    #  src_band = src_ds.GetRasterBand(1)
379    #  if nodata:
380    #      src_band.SetNoDataValue(nodata)
381    #  src_band.WriteArray(label_map)
382
383    # create polygon output
384    drv = ogr.GetDriverByName(poly_driver)
385    dst_ds = drv.CreateDataSource(output_poly)
386    sp_ref = osr.SpatialReference()
387    sp_ref.SetFromUserInput(authid)  # != 0 ?
388    dst_lyr = dst_ds.CreateLayer("clusters", srs=sp_ref, geom_type=ogr.wkbPolygon)
389    dst_lyr.CreateField(ogr.FieldDefn("DN", OFT))  # != 0 ?
390    dst_lyr.CreateField(ogr.FieldDefn("pixel_count", OFT))
391    # dst_lyr.CreateField(ogr.FieldDefn("area", OFT))
392    # dst_lyr.CreateField(ogr.FieldDefn("perimeter", OFT))
393
394    # 0 != gdal.Polygonize( srcband, maskband, dst_layer, dst_field, options, callback = gdal.TermProgress)
395
396    # FAIL: All together it merges labels into a single polygon
397    #  src_band = src_ds.GetRasterBand(1)
398    #  if nodata:
399    #      src_band.SetNoDataValue(nodata)
400    #  src_band.WriteArray(label_map)
401    # gdal.Polygonize(src_band, None, dst_lyr, 0, callback=gdal.TermProgress)  # , ["8CONNECTED=8"])
402
403    # B separado
404    # for loop for creating each label_map value into a different polygonized feature
405    mem_drv = ogr.GetDriverByName("Memory")
406    tmp_ds = mem_drv.CreateDataSource("tmp_ds")
407    # itera = iter(np.unique(label_map))
408    # cluster_id = next(itera)
409    areas = []
410    pixels = []
411    data = np.zeros_like(label_map)
412    for cluster_id, px_count in zip(*np.unique(label_map, return_counts=True)):
413        # temporarily write band
414        src_band = src_ds.GetRasterBand(1)
415        src_band.SetNoDataValue(0)
416        data[label_map == cluster_id] = 1
417        src_band.WriteArray(data)
418        # create feature
419        tmp_lyr = tmp_ds.CreateLayer("", srs=sp_ref)
420        gdal.Polygonize(src_band, src_band.GetMaskBand(), tmp_lyr, -1)
421        # unset tmp data
422        data[label_map == cluster_id] = 0
423        # set polygon feat
424        feat = tmp_lyr.GetNextFeature()
425        geom = feat.GetGeometryRef()
426        featureDefn = dst_lyr.GetLayerDefn()
427        feature = ogr.Feature(featureDefn)
428        feature.SetGeometry(geom)
429        feature.SetField("DN", float(cluster_id))
430        areas += [geom.GetArea()]
431        pixels += [px_count]
432        feature.SetField("pixel_count", float(px_count))
433        # feature.SetField("area", int(geom.GetArea()))
434        # feature.SetField("perimeter", int(geom.Boundary().Length()))
435        dst_lyr.CreateFeature(feature)
436
437    fprint(f"Polygon Areas: {min(areas)=} {max(areas)=}", level="info", feedback=feedback, logger=logger)
438    fprint(f"Cluster PixelCounts: {min(pixels)=} {max(pixels)=}", level="info", feedback=feedback, logger=logger)
439    # RESTART RASTER
440    # src_ds = None
441    # src_band = None
442    # src_ds = gdal.GetDriverByName(raster_driver).Create(output_raster, width, height, 1, GDT)
443    # src_ds.SetGeoTransform(geotransform)  # != 0 ?
444    # src_ds.SetProjection(authid)  # != 0 ?
445    src_band = src_ds.GetRasterBand(1)
446    if nodata:
447        src_band.SetNoDataValue(nodata)
448    else:
449        # useless paranoia ?
450        src_band.SetNoDataValue(-1)
451    src_band.WriteArray(label_map)
452    # close datasets
453    src_ds.FlushCache()
454    src_ds = None
455    dst_ds.FlushCache()
456    dst_ds = None
457    return True
458
459
460def plot(labels_reshaped, pipe1, pipe2, info_list, **kwargs):
461    """Plot the observed values of the input data, the rescaled data, and the cluster size history and histogram.
462    Args:
463        labels_reshaped (np.ndarray): The reshaped labels of the clusters
464        pipe1 (Pipeline): The first pipeline object containing imputer and feature scaling steps
465        pipe2 (Pipeline): The second pipeline object containing the rescaling and clustering steps
466        info_list (list): A list of dictionaries containing information about each feature
467        **kargs: Additional keyword arguments
468            n_clusters (int): The number of clusters
469            distance_threshold (float): The linkage distance threshold
470            sieve (int): The number of pixels to use as a sieve filter
471            block (bool): Block the execution until the plot window is closed
472            filename (str): The filename to save the plot
473    """
474    from matplotlib import pyplot as plt
475
476    no_data_imputed = pipe1.named_steps["no_data_imputer"].output_data
477    pre_aggclu_data = pipe2.named_steps["agglomerative_clustering"].input_data
478
479    # filtrar onehot
480    num_onehots = sum([1 for i in info_list if i["scaling_strategy"] == "onehot"])
481    num_no_onehots = len(info_list) - num_onehots
482    pre_aggclu_data = pre_aggclu_data[:, :num_no_onehots]
483
484    # indices sin onehots
485    nohots_idxs = [i for i, info in enumerate(info_list) if info["scaling_strategy"] != "onehot"]
486
487    # filtrar onehot de no_data_treated
488    no_data_imputed = no_data_imputed[:, nohots_idxs]
489
490    # reordenados en robust y despues standard
491    rob_std_idxs = [i for i, j in enumerate(nohots_idxs) if info_list[j]["scaling_strategy"] == "robust"]
492    rob_std_idxs += [i for i, j in enumerate(nohots_idxs) if info_list[j]["scaling_strategy"] == "standard"]
493
494    # reordenar rob then std
495    pre_aggclu_data = pre_aggclu_data[:, rob_std_idxs]
496
497    names = [info_list[i]["fname"] for i in nohots_idxs]
498
499    fgs = np.array(plt.rcParams["figure.figsize"]) * 5
500    fig, axs = plt.subplots(3, 2, figsize=fgs)
501    suptitle = ""
502    if n_clusters := kwargs.get("n_clusters"):
503        suptitle = f"n_clusters: {n_clusters}"
504    if distance_threshold := kwargs.get("distance_threshold"):
505        suptitle = f"distance_threshold: {distance_threshold}"
506    if sieve := kwargs.get("sieve"):
507        suptitle += f", sieve: {sieve}"
508    if n_clusters or distance_threshold or sieve:
509        suptitle += f", resulting clusters: {len(np.unique(labels_reshaped))}"
510    suptitle += "\n(Not showing categorical data)"
511    fig.suptitle(suptitle)
512
513    # plot violin plot
514    axs[0, 0].violinplot(no_data_imputed, showmeans=False, showmedians=True, showextrema=True)
515    axs[0, 0].set_title("Violin Plot of NoData Imputed")
516    axs[0, 0].yaxis.grid(True)
517    axs[0, 0].set_xticks([y + 1 for y in range(num_no_onehots)], labels=names)
518    axs[0, 0].set_ylabel("Observed values")
519
520    # plot boxplot
521    axs[0, 1].boxplot(no_data_imputed)
522    axs[0, 1].set_title("Box Plot of NoData Imputed")
523    axs[0, 1].yaxis.grid(True)
524    axs[0, 1].set_xticks([y + 1 for y in range(num_no_onehots)], labels=names)
525    axs[0, 1].set_ylabel("Observed values")
526
527    # plot violin plot
528    axs[1, 0].violinplot(pre_aggclu_data, showmeans=False, showmedians=True, showextrema=True)
529    axs[1, 0].set_title("Violin Plot of Common Rescaled")
530    axs[1, 0].yaxis.grid(True)
531    axs[1, 0].set_xticks([y + 1 for y in range(num_no_onehots)], labels=names)
532    axs[0, 1].set_ylabel("Adjusted range")
533
534    # plot boxplot
535    axs[1, 1].boxplot(pre_aggclu_data)
536    axs[1, 1].set_title("Box Plot of Common Rescaled")
537    axs[1, 1].yaxis.grid(True)
538    axs[1, 1].set_xticks([y + 1 for y in range(num_no_onehots)], labels=names)
539    axs[0, 1].set_ylabel("Adjusted range")
540
541    # cluster history
542    unique_labels, counts = np.unique(labels_reshaped, return_counts=True)
543    axs[2, 0].plot(unique_labels, counts, marker="o", color="blue")
544    axs[2, 0].set_title("Cluster history size (in pixels)")
545    axs[2, 0].set_xlabel("Algorithm Step")
546    axs[2, 0].set_ylabel("Size (in pixels)")
547
548    # cluster histogram
549    axs[2, 1].hist(counts, log=True)
550    axs[2, 1].set_xlabel("Cluster Size (in pixels)")
551    axs[2, 1].set_ylabel("Number of Clusters")
552    axs[2, 1].set_title("Histogram of Cluster Sizes")
553
554    plt.tight_layout()
555    if filename := kwargs.get("filename"):
556        logger.info(f"Saving plot to {filename}")
557        plt.savefig(filename)
558    else:
559        if block := kwargs.get("block"):
560            plt.show(block=block)
561        else:
562            plt.show()
563
564
565def sieve_filter(data, threshold=2, connectedness=4, feedback=None):
566    """Apply a sieve filter to the data to remove small clusters. The sieve filter is applied using the GDAL library. https://gdal.org/en/latest/programs/gdal_sieve.html#gdal-sieve
567    Args:
568        data (np.ndarray): The input data to filter
569        threshold (int): The maximum number of pixels in a cluster to keep
570        connectedness (int): The number of connected pixels to consider when filtering 4 or 8
571        feedback (QgsTaskFeedback): A feedback object to report progress to use inside QGIS plugins
572    Returns:
573        np.ndarray: The filtered data
574    """
575    logger.info("Applying sieve filter")
576
577    height, width = data.shape
578    # fprint("antes", np.sort(np.unique(data, return_counts=True)), len(np.unique(data)), level="info", feedback=feedback, logger=logger)
579    num_clusters = len(np.unique(data))
580    src_ds = gdal.GetDriverByName("MEM").Create("sieve", width, height, 1, GDT)
581    src_band = src_ds.GetRasterBand(1)
582    src_band.WriteArray(data)
583    if 0 != gdal.SieveFilter(src_band, None, src_band, threshold, connectedness):
584        fprint("Error applying sieve filter", level="error", feedback=feedback, logger=logger)
585    else:
586        sieved = src_band.ReadAsArray()
587        src_band = None
588        src_ds = None
589        num_sieved = len(np.unique(sieved))
590        # fprint("despues", np.sort(np.unique(sieved, return_counts=True)), len(np.unique(sieved)), level="info", feedback=feedback, logger=logger)
591        fprint(
592            f"Reduced from {num_clusters} to {num_sieved} clusters, {num_clusters-num_sieved} less",
593            level="info",
594            feedback=feedback,
595            logger=logger,
596        )
597        fprint(
598            "Please try again increasing distance_threshold or reducing n_clusters instead...",
599            level="info",
600            feedback=feedback,
601            logger=logger,
602        )
603        # from matplotlib import pyplot as plt
604        # fig, (ax1, ax2) = plt.subplots(1, 2)
605        # ax1.imshow(data)
606        # ax1.set_title("before sieve" + str(len(np.unique(data))))
607        # ax2.imshow(sieved)
608        # ax2.set_title("after sieve" + str(len(np.unique(sieved))))
609        # plt.show()
610        # data = sieved
611        return sieved
612
613
614def arg_parser(argv=None):
615    """Parse command line arguments."""
616    from argparse import ArgumentDefaultsHelpFormatter, ArgumentParser
617
618    parser = ArgumentParser(
619        description="Agglomerative Clustering with Connectivity for raster data",
620        formatter_class=ArgumentDefaultsHelpFormatter,
621        epilog="More at https://fire2a.github.io/fire2a-lib",
622    )
623    parser.add_argument(
624        "config_file",
625        nargs="?",
626        type=Path,
627        help="For each raster file, configure its preprocess: nodata & scaling methods",
628        default="config.toml",
629    )
630
631    aggclu = parser.add_mutually_exclusive_group(required=True)
632    aggclu.add_argument(
633        "-d",
634        "--distance_threshold",
635        type=float,
636        help="Distance threshold (a good starting point when scaling is 10, higher means less clusters, 0 could take a long time)",
637    )
638    aggclu.add_argument("-n", "--n_clusters", type=int, help="Number of clusters")
639
640    parser.add_argument("-or", "--output_raster", help="Output raster file, warning overwrites!", default="")
641    parser.add_argument("-op", "--output_poly", help="Output polygons file, warning overwrites!", default="output.gpkg")
642    parser.add_argument("-a", "--authid", type=str, help="Output raster authid", default="EPSG:3857")
643    parser.add_argument(
644        "-g", "--geotransform", type=str, help="Output raster geotransform", default="(0, 1, 0, 0, 0, 1)"
645    )
646    parser.add_argument(
647        "-nw",
648        "--no_write",
649        action="store_true",
650        help="Do not write outputs raster nor polygons",
651        default=False,
652    )
653    parser.add_argument(
654        "-s",
655        "--script",
656        action="store_true",
657        help="Run in script mode, returning the label_map and the pipeline object",
658        default=False,
659    )
660    parser.add_argument(
661        "--sieve",
662        type=int,
663        help="Use GDAL sieve filter to merge small clusters (number of pixels) into the biggest neighbor",
664    )
665    parser.add_argument("--verbose", "-v", action="count", default=0, help="WARNING:1, INFO:2, DEBUG:3")
666
667    plot = parser.add_argument_group(
668        "Plotting, Visually inspect input distributions: NoData treated observations, rescaled data, with violing plots and boxplots. Also check output clustering size history and histograms"
669    )
670    plot.add_argument(
671        "-p",
672        "--plots",
673        action="store_true",
674        help="Activate the plotting routines",
675    )
676    plot.add_argument(
677        "-b",
678        "--block",
679        action="store_false",
680        default=True,
681        help="Block the execution until the plot window is closed. Use False for interactive ipykernels or QGIS",
682    )
683    plot.add_argument(
684        "-f",
685        "--filename",
686        type=str,
687        help="Filename to save the plot. If not provided, matplotlib will raise a window",
688    )
689    args = parser.parse_args(argv)
690    args.geotransform = tuple(map(float, args.geotransform[1:-1].split(",")))
691    if Path(args.config_file).is_file() is False:
692        parser.error(f"File {args.config_file} not found")
693    return args
694
695
696def main(argv=None):
697    """
698
699    args = arg_parser(["-d","10.0", "-g","(0, 10, 0, 0, 0, 10)", "config2.toml"])
700    args = arg_parser(["-d","10.0"]])
701    args = arg_parser(["-d","10.0", "config2.toml"])
702    args = arg_parser(["-n","10"])
703    """
704    if argv is sys.argv:
705        argv = sys.argv[1:]
706    args = arg_parser(argv)
707
708    if args.verbose != 0:
709        global logger
710        from fire2a import setup_logger
711
712        logger = setup_logger(verbosity=args.verbose)
713
714    logger.info("args %s", args)
715
716    # 2 LEE CONFIG
717    config = read_toml(args.config_file)
718    # logger.debug(config)
719
720    # 2.1 ADD DEFAULTS
721    for filename, file_config in config.items():
722        if "no_data_strategy" not in file_config:
723            config[filename]["no_data_strategy"] = "mean"
724        if "scaling_strategy" not in file_config:
725            config[filename]["scaling_strategy"] = "robust"
726        if "fill_value" not in file_config:
727            config[filename]["fill_value"] = 0
728        if "weight" not in file_config:
729            config[filename]["weight"] = 1
730    logger.debug(config)
731
732    # 3. LEE DATA
733    from fire2a.raster import read_raster
734
735    data_list, info_list = [], []
736    for filename, file_config in config.items():
737        data, info = read_raster(filename)
738        info["fname"] = Path(filename).name
739        info["no_data_strategy"] = file_config["no_data_strategy"]
740        info["scaling_strategy"] = file_config["scaling_strategy"]
741        info["fill_value"] = file_config["fill_value"]
742        info["weight"] = file_config["weight"]
743        data_list += [data]
744        info_list += [info]
745        logger.debug("%s", data[:2, :2])
746        logger.debug("%s", info)
747
748    # 4. VALIDAR 2d todos mismo shape
749    height, width = check_shapes(data_list)
750
751    # 5. lista[mapas] -> OBSERVACIONES
752    observations = np.column_stack([data.ravel() for data in data_list])
753
754    # 6. nodata -> feature scaling -> all scaling -> clustering
755    labels_reshaped, pipe1, pipe2 = pipelie(
756        observations,
757        info_list,
758        height,
759        width,
760        n_clusters=args.n_clusters,
761        distance_threshold=args.distance_threshold,
762    )  # insert more keyworded arguments to the clustering algorithm here!
763
764    # SIEVE
765    if args.sieve:
766        logger.info(f"Number of clusters before sieving: {len(np.unique(labels_reshaped))}")
767        labels_reshaped = sieve_filter(labels_reshaped, args.sieve)
768
769    logger.info(f"Final number of clusters: {len(np.unique(labels_reshaped))}")
770
771    # 7 debbuging plots
772    if args.plots:
773        plot(labels_reshaped, pipe1, pipe2, info_list, **vars(args))
774
775    # 8. ESCRIBIR RASTER
776    if not args.no_write:
777        if not write(
778            labels_reshaped,
779            width,
780            height,
781            output_raster=args.output_raster,
782            output_poly=args.output_poly,
783            authid=args.authid,
784            geotransform=args.geotransform,
785        ):
786            logger.error("Error writing output raster")
787
788    # 9. SCRIPT MODE
789    if args.script:
790        return labels_reshaped, pipe1, pipe2
791
792    return 0
793
794
795if __name__ == "__main__":
796    sys.exit(main(sys.argv))