From f842af4a56a59714ea44408312155996122533ae Mon Sep 17 00:00:00 2001
From: Claudio Zandonella <claudiozandonella@gmail.com>
Date: Tue, 21 Feb 2023 20:13:34 +0100
Subject: [PATCH] add quarto report

- automatically created in the report folder with plots and tables
---
 justclust/analysis/bolzano.py |    7 +
 justclust/data/bolzano.py     |   38 +
 justclust/quarto.py           | 1452 +++++++++++++++++++++++++++++++++
 3 files changed, 1497 insertions(+)
 create mode 100644 justclust/quarto.py

diff --git a/justclust/analysis/bolzano.py b/justclust/analysis/bolzano.py
index 5981cca..733047e 100644
--- a/justclust/analysis/bolzano.py
+++ b/justclust/analysis/bolzano.py
@@ -19,6 +19,7 @@ from justclust.paths.bolzano import (
     selclfile,
 )
 from justclust.plots import reporting
+import justclust.quarto as qmd
 
 # %%
 # read and filter/select data
@@ -237,4 +238,10 @@ for clname in sel_clsts:
         preprocs=pre,
     )
 
+    qmd.create_qmd_report(
+        city = 'bolzano',
+        clname= clname,
+        cdir = report_dir / clname
+        )
+
 # %%
diff --git a/justclust/data/bolzano.py b/justclust/data/bolzano.py
index 1edfdf9..e9de98a 100644
--- a/justclust/data/bolzano.py
+++ b/justclust/data/bolzano.py
@@ -174,6 +174,44 @@ conv = {
     "%Single": "Families with one component [%]",
 }
 
+# %%
+# define subset of columns according to categories
+
+selected_col = {
+    'col_aq':[
+        'area very-low AQ-risk [%]',
+        'area low AQ-risk [%]',
+        'area medium-low AQ-risk [%]',
+        'area medium-high AQ-risk [%]',
+        'area high AQ-risk [%]',
+        'area very-high AQ-risk [%]'
+        ],
+    'col_carbon':[
+        "Carbon emission building [ton CO2/m²]",
+        "Carbon absorption vegetation [kg CO2/m²]"
+        ],
+    'col_char':[
+        "Protected area [%]",
+        "Accessibility urban green areas (<5 min) [%]",
+        "Walkability Index",
+        "Soil sealing between 2022 and 2018 [%]",
+        "Surface urban heat island"
+        ],
+    'col_socio':[
+        "Age Index",
+        "Foreign population [%]",
+        "Families with children [%]",
+        "Families with one component [%]"
+        ]
+    }
+
+
+
+
+
+
+
+
 # %%
 # Define function to read files
 
