
import sys
fo = sys.stdout

import numbers
from functools import singledispatch
import collections

import os.path
from warnings import warn

import imageio_embedded as imageio

import numpy as np
from matplotlib import pyplot as plt, gridspec
from matplotlib.figure import Figure
import matplotlib.patches as mpatches
from matplotlib.ticker import MaxNLocator

from geogst.core.geometries.grids.rasters import *
from geogst.core.geometries.points import *
from geogst.core.profiles.geoprofiles import *
from geogst.core.profiles.profilers import *
from geogst.core.profiles.profiletraces import *
from geogst.core.profiles.profiletraces import *
from geogst.plots.maps import *
from geogst.plots.parameters import *


topoprofile_default_linecolor = "peru"
INTERSECTION_LINES_LABEL_FONT_SIZE = 7
INTERSECTION_LINES_LABEL_X_OFFSET = 5
INTERSECTION_LINES_LABEL_Y_OFFSET = 45


def plot_points_as_pointtraces(
    ax,
    points_projections: Dict,
    **kwargs,
):

    point_projections_params = kwargs.get("point_projections_params", None)

    if point_projections_params is None:
        point_projections_params = PointPlotParams()

    projected_ids = []
    projected_s = []
    projected_z = []
    projected_dist = []

    for rec_id, point_projection in points_projections.items():
        projected_ids.append(rec_id)
        projected_s.append(point_projection.s)
        projected_z.append(point_projection.z)
        projected_dist.append(point_projection.dist)

    ax.plot(
        projected_s,
        projected_z,
        marker=point_projections_params.marker,
        color=point_projections_params.color,
        markersize=point_projections_params.markersize,
        alpha=point_projections_params.alpha,
        linestyle='None',
    )

    return


def ottieni_rapporto_aspetto(ax=None):
    """
    Phind, 2025-03-23

    Calcola il rapporto di aspetto effettivo dell'asse.

    Parametri:
        ax: matplotlib.axes.Axes (opzionale)
            L'asse da cui calcolare il rapporto. Se None, usa l'asse corrente.

    Restituisce:
        float: Il rapporto di aspetto numerico
    """
    if ax is None:
        ax = plt.gca()

    # Ottieni dimensioni figura e posizione dell'asse
    fig = ax.get_figure()
    fig_w, fig_h = fig.get_size_inches()
    _, _, w, h = ax.get_position().bounds

    # Calcola il rapporto delle unità di display
    disp_ratio = (fig_h * h) / (fig_w * w)

    # Calcola il rapporto dei dati
    data_ratio = (ax.get_ylim()[1] - ax.get_ylim()[0]) / \
                 (ax.get_xlim()[1] - ax.get_xlim()[0])

    return disp_ratio / data_ratio


def plot_point_projections_from_dicts(
    ax,
    aspect,
    points_projections: Dict,
    section_length: numbers.Real,
    **kwargs,
):

    for category, dataset in points_projections.items():

        # parameters

        points_projections_params = kwargs.pop(category, GenericPlotParams())

        if points_projections_params.plot_type == PointLikePlotTypes.POINTS:

            projected_categories = []
            projected_s = []
            projected_z = []
            projected_dist = []

            for rec_id, point_projection in dataset.items():
                projected_categories.append(rec_id)
                projected_s.append(point_projection.s)
                projected_z.append(point_projection.z)
                projected_dist.append(point_projection.dist)

            color = points_projections_params.color

            ax.plot(
                projected_s,
                projected_z,
                marker=points_projections_params.marker,
                color=color,
                markersize=points_projections_params.markersize,
                alpha=points_projections_params.alpha,
                linestyle='None',
            )

            if points_projections_params.labels:

                for category, s, z in zip(
                    projected_categories,
                    projected_s,
                    projected_z):

                    label = f"{category}"

                    ax.annotate(
                        label,
                        (s + 75, z + 75),
                        color=color,
                        fontsize=6,
                    )

        elif points_projections_params.plot_type == PointLikePlotTypes.ATTITUDES:

            projected_categories = []
            projected_s = []
            projected_z = []
            src_dip_dirs = []
            src_dip_angs = []

            for rec_id, profile_attitude in dataset.items():
                projected_categories.append(rec_id)
                projected_s.append(profile_attitude.s)
                projected_z.append(profile_attitude.z)
                src_dip_dirs.append(profile_attitude.src_dip_dir)
                src_dip_angs.append(profile_attitude.src_dip_ang)

            vertical_exaggeration = aspect

            ax.plot(
                projected_s,
                projected_z,
                marker=points_projections_params.marker,
                color=points_projections_params.color,
                markersize=points_projections_params.markersize,
                alpha=points_projections_params.alpha,
                linestyle='None',
            )

            # plot segments representing structural data

            for _, structural_attitude in dataset.items():
                structural_segment_s, structural_segment_z = structural_attitude.create_segment_for_plot(
                    section_length,
                    vertical_exaggeration)

                ax.plot(
                    structural_segment_s,
                    structural_segment_z,
                    '-',
                    color=points_projections_params.color,
                    alpha=points_projections_params.alpha,
                )

            if points_projections_params.label_orientations or points_projections_params.label_ids:

                for rec_id, src_dip_dir, src_dip_ang, s, z in zip(
                        projected_categories,
                        src_dip_dirs,
                        src_dip_angs,
                        projected_s,
                        projected_z):

                    if points_projections_params.label_orientations and points_projections_params.label_ids:
                        label = f"{rec_id}-{src_dip_dir:05.01F}/{src_dip_ang:04.01F}"
                    elif points_projections_params.label_orientations:
                        label = f"{src_dip_dir:05.01F}/{src_dip_ang:04.01F}"
                    else:
                        label = f"{rec_id}"

                    ax.annotate(
                        label,
                        (s + 15, z + 15),
                        size=7,  # Dimensione del testo
                        color=points_projections_params.color,  # Colore del testo
                    )

        else:
            raise NotImplementedError


