
import logging
from typing import Callable, Dict
from collections import defaultdict, namedtuple

from geogst.core.geometries.general.planes import versor_from_plane_attitude_inters
from geogst.core.geology.profiles.profiletraces import *
from geogst.core.geometries.general.planes import CPlane3D, vertical_plane_from_segment
from geogst.core.geometries.shapes.multilines import MultiLine
from geogst.core.geometries.shapes.segments import Segment
from geogst.core.geometries.shapes.polygons import Polygon
from geogst.core.geometries.grids.rasters import Grid
from geogst.core.geometries.calculations.intersections import (
    intersect,
    PreparedLine,
    prepare_line_for_intersection,
    intersect_segment_with_prepared_line,
    PreparedPolygon,
    PreparedPolygons,
    prepare_polygon_for_intersection,
    intersect_segment_with_prepared_polygon,
    select_segments_inside_prepared_polygon,
)
from geogst.core.geology.profiles.methods import ProjectionMethod
from geogst.core.geometries.projections.crs import Crs
from geogst.core.geometries.projections.project import project, project_xy_arrays
from geogst.core.utils.arrays import *

logger = logging.getLogger(__name__)


# How far, in metres, a vertex may lie from the straight chord of its run and
# still be treated as lying on it, when collapsing a densified line for point
# projection.
#
# This is deliberately a distance and not an angle. Profile vertices do not carry
# only the rounding of the densification that produced them: they are typically
# also draped on a DEM and reprojected into the project's CRS, and a reprojection
# does not preserve exact collinearity. An angular criterion tight enough to be
# safe against a genuine bend is then broken by sub-micron coordinate noise, and
# the collapse silently stops happening — measured, a lateral noise of a tenth of
# a micron was already enough to leave a straight 500-segment line uncollapsed.
#
# The error this admits is NOT bounded by the tolerance itself. Replacing a run by
# its chord tilts the projection plane by at most tolerance/segment_length, and a
# point is displaced along the profile by that angle times its own distance from
# the profile, so:
#
#     s displacement <= (tolerance / segment length) * max_profile_distance
#
# At a millimetre, with segments of tens of metres and points accepted up to a few
# kilometres away, that is a few centimetres at worst — negligible against the
# metre-scale resolution of a profile, and orders of magnitude above the noise
# reprojection leaves behind.
COLLINEARITY_DEVIATION_TOLERANCE = 1.0e-3

# Coordinate gap above which two consecutive segments are taken not to be joined,
# so that a run is never merged across a discontinuity (as in a multi-part line).
SEGMENT_CONTINUITY_TOLERANCE = 1.0e-6


def _max_deviation_from_chord(
    chord_start: np.ndarray,
    chord_end: np.ndarray,
    vertices: np.ndarray,
) -> float:
    """
    Returns the largest perpendicular distance of the given 2D vertices from the
    chord running between two positions, or 0.0 when there are no vertices.

    Used to decide whether a run of profile segments is straight enough to be
    replaced by that chord: the value returned is exactly the geometric error the
    replacement would introduce.

    :param chord_start: the chord start, as an (x, y) array.
    :param chord_end: the chord end, as an (x, y) array.
    :param vertices: the vertices to measure, as an (n, 2) array.
    :return: the largest perpendicular distance.
    """

    if len(vertices) == 0:
        return 0.0

    chord = chord_end - chord_start
    chord_length = float(np.hypot(chord[0], chord[1]))

    offsets = vertices - chord_start

    if chord_length == 0.0:
        # degenerate chord: fall back to the distance from its single position
        return float(np.max(np.hypot(offsets[:, 0], offsets[:, 1])))

    # magnitude of the 2D cross product, i.e. the perpendicular distance times
    # the chord length
    cross = offsets[:, 0] * chord[1] - offsets[:, 1] * chord[0]

    return float(np.max(np.abs(cross)) / chord_length)


ProfileBaseLineParameters = namedtuple(
    
    "ProfileBaseLineParameters",
    [
        "name", # str
        "source_type", # Enum
        "baseline", # Ln
        "wkt_crs", # str
    ]
)