diff --git a/justclust/quarto.py b/justclust/quarto.py
new file mode 100644
index 0000000..50a83a5
--- /dev/null
+++ b/justclust/quarto.py
@@ -0,0 +1,1452 @@
+#----    Settings    ----
+
+import os
+import sys
+from pathlib import Path
+
+import geopandas as gpd
+import numpy as np
+import pandas as pd
+import string
+from math import ceil, floor
+import re
+
+from sklearn import ensemble as ens
+
+import matplotlib.pyplot as plt
+import seaborn as sns
+from matplotlib.colors import LinearSegmentedColormap
+
+#----    get_dict_stat
+def get_dict_stat(
+        data,
+        selected_col,
+        var_cluster = 'cluster_lab',
+        exclude_Z = True
+        ):
+    data = data[['cluster_lab'] + selected_col].copy()
+
+    if exclude_Z:
+        data = data[data['cluster_lab'] != 'Z']
+    
+    selected_grouped = data.groupby(['cluster_lab'])
+    cluster_means = selected_grouped.mean()
+    cluster_medians = selected_grouped.median()
+    overall_mean = cluster_means.mean()
+    overall_median = cluster_medians.mean()
+
+    res = {
+        'cluster_means':cluster_means.T,
+        'cluster_means_cent':cluster_means.T\
+            .sub(overall_mean, axis = 0)\
+                .divide(np.abs(overall_mean), axis = 0),
+
+        'cluster_medians':cluster_medians.T,
+        'cluster_medians_cent':cluster_medians.T\
+            .sub(overall_median, axis = 0)\
+                .divide(np.abs(overall_median), axis = 0)
+    }
+
+    return res
+
+#----    summary_stat    ----
+
+def summary_stat(data, round = 2):
+    data_summary = data.describe().T
+
+    res = data_summary.round(decimals = round)
+    return res
+
+#----    table_latex    ----
+
+def table_latex(
+        data,
+        caption = None,
+        label = None,
+        size = None,
+        position = 'H'
+        ):
+    
+    res = data.to_latex(
+        position= position,
+        caption = caption,
+        label = label
+        )
+    
+    if size is not None:
+        res = re.sub(
+                pattern = 'centering',
+                repl = 'centering\\n\\\\'+size,
+                string = res
+                )
+        
+    return res
+
+#----    cluster_color_dict    ----
+
+def cluster_color_dict(cluster_labels):
+    n_clusters = len(cluster_labels)
+
+    palette = list(sns.color_palette("colorblind", n_colors=max(8, n_clusters)))
+    gray_color = palette[7]
+    palette.remove(gray_color)
+    palette =  palette[0: n_clusters-1] + [gray_color]
+    res = {
+        label:rgb for label, rgb in zip(cluster_labels.values(), palette)
+    }
+
+    return res
+
+#----    get_cluster_lab    ----
+
+def get_cluster_labels(cluster_var) :
+    
+    unique_values = sorted(cluster_var.unique().astype(int))
+
+    res = {
+        i:string.ascii_uppercase[i] if i >= 0 else 'Z' for i in unique_values
+        }
+    
+    # sort by value
+    res = dict(sorted(res.items(), key=lambda item: item[1]))
+    return res
+
+#----    get_clst_feature_importance    ----
+
+def get_clst_feature_importance(
+        data, 
+        labels, 
+        max_depth = 10
+        ):
+    
+    cl_cat = sorted(set(labels))
+    
+    res = pd.DataFrame()
+
+    clf = ens.RandomForestClassifier(
+        max_depth=max_depth,
+        n_estimators=500,
+        random_state=1,
+        n_jobs = -1
+        )
+
+    for cl in cl_cat:
+
+        clf.fit(
+            data, 
+            labels == cl)
+        fi = pd.DataFrame(
+            clf.feature_importances_,
+            index=pd.Series(data.columns, name="features"),
+            columns=[cl]
+            )
+       
+        res = pd.concat([res, fi], axis='columns')
+
+    return res
+
+#----    get_summary_clst_frequency    ----
+
+def get_summary_clst_frequency(
+        data,
+        var_clust = 'cluster_lab'
+        ):
+    res = pd.DataFrame({
+        'Frequency':data[var_clust].value_counts(sort = False).sort_index()
+        })
+    res['Percent'] = res['Frequency'].div(res['Frequency'].sum())
+    return res
+
+#----    is_light    ----
+
+def is_light(r, g, b, hsp_threshold=0.5):
+    # HSP (Highly Sensitive Poo) equation is from http://alienryderflex.com/hsp.html
+    hsp = np.sqrt(0.299 * (r * r) + 0.587 * (g * g) + 0.114 * (b * b))
+    return True if hsp > hsp_threshold else False
+
+#----    plot_map_included    ----
+
+def plot_map_included(
+        geo_out,
+        clname,
+        figsize = (12,12)
+        ):
+    
+    data_plot = geo_out[[clname, 'geometry']].copy()
+    data_plot['Inclusion'] = data_plot[clname].map(lambda x: 'Excluded' if np.isnan(x) else 'Included' )
+
+    sns_colors = list(sns.color_palette())
+    dict_colors = {
+        'Included':sns_colors[0],
+        'Excluded':sns_colors[1]}
+
+    fig, ax = plt.subplots(1, 1, figsize=figsize)
+    data_plot.plot(
+        ax=ax, alpha = .7, column = 'Inclusion', legend = True, categorical=True, 
+        categories = list(dict_colors.keys()),
+        cmap = LinearSegmentedColormap.from_list("a", list(dict_colors.values())))
+    data_plot.boundary.plot(ax=ax, color='white', linewidth = .35)
+    ax.axis('off')
+    plt.show()
+
+#----    plot_summary    ----
+
+def plot_summary(
+        data, 
+        bins = 25, 
+        text_size = 6, 
+        title_size = None, 
+        color = None):
+    n_var = data.shape[1]
+
+    if n_var <= 3:
+        n_col = n_var
+        n_row = 1
+    elif n_var == 4:
+        n_col = 2
+        n_row = 2
+    else:
+        n_col = 3
+        n_row = ceil(n_var/3)
+    
+    
+    
+    fig, ax = plt.subplots(1, 1, figsize=(n_col*2.5, n_row*2))
+    fig = pd.DataFrame(data).hist(
+        ax = ax,
+        layout = (n_row, n_col),
+        alpha = .7,
+        bins = bins)
+    
+    # Reduce font
+    if title_size is None:
+        title_size = text_size + 2
+
+    for x in fig.ravel():
+        for item in ([x.xaxis.label, 
+                      x.yaxis.label] +
+                      x.get_xticklabels() + 
+                      x.get_yticklabels()):
+            item.set_size(text_size)
+        x.title.set_size(title_size)
+
+        if color is not None:
+            for rect in x.patches:
+                rect.set_color(color)
+    
+
+    plt.show()
+
+#----    plot_heatmap_corr    ----
+
+def plot_heatmap_corr(
+        data_plot,
+        val_filter = 0.3,
+        figsize = (15, 15),
+        ):
+    fig, ax = plt.subplots(figsize=figsize)
+    data_cor = data_plot.corr()
+    
+    ax = sns.heatmap(
+        data_cor,
+        vmin = -1,
+        vmax = 1,
+        center = 0,
+        annot=True,
+        cmap=sns.color_palette("vlag", as_cmap=True),
+        )
+    
+    for t in ax.texts:
+        if np.abs(float(t.get_text()))>=val_filter:
+            t.set_text(t.get_text()) #if the value is greater than 0.3 then I set the text 
+        else:
+            t.set_text("")
+
+    plt.show()
+
+#----    plot_map_all_cluster    ----
+
+def plot_map_all_cluster(
+        data_plot,
+        cluster_colors,
+        figsize = (12, 12)
+        ):
+    fig, ax = plt.subplots(1, 1, figsize=figsize)
+    data_plot.plot(
+        column = 'cluster_lab', 
+        cmap = LinearSegmentedColormap.from_list("a", list(cluster_colors.values())),
+        ax=ax, 
+        categorical=True, 
+        legend=True
+    )
+    data_plot.boundary.plot(ax=ax, color='white', linewidth = .25)
+    ax.axis('off')
+    plt.show()
+
+#----    plot_clst_frequency    ----
+
+def plot_clst_frequency(
+        data,
+        cluster_colors,
+        var_cluster = 'cluster_lab',
+        figsize = (7, 4)
+        ):
+    gs_kw = dict(
+        width_ratios=[1, 2, 1], 
+        height_ratios=[1]
+        )
+    
+    fig, ax = plt.subplots(
+        ncols = 3,
+        nrows = 1,
+        figsize=figsize,
+        constrained_layout=True,
+        gridspec_kw=gs_kw
+        )
+    
+    ax[0].axis('off')
+    ax[2].axis('off')
+
+    sns.countplot(
+        ax = ax[1],
+        data=data, 
+        x=var_cluster, 
+        palette=cluster_colors,
+        order=cluster_colors)
+    ax[1].set(
+        ylabel="Count",
+        xlabel="Cluster",
+    )
+    plt.show()
+
+#----    plot_clst_boxplot_single    ----
+
+def plot_clst_boxplot_single(
+        data_plot,
+        data_plot_means,
+        data_plot_median,
+        data_fi,
+        sel_var,
+        list_clust,
+        cluster_colors,
+        ax
+        ):
+    fig_sns = sns.boxplot(
+        x=sel_var, 
+        y="cluster_lab",
+        data=data_plot, 
+        palette=cluster_colors,
+        order= list_clust,
+        width=0.6,
+        ax = ax,
+        )
+    fig_sns.set(ylabel=None)
+
+    # set alpha according to feature importance
+    max_fi = max(data_fi.max())
+    min_fi = min(data_fi.min())
+    delta_fi = max_fi - min_fi
+
+    for i, patch in enumerate(fig_sns.patches):
+        alpha = .05 + .95 * (data_fi.loc[sel_var][i] - min_fi)/ delta_fi
+        r, g, b, a = patch.get_facecolor()
+        patch.set_facecolor((r, g, b, alpha))
+
+    ax.plot(
+        data_plot_means.loc[sel_var,:],
+        data_plot_means.columns,
+        's',
+        color="k",
+        marker="*",
+        markeredgecolor="k"
+        )
+    
+    overall_mean =  data_plot_means.loc[sel_var,:].mean()
+    ax.axvline(overall_mean, color='r', linestyle='--', alpha=0.5)
+    overall_median = data_plot_median.loc[sel_var,:].mean()
+    ax.axvline(overall_median, color='r', linestyle='-', alpha=0.5)
+    
+    ax.xaxis.grid(True, linestyle="-", which="major", color="lightgrey", alpha=0.5)
+    max_x = max(data_plot[sel_var])
+    min_x = min(data_plot[sel_var])
+    delta = max_x - min_x
+    ax.set_xlim([min_x -.05*delta, max_x + .25*delta])
+
+    for i in range(0, len(list_clust)):
+        ax.text(max_x + .23*delta, i, f'{data_fi.loc[sel_var][i]*100:.1f}%', 
+                horizontalalignment='right', verticalalignment = 'center',
+                size='small', color='black')
+
+#----    plot_clst_boxplot_loop    ----
+
+def plot_clst_boxplot_loop(
+        data_plot,
+        data_plot_means,
+        data_plot_median,
+        data_fi,
+        list_vars,
+        list_clust,
+        cluster_colors,
+        figsize
+        ):
+    n_plots = len(list_vars)
+
+    if n_plots == 1:
+        widths = [1, 1, 1]
+        ncols = 3
+        grid = [0, 1, 0]
+    elif n_plots in [2, 4]:
+        widths = [1, 2, 2, 1]
+        ncols = 4
+        grid = [0, 1, 1, 0, 0, 1, 1, 0]
+    else:
+        widths = [1, 1, 1]
+        ncols = 3
+        grid = [1] * n_plots
+
+    nrows = ceil(n_plots / 3)
+
+
+    fig, axes = plt.subplots(
+        nrows=nrows,
+        ncols=ncols,
+        sharex=False,
+        sharey=False,
+        figsize=figsize,
+        constrained_layout=True,
+        gridspec_kw=dict(
+            width_ratios=widths, 
+            height_ratios= [1] * nrows
+        )
+    )
+
+    i = 0
+    for ax, to_plot in zip(axes.flat, grid + [0] * 4):
+        if to_plot == 0:
+            # removed unused plots in the grid
+            ax.axis('off')
+            continue
+
+        sel_var = list_vars[i]
+
+        plot_clst_boxplot_single(
+            data_plot = data_plot,
+            data_plot_means = data_plot_means,
+            data_plot_median = data_plot_median,
+            data_fi = data_fi.loc[list_vars],
+            sel_var = sel_var,
+            list_clust = list_clust,
+            cluster_colors = cluster_colors,
+            ax = ax
+        )
+        
+        i += 1
+                
+#----    plot_clst_heatmap    ----
+
+def plot_clst_heatmap(
+        data,
+        col_avg,
+        fig_size = (10, 8)
+        ):
+    fig, ax = plt.subplots(figsize=fig_size)
+    data_plot = data.sub(data[col_avg], axis = 0)\
+        .div(np.abs(data[col_avg]), axis = 0)
+    data_text = data.round(2)
+
+    ax = sns.heatmap(
+        data_plot,
+        # vmin = -1,
+        # vmax = 1,
+        center = 0,
+        annot=data_text,
+        cmap=sns.color_palette("vlag", as_cmap=True),
+        cbar=False,
+        fmt='.3g'
+        )
+    ax.set(xlabel = None)
+    plt.show()
+
+#----    plot_cluster_comp    ----
+
+def plot_cluster_comp(
+    serie_plot, 
+    ax,
+    cluster,
+    color,
+    serie_fi = None,
+    text_fi = True,
+    xmin=-2.5,
+    xmax=2.5
+    ):
+
+    fig_bar = ax.barh(serie_plot.index, serie_plot, label=cluster, color=color)
+
+    # set alpha according to feature importance
+    if serie_fi is not None:
+        max_fi = serie_fi.max()
+        min_fi = serie_fi.min()
+        delta_fi = max_fi - min_fi
+
+        for i, patch in enumerate(fig_bar.patches):
+            alpha = .19 + .8 * (serie_fi.iloc[i] - min_fi)/ delta_fi
+            r, g, b, a = patch.get_facecolor()
+            patch.set_facecolor((r, g, b, alpha))
+
+
+    light_bg = is_light(*color, hsp_threshold=0.6)
+    ax.set(title=f"Cluster {cluster}")
+    if xmin is not None or xmax is not None:
+        ax.set_xlim(xmin, xmax)
+        for idx, val in enumerate(serie_plot):
+            if val < xmin:
+                ax.text(
+                    xmin + 0.1 * abs(xmin),
+                    idx,
+                    f"{val:.2f}",
+                    horizontalalignment="center",
+                    verticalalignment="center",
+                    size=7,
+                    weight="semibold",
+                    color="black" if light_bg else "white",
+                )
+            elif val > xmax:
+                ax.text(
+                    xmax - 0.1 * xmax,
+                    idx,
+                    f"{val:.2f}",
+                    horizontalalignment="center",
+                    verticalalignment="center",
+                    size=7,
+                    weight="semibold",
+                    color="black" if light_bg else "white",
+                )
+            
+            # add fi text
+            if serie_fi is not None and text_fi is True:
+                ax.text(xmax * 1.3, idx, f'{serie_fi.iloc[idx]*100:.1f}%', 
+                    horizontalalignment='right', verticalalignment = 'center',
+                    size='small', color='black')
+
+#----    plot_cluster_comp_loop    ----
+
+def plot_cluster_comp_loop(
+    data_plot,
+    data_fi,
+    colors,
+    ncols = 3,
+    xmin=-2.5,
+    xmax=2.5,
+    xlabel=None,
+    ylabel=None
+    ):
+    
+    nvar, ncl = data_plot.shape
+    nrows = ceil(ncl / ncols)
+
+    fig, axes = plt.subplots(
+        nrows=nrows,
+        ncols=ncols,
+        sharex=False,
+        sharey=True,
+        figsize=(4 * 3, 4.5 * nrows),
+    )
+    for c, cl in enumerate(data_plot.columns):
+        irow = floor(c / ncols)
+        icol = floor(c - (irow * ncols))
+        ax = axes[irow][icol]
+    
+        plot_cluster_comp(
+            serie_plot = data_plot[cl],
+            serie_fi = data_fi[cl],
+            ax=ax,
+            cluster = cl,
+            color = colors[cl],
+            xmin=xmin,
+            xmax=xmax,
+        )
+        if icol == 0:
+            ax.set(ylabel=ylabel)
+        if irow == (nrows - 1):
+            ax.set(xlabel=xlabel)
+        plt.subplots_adjust(wspace=0.25)
+    
+    # removed unused plots in the grid
+    nplots = nrows * ncols
+    if nplots > ncl:
+        
+        for i in range(nplots - ncl, 0, -1):
+            axes[-1, -i].axis('off')
+
+    plt.show()
+
+#----    plot_cluster_comp_single    ----
+
+def plot_cluster_comp_single(
+    dict_stat,
+    data_fi,
+    cluster_lab,
+    colors,
+    xmin=-2.5,
+    xmax=2.5,
+    xlabel=None,
+    ylabel=None,
+    figsize=(10, 5),
+    sorted = False
+    ):
+    data_plot_mean = dict_stat.get('cluster_means_cent')[cluster_lab].copy()
+    data_plot_median = dict_stat.get('cluster_medians_cent')[cluster_lab].copy()
+
+    serie_fi = data_fi[cluster_lab].copy()
+    if sorted:
+        serie_fi = serie_fi.sort_values(ascending=True)
+        data_plot_mean = data_plot_mean.reindex(serie_fi.index)
+        data_plot_median = data_plot_median.reindex(serie_fi.index)
+        
+
+    fig, axes = plt.subplots(
+        nrows=1,
+        ncols=2,
+        sharex=True,
+        sharey=True,
+        figsize=figsize,
+    )
+    plot_cluster_comp(
+            serie_plot=data_plot_mean,
+            serie_fi= serie_fi,
+            text_fi = False,
+            ax=axes[0],
+            cluster = cluster_lab,
+            color = colors,
+            xmin=xmin,
+            xmax=xmax,
+        )
+    axes[0].set(title='Mean')
+
+    plot_cluster_comp(
+                serie_plot=data_plot_median,
+                serie_fi= serie_fi,
+                ax=axes[1],
+                cluster = cluster_lab,
+                color = colors,
+                xmin=xmin,
+                xmax=xmax,
+                )
+    axes[1].set(title='Median')
+    plt.show()
+
+#----    plot_map_single_cluster    ----
+
+def plot_map_single_cluster(
+    data_plot, 
+    mask,
+    cluster_lab,
+    cluster_colors,
+    figsize=(7,7)
+    ):
+    fig, ax = plt.subplots(1, 1, figsize=figsize)
+    data_plot.boundary.plot(ax=ax, color = '#bdbdbd', alpha = .4, linewidth = .7)
+    data_plot.loc[mask].plot(
+        column = 'cluster_lab', 
+        color = cluster_colors[cluster_lab],
+        ax=ax, 
+        categorical=True, 
+        legend=True
+        )
+    ax.set(title=f'Map cluster {cluster_lab}')
+    ax.axis('off')
+    plt.show()
+
+#----    plot_boxplot_comp_cluster_single    ----
+
+def plot_boxplot_comp_cluster_single(
+        ax,
+        data_cluster,
+        cluster_lab,
+        sel_var,
+        data_fi,
+        data_plot_means,
+        data_plot_median,
+        cluster_colors
+    ):
+    list_cluster = sorted(set(data_cluster['cluster_lab']))
+    list_cluster.remove(cluster_lab)
+    list_cluster = [cluster_lab] + list_cluster
+    list_cluster = list_cluster[::-1]
+
+    positions = [.15 * x for x in list(range(0,len(list_cluster)-1))]
+    positions.append(positions[-1] +.5)
+    widths = [.1] * (len(list_cluster) - 1) + [.2]
+
+    # set alpha according to feature importance
+    max_fi = data_fi[cluster_lab].max()
+    min_fi = data_fi[cluster_lab].min()
+    delta_fi = max_fi - min_fi
+    alpha = .05 + .95 * (data_fi.loc[sel_var, cluster_lab] - min_fi)/ delta_fi
+
+    for i, letter in enumerate(list_cluster):
+        mask_iter = data_cluster['cluster_lab'] == letter
+        color_alpha = (*cluster_colors[letter], alpha)
+        ax.boxplot(
+            data_cluster.loc[mask_iter,sel_var], positions=[positions[i]],  vert=False, widths = widths[i],
+            patch_artist=True,
+            boxprops=dict(facecolor= color_alpha)
+        )
+        ax.plot(
+            data_plot_means.loc[sel_var,letter],
+            positions[i], 
+            's', color="k", marker="*", markeredgecolor="k")
+
+    # add overall mean and median
+    overall_mean =  data_plot_means.loc[sel_var,:].mean()
+    ax.axvline(overall_mean, color='r', linestyle='--', alpha=0.5)
+    overall_median = data_plot_median.loc[sel_var,:].mean()
+    ax.axvline(overall_median, color='r', linestyle='-', alpha=0.5)
+
+    # add fi 
+    ax.text(
+        data_plot_means.loc[sel_var,cluster_lab], positions[-1] + .15, 
+        f'fi: {data_fi.loc[sel_var, cluster_lab]*100:.1f}%', 
+        horizontalalignment='center', verticalalignment = 'center',
+        size='small', color='black')
+
+    ax.set_yticks(positions, list_cluster)
+    ax.set_ylim(-.1, positions[-1] +.2 )
+    ax.set_xlabel(sel_var)
+
+#----    plot_boxplot_single_cluster_loop    ----
+
+def plot_boxplot_comp_cluster_loop(
+    data_cluster,
+    cluster_lab,
+    list_vars,
+    data_fi,
+    data_plot_means,
+    data_plot_median,
+    cluster_colors,
+    figsize,
+    sorted = True
+    ):
+
+    if sorted is True:
+        list_vars = list(data_fi.loc[list_vars]\
+                         .sort_values(by = [cluster_lab], ascending=False).index)
+
+    n_plots = len(list_vars)
+
+    if n_plots <= 3:
+        ncols = n_plots
+    elif n_plots == 4:
+        ncols = 2
+    else:
+        ncols = 3
+
+    nrows = ceil(n_plots / ncols)
+
+    fig, axes = plt.subplots(
+        nrows=nrows,
+        ncols=ncols,
+        sharex=False,
+        sharey=False,
+        figsize=figsize
+    )
+    for i, sel_var in enumerate(list_vars):
+        irow = floor(i / ncols)
+        icol = floor(i - (irow * ncols))
+        
+        if nrows > 1:
+            ax = axes[irow][icol]
+        else:
+            ax = axes[icol]
+
+        plot_boxplot_comp_cluster_single(
+            ax = ax,
+            data_cluster = data_cluster,
+            cluster_lab = cluster_lab,
+            sel_var = sel_var,
+            data_fi = data_fi,
+            data_plot_means = data_plot_means,
+            data_plot_median = data_plot_median,
+            cluster_colors = cluster_colors
+            )
+        # removed unused plots in the grid
+        n_grid = nrows * ncols
+        if n_grid > n_plots:
+
+            for i in range(n_grid - n_plots, 0, -1):
+                if nrows > 1:
+                    axes[-1, -i].axis('off')
+                else:
+                    axes[-i].axis('off')
+    plt.show()
+    
+#----    create_qmd_report    ----
+
+def create_qmd_report(
+        city,
+        clname,
+        cdir: Path,
+    ):
+    graph_open = "{"
+    graph_close = "}"
+
+    header =f"""---
+title: 'JUSTNature - {city}'
+format:
+  pdf:
+    toc: true
+    number-sections: true
+    colorlinks: true
+    keep-tex: false
+    fig-pos: 'H'
+execute: 
+  echo: false
+  warning: false
+---
+\pagebreak
+    """
+    settings =r"""
+```{python}
+import os
+import sys
+import geopandas as gpd
+import numpy as np
+import pandas as pd
+import re
+import matplotlib.pyplot as plt
+
+# Custom modules
+from justclust.data.bolzano import cols, conv, read_data, selected_col
+from justclust.paths.bolzano import clfile
+import justclust.quarto as qmd
+
+plt.rcParams['figure.dpi'] = 600
+plt.rcParams['savefig.dpi'] = 600
+
+```
+    """
+    data_load = f"""
+```{graph_open}python{graph_close}
+
+# Get selected data with the required cluster
+raw = read_data()
+selected = raw.loc[:, cols].dropna()
+selected.rename(columns=conv, inplace=True)
+
+# Remove error observation
+# mask = selected['Carbon emission building [ton CO2/m²]'] > 8
+# selected = selected[~ mask]
+
+clname = '{clname}'
+
+geo_out = gpd.read_file(clfile)
+geo_out.set_index('SEZ', drop=False, verify_integrity=True, inplace=True)
+clsts = geo_out.loc[selected.index,[clname, 'geometry']]
+selected = clsts.join(selected, validate='1:1')
+
+# Subset fo columns
+selected_col_aq = selected_col.get('col_aq')
+selected_col_carbon = selected_col.get('col_carbon')
+selected_col_char = selected_col.get('col_char')
+selected_col_socio = selected_col.get('col_socio')
+
+selected_col_all = []
+for x in [selected_col_aq, selected_col_carbon, selected_col_char, selected_col_socio]:
+    if x is not None:
+        selected_col_all = selected_col_all + x
+
+
+# Clustering
+cluster_labels = qmd.get_cluster_labels(selected[clname])
+cluster_colors = qmd.cluster_color_dict(cluster_labels)
+selected['cluster_lab'] = selected[clname].map(cluster_labels)
+
+data_cluster = selected[selected['cluster_lab'] != 'Z'][['cluster_lab'] + selected_col_all]
+list_clust =sorted(set(data_cluster['cluster_lab']))
+
+# Compute cluster statistics
+dict_stat = qmd.get_dict_stat(
+    data=selected,
+    selected_col=selected_col_all,
+    var_cluster = 'cluster_lab',
+    exclude_Z = True
+)
+
+# Compute feature importance
+X_fi = data_cluster[selected_col_all]
+y_fi = data_cluster['cluster_lab']
+
+data_fi = qmd.get_clst_feature_importance(X_fi, y_fi, max_depth=None)
+
+# Orther info
+n_units_total = geo_out.shape[0]
+n_units = selected.shape[0]
+
+# summary cluster frequencyz
+summary_cls_frequency = qmd.get_summary_clst_frequency(selected)
+
+```
+    """
+
+    descriptive_intro = f"""
+# Descriptive Statistics
+
+```{graph_open}python{graph_close}
+#| output: asis
+print(
+    f'The urban area of {city} is formed by {graph_open}n_units_total{graph_close} territorial units.',
+    'However, some territorial units are not included in the analysis due to the presence of missing data.',
+    f'The analysis covers {graph_open}n_units{graph_close} territorial units ({graph_open}n_units/n_units_total*100:.2f{graph_close}% of the total; see @fig-territtorial-units).')
+```
+
+```{graph_open}python{graph_close}
+#| fig-align: center
+#| fig-pos: 'H'
+#| label: fig-territtorial-units
+#| fig-cap: 'Territorial units included or excluded from the analysis' 
+
+qmd.plot_map_included(
+    geo_out=geo_out,
+    clname=clname,
+    figsize=(9,9)
+)
+```
+
+    """
+    descriptive_feature = f"""
+## Features
+
+In the next sections, descriptive statistics of the territorial units charateristics are presented.
+
+###  Air Quality
+
+Values regarding air quality are presented in @fig-air-quality and summarized in Table \ref{graph_open}tbl-air-quality{graph_close}.
+```{graph_open}python{graph_close}
+#| output: asis
+#| warning: false
+
+print(
+    qmd.table_latex(
+        data = qmd.summary_stat(selected[selected_col_aq]),
+        caption = 'Summary statistics air quality',
+        label = 'tbl-air-quality'
+    )
+)
+```
+```{graph_open}python{graph_close}
+#| fig-align: center
+#| fig-pos: 'H'
+#| label: fig-air-quality
+#| fig-cap: 'Distribution air quality'
+
+qmd.plot_summary(selected[selected_col_aq])
+```
+
+###  Carbon
+
+Values regarding carbon emission and absorption are presented in @fig-carbon and summarized in Table \ref{graph_open}tbl-carbon{graph_close}.
+```{graph_open}python{graph_close}
+#| output: asis
+#| warning: false
+
+print(
+    qmd.table_latex(
+        data = qmd.summary_stat(selected[selected_col_carbon]),
+        caption = 'Summary statistics carbon',
+        label = 'tbl-carbon'
+    )
+)
+
+```
+```{graph_open}python{graph_close}
+#| fig-align: center
+#| fig-pos: 'H'
+#| label: fig-carbon
+#| fig-cap: 'Distribution carbon'
+
+qmd.plot_summary(
+    selected[selected_col_carbon],
+    title_size=7)
+```
+
+###  Unit Charateristics
+
+Values regarding other territorial unit charateristics are presented in @fig-unit and summarized in Table \ref{graph_open}tbl-unit{graph_close}.
+```{graph_open}python{graph_close}
+#| output: asis
+#| warning: false
+
+print(
+    qmd.table_latex(
+        data = qmd.summary_stat(selected[selected_col_char]),
+        caption = 'Summary statistics unit charateristics',
+        label = 'tbl-unit',
+        size='small'
+    )
+)
+```
+```{graph_open}python{graph_close}
+#| fig-align: center
+#| fig-pos: 'H'
+#| label: fig-unit
+#| fig-cap: 'Distribution unit charateristics'
+
+qmd.plot_summary(selected[selected_col_char])
+```
+
+### Socio Demographic
+Values regarding socio demographic charateristics are presented in @fig-socio and summarized in Table \ref{graph_open}tbl-socio{graph_close}.
+
+```{graph_open}python{graph_close}
+#| output: asis
+#| warning: false
+
+print(
+    qmd.table_latex(
+        data = qmd.summary_stat(selected[selected_col_socio]),
+        caption = 'Summary statistics socio demographic',
+    label = 'tbl-socio'
+    )
+)
+
+```
+```{graph_open}python{graph_close}
+#| fig-align: center
+#| fig-pos: 'H'
+#| label: fig-socio
+#| fig-cap: 'Distribution socio demographic'
+
+qmd.plot_summary(selected[selected_col_socio])
+```
+
+    """
+    descriptive_overview = f"""
+## Overview
+
+The correlation between features is reportd in @fig-corr.
+
+
+
+```{graph_open}python{graph_close}
+#| fig-align: center
+#| fig-pos: 'H'
+#| label: fig-corr
+#| fig-cap: 'Corrrelation between features'
+
+qmd.plot_heatmap_corr(
+    data_plot = selected[selected_col_all],
+    val_filter = 0.3,
+    figsize = (15, 12),
+)
+```
+
+\clearpage
+
+    """
+    
+    cluster_intro = f"""
+# Cluster Analysis 
+
+In this section, we provide all the detail regarding the cluster analysis. In @sec-cluster-approach, the adopted clustering approach is described. In @sec-cluster-results, the results of the clustering are presented. Finally, in @sec-cluster-profiles, a profile with summary information is provided for each cluster.
+
+## Cluster Approach {graph_open}#sec-cluster-approach{graph_close}
+
+To conduct the clustering analysis, we explored different approaches for the data preprocessing and the differen clustering algorithm settings. 
+
+### Preprocessing
+
+Before applying cluster algorithms, data must be preprocessed. This process includ the following (possible) steps:
+
+1. **Scaling the data**. Scaling the data is important for clustering because it helps to ensure that all variables in the data set are on the same scale. In this way we avoid that the unit of measure of different features affects the clusteriing results. We considered the following scaler possibilities from `sklearn.preprocessing`:
+
+    - `None`
+    - `MinMaxScaler()`
+    - `MaxAbsScaler()`
+    - `StandardScaler()`
+    - `Normalizer()` with different values of `norm = ['l1', 'l2', 'max']` 
+    - `QuantileTransformer()` conosiderig both  `'uniform'` and` 'normal'` output distirbutions
+    - `PowerTransformer()`
+    - `RobustScaler()` with different quantile ranges `(0, 100)`, `(1, 99)`, `(2, 98)`, `(5, 95)`, `(10, 90)`, `(20, 80)`, and `(25, 75)`
+
+1. **Weighted features**. Weighting features allows to avoid that categories with multiple indicators overweigths other categories with less indicators. We considerd both the effect of weighting on not the features.
+
+1. **Decomposing techniques**. Decomposing techniques are used to define a new set of features that could properly represent the underlying structure of the data and facilitate the clusterability. We considered the following decomposition possibilities from `sklearn.decomposition`:
+
+    - `None`
+    - `PCA()` with different possible settings
+    - `FactorAnalysis()` 
+    - `FastICA()` with different possible settings
+    - `DictionaryLearning()` with different possible settings
+
+We elauated all the  combinations of the preprocessing steps. To select the preferred data preprocessing, we evaluated the *data clusterability* of the different procedures. That is the tendency of data to form clusters after the application of the preprocessing. To quantify the data clusterability we used the following metric: 
+
+- **Hopkins statistic**. This metric measures the cluster tendency of a data set by comparing the observed data distribution with a random uniform distribution. 
+
+The selected preprocessing method was:
+
+
+> Scaling: [TODO]
+>
+> Weighted features: [TODO]
+>
+> Decomposing: [TODO]
+
+
+### Clustering Algorithm
+
+To compute the clusters, we used the `HDBSCAN` (Hierarchical Density-based Spatial Clustering of Applications with Noise) algorithm. HDBSCAN  combines aspects of DBSCAN and hierarchical clustering. For an introduction on HDBSCAN, see the [official documentation](https://hdbscan.readthedocs.io/en/latest/how_hdbscan_works.html#:~:text=HDBSCAN%20is%20a%20clustering%20algorithm,in%20the%20stability%20of%20clusters.). Summarizing, advantages of HDBSCAN are:
+
+- Handles irregularly shaped and sized clusters. 
+- Robust to outliers. 
+- Does not require the number of clusters to be specified.
+
+We considered different parameter settings. In particular, we tried different values for
+
+- `min_cluster_size`, the minimum size of observations to be considered as a clusters.  Values considered `[2, 3, 5, 7, 10, 12, 15]`.
+- `min_samples`, the number of samples in a neighbourhood for a point to be considered a core point. Values considered `[None, 1, 2, 3, 5]`.
+- `metric`, the metric to use when calculating distance. Metric considered: euclidean, manhattan, braycurtis, canberra, chebyshev, correlation, minkowski, and sqeuclidean.
+- `cluster_selection_method`, the method used to select clusters from the condensed tree. Values considered `["eom", "leaf"]`.
+
+To evaluate the quality of the obtained clusters using all combinations of the different parameters settings, we considere the following metrics: 
+
+- `% covered by clusters`, percentage of observatiosn that have been assigned to a cluster
+- `silhouette`, measure of how similar an object is to its own cluster (cohesion) compared to other clusters (separation). The silhouette ranges from −1 to +1, where a high value indicates that the object is well matched to its own cluster and poorly matched to neighboring clusters. 
+- `davies bouldin`, the average similarity measure of each cluster with its most similar cluster, where similarity is the ratio of within-cluster distances to between-cluster distances. Thus, clusters which are farther apart and less dispersed will result in a better score. The minimum score is zero, with lower values indicating better clustering.
+- `calinski harabasz`, the score is defined as ratio of the sum of between-cluster dispersion and of within-cluster dispersion. Higher values indicates better clustering.
+
+
+The selected clusters were obtaied using the following settings:
+
+> `min_cluster_size = [TODO]`
+>
+> `min_samples = [TODO]`
+>
+> `metric = [TODO]
+>
+> `cluster_selection_method = [TODO]`
+ 
+\clearpage 
+
+    """
+    cluster_results = f"""
+## Cluster Results {graph_open}#sec-cluster-results{graph_close}
+
+In this section, the results of the cluster analysis are presented. In @fig-cluster-map-all, the territorial units are colored according to their assigned cluster. Note that `Z` is used to indicate the *outliers* that were not assigned to any cluster adn it should not be considered a cluster group per se. 
+
+```{graph_open}python{graph_close}
+#| fig-align: center
+#| fig-pos: 'H'
+#| label: fig-cluster-map-all
+#| fig-cap: 'Map cluster results'
+
+qmd.plot_map_all_cluster(
+    data_plot=selected,
+    cluster_colors=cluster_colors,
+    figsize = (12, 12)
+)
+```
+
+Number of territorial units for each cluster are presented in Table \ref{graph_open}tbl-clst-frequencies{graph_close} and @fig-cluster-cardinality.
+
+```{graph_open}python{graph_close}
+#| output: asis
+#| warning: false
+
+summary_cls_frequency['Percent'] = summary_cls_frequency['Percent']\
+    .map(lambda x: f'{graph_open}round(x*100, 2){graph_close}%')
+
+print(
+    qmd.table_latex(
+        data = summary_cls_frequency.T,
+        caption = 'Cluster frequencies',
+        label = 'tbl-clst-frequencies',
+        size = 'small'
+    )
+)
+```
+
+```{graph_open}python{graph_close}
+#| fig-align: center
+#| fig-pos: 'H'
+#| out-width: 50%
+#| label: fig-cluster-cardinality
+#| fig-cap: 'Cluster cardinality' 
+
+qmd.plot_clst_frequency(
+    data = selected,
+    cluster_colors = cluster_colors,
+    var_cluster = 'cluster_lab',
+    figsize = (10, 3.5)
+)
+```
+
+### Cluster Statistics
+
+In this section, descriptive statistics of clusters charateristics are presented. For each variable, the distribution of values for each cluster is presented using box-plots and the feature importance is indicated as percentage on the right for each cluster. This indicates the relevance that each variable has when classifing an observation as belonging to that specific cluster.
+
+
+- **Air Quality**
+
+```{graph_open}python{graph_close}
+#| fig-align: center
+#| fig-pos: 'H'
+#| label: fig-clst-air-quality
+#| fig-cap: 'Distribution air quality per cluster (red solid vertical line - Avg median; red dashed vertical line - Avg mean )'
+qmd.plot_clst_boxplot_loop(
+    data_plot = data_cluster,
+    data_plot_means = dict_stat.get('cluster_means'),
+    data_plot_median = dict_stat.get('cluster_medians'),
+    data_fi = data_fi,
+    list_vars = selected_col_aq,
+    list_clust = list_clust,
+    cluster_colors = cluster_colors,
+    figsize = (10,6.5)
+)
+```
+
+- **Carbon**
+
+```{graph_open}python{graph_close}
+#| fig-align: center
+#| fig-pos: 'H'
+#| label: fig-clst-carbono
+#| fig-cap: 'Distribution carbon per cluster (red solid vertical line - Avg median; red dashed vertical line - Avg mean )'
+qmd.plot_clst_boxplot_loop(
+    data_plot = data_cluster,
+    data_plot_means = dict_stat.get('cluster_means'),
+    data_plot_median = dict_stat.get('cluster_medians'),
+    data_fi = data_fi,
+    list_vars = selected_col_carbon,
+    list_clust = list_clust,
+    cluster_colors = cluster_colors,
+    figsize = (10, 3.5)
+)
+```
+
+- **Unit Charaterisics**
+
+```{graph_open}python{graph_close}
+#| fig-align: center
+#| fig-pos: 'H'
+#| label: fig-clst-unit
+#| fig-cap: 'Distribution unit charateristics per cluster (red solid vertical line - Avg median; red dashed vertical line - Avg mean )'
+qmd.plot_clst_boxplot_loop(
+    data_plot = data_cluster,
+    data_plot_means = dict_stat.get('cluster_means'),
+    data_plot_median = dict_stat.get('cluster_medians'),
+    data_fi = data_fi,
+    list_vars = selected_col_char,
+    list_clust = list_clust,
+    cluster_colors = cluster_colors,
+    figsize = (10,7)
+)
+```
+
+- **Socio Demographic**
+
+```{graph_open}python{graph_close}
+#| fig-align: center
+#| fig-pos: 'H'
+#| label: fig-clst-socio
+#| fig-cap: 'Distribution socio demographic per cluster (red solid vertical line - Avg median; red dashed vertical line - Avg mean )'
+qmd.plot_clst_boxplot_loop(
+    data_plot = data_cluster,
+    data_plot_means = dict_stat.get('cluster_means'),
+    data_plot_median = dict_stat.get('cluster_medians'),
+    data_fi = data_fi,
+    list_vars = selected_col_socio,
+    list_clust = list_clust,
+    cluster_colors = cluster_colors,
+    figsize = (10,7)
+)
+```
+
+### Overwie Clusters
+
+Too provide an overwie of the cluster charateristics, we consider cluster mean and median values. In @fig-mean-heatmap, cluster mean values are reprted for each variable and colored according to the difference with the average mean among clusters (`Avg mean`). In @fig-cluster-comp-mean, the cluster mean value is standardized with resepct to the average mean among clusters for each variable). That is,
+
+$$
+\\frac{graph_open}\\bar{graph_open}x{graph_close}_i  - \\bar{graph_open}x{graph_close}_{graph_open}Avg{graph_close}{graph_close}{graph_open}\\bar{graph_open}x{graph_close}_{graph_open}Avg{graph_close}{graph_close}
+$$
+
+where $\\bar{graph_open}x{graph_close}_i$ is the mean value of a given variable for cluster $i$ and $\\bar{graph_open}x{graph_close}_{graph_open}Avg{graph_close}$ is the average mean among clusters (i.e, $\\frac{graph_open}\\sum_i{graph_open}\\bar{graph_open}x{graph_close}_i{graph_close}{graph_close}{graph_open}\\#i{graph_close}$).
+
+```{graph_open}python{graph_close}
+#| fig-align: center
+#| fig-pos: 'H'
+#| label: fig-mean-heatmap
+#| fig-cap: 'Cluster feature mean values'
+#| 
+data_means = dict_stat['cluster_means']
+data_means['Avg mean'] = data_means.mean(axis = 1)
+
+qmd.plot_clst_heatmap(
+    data = data_means,
+    col_avg= 'Avg mean',
+    fig_size= (10,6)
+)
+```
+
+```{graph_open}python{graph_close}
+#| fig-width: 10
+#| fig-height: 12
+#| fig-align: center
+#| fig-pos: 'H'
+#| label: fig-cluster-comp-mean
+#| fig-cap: 'Cluster comparison Mean (the value is the statndardized cluster mean compared to othe Avg mean among clusters). Feature importance is repoorted on the right.' 
+
+
+qmd.plot_cluster_comp_loop(
+    data_plot = dict_stat['cluster_means_cent'],
+    data_fi=data_fi,
+    colors = cluster_colors,
+    ncols = 3,
+    xmin=-2.5,
+    xmax=2.5,
+    xlabel="Value",
+    ylabel=None
+)
+```
+
+In @fig-median-heatmap, cluster median values are reprted for each variable and colored according to the difference with the average median among clusters (`Avg median`). In @fig-cluster-comp-median, the cluster median value is standardized with resepct to the average median among clusters for each variable). That is,
+
+$$
+\\frac{graph_open}\\tilde{graph_open}x{graph_close}_i  - \\tilde{graph_open}x{graph_close}_{graph_open}Avg{graph_close}{graph_close}{graph_open}\\tilde{graph_open}x{graph_close}_{graph_open}Avg{graph_close}{graph_close}
+$$
+
+where $\\tilde{graph_open}x{graph_close}_i$ is the median value of a given variable for cluster $i$ and $\\tilde{graph_open}x{graph_close}_{graph_open}Avg{graph_close}$ is the average median among clusters (i.e, $\\frac{graph_open}\\sum_i{graph_open}\\tilde{graph_open}x{graph_close}_i{graph_close}{graph_close}{graph_open}\\#i{graph_close}$).
+
+
+```{graph_open}python{graph_close}
+#| fig-align: center
+#| fig-pos: 'H'
+#| label: fig-median-heatmap
+#| fig-cap: 'Cluster feature median values'
+#| 
+data_medians = dict_stat['cluster_medians']
+data_medians['Avg median'] = data_medians.mean(axis = 1)
+
+qmd.plot_clst_heatmap(
+    data = data_medians,
+    col_avg= 'Avg median',
+    fig_size= (10,6)
+)
+```
+
+```{graph_open}python{graph_close}
+#| fig-width: 10
+#| fig-height: 12
+#| fig-align: center
+#| fig-pos: 'H'
+#| label: fig-cluster-comp-median
+#| fig-cap: 'Map cluster comparison Median (the value is the statndardized cluster median compared to othe Avg median among clusters). Feature importance is repoorted on the right.' 
+
+data_plot = dict_stat.get('cluster_medians_cent')
+qmd.plot_cluster_comp_loop(
+    data_plot = data_plot,
+    data_fi=data_fi,
+    colors = cluster_colors,
+    ncols = 3,
+    xmin=-2.5,
+    xmax=2.5,
+    xlabel="Value",
+    ylabel=None
+)
+```
+
+\clearpage
+
+    """
+    cluster_profiles = f"""
+## Cluster Profiles {graph_open}#sec-cluster-profiles{graph_close}
+
+```{graph_open}python{graph_close}
+#| output: asis
+#| fig-align: center
+#| fig-pos: 'H'
+#| warning: false
+
+
+for cluster_lab in cluster_labels.values():
+    if(cluster_lab == 'Z'): continue
+
+    mask = selected['cluster_lab'] == cluster_lab
+
+    graph_par_open = "{graph_open}"
+    graph_par_close = "{graph_close}"
+    print(f'\\subsubsection{graph_open}graph_par_open{graph_close}Cluster {graph_open}cluster_lab{graph_close} (n = {graph_open}selected.loc[mask].shape[0]{graph_close}){graph_open}graph_par_close{graph_close}')
+    
+    qmd.plot_map_single_cluster(
+        data_plot=selected,
+        mask = mask,
+        cluster_lab=cluster_lab,
+        cluster_colors=cluster_colors
+    )
+
+    qmd.plot_cluster_comp_single(
+        dict_stat=dict_stat,
+        data_fi = data_fi,
+        cluster_lab=cluster_lab,
+        colors=cluster_colors[cluster_lab],
+        sorted = True
+    )
+
+    qmd.plot_boxplot_comp_cluster_loop(
+        data_cluster = data_cluster,
+        cluster_lab = cluster_lab,
+        list_vars = selected_col_all,
+        data_fi = data_fi,
+        data_plot_means = dict_stat['cluster_means'],
+        data_plot_median = dict_stat['cluster_medians'],
+        cluster_colors = cluster_colors,
+        figsize = (20,25)
+        )
+    
+    print('\\clearpage')
+  
+```
+
+```{graph_open}python{graph_close}
+#| output: asis
+#| fig-align: center
+#| fig-pos: 'H'
+#| warning: false
+cluster_lab == 'Z'
+
+mask = selected['cluster_lab'] == cluster_lab
+print(f'\\subsubsection{graph_open}graph_par_open{graph_close}Cluster {graph_open}cluster_lab{graph_close} (n = {graph_open}selected.loc[mask].shape[0]{graph_close}){graph_open}graph_par_close{graph_close}')
+    
+qmd.plot_map_single_cluster(
+        data_plot=selected,
+        mask = mask,
+        cluster_lab=cluster_lab,
+        cluster_colors=cluster_colors
+    )
+
+```
+```{graph_open}python{graph_close}
+#| fig-align: center
+#| fig-pos: 'H'
+#| fig-width: 20
+#| fig-height: 25
+#| warning: false
+qmd.plot_summary(
+    data = selected.loc[mask, selected_col_all],
+    color = cluster_colors[cluster_lab]
+)
+
+```
+
+\clearpage
+
+
+    """
+    
+    file = open(cdir / f'report-{city}-template.qmd', 'w') 
+    
+    
+    file.write(header + settings + data_load) 
+    file.write(descriptive_intro + descriptive_feature + descriptive_overview) 
+    file.write(cluster_intro + cluster_results + cluster_profiles) 
+                
+    file.close() 
-- 
GitLab