def plot_line_intersections(
    ax,
    topo_profile,
    line_intersections: list,
    line_intersections_style,
    **kwargs,
):

    parsed_s_vals, parsed_z_vals, parsed_labels = [], [], []

    for ndx_profile_line, profile_line_intersection_elements in enumerate(line_intersections, start=1):

        for ndx_intersection_element, intersected_element in enumerate(profile_line_intersection_elements, start=1):

            intersected_element_id = intersected_element.id
            intersected_element_parts = intersected_element.arrays

            for s_range in intersected_element_parts:

                s_start = s_range[0]
                s_end = s_range[1] if len(s_range) > 1 else None
                plot_symbol = '-' + line_intersections_style.marker if len(
                    s_range) > 1 else line_intersections_style.marker

                s_vals = topo_profile.s_subset(
                    s_start,
                    s_end
                )

                z_vals = topo_profile.zs_from_s_range(
                    s_start,
                    s_end
                )

                if s_vals is None or z_vals is None:
                    continue

                if len(s_vals) != len(z_vals):
                    print(
                        f"Error with numerosity of line intersection data to plot for profile {ndx_profile_line}: s values are {len(s_vals)} while z_values are {len(z_vals)}")
                    continue

                for s_val, z_val in zip(s_vals, z_vals):
                    if not isinstance(s_val, (numbers.Real, numbers.Integral)) or not isinstance(z_val, (
                    numbers.Real, numbers.Integral)):
                        print(
                            f"Discarding s-z couple for profile {ndx_profile_line} line intersections plot: s -> {s_val} z -> {z_val}")
                        continue
                    if not isfinite(s_val) or not isfinite(z_val):
                        print(
                            f"Discarding s-z couple for profile {ndx_profile_line} line intersections plot: s -> {s_val} z -> {z_val}")
                        continue
                    parsed_s_vals.append(s_val)
                    parsed_z_vals.append(z_val)
                    parsed_labels.append(intersected_element_id)

    if parsed_s_vals and parsed_z_vals:
        
        ax.plot(
            parsed_s_vals,
            parsed_z_vals,
            plot_symbol,
            color=line_intersections_style.color,
            markersize=line_intersections_style.markersize,
            alpha=line_intersections_style.alpha,
            linestyle='None',
        )

        if line_intersections_style.labels:

            for s_val, z_val, label in zip(parsed_s_vals, parsed_z_vals, parsed_labels):

                ax.annotate(
                    f"{label}",
                    (s_val+INTERSECTION_LINES_LABEL_X_OFFSET, z_val + INTERSECTION_LINES_LABEL_Y_OFFSET),
                    fontsize=INTERSECTION_LINES_LABEL_FONT_SIZE,
                    color=line_intersections_style.color,
                    alpha=line_intersections_style.alpha,
                )


def plot_points_projections(
    ax,
    aspect: numbers.Real,
    points_projections: Dict,
    section_length: numbers.Real,
    **kwargs,
):

    if not isinstance(points_projections, dict):
        print(
            f"Warning: Points projections are not of type 'Dict' but {type(points_projections)}")
        return

    values = points_projections.values()
    for value in values:
        print(f"Debug: point value: {type(value)} -> {value}")
    value_types = list(set(map(lambda value: type(value), values)))
    if len(value_types) != 1:
        print(
            f"Warning: Points projections values must be of a single type but {len(value_types)} types found"
        )
        return

    point_projection_type = value_types[0]

    if point_projection_type == PointTrace:

        plot_points_as_pointtraces(
            ax,
            points_projections,
            **kwargs,
        )

    elif point_projection_type == collections.defaultdict:

        plot_point_projections_from_dicts(
            ax,
            aspect,
            points_projections,
            section_length,
            **kwargs,
        )

    else:

        print(f"Got non-managed point projection type: {point_projection_type}")

    return


