import collections
from functools import singledispatch
from math import isfinite, radians, cos, sin
from warnings import warn

from matplotlib.figure import Figure
import matplotlib.patches as mpatches
from matplotlib.gridspec import GridSpec
from matplotlib.ticker import FuncFormatter, MaxNLocator

from geogst.core.geology.profiles.geoprofiles import GeoProfile, GeoProfiles, MultiGridsProfiles
from geogst.core.geology.profiles.profiletraces import PointTrace, ZTrace
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

import colorsys
import re
import random


def _get_cat_id(arrlist) -> str:
    """
    Prova in ordine:
    - .id (tipico di CategorizedArrays nel tuo arrays.py)
    - ._id
    - parsing dal repr: 'rec_id=XYZ'
    """
    for attr in ("id", "_id"):
        if hasattr(arrlist, attr):
            v = getattr(arrlist, attr)
            if v is not None:
                return str(v).strip()

    m = re.search(r"rec_id=([^\n\r\)]+)", repr(arrlist))
    return m.group(1).strip() if m else ""

def make_polygon_colors(
    polygons_intersections,
    seed: int | None = None,
) -> dict[str, tuple[float, float, float]]:
    """
      Ritorna: { "PL": (r, g, b), "Qt": (r, g, b), ... }
      Valori in [0, 1]. Compatibile con matplotlib e convertibile a QColor.
    """
    if seed is not None:
        random.seed(seed)

    cats = []
    def walk(x):
        if x is None:
            return
        if isinstance(x, (list, tuple)):
            for e in x:
                walk(e)
        else:
            # x dovrebbe essere ArrayList/CategorizedArrays
            cat = _get_cat_id(x)
            if cat:
                cats.append(cat)

    walk(polygons_intersections)

    uniq = sorted(set(cats))
    if not uniq:
        return {}

    # colori “ben separati”: hue diverso per categoria, sat/val random in range buono
    colors = {}
    n = len(uniq)
    for i, cat in enumerate(uniq):
        # Hue distribuito uniformemente
        h = (i / max(1, n)) % 1.0
        # not security/cryptography-relevant: only used to vary plot colors
        s = random.uniform(0.6, 0.9)  # nosec B311
        v = random.uniform(0.7, 0.95)  # nosec B311

        r, g, b = colorsys.hsv_to_rgb(h, s, v)
        colors[cat] = (r, g, b)

    return colors


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

    if not polygon_intersections_colors:

        polygon_intersections_colors = make_polygon_colors(polygons_intersections)

    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
    )

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)

            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,
                    aspect)

                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
                    )

        elif points_projections_params.plot_type == PointLikePlotTypes.FOCAL_MECHANISMS:

            plot_focal_mechanisms(
                ax,
                dataset,
                points_projections_params,
            )

        else:
            raise NotImplementedError


def _section_view_moment_tensor(
    strike: numbers.Real,
    dip: numbers.Real,
    rake: numbers.Real,
    profile_azimuth: numbers.Real,
):
    """
    Rotates a double-couple focal mechanism (given as strike/dip/rake, in the true
    North-East-Down frame) into a "fake NED" frame aligned with the vertical section:
    fake-North = true vertical up, fake-East = local along-profile direction, fake-Down =
    horizontal direction perpendicular to the profile.

    Feeding the resulting 6 moment tensor components back in as if they were a regular
    NED mechanism yields the standard map-view rendering of this fake frame, which is
    exactly the true cross-section (side) view of the original mechanism.

    :param strike: the RHR strike of the focal plane, in degrees.
    :param dip: the dip angle of the focal plane, in degrees.
    :param rake: the rake of the focal mechanism, in degrees.
    :param profile_azimuth: the azimuth (degrees from North) of the local along-profile direction.
    :return: the 6 moment tensor components (M_nn, M_ee, M_dd, M_ne, M_nd, M_ed) in the fake frame.
    """

    from obspy.imaging.scripts.mopad import strikediprake_2_moments

    m_nn, m_ee, m_dd, m_ne, m_nd, m_ed = strikediprake_2_moments(strike, dip, rake)

    moment_tensor = np.array([
        [m_nn, m_ne, m_nd],
        [m_ne, m_ee, m_ed],
        [m_nd, m_ed, m_dd],
    ])

    az_rad = radians(profile_azimuth)
    profile_n, profile_e = cos(az_rad), sin(az_rad)

    fake_north = np.array([0.0, 0.0, -1.0])              # true vertical up
    fake_east = np.array([profile_n, profile_e, 0.0])    # local along-profile direction
    fake_down = np.array([profile_e, -profile_n, 0.0])   # perpendicular-to-profile direction

    rotation = np.vstack([fake_north, fake_east, fake_down])
    rotated = rotation @ moment_tensor @ rotation.T

    return (
        rotated[0, 0],
        rotated[1, 1],
        rotated[2, 2],
        rotated[0, 1],
        rotated[0, 2],
        rotated[1, 2],
    )