class SegmentProfiler:
    """
    Class storing a linear (straight) profile.
    It intends to represent a segmented vertical profile.
    In a possible future implementations, it would be superseded by a
    plane profilers, not necessarily vertical.
    """

    def __init__(self,
                 segment: Segment,
                 crs_code: Union[numbers.Integral, str] = 4326,
                 ):
        """
        Instantiates a 2D segment profile object.


        :param segment: the profile segment.
        """

        self._segment = Segment(
            start_pt=segment.start_pt.as_point2d(),
            end_pt=segment.end_pt.as_point2d()
        )
        self._crs_code = crs_code

    @classmethod
    def from_points(cls,
                    start_pt: Point,
                    end_pt: Point,
                    crs_code: Union[numbers.Integral, str] = 4326,
                    ):
        """
        Instantiates a 2D segment profile object.
        It is represented by two points and an EPSG code.

        :param start_pt: the profile start point.
        :param end_pt: the profile end point.
        """

        _start_pt = Point(
            start_pt.x,
            start_pt.y
        )

        _end_pt = Point(
            end_pt.x,
            end_pt.y
        )

        return cls(
            Segment(_start_pt, _end_pt),
            crs_code=crs_code,
        )

    @property
    def segment(self) -> Segment:
        """
        Returns the horizontal segment representing the profile.

        :return: segment representing the profile.
        """

        return self._segment.clone()

    @property
    def crs_code(self) -> numbers.Integral:
        return self._crs_code

    def start_pt(self) -> Point:
        """
        Returns a copy of the segment start point.

        :return: segment start point copy.
        """

        return Point(
            self.segment.start_pt.x,
            self.segment.start_pt.y
        )

    def end_pt(self) -> Point:
        """
        Returns a copy of the segment end point.

        :return: segment end point copy.
        """

        return Point(
            self.segment.end_pt.x,
            self.segment.end_pt.y
        )

    def to_line(self) -> Ln:
        """
        Convert to a line.

        """

        return Ln.from_points(
                self.start_pt(),
                self.end_pt()
        )

    def __repr__(self):
        """
        Representation of a profile instance.

        :return: the textual representation of the instance.
        """

        return f"SegmentProfiler(\n\tstart_pt = {self.start_pt()},\n\tend_pt = {self.end_pt()},\n\tepsg_code = {self.crs_code})"

    def clone(self) -> 'SegmentProfiler':
        """
        Returns a deep copy of the current profilers.

        :return: a deep copy of the current profilers.
        """

        return SegmentProfiler(
            segment=self.segment,
            crs_code=self.crs_code,
        )

    def length(self) -> numbers.Real:
        """
        Returns the length of the profilers section.

        :return: length of the profilers section.
        """

        return self.segment.length

    def vector(self) -> Vect3D:
        """
        Returns the horizontal 3D vector representing the profile.

        :return: vector representing the profile.
        """

        return Vect3D(
            x=self.segment.delta_x(),
            y=self.segment.delta_y(),
            z=0.0
        )

    def versor(self) -> Vect3D:
        """
        Returns the horizontal 3D versor (unit vector) representing the profile.

        :return: versor representing the profile.
        """

        return self.vector().versor()

    def densify_as_steps_array(self,
                               sampling_distance: numbers.Real
                               ) -> Tuple[bool, Union[str, array]]:
        """
        Returns an array made up by the incremental steps (2D distances) along the profile.

        :param sampling_distance: the segment profilers sampling distance
        :return: optional array storing the incremental steps, with the last step being equal to the segment length.
        """

        return densify_as_array1d_(
            segment_length=self.segment.length_2d(),
            densify_distance=sampling_distance
        )

    def densify_as_2D_points(self,
                             densify_distance: numbers.Real
                             ) -> Tuple[Union[None, List[Point]], Error]:
        """
        Returns a list of densified points.

        :param densify_distance: the profilers densify distance.
        :return: optional list of densified points and error status.
        """

        return self.segment.densify_as_pts2d(
            densify_distance=densify_distance
        )

    def segment_as_vertical_plane(self) -> Optional[CPlane3D]:
        """
        Returns the vertical plane of the segment, as a Cartesian plane.

        :return: the vertical plane of the segment, as a Cartesian plane.
        """

        return vertical_plane_from_segment(
            segment=self.segment
        )

    def normal_versor(self) -> Tuple[Union[None, Vect3D], Error]:
        """
        Returns the perpendicular (horizontal) versor to the profile (vertical) plane.

        :return: the perpendicular (horizontal) versor to the profile (vertical) plane.
        """

        vert_plane = self.segment_as_vertical_plane()
        if vert_plane is None:
            return None, Error(
                True,
                caller_name(),
                Exception("Vertical plane is None"),
                traceback.format_exc())

        return vert_plane.normal_versor()

    def left_norm_vers(self) -> Vect2D:
        """
        Returns the left horizontal normal versor.

        :return: the left horizontal normal versor.
        """

        return Vect2D(
            x=-self.versor().y,
            y=self.versor().x
        )

    def right_norm_vers(self) -> Vect2D:
        """
        Returns the right horizontal normal versor.

        :return: the right horizontal normal versor.
        """

        return Vect2D(
            x=self.versor().y,
            y=-self.versor().x
        )

    def shift(self,
              dx: numbers.Real,
              dy: numbers.Real
    ) -> 'SegmentProfiler':
        """
        Returns a new LinearProfiler instance, horizontally offset by the
        provided horizontal components.
        """

        return SegmentProfiler(
            segment=self.segment.shift(dx, dy),
            crs_code=self.crs_code,
        )

    def vector_offset(self,
                      vect: Vect2D
                      ) -> 'SegmentProfiler':
        """
        Returns a new SegmentProfiler instance, horizontally offset by the
        provided vector horizontal components.
        """

        return SegmentProfiler(
            segment=self.segment.shift(vect.x, vect.y),
            crs_code=self.crs_code,
        )

    def right_parallel_offset(self,
                              offset: numbers.Real
                              ) -> 'SegmentProfiler':
        """
        Returns a copy of the current profilers, offset to the right by the provided offset distance.

        :param offset: the lateral offset to apply to create the new profilers.
        :return: the offset profilers.
        """

        return self.vector_offset(
            vect=self.right_norm_vers().scale(offset),
        )

    def left_parallel_offset(self,
                             offset: numbers.Real
                             ) -> 'SegmentProfiler':
        """
        Returns a copy of the current profilers, offset to the left by the provided offset distance.

        :param offset: the lateral offset to apply to create the new profilers.
        :return: the offset profilers.
        """

        return self.vector_offset(vect=self.left_norm_vers().scale(offset))

    def point_in_vertical_profile(self,
                                  pt: Point
                                  ) -> bool:
        """
        Checks whether a point lies in the profilers plane.

        :param pt: the point to check.
        :return: whether the point lies in the profilers plane.
        :raise: Exception.
        """

        pt_use = Point(
                pt.x,
                pt.y,
                0.0
            )

        is_point_in_vertical_plane = self.segment_as_vertical_plane().contains_point(pt_use)

        return is_point_in_vertical_plane

    def point_horizontal_distance(self,
                                  pt: Point
                                  ) -> numbers.Real:
        """
        Calculates the point distance from the profilers plane.

        :param pt: the point to check.
        :return: the point distance from the profilers plane.
        :raise: Exception.
        """

        pt_use = Point(
                pt.x,
                pt.y,
                0.0
            )

        return self.segment_as_vertical_plane().absolute_distance_to_point(pt_use)

    def sample_grid(
        self,
        grid: Grid,
        sampling_distance: numbers.Real,
    ) -> Tuple[Union[None, array], Error]:
        """
        Sample grid values along the profilers points.

        :param sampling_distance: the grid sampling distance along the trace.
        :param grid: the input grid.
        :return: array storing the z values sampled from the grid and an error status.
        """

        try:

            result, error = self.densify_as_2D_points(sampling_distance)

            if error:
                return None, error

            pts_2d = result

            values_arr = array(
                'd',
                [
                    grid.interpolate_bilinear(pt_2d.x, pt_2d.y) for pt_2d in pts_2d
                ]
            )

            return values_arr, Error()

        except Exception as e:

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

    def intersect_mline(self,
                        mline: MultiLine,
                        ) -> List[UnionPtSegment2D]:
        """
        Calculates the intersection with a line/multiline.
        Note: the intersections are intended flat (in a 2D plane, not 3D).

        :param mline: the line/multiline to intersect profile with
        :return: the possible intersections
        """

        if isinstance(mline, (Ln, MultiLine)):
            mline = mline.as_2d()

        return intersect(self.segment, mline)

    def intersect_line(self,
        line: Ln,
    ) -> List[UnionPtSegment2D]:
        """
        Calculates the intersection with a line.
        Note: the intersections are intended flat (in a 2D plane, not 3D).

        :param line: the line to intersect profile with
        :return: the possible intersections
        """

        line = line.as_2d()

        return intersect(
            self.segment,
            line,
        )

    def intersect_lines(self,
        lines: Dict[Category, List[PreparedLine]],
    ) -> Dict[Category, List[UnionPtSegment2D]]:
        """
        Calculates the intersection with a set of lines.
        Note: the intersections are calculated in a 2D plane (not 3D).

        :param lines: a dict of category -> already-prepared lines
            (see prepare_line_for_intersection) to intersect profile with.
            Preparing lines once, outside the per-profile/per-segment loop
            that calls this method, avoids recomputing the same bounding
            boxes and segment lists over and over.
        :return: the optional line-identified intersections
        :raise: Exception
        """

        segment_intersections = defaultdict(list)

        profile_segment = self.segment
        profile_segment_bbox = profile_segment.bounding_box()

        for category, prepared_lines in lines.items():

            for prepared_line in prepared_lines:
                result = intersect_segment_with_prepared_line(
                    profile_segment,
                    profile_segment_bbox,
                    prepared_line,
                )
                if result is not None:
                    segment_intersections[category].append(result)

        return segment_intersections

    def intersect_line_with_attitude(self,
        line: Ln,
        attitude: Plane,
    ) -> Tuple[Union[None, List[Tuple[Vect3D, Point]]], Error]:
        """
        Calculates the intersection with a line with a geological attitude.
        Note: the intersections are intended flat (in a 2D plane, not 3D).

        :param line: the line to intersect profile with
        :param attitude: the surface attitude, property of the line to intersect.
        :return: the possible intersections
        """

        try:

            line = line.as_2d()

            intersection_geometries = intersect(self.segment, line)  # List[Union[None, Point, 'Segment']]

            if intersection_geometries is None:
                logger.warning("intersecting geometries are None")
                return None, Error()

            # get vertical section plane

            section_plane = self.segment_as_vertical_plane()
            if section_plane is None:
                return None, Error(
                    True,
                    caller_name(),
                    Exception("Returned vertical plane from segment profile is None"),
                    traceback.format_exc())

            intersection_results = []

            for ndx, intersection_geometry in enumerate(intersection_geometries):

                if intersection_geometry is None:
                    logger.warning("intersecting geometry %d is None", ndx)
                    continue

                if isinstance(intersection_geometry, Point):

                    # calculate intersection geometry between attitude and plane
                    # note that intersection point must be transformed
                    # to a 3D point.

                    # Tuple[Union[None, Vect3D, CPlane3D], Error]
                    intersection_attitude, err = versor_from_plane_attitude_inters(
                        cartes_plane=section_plane,
                        attitude=attitude,
                        point=intersection_geometry
                    )

                    if err:
                        logger.error("error with intersection geometry %s: %r", intersection_geometry, err)
                        continue

                    if isinstance(intersection_attitude, CPlane3D):
                        logger.warning("intersected geometry is a CPlane3D — likely parallel to the section plane")
                        continue
                    elif isinstance(intersection_attitude, Vect3D):
                        intersection_results.append((intersection_attitude, intersection_geometry))
                    else:
                        logger.warning("intersection geometry is a %s, not a Vect3D", type(intersection_attitude))
                        continue

                else:

                    logger.warning("intersection geometry is a %s, not a Point", type(intersection_geometry))
                    continue

            return intersection_results, Error()

        except Exception as e:

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

    def intersect_lines_with_attitude(self,
        lines: List[Ln],
        attitude: Plane
    ) -> Tuple[Union[None, List[Tuple[Vect3D, Point]]], Error]:
        """
        Calculates the intersection with a set of lines.
        Note: the intersections are calculated in a 2D plane (not 3D).

        :param lines: an iterable of georeferenced Lines or MultiLines to intersect profile with
        :param attitude: the plane expressing the local orientation of the surface traces.
        :return: the optional line-identified intersections
        :raise: Exception
        """

        try:

            point_traces = []

            for line in lines:

                result, err = self.intersect_line_with_attitude(
                    line,
                    attitude)

                if err:
                    logger.error("%r", err)
                    continue

                if result is None:
                    logger.warning("segment intersection with line/attitude is None")
                    continue

                point_traces.extend(result)

            return point_traces, Error()

        except Exception as e:

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

    def intersect_polygon(self,
                          polygon: Polygon) -> List[Segment]:
        """
        Calculates the intersection with a polygon.

        """

        return self.intersect_prepared_polygon(
            prepare_polygon_for_intersection(polygon)
        )

    def intersect_prepared_polygon(self,
                                   prepared_polygon: PreparedPolygon,
                                   segment_bbox: Optional[tuple] = None,
                                   ) -> List[Segment]:
        """
        Same result as intersect_polygon, but on an already-prepared polygon
        (see prepare_polygon_for_intersection), whose boundary segments,
        bounding boxes and ring arrays are computed once for all the profile
        segments (and all the parallel profiles) tested against it.

        :param prepared_polygon: the polygon to intersect, already prepared.
        :param segment_bbox: the optional precomputed bounding box of this
            profiler segment, to avoid recomputing it on every call.
        :return: the parts of the profiler segment lying inside the polygon.
        """

        if segment_bbox is None:
            segment_bbox = self._segment.bounding_box()

        full_intersections = intersect_segment_with_prepared_polygon(
            self._segment,
            segment_bbox,
            prepared_polygon,
        )

        if full_intersections is None:
            return []

        if not full_intersections:
            # the segment does not cross the polygon boundaries: it either
            # lies wholly inside the polygon or wholly outside it, with no
            # need to split it at the (absent) boundary crossings
            if self._segment.start_pt.is_coincident(self._segment.end_pt):
                return []
            return select_segments_inside_prepared_polygon([self._segment], prepared_polygon)

        full_sorted_pts = self.sort_intersections_by_start_distance_as_pts(full_intersections)

        segments = Ln.from_points(self._segment.start_pt, *full_sorted_pts, self._segment.end_pt).remove_coincident_points().segments()

        return select_segments_inside_prepared_polygon(segments or [], prepared_polygon)

    def sort_intersections_by_start_distance_as_pts(self,
                                                    intersections: List[Union[None, Point, Segment]]
                                                    ) -> List[Point]:
        """
        Sort the intersections based on the distance from the profile start.
        """

        inters_points = list(filter(lambda el: isinstance(el, Point), intersections))
        inters_segments = list(filter(lambda el: isinstance(el, Segment), intersections))

        inters_segments_points = []
        for segm in inters_segments:
            inters_segments_points.extend([segm.start_pt, segm.end_pt])

        total_points = inters_points + inters_segments_points
        total_points.sort(key=self._segment.point_distance_from_start)

        return total_points

    def point_along_profile_signed_s(
            self,
            pt: Point,
            offset: numbers.Real = 0.0
    ) -> Union[None, numbers.Real]:
        """
        Calculates the point along-profile signed distance (positive in the segment direction, negative otherwise)
        from the profile start.
        The projected point must already lay in the profile vertical plane, otherwise an exception is raised.

        The implementation assumes (and also verifies) that the point lies in the profile vertical plane.
        Given that case, it calculates the signed distance from the section start point,
        by using the triangle law of sines.

        :param pt: the point on the section.
        :param offset: the offset to apply to the distance from segment start (to comply with line profiles).
        :return: the signed distance along the profile or None if outside the segment.
        """

        if pt.is_2d():
            use_pt = Point(
                pt.x,
                pt.y,
                0.0
            )
        else:
            use_pt = pt.clone()

        if not self.point_in_vertical_profile(use_pt):
            raise Exception(f"Projected point should lie in the profile plane but there is a distance of {self.point_horizontal_distance(use_pt)} units")

        if use_pt.as_point2d().is_coincident(self.start_pt()):
            return 0.0 + offset

        # the vector starting at the profile start and ending at the given point

        projected_vector = Segment(
            Point(
                self.start_pt().x,
                self.start_pt().y,
                0.0),
            use_pt
        ).as_vector3d()

        # the rot_angle between the profile vector and the previous vector
        cos_alpha = self.vector().cosine_of_angle(projected_vector)
        if cos_alpha is None:
            return None

        signed_distance = projected_vector.length * cos_alpha

        return signed_distance + offset

    def segment_along_profile(self,
                              segment: Segment,
                              offset: numbers.Real = 0.0
                              ) -> Tuple[Optional[numbers.Real], Optional[numbers.Real]]:
        """
        Calculates the segment distances from the profiles start.
        The segment must already lay in the profile vertical plane, otherwise None is returned.

        :param segment: the analysed segment.
        :param offset: the offset to apply to the distance from segment start (to comply with line profiles).
        :return: the segment vertices distances from the profile start.
        """

        segment_start_distance = self.point_along_profile_signed_s(segment.start_pt, offset)
        segment_end_distance = self.point_along_profile_signed_s(segment.end_pt, offset)

        return segment_start_distance, segment_end_distance

    def pt_segm_along_profile_signed_s(self,
                                       ptsegm2d: UnionPtSegment2D,
                                       offset: numbers.Real = 0.0
                                       ) -> array:
        """
        Calculates the point or segment signed distances from the segment start.

        :param ptsegm2d: point or segment.
        :param offset: the offset to apply to the distance from segment start (to comply with line profiles).
        :return: the distance(s) from the segment start.
        """

        if isinstance(ptsegm2d, Point):
            arr = array('d', [self.point_along_profile_signed_s(ptsegm2d, offset)])
        elif isinstance(ptsegm2d, Segment):
            arr = array('d', [*self.segment_along_profile(ptsegm2d, offset)])
        else:
            raise Exception(f"Got type {type(ptsegm2d)} for point/segment")

        return arr

    def get_intersection_slope(self,
        intersection_vector: Vect3D
    ) -> Tuple[Union[Tuple[numbers.Real, str], None], Error]:
        """
        Calculates the slope (in radians) and the downward sense ('left', 'right' or 'vertical')
        for a profile-laying vector.

        :param intersection_vector: the profile-plane lying vector.
        :return: the slope (in radians) and the downward sense.
        """

        try:

            norm_versor, error = self.normal_versor()
            if error:
                return None, error

            angle = degrees(acos(norm_versor.cosine_of_angle(intersection_vector)))
            if abs(90.0 - angle) > 1.0e-4:
                return None, Error(
                    True,
                    caller_name(),
                    Exception("Input argument should lay in the profile plane"),
                    traceback.format_exc())

            slope_radians = abs(radians(intersection_vector.dip_angle()))

            scalar_product_for_downward_sense = self.vector().dot_product(intersection_vector.downward())
            if scalar_product_for_downward_sense > 0.0:
                intersection_downward_sense = "right"
            elif scalar_product_for_downward_sense == 0.0:
                intersection_downward_sense = "vertical"
            else:
                intersection_downward_sense = "left"

            return (slope_radians, intersection_downward_sense), Error()

        except Exception as e:

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

    def calculate_pt_intersection_along_axis(
        self,
        project_axis: Axis,
        pt: Point
    ) -> Tuple[Union[Point, 'ParamLine3D', None], Error]:
        """
        Calculates the optional intersection point between an axis passing through a point
        and the profilers plane.

        :param project_axis: the projection axis.
        :param pt: the point through which the axis passes.
        :return: the optional intersection point.
        :raise: Exception.
        """

        axis_versor = project_axis.as_direction().as_versor()

        l, m, n = axis_versor.x, axis_versor.y, axis_versor.z

        axis_param_line = ParamLine3D(pt, l, m, n)

        return intersect_plane_with_line(axis_param_line, self.segment_as_vertical_plane())

    def intersect_with_plane(
            self,
            plane: Plane,
            pt: Point
    ) -> Tuple[Union[None, Vect3D], Error]:
        """
        Calculate the intersection versor between the segment profilers and
        a plane passing through a point.

        :param plane: the plane orientation.
        :param pt: the record location.
        :return: an optional intersection versor and the error status.
        """

        try:

            # input type checks

            check_type(plane, "Attitude plane", Plane)
            check_type(pt, "Attitude point", Point)

            # calculate intersection versor

            vert_pln = self.segment_as_vertical_plane()
            if vert_pln is None:
                return None, Error(
                    True,
                    caller_name(),
                    Exception("the vertical plane is None"),
                    traceback.format_exc())

            intersection, err = vert_pln.intersects_other(
                CPlane3D.from_geol_plane(plane, pt))

            if err:
                return None, err

            if isinstance(intersection, CPlane3D):
                return None, Error(
                    True,
                    caller_name(),
                    Exception("the section and the attitude plane coincide"),
                    traceback.format_exc())

            if not intersection.is_valid():
                return None, Error(
                    True,
                    caller_name(),
                    Exception("intersection is not valid"),
                    traceback.format_exc())

            return intersection, Error()

        except Exception as e:

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

    def nearest_coplanar_projection(
            self,
            plane_attitude: Plane,
            position: Point
    ) -> Tuple[Union[None, Point], Error]:
        """
        Calculates the nearest projection of a given attitude on a vertical plane.

        :param plane_attitude: plane attitude.
        :param  position: attitude position.
        :return: the nearest projected point on the vertical section.
        :raise: Exception.
        """

        try:

            check_type(plane_attitude, "Plane attitude", Plane)
            check_type(position, "Position", Point)

            attitude_cplane = CPlane3D.from_geol_plane(
                geol_plane=plane_attitude,
                pt=position)

            intersection, err = self.segment_as_vertical_plane().intersects_other(attitude_cplane)

            if err:
                return None, err

            if isinstance(intersection, CPlane3D):
                return None, Error(
                    True,
                    caller_name(),
                    Exception("unhandled degenerate case where section and attitude coincides."),
                    traceback.format_exc()
                )

            if not intersection.is_valid():
                return None, Error(
                    True,
                    caller_name(),
                    Exception(f"intersection vector is not valid: {intersection}"),
                    traceback.format_exc()
                )

            vert_plane = self.segment_as_vertical_plane()
            if vert_plane is None:
                return None, Error(
                    True,
                    caller_name(),
                    Exception("vertical plane is None"),
                    traceback.format_exc()
                )

            dummy_inters_pt, err = vert_plane.intersects_other_as_pt(attitude_cplane)
            if err:
                return None, err

            dummy_structural_vect = Segment(dummy_inters_pt, position).as_vector3d()
            dummy_distance = dummy_structural_vect.dot_product(intersection)
            offset_vector = intersection.scale(dummy_distance)

            projected_pt = Point(
                dummy_inters_pt.x + offset_vector.x,
                dummy_inters_pt.y + offset_vector.y,
                dummy_inters_pt.z + offset_vector.z
            )

            return projected_pt, Error()

        except Exception as e:

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

    def project_point(
            self,
            point: Point,
            projection_method: ProjectionMethod = ProjectionMethod.NEAREST,
            min_profile_distance: numbers.Real = MIN_DISTANCE_TOLERANCE,
            map_axis: Optional[Axis] = None,
            attitude: Optional[Plane] = None,
            axis_angular_tolerance: numbers.Real = MIN_DISORIENTATION_TOLERANCE,
            precomputed_plane: Optional[CPlane3D] = None,
    ) -> Tuple[Union[None, Point], Error]:
        """
        Project a point to the segment according to various methods.

        :param point: the point to project.
        :param projection_method: the projection method.
        :param min_profile_distance: the minimum distance of the point from the section to be considered not lying on the section plane.
        :param map_axis: the optional axis to project the point to the section.
        :param attitude: the optional attitude along which to profiles the point to the plane.
        :param axis_angular_tolerance: the minimum angular tolerance (in degrees) between axis and plane.
        :param precomputed_plane: the segment's vertical plane, if already computed by the caller
            (e.g. once per segment before projecting many points), to avoid rebuilding it on every call.
        :return: the projected point and the error status.
        """

        try:

            # input type checks

            check_type(point, "Position", Point)
            check_type(min_profile_distance, "Min point distance from profile", numbers.Real)

            # get point distance from section plane and possibly return result

            section_plane = precomputed_plane if precomputed_plane is not None else self.segment_as_vertical_plane()
            if section_plane is None:
                return None, Error(
                    True,
                    caller_name(),
                    Exception("Returned vertical plane from segment profile is None"),
                    traceback.format_exc()
                )

            if projection_method == ProjectionMethod.NEAREST:

                return section_plane.project_point(
                    pt=point,
                    distance_atol=min_profile_distance)

            elif projection_method in (ProjectionMethod.COMMON_AXIS, ProjectionMethod.INDIVIDUAL_AXES):

                return section_plane.project_point(
                    pt=point,
                    axis=map_axis,
                    distance_atol=min_profile_distance,
                    angular_atol=axis_angular_tolerance)

            elif projection_method == ProjectionMethod.ALONG_PLANE:  # project along attitude

                result, err = intersect_plane_with_attitude(
                    section_plane,
                    attitude,
                    point)

                if err:
                    return None, err
                if result is None:
                    return None, Error()

                _, inters_point = result
                return inters_point, Error()

            else:

                return None, Error(
                    True,
                    caller_name(),
                    Exception(f"Projection method is {projection_method}"),
                    traceback.format_exc()
                )

        except Exception as e:

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

    def project_attitude(
            self,
            attitude: Plane,
            point: Point,
            projection_method: ProjectionMethod = ProjectionMethod.NEAREST,
            min_profile_distance: numbers.Real = MIN_DISTANCE_TOLERANCE,
            map_axis: Optional[Axis] = None,
            axis_angular_tolerance: numbers.Real = MIN_DISORIENTATION_TOLERANCE,
            precomputed_plane: Optional[CPlane3D] = None,
    ) -> Tuple[Union[None, Tuple[Vect3D, Point]], Error]:
        """
        Project an attitude point to the segment according to various methods.

        :param attitude: the attitude to plot to the plane section.
        :param point: the point defining the attitude location.
        :param projection_method: the projection method.
        :param min_profile_distance: the minimum distance of the point from the section to be considered not lying on the section plane.
        :param map_axis: the optional axis to project the point to the section.
        :param axis_angular_tolerance: the minimum angular tolerance (in degrees) between axis and plane.
        :param precomputed_plane: the segment's vertical plane, if already computed by the caller
            (e.g. once per segment before projecting many attitudes), to avoid rebuilding it on every call.
        :return: the intersection and an error status.
        """

        try:

            # get vertical section plane

            section_plane = precomputed_plane if precomputed_plane is not None else self.segment_as_vertical_plane()

            if section_plane is None:
                return None, Error(
                    True,
                    caller_name(),
                    Exception("Returned vertical plane from segment profile is None"),
                    traceback.format_exc())

            # define projected attitude location

            if projection_method == ProjectionMethod.NEAREST:

                result, err = section_plane.project_point(
                    pt=point,
                    distance_atol=min_profile_distance)

            elif projection_method in (ProjectionMethod.COMMON_AXIS, ProjectionMethod.INDIVIDUAL_AXES):

                result, err = section_plane.project_point(
                    pt=point,
                    axis=map_axis,
                    distance_atol=min_profile_distance,
                    angular_atol=axis_angular_tolerance
                )

            else:  # project along attitude

                result, err = intersect_plane_with_attitude(
                    cartes_plane=section_plane,
                    attitude=attitude,
                    point=point
                )

            if err:
                return None, err

            if result is None:
                return None, Error()

            if isinstance(result, Point):
                projected_pt = result
            else:
                _, projected_pt = result

            # calculate intersection geometry between attitude and plane

            intersection, err = versor_from_plane_attitude_inters(
                cartes_plane=section_plane,
                attitude=attitude,
                point=point
            )

            if err:
                return None, err

            return (intersection, projected_pt), Error()

        except Exception as e:

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

    def project_line(
        self,
        line: Ln,
        projection_method: ProjectionMethod = ProjectionMethod.NEAREST,
        min_profile_distance: numbers.Real = MIN_DISTANCE_TOLERANCE,
        map_axis: Optional[Axis] = None,
        attitude: Optional[Plane] = None,
        axis_angular_tolerance: numbers.Real = MIN_DISORIENTATION_TOLERANCE
    ) -> Tuple[Union[None, Ln], Error]:
        """
        Project a line to the segment according to various methods.

        :param line: the line to project.
        :param projection_method: the projection method.
        :param min_profile_distance: the minimum distance of the point from the section to be considered not lying on the section plane.
        :param map_axis: the optional axis to project the point to the section.
        :param attitude: the optional attitude along which to plot the point to the plane.
        :param axis_angular_tolerance: the minimum angular tolerance (in degrees) between axis and plane.
        """

        try:

            projected_line = Ln()

            for pt in line.pts():

                projected_pt, err = self.project_point(
                    point=pt,
                    projection_method=projection_method,
                    min_profile_distance=min_profile_distance,
                    map_axis=map_axis,
                    attitude=attitude,
                    axis_angular_tolerance=axis_angular_tolerance
                )

                if err:
                    return None, err

                if projected_pt is not None:
                    projected_line.add_pt(projected_pt)

            return projected_line, Error()

        except Exception as e:

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

    def parse_multipoints_intersections(
        self,
        intersections: Dict[Category, List[UnionPtSegment2D]],
        offset: numbers.Real = 0.0,
    ) -> List[CategorizedArrays]:
        """
        Parse the profile intersections for incorporation
        as elements in a geoprofile.

        :param intersections: the intersections
        :param offset: the offset to apply to the calculated values, to comply with LineProfilers
        :return: the list of resulting along-profile distances.
        """

        parsed_intersections = []

        for rec_id, ptsegm2d_list in intersections.items():
            intersections_arrays = []
            for ptsegm2d in ptsegm2d_list:
                if isinstance(ptsegm2d, (Point, Segment)):
                    value = self.pt_segm_along_profile_signed_s(ptsegm2d, offset)
                    intersections_arrays.append(value)
                else:
                    for subelem in ptsegm2d:
                        value = self.pt_segm_along_profile_signed_s(subelem, offset)
                        intersections_arrays.append(value)
            parsed_intersections.append(CategorizedArrays(rec_id, intersections_arrays))

        return parsed_intersections

    def parse_single_point_projection_result(self,
                                             source_pt: Point,
                                             intersection_point: Point,
                                             max_profile_distance,
                                             offset: numbers.Real = 0.0,
                                             ) -> Tuple[Union[None, PointTrace], Error]:
        """
        Parse a point intersection for plotting in a section.

        :param source_pt: the source point.
        :param intersection_point: the projected point in the section plane.
        :param max_profile_distance: the maximum allowed distance of the source point form the section plane.
        :param offset: the s distance offset to apply to the segment-calculated one.
        :return: the potential point trace for plotting, and the error status.
        """

        try:

            # distance between source and projected points

            dist = source_pt.distance(intersection_point)

            if dist > max_profile_distance:
                return None, Error()

            # horizontal spat_distance between projected point and profile start

            signed_distance_from_section_start = self.point_along_profile_signed_s(intersection_point.to2d())
            if signed_distance_from_section_start is None:
                return None, Error(
                    True,
                    caller_name(),
                    Exception("Error: <signed_distance_from_section_start> is None"),
                    traceback.format_exc()
                )
            elif not (0.0 <= signed_distance_from_section_start <= self.length()):
                return None, Error()

            # solution for current point

            parsed_point_intersection = PointTrace(
                s=signed_distance_from_section_start + offset,
                z=intersection_point.z,
                dist=dist,
            )

            return parsed_point_intersection, Error()

        except Exception as e:

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

    def parse_single_attitude_projection_result(
        self,
        attitude: Plane,
        attitude_intersection: Vect3D,
        source_pt: Point,
        intersection_point: Point,
        max_profile_distance,
        offset: numbers.Real = 0.0,
    ) -> Union[None, Error, PlaneTrace]:
        """
        Parse an attitude intersection for plotting in a section.

        :param attitude: the attitude plane.
        :param attitude_intersection: the section-attitude intersection vector.
        :param source_pt: the attitude source point.
        :param intersection_point: the projected point in the section plane.
        :param max_profile_distance: the maximum allowed distance of the source point form the section plane.
        :param offset: the s distance offset to apply to the segment-calculated one.
        :return: the potential plane trace for plotting, and the error status.
        """

        try:

            # calculate slope of geological plane onto section plane

            result, error = self.get_intersection_slope(attitude_intersection)
            if error:
                return error

            slope_radians, intersection_downward_sense = result

            # distance between source and projected points

            dist = source_pt.distance(intersection_point)

            if dist > max_profile_distance:
                return None

            # horizontal spat_distance between projected structural point and profile start

            signed_distance_from_section_start = self.point_along_profile_signed_s(intersection_point.to2d())
            if signed_distance_from_section_start is None:
                return Error(
                    True,
                    caller_name(),
                    Exception("Error: <signed_distance_from_section_start> is None"),
                    traceback.format_exc()
                )

            logger.debug("signed distance from section start = %s (length = %s)", signed_distance_from_section_start, self.length())

            if not (0.0 <= signed_distance_from_section_start <= self.length()):
                return None

            # solution for current structural point

            parsed_plane_trace = PlaneTrace(
                s=signed_distance_from_section_start + offset,
                z=intersection_point.z,
                slope_degr=degrees(slope_radians),
                down_sense=intersection_downward_sense,
                dist=dist,
                src_dip_dir=attitude.dipazim,
                src_dip_ang=attitude.dipang
            )

            return parsed_plane_trace

        except Exception as e:

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

    def parse_traces_with_attitudes_intersections(self,
        vect_points_intersections: List[Tuple[Vect3D, Point]],
        plane_orientation: Plane,
        offset: numbers.Real = 0.0,
    ) -> Tuple[Union[None, List[PlaneTrace]], Error]:
        """
        Parse an attitude intersection for plotting in a section.

        :param vect_points_intersections: the attitude plane.
        :param plane_orientation: the section-attitude intersection vector.
        :param offset: the s distance offset to apply to the segment-calculated one.
        :return: the potential plane trace for plotting, and the error status.
        """

        try:

            parsed_traces = []

            for attitude_intersection, intersection_point_2d in vect_points_intersections:

                # calculate slope of geological plane onto section plane

                result, error = self.get_intersection_slope(attitude_intersection)

                if error:
                    return None, error

                slope_radians, intersection_downward_sense = result

                # horizontal distance between structural point and profile start

                signed_distance_from_section_start = self.point_along_profile_signed_s(intersection_point_2d.to2d())

                if signed_distance_from_section_start is None:
                    return None, Error(
                        True,
                        caller_name(),
                        Exception("Error: <signed_distance_from_section_start> is None"),
                        traceback.format_exc()
                    )

                if not (0.0 <= signed_distance_from_section_start <= self.length()):
                    return None, Error()

                # solution for current structural point

                parsed_plane_trace = PlaneTrace2D(
                    s=signed_distance_from_section_start + offset,
                    slope_degr=degrees(slope_radians),
                    down_sense=intersection_downward_sense,
                    src_dip_dir=plane_orientation.dipazim,
                    src_dip_ang=plane_orientation.dipang
                )

                parsed_traces.append(parsed_plane_trace)

            return parsed_traces, Error()

        except Exception as e:

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