def plot_polygons_intersections(
    ax,
    topo_profile,
    polygons_intersections,
    polygon_intersections_colors,
    polygon_intersections_linewidth,
    polygon_inters_marker_label,
    z_min,
    z_max,
    polygon_inters_legend_on,
    **kwargs
):

    for ndx_profile_line, profile_intersection_elements in enumerate(polygons_intersections):

        for polygon_intersection_element in profile_intersection_elements:

            polygon_intersection_id = polygon_intersection_element.id
            polygon_intersection_subparts = polygon_intersection_element.arrays

            for s_range in polygon_intersection_subparts:

                s_start = s_range[0]
                s_end = s_range[1] if len(s_range) > 1 else None
                s_mid = s_start if s_end is None else (s_start + s_end) / 2

                plot_symbol = '-' if len(s_range) > 1 else 'o'

                s_vals = topo_profile.s_subset(
                    s_start,
                    s_end
                )

                if s_vals is None:
                    continue

                z_vals = topo_profile.zs_from_s_range(
                    s_start,
                    s_end
                )

                if z_vals is None:
                    continue

                ax.plot(
                    s_vals,
                    z_vals,
                    plot_symbol,
                    color=polygon_intersections_colors[polygon_intersection_id],
                    linewidth=polygon_intersections_linewidth
                )

                '''
                if polygon_inters_marker_label:
                    ax.annotate(
                        f"{polygon_intersection_id}",
                        (s_mid, z_min + int((z_max - z_min) / 20)),
                        color=polygon_intersections_colors[polygon_intersection_id]
                    )
                '''

    if polygon_inters_legend_on:

        legend_patches = []
        for polygon_code in polygon_intersections_colors:
            legend_patches.append(mpatches.Patch(color=polygon_intersections_colors[polygon_code], label=str(polygon_code)))

        # from: https://stackoverflow.com/questions/4700614/how-to-put-the-legend-out-of-the-plot
        # Shrink current axis by 20%
        box = ax.get_position()
        ax.set_position([box.x0, box.y0, box.width * 0.8, box.height])

        # Put a legend to the right of the current axis
        ax.legend(
            handles=legend_patches,
            loc='center left',
            bbox_to_anchor=(1, 0.5)
        )


@singledispatch
def profiles(
    obj,
    **kargs
) -> Figure:
    """

    :param obj:
    :param kargs:
    :return:
    """

    fig = kargs.pop("fig", None)
    aspect = kargs.pop("aspect", 1)
    width = kargs.pop("width", FIG_WIDTH_INCHES_DEFAULT)
    height = kargs.pop("height", FIG_HEIGHT_INCHES_DEFAULT)

    if fig is None:

        fig, ax = plt.subplots()
        fig.set_size_inches(width, height)

        ax.set_aspect(aspect)

    else:

        plt.gca()

    return fig


@profiles.register(ZTrace)
def _(
    ztrace: ZTrace,
    **kargs
) -> Figure:

    axis_parameters = kargs.pop("axis_parameters", None)
    if axis_parameters is not None:
        z_min = axis_parameters.z_min
        z_max = axis_parameters.z_max
        aspect = axis_parameters.vertical_exaggeration
        grid = axis_parameters.grid
        grid_color = axis_parameters.grid_color
        grid_linestyle = axis_parameters.grid_linestyle
        grid_width = axis_parameters.grid_linewidth
        breaklines = axis_parameters.breaklines
        breaklines_color = axis_parameters.breaklines_color
        breaklines_width = axis_parameters.breaklines_width
        breaklines_style = axis_parameters.breaklines_style
    else:
        z_min = None
        z_max = None
        aspect = None
        grid = False
        grid_color = 'tan'
        grid_linestyle = None
        grid_width = 0.2
        breaklines = True
        breaklines_color = breakline_color
        breaklines_width = breakline_width
        breaklines_style = 'dotted'

    fig = kargs.pop("fig", None)
    width = kargs.pop("width", FIG_WIDTH_INCHES_DEFAULT)
    height = kargs.pop("height", FIG_HEIGHT_INCHES_DEFAULT)

    file_path = kargs.pop("file_path", None)

    if z_min is None or z_max is None:
        z_range = ztrace.z_max() - ztrace.z_min()
        z_min = ztrace.z_min() - FIG_Z_PADDING_DEFAULT * z_range
        z_max = ztrace.z_max() + FIG_Z_PADDING_DEFAULT * z_range

    if np.isnan(z_min) or np.isnan(z_max):
        return

    if fig is None:

        fig = plt.figure()
        fig.set_size_inches(width, height)

    ax = fig.add_subplot()

    if aspect is not None:
        ax.set_aspect(aspect)

    if z_min is not None or z_max is not None:
        ax.set_ylim([z_min, z_max])

    if grid:
        ax.grid(
            True,
            linestyle='-',
            color=grid_color,
            linewidth=grid_width)

    ax.plot(
        ztrace.s_arr(),
        ztrace.z_arr(),
        **kargs
    )

    if breaklines:
        bottom, top = ax.get_ylim()
        ax.vlines(
            ztrace.s_breaks(),
            bottom,
            top,
            color=breaklines_color,
            linewidth=breaklines_width,
            linestyles=breaklines_style
        )

    fig.suptitle(f"{ztrace.rec_id()}")

    if file_path is not None:
        plt.savefig(file_path)

    return fig