def _format_axis_tick(
    metres: numbers.Real,
    metres_per_unit: numbers.Real,
    negate: bool = False,
) -> str:
    """
    Renders one tick value, given in metres, as a label in the chosen unit.

    The plotted values are always metres — elevations upwards on the vertical axis,
    distance along the profile on the horizontal one. Only their presentation
    changes here, so that a section of hypocentres can be read in the km-below-sea-
    level convention it is catalogued in without altering the data.

    :param metres: the tick's value, in metres.
    :param metres_per_unit: metres per label unit, e.g. 1000 to label in km.
    :param negate: whether to flip the sign, i.e. to label an elevation as a depth.
    :return: the tick label.
    """

    value = metres / metres_per_unit

    if negate:
        value = -value

    # a whole number reads as one: dropping the always-zero decimal frees real
    # width on the left of every stacked profile
    if value == int(value):
        return f"{int(value)}"

    return f"{value:g}"


def _label_vertical_axis(
    ax,
    axis_params,
):
    """
    Sets how a profile's vertical axis is ticked and labelled.

    Done with a locator and a formatter, and deliberately not by writing the ticks
    out with set_yticks: that widens the axis to cover whatever range the ticks
    span, so the elevation range the user asked for was silently overshot — the
    wider the tick spacing, the worse. A formatter leaves the limits alone.

    A tick landing exactly on the top of the view is left unlabelled, since
    stacked profiles share that edge and the label would sit against the bottom
    one of the profile above.

    :param ax: the matplotlib axis to label.
    :param axis_params: the AxisPlotParams holding the labelling choices.
    """

    max_ticks = getattr(axis_params, "z_axis_max_ticks", ZAXIS_MAX_TICKS_DEFAULT)
    units = getattr(axis_params, "z_axis_units", ZAXIS_UNITS_DEFAULT)
    as_depth = getattr(axis_params, "z_axis_as_depth", ZAXIS_AS_DEPTH_DEFAULT)

    metres_per_unit = axis_units.get(units, 1.0)

    # MaxNLocator counts intervals, which is one fewer than the labels they carry:
    # asking it for max_ticks would show one label more than requested
    ax.yaxis.set_major_locator(MaxNLocator(nbins=max(1, max_ticks - 1)))

    def format_tick(elevation, _position):

        top_of_view = ax.get_ylim()[1]

        if abs(elevation - top_of_view) <= ZAXIS_TOP_LABEL_TOLERANCE * max(1.0, abs(top_of_view)):
            return ""

        return _format_axis_tick(elevation, metres_per_unit, as_depth)

    ax.yaxis.set_major_formatter(FuncFormatter(format_tick))

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


def _label_horizontal_axis(
    ax,
    axis_params,
):
    """
    Sets how a profile's along-profile distance axis is labelled.

    Only the unit is configurable here: the tick positions are left to matplotlib,
    since stacked profiles share this axis and only the bottom one shows its
    labels, so they are not competing for room the way the vertical ones are.

    :param ax: the matplotlib axis to label.
    :param axis_params: the AxisPlotParams holding the labelling choices.
    """

    units = getattr(axis_params, "s_axis_units", SAXIS_UNITS_DEFAULT)

    metres_per_unit = axis_units.get(units, 1.0)

    ax.xaxis.set_major_formatter(
        FuncFormatter(lambda distance, _position: _format_axis_tick(distance, metres_per_unit))
    )