class LineProfiler:
    """
    Ln profilers.
    """

    def __init__(self,
                 line: Ln,
                 crs: Union[str, Crs],
                 s_breaks: Optional[np.ndarray] = None,
                 ):
        """
        Initialize the linear profiler.

        :param line: the source line for profiler creation
        :param crs: the coordinate reference system
        :param s_breaks: the optional array of distances breaks
        """

        self._line = line.clone()

        self._crs = crs if isinstance(crs, Crs) else Crs(crs)

        if s_breaks is None:
            self._breaks = np.array(line.accumulated_length_2d())
        else:
            self._breaks = s_breaks


    def clone(self) -> 'LineProfiler':
        """Clone the LineProfiler instance"""

        return LineProfiler(
            line=self._line,
            crs=self.crs,
            s_breaks=self.breaks,
        )

    def __iter__(self):
        """
        Return the segment profilers of the current line profilers.
        """

        return (SegmentProfiler(segment, self.wkt_crs) for segment in self.line)

    def direction_change_segment_profilers(
            self,
            deviation_tolerance: numbers.Real = COLLINEARITY_DEVIATION_TOLERANCE,
    ) -> List[Tuple['SegmentProfiler', numbers.Real]]:
        """
        Returns the segment profilers with consecutive collinear segments merged,
        each paired with its along-profile offset (the 2D distance from the line start).

        Profile lines are normally densified — one vertex per sampling step, to drape
        them on a DEM — but densification adds no direction changes, and projecting a
        point only depends on the direction of the segment it falls against. Iterating
        the line vertex by vertex therefore repeats the same projection against
        hundreds of collinear sub-segments, all but one of which are then discarded by
        the segment extent check. Collapsing each straight run to a single segment
        makes the projection cost depend on the number of direction changes rather
        than on the sampling step. Offsets are taken from the original per-segment
        lengths rather than from the merged chord, so the merge cannot accumulate
        drift along the profile; the only displacement it introduces comes from the
        chord's slightly different heading, bounded as described on
        COLLINEARITY_DEVIATION_TOLERANCE (centimetres at worst, and exactly zero on
        a run that is straight to the bit).

        Zero-length segments carry no direction and are skipped; they would otherwise
        yield an undefined vertical plane and abort the whole projection. A run is
        also never merged across a gap between two segments that do not join, so a
        multi-part line keeps its parts separate.

        :param deviation_tolerance: how far, in metres, a vertex may lie from the
            straight chord of its run and still be merged into it. Raise it if real
            profiles are not being collapsed (their vertices deviate more than this
            from a straight line); the resulting error scales as described on
            COLLINEARITY_DEVIATION_TOLERANCE.
        :return: the merged segment profilers, each with its along-profile offset.
        """

        # taken from the segments, i.e. the same source __iter__ reads, rather than
        # from a vertex accessor the line type may not expose
        segments = list(self.line)

        if not segments:
            return []

        starts = np.array([[seg.start_pt.x, seg.start_pt.y] for seg in segments], dtype=float)
        ends = np.array([[seg.end_pt.x, seg.end_pt.y] for seg in segments], dtype=float)

        lengths = np.hypot(ends[:, 0] - starts[:, 0], ends[:, 1] - starts[:, 1])

        # offset of each segment start along the profile, from the original
        # per-segment lengths, so merging cannot shift the s of a projected point
        offsets = np.concatenate(([0.0], np.cumsum(lengths)))

        merged = []
        run_first = None  # first segment of the run being built
        run_last = None  # last segment absorbed into it so far

        def _close_run():
            merged.append(
                (
                    SegmentProfiler(
                        Segment(
                            start_pt=Point(*starts[run_first]),
                            end_pt=Point(*ends[run_last]),
                        ),
                        self.wkt_crs,
                    ),
                    float(offsets[run_first]),
                )
            )

        for ndx in range(len(segments)):

            if lengths[ndx] == 0.0:
                continue

            if run_first is None:
                run_first, run_last = ndx, ndx
                continue

            joins_run = (
                float(np.hypot(*(starts[ndx] - ends[run_last]))) <= SEGMENT_CONTINUITY_TOLERANCE
            )

            # the run may only absorb this segment if every vertex it already
            # spans still lies on the chord that would then close it
            if joins_run and _max_deviation_from_chord(
                starts[run_first], ends[ndx], ends[run_first:ndx]
            ) <= deviation_tolerance:
                run_last = ndx
                continue

            _close_run()
            run_first, run_last = ndx, ndx

        if run_first is not None:
            _close_run()

        return merged

    def __repr__(self):
        """
        Representation of a profile instance.

        :return: the textual representation of the instance.
        """

        return f"LineProfiler(\n\tsource line = {self._line} - CRS = {self._crs})"

    @property
    def line(self) -> Ln:
        """
        Returns the profiler line.

        :return: the profiler line.
        """

        return self._line

    @property
    def breaks(self) -> np.ndarray:
        """
        Returns the along-profile cumulated segment lengths.

        :return: the along-profile cumulated segment lengths.
        """

        return self._breaks

    @property
    def crs(self) -> Crs:
        """Returns the coordinate reference system."""

        return self._crs

    @property
    def wkt_crs(self) -> Optional[str]:
        """Returns the WKT string of the CRS, or None if undefined."""

        return self._crs.to_wkt()

    def length(self) -> numbers.Real:
        """
        Returns the length of the profilers section.

        :return: length of the profilers section.
        """

        return self._line.length_2d()

    def densified_points(self,
                         sampling_distance) -> List[Point]:
        """
        Returns the list of densified 2D points when respecting source line original vertices.
        """

        segments = list(self.line.segments())
        pts = []

        for ndx, segment in enumerate(segments):

            densified_points = segment.densify_as_points2d(
                densify_distance=sampling_distance,
            )

            if ndx == len(segments) - 1:
                pts.extend(densified_points)
            else:
                pts.extend(densified_points[:-1])

        return pts

    def sample_grid(
            self,
            grid: Grid,
            sampling_distance: Optional[numbers.Real] = None,
            enforce_segment_vertices: bool = False,
    ) -> Tuple[Union[None, ZTrace], Error]:
        """
        Create profile from a grid, given a sampling distance.

        :param grid: the source grid.
        :param sampling_distance: the sampling distance
        :param enforce_segment_vertices: whether to constrain the profile to use segment vertices
        :return: the profile of the scalar variable stored in the grid.
        """

        try:

            if enforce_segment_vertices:

                # densify each segment independently (so that every original
                # vertex is guaranteed to be an exact sample point), then sample
                # the grid for the whole line in a single vectorized call
                if sampling_distance is None:
                    sampling_distance = grid.cellsize_mean

                segments = list(self._line.segments())

                xs_parts = []
                ys_parts = []
                s_parts = []
                s_offset = 0.0

                for ndx, segment in enumerate(segments):

                    seg_length = segment.length_2d()
                    start_pt = segment.start_pt.as_point2d()
                    versor = segment.as_versor2d()

                    n_steps = max(int(np.ceil(seg_length / sampling_distance)) - 1, 0)
                    local_s = np.arange(n_steps + 1) * sampling_distance
                    local_s = np.append(local_s, seg_length)

                    if ndx < len(segments) - 1:
                        # drop the segment-end point: it coincides with the
                        # next segment's start point and would else be duplicated
                        local_s = local_s[:-1]

                    xs_parts.append(start_pt.x + versor.x * local_s)
                    ys_parts.append(start_pt.y + versor.y * local_s)
                    s_parts.append(local_s + s_offset)

                    s_offset += seg_length

                x_values = np.concatenate(s_parts)
                pts_x = np.concatenate(xs_parts)
                pts_y = np.concatenate(ys_parts)
                y_values = grid.interpolate_bilinear_arrays(pts_x, pts_y)

                joined_ztrace = ZTrace(
                    s_array=x_values,
                    z_array=y_values,
                    s_breaks=self._breaks
                )

            else:

                # convert line to equally spaced points

                if sampling_distance is None:
                    sampling_distance = grid.cellsize_mean
                equally_spaced_points = self._line.densify_as_equally_spaced_points2d(
                    sample_distance=sampling_distance
                )

                # sample grid using points (vectorized: one batched call to the
                # grid instead of one grid.interpolate_bilinear call per point)
                num_points = len(equally_spaced_points)
                x_values = np.arange(num_points) * sampling_distance
                pts_x = np.fromiter((pt.x for pt in equally_spaced_points), dtype=np.float64, count=num_points)
                pts_y = np.fromiter((pt.y for pt in equally_spaced_points), dtype=np.float64, count=num_points)
                y_values = grid.interpolate_bilinear_arrays(pts_x, pts_y)

                joined_ztrace = ZTrace(
                    s_array=np.array(x_values),
                    z_array=np.array(y_values),
                    s_breaks=self._breaks
                )

            return joined_ztrace, Error()

        except Exception as e:

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

    def intersect_lines(
        self,
        lines: Dict[Category, List[PreparedLine]],
    ) -> List[Dict[Category, List[UnionPtSegment2D]]]:
        """
        Calculates the intersection with a set of lines.
        Note: the intersections are intended flat (in a 2D plane, not 3D).

        :param lines: a dict of category -> already-prepared lines
            (see prepare_line_for_intersection) to intersect profile with.
        :return: the possible intersections
        """

        line_profiler_intersections = []
        for segment_profiler in self:
            result = segment_profiler.intersect_lines(
                lines,
            )
            line_profiler_intersections.append(result)
        return line_profiler_intersections

    def intersect_lines_with_attitude(
        self,
        lines: List[Ln],
        plane_orientation: Plane
    ) -> Tuple[Union[None, List[PlaneTrace]], Error]:

        try:

            plane_traces = []
            segment_offset = 0.0

            for segment_profiler in self:

                result, err = segment_profiler.intersect_lines_with_attitude(
                    lines,
                    plane_orientation
                )

                if err:
                    logger.error("%r", err)
                    continue

                if result is None:
                    logger.warning("segment-profilers / line-attitude intersections are None")
                    continue

                vect_points_intersections = result  #List[Tuple[Vect3D, Point]]]

                result, err = segment_profiler.parse_traces_with_attitudes_intersections(
                    vect_points_intersections,
                    plane_orientation,
                    segment_offset,
                )

                if err:
                    logger.error("%r", err)
                    continue

                if result is None:
                    logger.warning("segment-profilers / line-attitude intersections are None")
                    continue

                plane_traces.extend(result)
                segment_offset += segment_profiler.length()

            return plane_traces, Error()

        except Exception as e:

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

    def intersect_polygon(self,
                          polygon: Polygon,
                          ) -> List[Segment]:
        """
        Calculates the attitude_intersection with a polygon.
        Note: the intersections are considered flat, i.e., in a 2D plane, not 3D.

        :param polygon: the polygon/multipolygon to intersect profile with
        :return: the possible intersections
        """

        segments = []

        prepared_polygon = prepare_polygon_for_intersection(polygon)

        for segment_profiler in self:
            segments.extend(segment_profiler.intersect_prepared_polygon(prepared_polygon))

        return segments

    def intersect_polygons(self,
                           polygons: Dict[Category, PreparedPolygons]
                           ) -> List[Dict[Category, List[Segment]]]:
        """
        Calculates the attitude_intersection with a list of polygon.
        Note: the intersections are intended flat (in a 2D plane, not 3D).

        :param polygons: a dict of category -> already-prepared, bounding-box
            indexed polygons (see PreparedPolygons) to intersect profile with.
            Preparing them once, outside the per-profile/per-segment loop that
            calls this method, avoids rebuilding the same boundary segments
            and bounding boxes over and over.
        :return: the possible intersections
        """

        lineprofiler_intersections = []

        for segment_profiler in self:

            segment_intersections = defaultdict(list)

            segment_bbox = segment_profiler.segment.bounding_box()

            for polygons_id, prepared_polygons in polygons.items():

                if prepared_polygons.num_source_polygons == 0:
                    continue

                # the category is reported (possibly as empty) as it was when
                # every polygon of the category was tested one by one

                category_intersections = segment_intersections[polygons_id]

                for prepared_polygon in prepared_polygons.candidates(segment_bbox):
                    category_intersections.extend(
                        segment_profiler.intersect_prepared_polygon(
                            prepared_polygon,
                            segment_bbox,
                        )
                    )

            lineprofiler_intersections.append(segment_intersections)

        return lineprofiler_intersections

    def vector_offset(self,
                      vect: Vect2D
                      ) -> 'LineProfiler':
        """
        Returns a new LineProfiler instance, horizontally offset by the
        provided vector horizontal components.
        """

        return LineProfiler(
            line=self._line.shift_horiz(vect.x, vect.y),
            crs=self.crs,
            s_breaks=self.breaks
        )

    def right_parallel_offset(self,
                              offset: numbers.Real
                              ) -> Union[None, 'LineProfiler']:
        """
        Returns a copy of the current linear profilers, offset to the right by the provided offset distance.

        :param offset: the lateral offset to apply to create the new LinearProfiler.
        :return: the offset linear profilers.
        """

        try:

            segm = self._line.as_segment()
            if segm is None:
                raise Exception(f"Unable to convert line to segment")
            vect = segm.right_norm_vers2d().scale(offset)
            return self.vector_offset(
                vect=vect,
            )

        except Exception as e:

            logger.error(traceback.format_exc())
            return None

    def left_parallel_offset(self,
                             offset: numbers.Real
                             ) -> Union[None, 'LineProfiler']:
        """
        Returns a copy of the current linear profilers, offset to the left by the provided offset distance.

        :param offset: the lateral offset to apply to create the new LinearProfiler.
        :return: the offset linear profilers.
        """

        try:

            segm = self._line.as_segment()
            if segm is None:
                raise Exception(f"Unable to convert line to segment")
            vect = segm.left_norm_vers2d().scale(offset)
            return self.vector_offset(
                vect=vect,
            )

        except Exception as e:

            logger.error(traceback.format_exc())
            return None

    def parse_multipoints_intersections(
            self,
            intersections: List[Dict[Category, List[UnionPtSegment2D]]],
    ) -> List[CategorizedArrays]:
        """
        Parse the profile intersections for incorporation
        as elements in a geoprofile.

        :param intersections: the intersections
        :return: the list of resulting along-profile distances.
        """

        array_lists = []

        offset = 0.0

        for ndx, (segment_profiler, segment_intersection_points) in enumerate(zip(self, intersections), start=1):
            #print(f"\nSegment profiler level: segment_intersection_points -> {segment_intersection_points}")
            array_lists.append(segment_profiler.parse_multipoints_intersections(segment_intersection_points, offset))
            offset += segment_profiler.length()

        return array_lists


