
import sys
fo = sys.stdout

import numbers
from functools import singledispatch
import collections

import os.path
from warnings import warn

import 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