@profiles.register(list)
def _(
    z_profiles: List[ZTrace],
    **kargs
) -> Figure:

    axis_parameters = kargs.pop("axis_params", None)

    if axis_parameters is not None:

        z_min = axis_parameters.z_min
        z_max = axis_parameters.z_max
        aspect = axis_parameters.vertical_exaggeration
        grid = axis_parameters.grid
        grid_color = axis_parameters.grid_color
        grid_linestyle = axis_parameters.grid_linestyle
        grid_width = axis_parameters.grid_linewidth
        breaklines = axis_parameters.breaklines
        breaklines_color = axis_parameters.breaklines_color
        breaklines_width = axis_parameters.breaklines_width
        breaklines_style = axis_parameters.breaklines_style

    else:

        z_min = None
        z_max = None
        aspect = None
        grid = False
        grid_color = 'tan'
        grid_linestyle = None
        grid_width = 0.2
        breaklines = True
        breaklines_color = breakline_color
        breaklines_width = breakline_width
        breaklines_style = 'dotted'


    width = kargs.pop("width", FIG_WIDTH_INCHES_DEFAULT)
    height = kargs.pop("height", FIG_HEIGHT_INCHES_DEFAULT)

    elevation_params = kargs.pop("elevation_params", None)
    points = kargs.pop("points", None)
    attitudes = kargs.pop("attitudes", None)
    line_intersections = kargs.pop("line_intersections", None)
    polygon_intersections = kargs.pop("polygon_intersections", None)

    single_plot = kargs.pop("single_plot", False)
    file_path = kargs.pop("file_path", None)

    #

    num_profiles = len(z_profiles)

    if z_min is None or z_max is None:

        max_of_zs = np.nanmax([ztrace.z_max() for ztrace in z_profiles])
        min_of_zs = np.nanmin([ztrace.z_min() for ztrace in z_profiles])
        z_range = max_of_zs - min_of_zs
        z_min = min_of_zs - FIG_Z_PADDING_DEFAULT * z_range
        z_max = max_of_zs + FIG_Z_PADDING_DEFAULT * z_range

    if np.isnan(z_min) or np.isnan(z_max):
        return

    if not single_plot:

        single_relative_heights = 1.0 / (num_profiles + 2)
        fig = plt.figure(figsize=(width, height * num_profiles * 1.2))

    else:

        fig = plt.figure(figsize=(width, height))

    if aspect is None:
        aspect = 1

    if not single_plot:

        fig, axs = plt.subplots(
            num_profiles,
            ncols=1,
            sharex=True,
        )

        for ndx in range(num_profiles):

            z_prof = z_profiles[ndx]

            if not isinstance(axs, np.ndarray):
                ax = axs
            else:
                ax = axs[ndx]

            if grid:
                ax.grid(
                    True,
                    linestyle='-',
                    color=grid_color,
                    linewidth=grid_width)

            ax.set_aspect(aspect)

            if z_min is not None or z_max is not None:
                ax.set_ylim([z_min, z_max])

            ax.plot(
                z_prof.s_arr(),
                z_prof.z_arr(),
                color=elevation_params.color,
                ** kargs
            )

            ax.tick_params(axis='y', labelsize=8)
            if ndx < num_profiles - 1:
                ax.tick_params(labelbottom=False)
            else:
                ax.tick_params(axis='x', labelsize=8)

            # Etichetta interna in alto a sinistra

            ax.text(0.02, 0.95, str(z_prof.rec_id()), transform=ax.transAxes,
                    fontsize=7, va='top', ha='left'
            )

    else:

        axs = fig.add_subplot()

        axs.set_aspect(aspect)

        if z_min is not None or z_max is not None:
            axs.set_ylim([z_min, z_max])

        for ndx in range(num_profiles):

            z_prof = z_profiles[ndx]
            axs.plot(
                z_prof.s_arr(),
                z_prof.z_arr(),
                color=elevation_params.color,
                ** kargs
            )

        axs.tick_params(axis='x', labelsize=8)
        axs.tick_params(axis='y', labelsize=8)

    if breaklines:

        if not single_plot:

            for ndx in range(num_profiles):

                z_prof = z_profiles[ndx]

                if not isinstance(axs, np.ndarray):
                    ax = axs
                else:
                    ax = axs[ndx]

                bottom, top = ax.get_ylim()
                ax.vlines(
                    z_prof.s_breaks(),
                    bottom,
                    top,
                    color=breaklines_color,
                    linewidth=breaklines_width,
                    linestyles=breaklines_style
                )

        else:

            bottom, top = axs.get_ylim()
            axs.vlines(
                z_profiles.s_breaks(),
                bottom,
                top,
                color=breaklines_color,
                linewidth=breaklines_width,
                linestyles=breaklines_style
            )

    plt.tight_layout()

    if file_path is not None:
        plt.savefig(file_path)

    return fig