def _name_axes(
    ax,
    axis_params,
    panel_ndx: numbers.Integral,
    num_panels: numbers.Integral,
):
    """
    Names the two axes, with their units.

    Left to matplotlib to place, by way of the axes' own labels, rather than
    positioned by hand on the figure. A profile's box is only as wide as its
    vertical exaggeration allows, so it can end up in a narrow band of a much
    wider figure — and that narrowing happens when the figure is drawn, after
    this runs, so any position computed here from the box would be stale.

    Stacked profiles share both axes, so each name is written once: the vertical
    one against the middle profile, where it reads as belonging to the stack, and
    the horizontal one against the lowest, the only one showing its tick labels.

    Which panel gets each name is decided by the panel's place in the figure, not
    by the profile's own number: a profile exploded on its own is the middle and
    the lowest panel of its figure whatever number it carries in the stack it came
    from, and must be named all the same.

    :param ax: the profile's axis.
    :param axis_params: the AxisPlotParams holding the labelling choices.
    :param panel_ndx: the panel's place in the figure, from 0.
    :param num_panels: how many panels the figure holds.
    """

    if ax is None or not getattr(axis_params, "show_axis_labels", AXIS_LABELS_DEFAULT):
        return

    z_units = getattr(axis_params, "z_axis_units", ZAXIS_UNITS_DEFAULT)
    s_units = getattr(axis_params, "s_axis_units", SAXIS_UNITS_DEFAULT)
    as_depth = getattr(axis_params, "z_axis_as_depth", ZAXIS_AS_DEPTH_DEFAULT)

    if panel_ndx == num_panels // 2:

        vertical_quantity = "depth" if as_depth else "elevation"

        ax.set_ylabel(
            f"{vertical_quantity} ({z_units})",
            fontsize=AXIS_LABEL_FONT_SIZE,
            # clears the tick labels and the profile numbers alike, both of which
            # already sit to the left of the box
            labelpad=AXIS_LABEL_PAD_POINTS,
        )

    if panel_ndx == num_panels - 1:

        ax.set_xlabel(
            f"distance along profile ({s_units})",
            fontsize=AXIS_LABEL_FONT_SIZE,
        )


def plot_focal_mechanisms(
    ax,
    dataset: Dict,
    focal_mechanism_params,
):
    """
    Plots earthquake focal mechanisms as obspy beachballs, at their (s, z) projected location.
    Sized in points (screen units) so it stays circular regardless of the profile's vertical
    exaggeration.

    Two views are supported, via focal_mechanism_params.view:
    - FocalMechanismView.MAP: the original map-view beachball glyph (strike/dip/rake as given).
    - FocalMechanismView.SECTION: the mechanism rotated into the vertical section plane, i.e.,
      as it would appear viewed horizontally, perpendicular to the local profile direction.

    :param ax: the matplotlib axis to plot onto.
    :param dataset: mapping from record id to FocalMechanismTrace.
    :param focal_mechanism_params: the FocalMechanismPlotParams instance controlling the rendering.
    """

    try:
        from obspy.imaging.mopad_wrapper import beach
    except ImportError:
        raise ImportError(
            "Plotting focal mechanisms requires obspy. Install it with: pip install geogst[seismology]"
        )

    for rec_id, focal_mechanism in dataset.items():

        if focal_mechanism_params.view == FocalMechanismView.SECTION:
            fm = _section_view_moment_tensor(
                focal_mechanism.strike,
                focal_mechanism.dip,
                focal_mechanism.rake,
                focal_mechanism.profile_azimuth,
            )
            mopad_basis = 'NED'
        else:
            fm = (focal_mechanism.strike, focal_mechanism.dip, focal_mechanism.rake)
            mopad_basis = 'USE'

        beachball = beach(
            fm,
            mopad_basis=mopad_basis,
            xy=(focal_mechanism.s, focal_mechanism.z),
            width=focal_mechanism_params.width,
            facecolor=focal_mechanism_params.facecolor,
            bgcolor=focal_mechanism_params.bgcolor,
            edgecolor=focal_mechanism_params.edgecolor,
            linewidth=focal_mechanism_params.linewidth,
            axes=ax,
            alpha=focal_mechanism_params.alpha,
        )

        ax.add_collection(beachball)

        if focal_mechanism_params.labels:

            # the beachball is sized in points, so its label has to be offset in
            # points too: a data-unit offset is meaningless against a glyph whose
            # size does not follow the data, and at profile scale it left the
            # label on top of the beachball, unreadable against its own colours.
            # Offset past the radius so the label always clears the glyph.
            label_offset = focal_mechanism_params.width / 2 + FMEC_LABEL_GAP_POINTS

            ax.annotate(
                f"{rec_id}",
                (focal_mechanism.s, focal_mechanism.z),
                textcoords="offset points",
                xytext=(label_offset, label_offset),
                size=FMEC_LABEL_FONT_SIZE,
                color=focal_mechanism_params.edgecolor,
            )


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)):
                       continue
                    if not isfinite(s_val) or not isfinite(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):
        return

    values = points_projections.values()
    value_types = list(set(map(lambda value: type(value), values)))
    if len(value_types) != 1:
        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,
        )

    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
):

    if 'all features' in polygon_intersections_colors:
        default_color = polygon_intersections_colors['all features']
    else:
        default_color = 'red'

    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

                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.get(polygon_intersection_id, default_color),
                    linewidth=polygon_intersections_linewidth
                )