class Profilers(list):
    """
    ParalleProfilers.
    """

    def __init__(self,
                 profilers: Union[None, List[LineProfiler]] = None,
                 ):

        if profilers is None:

            super().__init__()

        else:

            super().__init__(profilers)

    @classmethod
    def make_parallel_from_line(cls,
                                src_trace: Ln,
                                src_crs: Crs,
                                num_profiles: numbers.Integral = 1,
                                offset: numbers.Real = 50,  # horizontal distance unit, e.g. meters
                                profiles_arrangement: str = "central",  # one of: "left", "central", "right"
                                ):

        if num_profiles < 1:
            raise Exception("Profilers number must be >= 1")

        if profiles_arrangement == "central" and num_profiles % 2 != 1:
            raise Exception("When profilers arrangement is 'central' profilers number must be odd")

        if profiles_arrangement == "central":

            side_profs_num = num_profiles // 2
            num_left_profs = num_right_profs = side_profs_num

        elif profiles_arrangement == "left":

            num_left_profs = num_profiles - 1
            num_right_profs = 0

        else:

            num_right_profs = num_profiles - 1
            num_left_profs = 0

        profilers = []

        base_profiler = LineProfiler(
            line=src_trace,
            crs=src_crs,
        )

        for i in range(num_left_profs, 0, -1):
            current_offset = offset * i

            profilers.append(base_profiler.left_parallel_offset(offset=current_offset))

        profilers.append(base_profiler.clone())

        for i in range(1, num_right_profs + 1):
            current_offset = offset * i

            profilers.append(base_profiler.right_parallel_offset(offset=current_offset))

        return cls(profilers)

    @classmethod
    def make_parallel_from_lines(cls,
        src_traces: List[ProfileBaseLineParameters],
        num_profiles: numbers.Integral = 1,
        offset: numbers.Real = 50,  # horizontal distance unit, e.g. meters
        profiles_arrangement: str = "central",  # one of: "left", "central", "right"
    ):

        profilers = Profilers()

        for trace_params in src_traces:

            src_trace = trace_params.baseline
            src_crs = Crs.from_wkt(trace_params.wkt_crs)

            profilers_to_add = Profilers.make_parallel_from_line(
                src_trace=src_trace,
                src_crs=src_crs,
                num_profiles=num_profiles,
                offset=offset,
                profiles_arrangement=profiles_arrangement,
            )

            profilers.extend(profilers_to_add)

        return profilers

    def __repr__(self) -> str:
        """
        Represents a linear profilers set.

        :return: the textual representation of the linear profilers set.
        """

        inner_profilers = "\n".join([repr(profiler) for profiler in self])
        return f"Profilers([\n{inner_profilers}]\n)"

    @property
    def lines(self) -> List[Ln]:
        """
        Returns the profiles lines.

        :return: the parallel profilers lines.
        """

        return [profiler.line for profiler in self]

    def sample_grid(
            self,
            grid: Grid,
            sampling_distance: Optional[numbers.Real] = None
    ) -> Tuple[Union[None, ZTraces], Error]:
        """
        Create profile from a grid.

        :param grid: the source grid.
        :param sampling_distance: the grid sampling distance.
        :return: list of profiles of the scalar variable stored in the grid and an error status.
        """

        try:

            z_traces = []

            for line_profiler in self:

                z_trace, err = line_profiler.sample_grid(
                    grid=grid,
                    sampling_distance=sampling_distance
                )

                if err:
                    return None, err
                z_traces.append(z_trace)

            return ZTraces.fromProfiles(z_traces), Error()

        except Exception as e:

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

    def sample_grid_as_zarrays(
            self,
            grid: Grid,
            sampling_distance: Optional[numbers.Real] = None
    ) -> Tuple[Union[None, np.ndarray], Error]:
        """
        Create profile from a grid as ZArray.

        :param grid: the source grid.
        :param sampling_distance: the grid sampling distance.
        :return: list of profiles of the scalar variable stored in the grid and an error status.
        """

        try:

            topo_profiles = []

            for line_profiler in self:

                z_array, err = line_profiler.sample_grid(
                    grid=grid,
                    sampling_distance=sampling_distance
                )

                if err:
                    return None, err
                topo_profiles.append(z_array)

            return ZTraces.fromProfiles(topo_profiles), Error()

        except Exception as e:

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

    def intersect_lines(self,
        lines: Dict[Category, List[Ln]],
        progress_callback: Optional[Callable[[numbers.Integral, numbers.Integral], None]] = None,
    ) -> List[List[Dict[Category, List[UnionPtSegment2D]]]]:
        """
        :param lines: the lines to intersect, by category.
        :param progress_callback: optional callable invoked as
            progress_callback(done, total) after each of the 'total' parallel
            profiles is intersected, e.g. to drive a progress bar. Errors
            raised by the callback itself are not caught.
        """

        try:

            # Prepare each line's bounding box and segment list once here,
            # instead of once per (profile x profile-segment x line) triple:
            # with several parallel profiles, each made of many densified
            # segments, the same lines would otherwise be re-scanned
            # (bounding boxes, segment lists) many thousands of times over.
            prepared_lines = {
                category: [prepare_line_for_intersection(line) for line in lines_geoms]
                for category, lines_geoms in lines.items()
            }

            num_profilers = len(self)
            profilers_intersections = []

            for ndx, line_profiler in enumerate(self):
                line_profiler_intersections = line_profiler.intersect_lines(
                    prepared_lines,
                )
                profilers_intersections.append(line_profiler_intersections)

                if progress_callback is not None:
                    progress_callback(ndx + 1, num_profilers)

            return profilers_intersections

        except Exception:

            raise

    def intersect_lines_with_attitudes(self,
        lines: Dict[Category, Tuple[Plane, List[Ln]]]
    ) -> Tuple[Union[None, List[Dict[Category, List[PlaneTrace]]]], Error]:

        try:

            profiles_plane_traces = []

            for line_profiler in self:

                line_profile_plane_traces = defaultdict(list)

                for lines_id, unpacked_data in lines.items():

                    if isinstance(unpacked_data, List):

                        for plane_orientation, rec_id_lines in unpacked_data:

                            plane_traces, err = line_profiler.intersect_lines_with_attitude(
                                lines=rec_id_lines,
                                plane_orientation=plane_orientation)
                            if err:
                                repr(err)
                                continue
                            if plane_traces:
                                line_profile_plane_traces[lines_id] = plane_traces
                            else:
                                logger.warning("plane traces for record id %s is %s", lines_id, plane_traces)

                    else:

                        plane_orientation, rec_id_lines = unpacked_data

                        plane_traces, err = line_profiler.intersect_lines_with_attitude(
                            lines=rec_id_lines,
                            plane_orientation=plane_orientation)
                        if err:
                            repr(err)
                            continue
                        if plane_traces:
                            line_profile_plane_traces[lines_id] = plane_traces
                        else:
                            logger.warning("plane traces for record id %s is %s", lines_id, plane_traces)

                profiles_plane_traces.append(line_profile_plane_traces)

            return profiles_plane_traces, Error()

        except Exception as e:

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

    def intersect_polygons(self,
                           polygons: Dict[Category, List[Polygon]],
                           progress_callback: Optional[Callable[[numbers.Integral, numbers.Integral], None]] = None,
                           ) -> List[Dict[Category, List[Segment]]]:
        """
        Calculates the profilers intersections with a set of polygon/multipolygon.
        Note: the intersections are intended flat (in a 2D plane, not 3D).

        :param polygons: the set of polygon/multipolygon to intersect profile with
        :param progress_callback: optional callable invoked as
            progress_callback(done, total) after each of the 'total' parallel
            profiles is intersected, e.g. to drive a progress bar. Errors
            raised by the callback itself are not caught.
        :return: the possible intersections
        """

        # Prepare each category's polygons (boundary segments, bounding boxes,
        # ring arrays, bounding-box index) once here, instead of once per
        # (profile x profile-segment x polygon) triple: with several parallel
        # profiles and an extensive polygon layer, the same polygons would
        # otherwise be entirely rebuilt and re-walked many thousands of times.
        prepared_polygons = {
            category: PreparedPolygons(polygons_geometries)
            for category, polygons_geometries in polygons.items()
        }

        num_profilers = len(self)
        intersecting_segments = []

        for ndx, line_profiler in enumerate(self):
            intersecting_segments.append(line_profiler.intersect_polygons(prepared_polygons))

            if progress_callback is not None:
                progress_callback(ndx + 1, num_profilers)

        return intersecting_segments

    def parse_multipoints_intersections(
            self,
            intersections_list: List[List[Dict[Category, List[UnionPtSegment2D]]]]
    ) -> Optional[List[List[CategorizedArrays]]]:
        """
        Parse the profile intersections for incorporation
        as elements in a geoprofile.

        :param intersections_list: the intersections for the various parallel profiles
        :return: the optional list of resulting along-profile distances.
        """

        if len(self) != len(intersections_list):
            raise Exception(f"Number of profilers (#{len(self)}) is different from number in intersections lists (#{len(intersections_list)}) ")

        array_lists = []

        for line_profiler, intersections in zip(self, intersections_list):
            #print(f"\nLine profiler level: intersections -> {intersections}")
            array_lists.append(line_profiler.parse_multipoints_intersections(intersections))

        return array_lists