@profiles.register(GeoProfile)
def _(
    geoprofile: GeoProfile,
    **kwargs
) -> Figure:
    """
    Plot a single geological profile.

    :param geoprofile: the geoprofile to plot
    :return: the figure.
    """

    if not geoprofile.has_topography():
        print("Warning: geoprofile has no topography defined")
        return

    # keyword parameters extraction

    fig = kwargs.pop("fig", None)
    ax = kwargs.pop("ax", None)

    width = kwargs.pop("width", None)
    height = kwargs.pop("height", None)

    profile_ndx = kwargs.pop("profile_ndx", 0)
    superposed = kwargs.pop("superposed", TOPOPROF_SUPERPOSED_CHOICE_DEFAULT)

    # axis params

    axis_params = kwargs.pop("axis_params", AxisPlotParams())
    z_min = axis_params.z_min
    z_max = axis_params.z_max
    aspect = axis_params.vertical_exaggeration

    # elevation parameters

    elevation_params = kwargs.pop("elevation_params", ElevationPlotParams())

    # line intersections parameters

    line_intersections_style = kwargs.pop("line_intersections", PointPlotParams())

    # polygons intersections parameters

    polygon_intersections = kwargs.pop("polygon_intersections", None)
    polygon_intersections_linewidth = PLINT_LINE_WIDTH_DEFAULT  #if polygon_intersections is None else polygon_intersections.get("line_width", PLINT_LINE_WIDTH_DEFAULT)
    polygon_intersections_colors = polygon_intersections  #None if polygon_intersections is None else polygon_intersections.get("colors", None)
    polygon_inters_marker_label = PLINT_LABELS_DEFAULT  #if polygon_intersections is None else polygon_intersections.get("linelabels", PLINT_LABELS_DEFAULT)
    polygon_inters_legend_on = PLINT_LEGEND_DEFAULT  #if polygon_intersections is None else polygon_intersections.get("legend", PLINT_LEGEND_DEFAULT)

    # figure definitions

    if fig is None:

        fig = plt.figure()
        fig.set_size_inches(width, height)

    if superposed:

        ax = fig.add_axes(
            [0.1, 0.1, 0.8, 0.8]
        )

    elif ax is None:

        ax = fig.add_subplot()

    # profile parameters

    section_length = geoprofile.length_2d()

    # definition of elevation range

    if z_min is None or z_max is None:
        z_range = geoprofile.z_max() - geoprofile.z_min()
        z_min = geoprofile.z_min() - FIG_Z_PADDING_DEFAULT * z_range
        z_max = geoprofile.z_max() + FIG_Z_PADDING_DEFAULT * z_range

    if np.isnan(z_min) or np.isnan(z_max):
        return

    # general parameters

    ax.set_ylim([z_min, z_max])

    #ax.yaxis.set_major_locator(MaxNLocator(nbins=5))

    yticks = ax.get_yticks()
    ymax = max(yticks)

    # Costruisci le etichette, vuota per il massimo
    yticklabels = ["" if tick == ymax else str(tick) for tick in yticks]

    # Applica i tick e le etichette
    ax.set_yticks(yticks)
    ax.set_yticklabels(yticklabels)

    ax.tick_params(axis='y', labelsize=8)

    ax.set_aspect(aspect)

    # plot of elevation profiles

    if geoprofile._topo_profile:

        if superposed:
            linecolor = LNINT_ADDITIONAL_COLORS[profile_ndx % len(LNINT_ADDITIONAL_COLORS)]
        else:
            if elevation_params.color is None:
                linecolor = topoprofile_default_linecolor
            else:
                linecolor = elevation_params.color

        if axis_params.grid:
            ax.grid(
                True,
                color=axis_params.grid_color,
                linestyle=axis_params.grid_linestyle,
                linewidth=axis_params.grid_linewidth)

        ax.plot(
            geoprofile._topo_profile.s_arr(),
            geoprofile._topo_profile.z_arr(),
            color=linecolor,
            linestyle=elevation_params.linestyle,
            linewidth=elevation_params.width
        )

        if axis_params.breaklines:
            bottom, top = ax.get_ylim()
            ax.vlines(
                geoprofile._topo_profile.s_breaks(),
                bottom,
                top,
                color=axis_params.breaklines_color,
                linewidth=axis_params.breaklines_width,
                linestyles=axis_params.breaklines_style
            )

    # plot of polygons intersections

    if geoprofile._polygons_intersections:

        if not geoprofile._topo_profile:

            print('Warning: topographic profile is not defined, so intersections cannot be plotted')

        elif not polygon_intersections:

            print('Warning: polygon intersection styles are not defined, so intersections cannot be plotted')

        else:

            plot_polygons_intersections(
                ax,
                geoprofile._topo_profile,
                geoprofile._polygons_intersections,
                polygon_intersections_colors,
                polygon_intersections_linewidth,
                polygon_inters_marker_label,
                z_min,
                z_max,
                polygon_inters_legend_on,
                **kwargs
            )

    # plot of line intersections

    if geoprofile._lines_intersections:

        if not geoprofile._topo_profile:

            warn('Topographic profile is not defined, so intersections cannot be plotted')

        else:

            plot_line_intersections(
                ax,
                geoprofile._topo_profile,
                geoprofile._lines_intersections,
                line_intersections_style,
                **kwargs
            )

    # plot of point projections

    if geoprofile._points_projections:

        plot_points_projections(
            ax,
            aspect,
            geoprofile._points_projections,
            section_length,
            **kwargs
        )

    # plot of traces with attitudes intersections

    if geoprofile._lines_intersections_with_attitudes:

        section_length = geoprofile.length_2d()

        projected_ids = []
        projected_s = []
        projected_z = []
        src_dip_dirs = []
        src_dip_angs = []

        for rec_id, profile_attitudes in geoprofile._lines_intersections_with_attitudes.items():

            for profile_attitude in profile_attitudes:
                projected_ids.append(rec_id)
                projected_s.append(profile_attitude.s)
                projected_z.append(geoprofile._topo_profile.z_linear_interpol(profile_attitude.s))
                src_dip_dirs.append(profile_attitude.src_dip_dir)
                src_dip_angs.append(profile_attitude.src_dip_ang)

        axes = fig.gca()
        vertical_exaggeration = axes.get_aspect()

        line_attitudes_projections_params = GenericPlotParams(
            color="yellow",
            markersize=8,
            alpha=0.4,
            label_orientations=False
        )

        axes.plot(
            projected_s,
            projected_z,
            marker="o",  #line_attitudes_projections_params.marker,
            color="yellow",  #qcolor2rgbmpl(line_attitudes_projections_params.color) if isinstance(line_attitudes_projections_params.color, QColor) else line_attitudes_projections_params.color,
            markersize=8,  #line_attitudes_projections_params.markersize,
            alpha=line_attitudes_projections_params.alpha,
            linestyle='None',
        )

        # plot segments representing structural data

        for _, profile_attitudes in geoprofile._lines_intersections_with_attitudes.items():
            for profile_attitude in profile_attitudes:
                structural_segment_s, structural_segment_z = profile_attitude.create_segment_for_plot(
                    section_length,
                    geoprofile._topo_profile.z_linear_interpol(profile_attitude.s),
                    vertical_exaggeration,
                    segment_scale_factor=3)

                fig.gca().plot(
                    structural_segment_s,
                    structural_segment_z,
                    '-',
                    linewidth=10,
                    color=line_attitudes_projections_params.color,
                    alpha=line_attitudes_projections_params.alpha,
                )

        if line_attitudes_projections_params.label_orientations or line_attitudes_projections_params.label_ids:

            for rec_id, src_dip_dir, src_dip_ang, s, z in zip(
                    projected_ids,
                    src_dip_dirs,
                    src_dip_angs,
                    projected_s,
                    projected_z):

                if line_attitudes_projections_params.label_orientations and line_attitudes_projections_params.label_ids:
                    label = f"{rec_id}-{src_dip_dir:05.01F}/{src_dip_ang:04.01F}"
                elif line_attitudes_projections_params.label_orientations:
                    label = f"{src_dip_dir:05.01F}/{src_dip_ang:04.01F}"
                else:
                    label = f"{rec_id}"

                axes.annotate(label, (s + 15, z + 15))

    # Etichetta interna al centro a sinistra

    ax.text(0.02, 0.55, str(geoprofile._topo_profile.rec_id), transform=ax.transAxes,
            fontsize=7, va='top', ha='left'
            )

    return fig