@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():
        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)
    # the panel's place in this figure, distinct from the profile's own number:
    # an exploded profile is the only panel of its figure whatever number it bears
    panel_ndx = kwargs.pop("panel_ndx", 0)
    num_panels = kwargs.pop("num_panels", 1)
    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_colors = 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_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])

    _label_vertical_axis(ax, axis_params)
    _label_horizontal_axis(ax, axis_params)
    _name_axes(ax, axis_params, panel_ndx, num_panels)

    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:

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

        else:

            manage_polygons_intersections_plotting(
                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 line 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:

        if not geoprofile._topo_profile:

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

        else:

            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)

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

            ax.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),
                        aspect,
                        segment_scale_factor=3)

                    ax.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}"

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

    # Profile id, in the top left corner inside the plot box. Left of the box is
    # the column of elevation tick labels, against which the number read as one of
    # them, and how wide that column is depends on the values shown, so no fixed
    # gap clears it for every profile. The corner does, and it is the part of the
    # box a topographic profile is least likely to reach — where it does, the patch
    # the number is written on keeps it legible.
    # Offset in points, which survives the layout being reworked when the figure
    # is drawn, and leaves the axis label itself free to name the axis.

    ax.annotate(
        str(profile_ndx + 1),
        xy=(0, 1),
        xycoords='axes fraction',
        xytext=(PROFILE_LABEL_INSET_POINTS, -PROFILE_LABEL_INSET_POINTS),
        textcoords='offset points',
        fontsize=PROFILE_LABEL_FONT_SIZE,
        ha='left',
        va='top',
        bbox=dict(
            boxstyle='square,pad=0.2',
            facecolor='white',
            edgecolor='none',
            alpha=PROFILE_LABEL_BACKGROUND_ALPHA,
        ),
    )

    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")
        kargs.pop("aspect", None)
    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_params.vertical_exaggeration = kargs.pop("aspect", 1)

    # others

    num_profiles = geoprofiles.num_profiles()
    polygon_intersections_colors = kargs.get("polygon_intersections", None)
    has_polygon_legend = bool(polygon_intersections_colors)

    if not superposed:

        fig = plt.figure()

        if has_polygon_legend:
            gs = GridSpec(
                nrows=num_profiles,
                ncols=2,
                figure=fig,
                width_ratios=[20, 4],
                height_ratios=[1] * num_profiles,
                wspace=0.05,
            )
        else:
            gs = GridSpec(
                nrows=num_profiles,
                ncols=1,
                figure=fig,
                height_ratios = [1]*num_profiles,
                wspace=0.05,
            )

        axes = []

        if has_polygon_legend:
            ax = fig.add_subplot(gs[0, 0])
        else:
            ax = fig.add_subplot(gs[0])


        if num_profiles > 1:
            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()):

            if has_polygon_legend:
                ax = fig.add_subplot(gs[ndx, 0], sharex=axes[0])
            else:
                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)

    polygon_intersections_colors = kargs.get("polygon_intersections", None)

    for ndx in range(num_profiles):

        geoprofile = geoprofiles[ndx]

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

    polygon_inters_legend_on = PLINT_LEGEND_DEFAULT #if polygon_intersections is None else polygon_intersections.get("legend", PLINT_LEGEND_DEFAULT)

    if has_polygon_legend and 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)))

        legend_ax = fig.add_subplot(gs[:, 1])
        legend_ax.axis("off")

        # Put a legend to the right of the current axis
        legend_ax.legend(
            handles=legend_patches,
            loc='center left',
            frameon=False,
        )

    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(
            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


# 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 = figure()

    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