@singledispatch
def sample_grid_as_pts3d(
    profilers,
    grid,
    sampling_distance,
):

    return NotImplemented


@sample_grid_as_pts3d.register(list)
def _(
    profilers: Profilers,
    grid: Grid,
    sampling_distance: Optional[numbers.Real] = None
) -> Tuple[Union[None, List[List[Point]]], Error]:
    """
    Create profile from a grid.

    :param: profilers: the profilers to use
    :param grid: the source grid.
    :param sampling_distance: the grid sampling distance.
    :return: list of profiles of the scalar variable stored in the grid and an error status.
    """

    try:

        profiles_3d = []

        for line_profiler in profilers:

            pts_3d, err = sample_grid_as_pts3d(
                line_profiler,
                grid=grid,
                sampling_distance=sampling_distance
            )

            if err:
                return None, err

            profiles_3d.append(pts_3d)

        logger.debug("there are %d line profilers in profiles_3d", len(profiles_3d))

        return profiles_3d, Error()

    except Exception as e:

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


@sample_grid_as_pts3d.register(LineProfiler)
def _(
    line_profiler: LineProfiler,
    grid: Grid,
    sampling_distance: Optional[numbers.Real] = None,
    enforce_segment_vertices: bool = True,
    atol: numbers.Real = MIN_POINT_POS_DIFF,
) -> Tuple[Union[None, List[Point]], Error]:
    """
    Create profile from a grid, given a sampling distance.

    :param: profiler: the profiler to use
    :param grid: the source grid.
    :param sampling_distance: the sampling distance
    :param enforce_segment_vertices: whether to constrain the profile to use segment vertices
    :param atol: the absolute tolerance between consecutive points distance.
    :return: the profile of the scalar variable stored in the grid.
    """

    try:

        if enforce_segment_vertices:

            total_points = []

            for segment_profiler in line_profiler:

                points_3d, err = sample_grid_as_pts3d(
                    segment_profiler,
                    grid,
                    sampling_distance
                )

                if err:
                    return None, err

                total_points.extend(points_3d)

        else:

            # convert line to equally spaced points

            if sampling_distance is None:
                sampling_distance = grid.cellsize_mean
            equally_spaced_points = line_profiler._line.densify_as_equally_spaced_points2d(
                sample_distance=sampling_distance
            )

            # sample grid values for all points in a single vectorized call,
            # instead of one grid.interpolate_bilinear_point_with_nan call
            # (and, when CRS differ, one Transformer construction) per point
            pts_crs = line_profiler.crs
            num_pts = len(equally_spaced_points)
            orig_x = np.fromiter((pt.x for pt in equally_spaced_points), dtype=np.float64, count=num_pts)
            orig_y = np.fromiter((pt.y for pt in equally_spaced_points), dtype=np.float64, count=num_pts)

            if grid.crs != pts_crs:
                grid_x, grid_y = project_xy_arrays(
                    orig_x, orig_y,
                    src_crs=pts_crs.to_wkt(),
                    dest_crs=grid.wkt_crs,
                )
                if grid_x is None:
                    raise Exception("Unable to project points to grid CRS")
            else:
                grid_x, grid_y = orig_x, orig_y

            z_values = grid.interpolate_bilinear_arrays(grid_x, grid_y)

            total_points = [
                Point(x, y, z) for x, y, z in zip(orig_x, orig_y, z_values)
            ]

        total_points = remove_coincident_points_2d_fast(total_points, atol=atol)

        return total_points, Error()

    except Exception as e:

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