@profiles.register(GeoProfiles)
def _(
    geoprofiles: GeoProfiles,
    **kargs
) -> Figure:
    """
    Plot a set of geological profiles.

    :param geoprofiles: the geoprofiles to plot
    :return: the figures.
    """

    if not geoprofiles.have_topographies():
        print("Geoprofiles have no topographic set defined")
        return

    # keyword parameters extraction

    width = kargs.pop("width", FIG_WIDTH_INCHES_DEFAULT)
    height = kargs.pop("height", FIG_WIDTH_INCHES_DEFAULT / 4) * geoprofiles.num_profiles()

    superposed = kargs.pop("superposed", False)

    # axis parameters

    if "axis_params" in kargs:
        axis_params = kargs.pop("axis_params")
    else:
        axis_params = AxisPlotParams()
        z_range = geoprofiles.z_max() - geoprofiles.z_min()
        z_min = geoprofiles.z_min() - FIG_Z_PADDING_DEFAULT * z_range
        z_max = geoprofiles.z_max() + FIG_Z_PADDING_DEFAULT * z_range
        axis_params.z_min = z_min
        axis_params.z_max = z_max

    '''
    axis_parameters = kargs.pop("axis_params", None)

    if axis_parameters is not None:

        z_min = axis_parameters.z_min
        z_max = axis_parameters.z_max
        aspect = axis_parameters.vertical_exaggeration
        grid = axis_parameters.grid
        grid_color = axis_parameters.grid_color
        grid_linestyle = axis_parameters.grid_linestyle
        grid_width = axis_parameters.grid_linewidth
        breaklines = axis_parameters.breaklines
        breaklines_color = axis_parameters.breaklines_color
        breaklines_width = axis_parameters.breaklines_width
        breaklines_style = axis_parameters.breaklines_style

    else:

        z_min = None
        z_max = None
        aspect = None
        grid = False
        grid_color = 'tan'
        grid_linestyle = None
        grid_width = 0.2
        breaklines = True
        breaklines_color = breakline_color
        breaklines_width = breakline_width
        breaklines_style = 'dotted'
    '''

    # others

    num_profiles = geoprofiles.num_profiles()

    if not superposed:

        fig = plt.figure()
        gs = gridspec.GridSpec(
            ncols=1,
            nrows=num_profiles,
            height_ratios = [1]*num_profiles,
        )
        axes = []

        ax = fig.add_subplot(gs[0])
        if num_profiles > 1:
            #ax.set_xticks([])
            #plt.setp(ax.get_xticks(), visible=False)
            plt.setp(ax.get_xticklabels(), visible=False)
        else:
            ax.tick_params(axis='x', labelsize=8)
        axes.append(ax)

        for ndx in range(1, geoprofiles.num_profiles()):
            ax = fig.add_subplot(gs[ndx], sharex=axes[0])
            if ndx < num_profiles - 1:
                plt.setp(ax.get_xticklabels(), visible=False)
            else:
                ax.tick_params(axis='x', labelsize=8)
            axes.append(ax)

        plt.subplots_adjust(hspace=0)

    else:
        fig = plt.figure()

    fig.set_size_inches(width, height)

    for ndx in range(num_profiles):

        geoprofile = geoprofiles[ndx]

        profiles(
            geoprofile,
            fig=fig,
            ax=axes[ndx],
            profile_ndx=ndx,
            axis_params=axis_params,
            **kargs
        )

    return fig


@profiles.register(MultiGridsProfiles)
def _(
    grid_profiles: MultiGridsProfiles,
    **kargs
): #-> Optional[List[Optional[Figure]]]:
    """
    Plot a set of geological profiles.

    :param geoprofiles: the geoprofiles to plot
    :return: the figures.
    """

    # keyword parameters extraction

    width = kargs.pop("width", FIG_WIDTH_INCHES_DEFAULT)
    height = kargs.pop("height", FIG_WIDTH_INCHES_DEFAULT / 4) * grid_profiles.num_profiles()

    superposed = kargs.pop("superposed", True)

    # axis parameters

    if "axis_params" in kargs:
        axis_params = kargs.pop("axis_params")
        aspect = axis_params.vertical_exaggeration
    else:
        axis_params = AxisPlotParams()
        z_range = grid_profiles.z_max() - grid_profiles.z_min()
        z_min = grid_profiles.z_min() - FIG_Z_PADDING_DEFAULT * z_range
        z_max = grid_profiles.z_max() + FIG_Z_PADDING_DEFAULT * z_range
        axis_params.z_min = z_min
        axis_params.z_max = z_max
        aspect = 1

    # elevation parameters

    elevation_params = kargs.pop("elevation_params", ElevationPlotParams())

    # others

    num_profiles = grid_profiles.num_profiles()

    if not superposed:
        fig = plt.figure(constrained_layout=True)
        spec = gridspec.GridSpec(
            ncols=1,
            nrows=num_profiles,
            figure=fig)
    else:
        fig = plt.figure()
        ax = fig.add_axes(
            [0.1, 0.1, 0.8, 0.8]
        )
        spec = None

    fig.set_size_inches(width, height)


    if axis_params.grid:
        ax.grid(
            True,
            color=axis_params.grid_color,
            linestyle=axis_params.grid_linestyle,
            linewidth=axis_params.grid_linewidth)

    x_arr = grid_profiles._zarray._s_array

    for profile_ndx in range(grid_profiles.num_profiles()):

        if superposed:
            linecolor = LNINT_ADDITIONAL_COLORS[profile_ndx % len(LNINT_ADDITIONAL_COLORS)]
        else:
            if elevation_params.color is None:
                linecolor = topoprofile_default_linecolor
            elif isinstance(elevation_params.color, str):
                linecolor = elevation_params.color
            else:
                linecolor = topoprofile_default_linecolor

        ax.plot(
            x_arr,
            grid_profiles._zarray._z_array[profile_ndx, :],
            color=linecolor,
            linestyle=elevation_params.linestyle,
            linewidth=elevation_params.width
        )

    ax.set_aspect(aspect)
    ax.set_ylim([axis_params.z_min, axis_params.z_max])

    if axis_params.breaklines:
        bottom, top = ax.get_ylim()
        ax.vlines(
            grid_profiles._zarray._s_breaks_array,
            bottom,
            top,
            color=axis_params.breaklines_color,
            linewidth=axis_params.breaklines_width,
            linestyles=axis_params.breaklines_style
        )

    return fig