@sample_grid_as_pts3d.register(SegmentProfiler)
def _(
    segment_profiler: SegmentProfiler,
    grid: Grid,
    sampling_distance: Optional[numbers.Real] = None
) -> Tuple[Union[None, List[Point]], Error]:
    """
    Create a point-based profile from one grid.

    :param grid: the source grid.
    :param sampling_distance: the distance along to which to sample the grid.
    :return: the optional point-like profile along the grid and an error status.
    """

    try:

        pts_crs_code = segment_profiler.crs_code
        logger.debug("pts_crs_code: %s", pts_crs_code)
        logger.debug("grid.wkt_crs: %s", grid.wkt_crs)

        if sampling_distance is None:
            sampling_distance = grid.cellsize_mean

        if segment_profiler.crs_code != grid.wkt_crs:
            projected_segment = project(
                segment_profiler.segment,
                source_crs=segment_profiler.crs_code,
                dest_crs=grid.wkt_crs,
            )
            if projected_segment is None:
                raise Exception("Unable to project segment profiler")
        else:
            projected_segment = segment_profiler.segment

        # densify the segment into equally spaced points (vectorized: a closed-form
        # array computation instead of the point-by-point 'densify_as_pts2d' loop)
        seg_length = projected_segment.length_2d()
        seg_start_pt = projected_segment.start_pt.as_point2d()
        seg_versor = projected_segment.as_versor2d()

        n_steps = max(int(np.ceil(seg_length / sampling_distance)) - 1, 0)
        local_s = np.arange(n_steps + 1) * sampling_distance
        local_s = np.append(local_s, seg_length)

        pts_x = seg_start_pt.x + seg_versor.x * local_s
        pts_y = seg_start_pt.y + seg_versor.y * local_s

        # sample grid values for all points in a single vectorized call,
        # instead of one grid.interpolate_bilinear_point_with_nan call per point
        z_values = grid.interpolate_bilinear_arrays(pts_x, pts_y)

        projected_pts_3d = [
            Point(x, y, z) for x, y, z in zip(pts_x, pts_y, z_values)
        ]

        if segment_profiler.crs_code != grid.wkt_crs:
            pts_3d = project(
                projected_pts_3d,
                source_crs=grid.wkt_crs,
                dest_crs=segment_profiler.crs_code,
            )
        else:
            pts_3d = projected_pts_3d

        return pts_3d, Error()

    except Exception as e:

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