def map_profile(
    grid: Grid,
    geoprofile: GeoProfile,
    width: numbers.Real = 5,  # inches
    height: numbers.Real = 2.5,  # inches
    width_ratios_map: numbers.Real = 1.0,
    width_ratios_profile: numbers.Real = 5.0,
    grid_colormap="gist_earth",
    lines: List[Ln] = None,
    linecolor: str = "red",
    linestyle: str = '-',
    linewidth: numbers.Real = 1.5,
    linelabels: bool = True,
    map_zoom: numbers.Real = 1,
    plot_colorbar: bool = False,
    hillshade: bool = False,
    hs_vert_exagg: numbers.Real = 1.0,
    hs_blend_mode: str = 'hillshade',  # one of 'hillshade', 'hsv', 'overlay', 'soft'
    hs_light_source_azim: numbers.Real = 315.0,
    hs_light_source_degr: numbers.Real = 45.0,
    file_path: Optional[str] = None,
    **kargs
) -> Tuple[Union[None, Figure], Error]:
    """

    """

    # plot grid

    if grid.has_rotation:
        return None, Error(
            True,
            caller_name(),
            Exception(f"Grids with rotations are not supported"),
            traceback.format_exc()
        )

    if not geoprofile.has_topography():
        return None, Error(
            True,
            caller_name(),
            Exception(f"Geoprofile has no topography defined"),
            traceback.format_exc()
        )

    fig, (ax1, ax2) = plt.subplots(
        1,
        2,
        gridspec_kw={
            'width_ratios': [
                width_ratios_map,
                width_ratios_profile
            ]
        }
    )

    fig.set_size_inches(
        width,
        height
    )

    err = subplot_map(
        ax=ax1,
        grid=grid,
        lines=lines,
        grid_colormap=grid_colormap,
        linecolor=linecolor,
        linestyle=linestyle,
        linewidth=linewidth,
        linelabels=linelabels,
        map_zoom=map_zoom,
        plot_colorbar=plot_colorbar,
        hillshade=hillshade,
        hs_vert_exagg=hs_vert_exagg,
        hs_blend_mode=hs_blend_mode,
        hs_light_source_azim=hs_light_source_azim,
        hs_light_source_degr=hs_light_source_degr,
    )

    if err:
        return None, err

    subplot(
        geoprofile,
        ax2,
        **kargs
    )

    if file_path is not None:
        plt.savefig(file_path)

    return fig, Error()


@singledispatch
def animated_profiles(

):
    """
    Create an animation of profiles with their map traces.
    """


@animated_profiles.register(GeoProfiles)
def _(
    geoprofiles: GeoProfiles,
    grid: Grid,
    traces: List[Ln],
    animation_flpth: str,
    dpi_resolution: numbers.Integral = 250,
    **kargs
) -> Error:
    """
    Base commands from: https://towardsdatascience.com/basics-of-gifs-with-pythons-matplotlib-54dd544b6f30
    (cons.2021-08-01)
    """

    print("Creating geoprofiles animation")

    try:

        animation_folder_path = os.path.dirname(animation_flpth)

        z_min = kargs.pop("z_min", None)
        z_max = kargs.pop("z_max", None)

        if z_min is None or z_max is None:
            z_range = geoprofiles.z_max() - geoprofiles.z_min()
            z_min = geoprofiles.z_min() - FIG_Z_PADDING_DEFAULT * z_range
            z_max = geoprofiles.z_max() + FIG_Z_PADDING_DEFAULT * z_range

        if np.isnan(z_min) or np.isnan(z_max):
            return Error(
                True,
                caller_name(),
                Exception("z min and/or z max are Nan"),
                traceback.format_exc()
            )

        filepaths = []

        for ndx in range(geoprofiles.num_profiles()):

            print(f"Creating profile {ndx}")

            curr_trace = traces[ndx]
            geoprofile = geoprofiles[ndx]

            fig, err = map_profile(
                grid,
                geoprofile,
                lines=[curr_trace],
                z_min=z_min,
                z_max=z_max,
                profile_ndx=ndx,
                **kargs
            )

            if err:
                return err

            filepath = os.path.join(
                animation_folder_path,
                f"map_profile_{ndx:02d}.png"
            )
            filepaths.append(filepath)

            fig.savefig(filepath, dpi=dpi_resolution)

        # build gif
        with imageio.get_writer(animation_flpth, mode='I', fps=3) as writer:
            for filepath in filepaths:
                image = imageio.imread(filepath)
                writer.append_data(image)

        print(f"Geoprofiles animation saved as {animation_flpth}")

        return Error()

    except Exception as e:

        return Error(
            True,
            caller_name(),
            e,
            traceback.format_exc())


# bokeh profiles

default_width = 18.5
default_height = 10.5


@singledispatch
def plot(
    obj,
    **kargs
) -> Optional[Figure]:
    """

    :param obj:
    :param kargs:
    :return:
    """

    fig = kargs.get("fig", None)
    aspect = kargs.get("aspect", 1)
    width = kargs.get("width", default_width)
    height = kargs.get("height", default_height)

    if fig is None:

        output_notebook()
        fig = figurefigure()

    show(fig)

    return fig


@plot.register(ZTrace)
def _(
    xyarrays: ZTrace,
    **kargs
) -> Optional[Figure]:

    fig = kargs.get("fig", None)

    if fig is None:

        output_notebook()
        fig = Figure()

    fig.match_aspect = True

    fig.line(
        xyarrays.s_arr(),
        xyarrays.z_arr(),
        line_width=0.75,
    )

    show(fig)

    return fig

