from array import array
from copy import deepcopy

from geogst.core.geometries.shapes.lines import *
from geogst.core.geometries.shapes.multilines import MultiLine
from geogst.core.utils.types import *

# from: https://stackoverflow.com/questions/47004506/check-if-a-numpy-array-is-sorted
# answer by: luca
is_sorted = lambda a: np.all(a[:-1] <= a[1:])


class PointTrace:
    """
    Represents a point projected onto a vertical section.
    """

    def __init__(
        self,
        s: numbers.Real,
        z: numbers.Real,
        dist: numbers.Real,
    ):
        """
        :param s: the signed horizontal plane location along the profile.
        :param z: the height of the plane point location in the profile.
        :param dist: the distance between the source point and the point projection on the profile.
        """

        self.s = s
        self.z = z
        self.dist = dist

    def __repr__(self) -> str:
        """
        Creates the representation of a PointTrace instance.

        :return: the representation of a PointTrace instance.
        """

        return f"PointTrace(s={self.s:.2f}, z={self.z:.2f}, dist={self.dist:.2f})"

    def shift_by_s(
        self,
        shift_value: numbers.Real) -> 'PointTrace':
        """
        Shifts a point trace by the given value along the profile direction.

        :param shift_value: the shift value to apply.
        :returns: the shifted point trace.
        """

        return PointTrace(
            self.s + shift_value,
            self.z,
            self.dist,
        )


class FocalMechanismTrace:
    """
    Represents an earthquake focal mechanism projected onto a vertical section.
    Stores both the original strike/dip/rake (for a map-view beachball glyph) and the
    local along-profile azimuth (needed to optionally rotate the mechanism into the true
    section plane, i.e., as it would appear viewed horizontally from the section side).
    """

    def __init__(
        self,
        s: numbers.Real,
        z: numbers.Real,
        dist: numbers.Real,
        strike: numbers.Real,
        dip: numbers.Real,
        rake: numbers.Real,
        profile_azimuth: numbers.Real,
    ):
        """
        :param s: the signed horizontal plane location along the profile.
        :param z: the height of the hypocenter location in the profile.
        :param dist: the distance between the source point and the point projection on the profile.
        :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 (Aki & Richards, 1980) of the focal mechanism, in degrees.
        :param profile_azimuth: the azimuth (degrees from North) of the local along-profile
            direction (the segment the point was projected onto), used for the section view.
        """

        self.s = s
        self.z = z
        self.dist = dist
        self.strike = strike
        self.dip = dip
        self.rake = rake
        self.profile_azimuth = profile_azimuth

    def __repr__(self) -> str:
        """
        Creates the representation of a FocalMechanismTrace instance.

        :return: the representation of a FocalMechanismTrace instance.
        """

        return f"FocalMechanismTrace(s={self.s:.2f}, z={self.z:.2f}, dist={self.dist:.2f}, " \
               f"strike={self.strike:.2f}, dip={self.dip:.2f}, rake={self.rake:.2f}, " \
               f"profile_azimuth={self.profile_azimuth:.2f})"

    def shift_by_s(
        self,
        shift_value: numbers.Real) -> 'FocalMechanismTrace':
        """
        Shifts a focal mechanism trace by the given value along the profile direction.

        :param shift_value: the shift value to apply.
        :returns: the shifted focal mechanism trace.
        """

        return FocalMechanismTrace(
            self.s + shift_value,
            self.z,
            self.dist,
            self.strike,
            self.dip,
            self.rake,
            self.profile_azimuth,
        )


class PlaneTrace:
    """
    Represents a plane projected onto/intersected by a vertical section,
    expressed by a selected location of the plane onto the section,
    as well as its intersecting and original attitude.
    """

    def __init__(
        self,
        s: numbers.Real,
        z: numbers.Real,
        slope_degr: numbers.Real,
        down_sense: str,
        dist: numbers.Real,
        src_dip_dir: numbers.Real,
        src_dip_ang: numbers.Real
    ):
        """
        :param s: the signed horizontal plane location along the profile.
        :param z: the height of the plane point location in the profile.
        :param slope_degr: the slope of the plane attitude in the profile. Unit is degrees.
        :param down_sense: downward sense, to the right or to the profile left.
        :param dist: the distance between the plane attitude point and the point projection on the profile.
        :param src_dip_dir: plane source dip direction.
        :param src_dip_ang: plane source dip angle.
        """

        self.s = s
        self.z = z
        self.slope_degr = slope_degr
        self.down_sense = down_sense
        self.dist = dist
        self.src_dip_dir = src_dip_dir
        self.src_dip_ang = src_dip_ang

    def __repr__(self) -> str:
        """
        Creates the representation of a PlaneTrace instance.

        :return: the representation of a PlaneTrace instance.
        """

        return f"PlaneTrace(s={self.s:.2f}, z={self.z:.2f}, " \
               f"dip_angle={self.slope_degr:.2f}, down_sense={self.down_sense}, " \
               f"dist={self.dist:.2f}, src_dip_dir={self.src_dip_dir:.2f}, src_dip_ang={self.src_dip_ang:.2f})"

    def shift_by_s(
        self,
        shift_value: numbers.Real) -> 'PlaneTrace':
        """
        Shifts a plane attitude by the given value along the profile direction.

        :param shift_value: the shift value to apply.
        :returns: the shifted plane attitude.
        """

        return PlaneTrace(
            self.s + shift_value,
            self.z,
            self.slope_degr,
            self.down_sense,
            self.dist,
            self.src_dip_dir,
            self.src_dip_ang,
        )

    def create_segment_for_plot(
            self,
            profile_length: numbers.Real,
            vertical_exaggeration: numbers.Real = 1,
            segment_scale_factor: numbers.Real = 50.0):
        """

        :param profile_length:
        :param vertical_exaggeration:
        :param segment_scale_factor: the (inversely proportional) scale factor controlling the plane_attitude segment length in the plot.
        :return:
        """

        ve = float(vertical_exaggeration)

        z0 = self.z

        h_dist = self.s
        slope_rad = radians(self.slope_degr)
        intersection_downward_sense = self.down_sense
        length = profile_length / segment_scale_factor

        s_slope = sin(float(slope_rad))
        c_slope = cos(float(slope_rad))

        if c_slope == 0.0:
            height_corr = length / ve
            structural_segment_s = [h_dist, h_dist]
            structural_segment_z = [z0 + height_corr, z0 - height_corr]
        else:
            t_slope = s_slope / c_slope
            width = length * c_slope

            length_exag = width * sqrt(1 + ve*ve * t_slope*t_slope)

            corr_width = width * length / length_exag
            corr_height = corr_width * t_slope

            structural_segment_s = [h_dist - corr_width, h_dist + corr_width]
            structural_segment_z = [z0 + corr_height, z0 - corr_height]

            if intersection_downward_sense == "left":
                structural_segment_z = [z0 - corr_height, z0 + corr_height]

        return structural_segment_s, structural_segment_z


class PlaneTrace2D:
    """
    Represents a plane intersected by a vertical section,
    expressed by a selected location of the plane onto the section,
    as well as its intersecting and original attitude.
    """

    def __init__(
        self,
        s: numbers.Real,
        slope_degr: numbers.Real,
        down_sense: str,
        src_dip_dir: numbers.Real,
        src_dip_ang: numbers.Real
    ):
        """
        :param s: the signed horizontal plane location along the profile.
        :param slope_degr: the slope of the plane attitude in the profile. Unit is degrees.
        :param down_sense: downward sense, to the right or to the profile left.
        :param src_dip_dir: plane source dip direction.
        :param src_dip_ang: plane source dip angle.
        """

        self.s = s
        self.slope_degr = slope_degr
        self.down_sense = down_sense
        self.src_dip_dir = src_dip_dir
        self.src_dip_ang = src_dip_ang

    def __repr__(self) -> str:
        """
        Creates the representation of a PlaneTrace2D instance.

        :return: the representation of a PlaneTrace2D instance.
        """

        return f"PlaneTrace2D(s={self.s:.2f}, " \
               f"dip_angle={self.slope_degr:.2f}, down_sense={self.down_sense}, " \
               f"src_dip_dir={self.src_dip_dir:.2f}, src_dip_ang={self.src_dip_ang:.2f})"

    def shift_by_s(
        self,
        shift_value: numbers.Real) -> 'PlaneTrace2D':
        """
        Shifts a plane attitude by the given value along the profile direction.

        :param shift_value: the shift value to apply.
        :returns: the shifted plane attitude.
        """

        return PlaneTrace2D(
            self.s + shift_value,
            self.slope_degr,
            self.down_sense,
            self.src_dip_dir,
            self.src_dip_ang,
        )

    def create_segment_for_plot(
            self,
            profile_length: numbers.Real,
            z0: numbers.Real,
            vertical_exaggeration: numbers.Real = 1,
            segment_scale_factor: numbers.Real = 50.0):
        """

        :param profile_length:
        :param z0: the elevation of the observation.
        :param vertical_exaggeration:
        :param segment_scale_factor: the (inversely proportional) scale factor controlling the plane_attitude segment length in the plot.
        :return:
        """

        ve = float(vertical_exaggeration)

        h_dist = self.s
        slope_rad = radians(self.slope_degr)
        intersection_downward_sense = self.down_sense
        length = profile_length / segment_scale_factor

        s_slope = sin(float(slope_rad))
        c_slope = cos(float(slope_rad))

        if c_slope == 0.0:
            height_corr = length / ve
            structural_segment_s = [h_dist, h_dist]
            structural_segment_z = [z0 + height_corr, z0 - height_corr]
        else:
            t_slope = s_slope / c_slope
            width = length * c_slope

            length_exag = width * sqrt(1 + ve*ve * t_slope*t_slope)

            corr_width = width * length / length_exag
            corr_height = corr_width * t_slope

            structural_segment_s = [h_dist - corr_width, h_dist + corr_width]
            structural_segment_z = [z0 + corr_height, z0 - corr_height]

            if intersection_downward_sense == "left":
                structural_segment_z = [z0 - corr_height, z0 + corr_height]

        return structural_segment_s, structural_segment_z


class FaultTrace:
    """
    Represents a fault projected onto a vertical section,
    expressed by a selected location of the plane onto the section,
    as well as the fault attitude.
    """

    def __init__(
        self,
        s: numbers.Real,
        z: numbers.Real,
        fault: Fault
    ):
        """
        :param s: the signed horizontal plane location along the profile.
        :param z: the height of the plane point location in the profile.
        :param fault: the projected fault.
        """

        self.s = s
        self.z = z
        self.fault = fault

    def __repr__(self) -> str:
        """
        Creates the representation of a FaultTrace instance.

        :return: the representation of a FaultTrace instance.
        :rtype: str.
        """

        return f"FaultTrace(s={self.s:.2f}, z={self.z:.2f}, fault={self.fault})"

    def shift_by_s(
        self,
        shift_value: numbers.Real) -> 'FaultTrace':
        """
        Shifts a fault projection by the given value along the profile direction.

        :param shift_value: the shift value to apply.
        :returns: the shifted fault projection.
        """

        return FaultTrace(
            self.s + shift_value,
            self.z,
            self.fault,
        )

class ZTrace:
    """
    Represent an s-z array pair (i.e., a single z value for each s value).
    S values should be sorted and should start from zero.
    """

    def __init__(self,
                 s_array: Union[array, np.ndarray],
                 z_array: Union[array, np.ndarray],
                 s_breaks: Union[None, array, np.ndarray] = None,
                 point_breaks: Optional[List[Point]] = None,
                 wkt_crs: Optional[str] = None,
                 rec_id: Optional[numbers.Integral] = None,
                 ):
        """
        Constructor.

        :param s_array: the s values array.
        :param z_array: the z values array.
        :param s_breaks: the internal s breaks.
        :param point_breaks: the optional real-world points (2D) at the s_breaks positions.
        :param wkt_crs: the optional WKT of the coordinate reference system of point_breaks.
        :param rec_id: an optional identifier for the profile (e.g. a database record id).
        """

        check_type(s_array, "S array", (array, np.ndarray))
        if not np.all(np.isfinite(s_array)):
            raise Exception("S array values must all be finite")
        if s_array[0] != 0.0:
            raise Exception("First value of S array should be zero")

        if not is_sorted(s_array):
            raise Exception("S array must be already sorted")

        check_type(z_array, "Z array", (array, np.ndarray))

        if len(z_array) != len(s_array):
            raise Exception("Z array must have the same length as S array")

        if s_breaks is not None:
            check_type(s_breaks, "S breaks", (array, np.ndarray))
            if s_breaks[0] != 0.0:
                raise Exception("First element of S breaks must always be zero")

        self._s = np.array(s_array, dtype=np.float64)
        self._z = np.array(z_array, dtype=np.float64)
        self._s_breaks = np.array(s_breaks, dtype=np.float64)  if s_breaks is not None else np.array([s_array[0], s_array[-1]])
        self._point_breaks = list(point_breaks) if point_breaks is not None else None
        self._wkt_crs = wkt_crs
        self._rec_id = rec_id

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

        return ZTrace(
            s_array=np.copy(self._s),
            z_array=np.copy(self._z),
            s_breaks=None if self._s_breaks is None else np.copy(self._s_breaks),
            point_breaks=None if self._point_breaks is None else [pt.clone() for pt in self._point_breaks],
            wkt_crs=self._wkt_crs,
            rec_id=self._rec_id,
        )

    def s_arr(self) -> np.ndarray:
        """
        Return the s array.

        :return: the s array.
        :rtype: array.
        """

        return self._s

    def z_arr(self) -> np.ndarray:
        """
        Return the z array.

        :return: the z array.
        """

        return self._z

    def s_breaks(self) -> Optional[np.ndarray]:

        return self._s_breaks

    def point_breaks(self) -> Optional[List[Point]]:
        """
        Return the optional real-world points (2D) at the s_breaks positions.
        """

        return self._point_breaks

    def wkt_crs(self) -> Optional[str]:
        """
        Return the optional WKT of the coordinate reference system of point_breaks.
        """

        return self._wkt_crs

    def rec_id(self) -> Optional[numbers.Integral]:
        """
        Return the optional profile identifier (e.g. a database record id).
        """

        return self._rec_id

    def __repr__(self) -> str:
        """
        Representation of a topographic profile instance.

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

        return f"ZTrace with {len(self.s_arr())} nodes\ns = {self._s},\nz = {self._z}"

    @classmethod
    def from_mline(cls,
                   mline: Union[Ln, MultiLine],
                   wkt_crs: Optional[str] = None,
                   rec_id: Optional[numbers.Integral] = None,
                   ):

        if isinstance(mline, MultiLine):
            line = mline.to_line()
        elif isinstance(mline, Ln):
            line = mline
        else:
            raise Exception(f"Invalid line type: {type(mline)}")

        s_array = np.array(line.accumulated_length_2d())
        z_array = line.z_array()

        breaks_result = line.breaks()
        if breaks_result is not None:
            break_points_arr, break_distances = breaks_result
            point_breaks = [Point(*coord[:2]) for coord in break_points_arr]
        else:
            break_distances = None
            point_breaks = None

        return cls(
            s_array=s_array,
            z_array=z_array,
            s_breaks=break_distances,
            point_breaks=point_breaks,
            wkt_crs=wkt_crs,
            rec_id=rec_id,
        )

    def num_steps(self) -> numbers.Integral:
        """
        Return the number of steps in the array pair.

        :return: number of steps in the array pair.
        """

        return len(self._s)

    def s(self,
          ndx: numbers.Integral
          ) -> numbers.Real:
        """
        Returns the x value with the index ndx.

        :param ndx: the index in the x array
        :return: the s value corresponding to the ndx index
        """

        return self._s[ndx]

    def z(self,
          ndx: numbers.Integral
          ) -> numbers.Real:
        """
        Returns the z value with the index ndx.

        :param ndx: the index in the y array
        :return: the z value corresponding to the ndx index
        """

        return self._z[ndx]

    def s_min(self) -> numbers.Real:
        """
        Returns the minimum s value.

        :return: the minimum s value.
        """

        return float(np.nanmin(self._s))

    def s_max(self) -> numbers.Real:
        """
        Returns the maximum s value.

        :return: the maximum s value.
        """

        return float(np.nanmax(self._s))

    def z_min(self) -> numbers.Real:
        """
        Returns the minimum z value.

        :return: the minimum z value.
        """

        return float(np.nanmin(self._z))

    def z_max(self) -> numbers.Real:
        """
        Returns the maximum z value.

        :return: the maximum z value.
        :rtype: numbers.Real.
        """

        return float(np.nanmax(self._z))

    def s_length(self) -> numbers.Real:
        """
        Returns the geographic length of the profile.

        :return: length of profile.
        :rtype: numbers.Real.
        """

        return self._s[-1]

    def find_index_ge(self,
                      s_val: numbers.Real) -> Union[None, numbers.Integral]:
        """

        Examples:
          >>> p = ZTrace(array('d', [ 0.0,  1.0,  2.0,  3.0, 3.14]), array('d', [10.0, 20.0, 0.0, 14.5, 17.9]))
          >>> p.find_index_ge(-1)
          0
          >>> p.find_index_ge(0.0)
          0
          >>> p.find_index_ge(0.5)
          1
          >>> p.find_index_ge(0.75)
          1
          >>> p.find_index_ge(1.0)
          1
          >>> p.find_index_ge(2.0)
          2
          >>> p.find_index_ge(2.5)
          3
          >>> p.find_index_ge(3.08)
          4
          >>> p.find_index_ge(3.14)
          4
          >>> p.find_index_ge(5) is None
          True
        """

        check_type(s_val, "S value", numbers.Real)
        if not np.isfinite(s_val):
            raise Exception(f"S value must be finite but {s_val} got")

        if s_val <= self.s_min():
            return 0
        elif s_val > self.s_max():
            return None
        else:
            return int( np.argmax(self._s >= s_val) )

    def s_upper_ndx(self,
                    s_val: numbers.Real
                    ) -> Optional[numbers.Integral]:
        """
        To be possibly deprecated.
        Returns the optional index in the s array of the provided value.

        :param s_val: the value to search the index for in the x array
        :return: the optional index in the s array of the provided value

        Examples:
          >>> p = ZTrace(array('d', [ 0.0,  1.0,  2.0,  3.0, 3.14]), array('d', [10.0, 20.0, 0.0, 14.5, 17.9]))
          >>> p.s_upper_ndx(-1) is None
          True
          >>> p.s_upper_ndx(5) is None
          True
          >>> p.s_upper_ndx(0.5)
          1
          >>> p.s_upper_ndx(0.75)
          1
          >>> p.s_upper_ndx(1.0)
          1
          >>> p.s_upper_ndx(2.0)
          2
          >>> p.s_upper_ndx(2.5)
          3
          >>> p.s_upper_ndx(3.11)
          4
          >>> p.s_upper_ndx(3.14)
          4
          >>> p.s_upper_ndx(0.0)
          0
        """

        check_type(s_val, "Input value", numbers.Real)

        if s_val < self.s_min() or s_val > self.s_max():
            return None

        return int( np.searchsorted(self._s, s_val) )

    def z_linear_interpol(self,
                          s_val: numbers.Real
                          ) -> Optional[numbers.Real]:
        """
        Returns the optional interpolated z value in the z array of the provided s value.

        :param s_val: the value to search the index for in the s array
        :return: the optional interpolated z value

        Examples:
          >>> p = ZTrace(array('d', [ 0.0,  1.0,  2.0,  3.0, 3.14]), array('d', [10.0, 20.0, 0.0, 14.5, 17.9]))
          >>> p.z_linear_interpol(-1) is None
          True
          >>> p.z_linear_interpol(5) is None
          True
          >>> p.z_linear_interpol(0.5)
          15.0
          >>> p.z_linear_interpol(0.75)
          17.5
          >>> p.z_linear_interpol(2.5)
          7.25
          >>> p.z_linear_interpol(3.14)
          17.9
          >>> p.z_linear_interpol(0.0)
          10.0
          >>> p.z_linear_interpol(1.0)
          20.0
        """

        check_type(s_val, "Input value", numbers.Real)

        ndx = self.s_upper_ndx(s_val)

        if ndx is not None:

            if ndx == 0:
                return float( self.z(0) )

            val_z_i = self.z(ndx - 1)
            val_z_i_next = self.z(ndx)
            delta_val_z = val_z_i_next - val_z_i

            if delta_val_z == 0.0:
                return float( val_z_i )

            val_s_i = self.s(ndx - 1)
            val_s_i_next = self.s(ndx)
            delta_val_s = val_s_i_next - val_s_i

            if delta_val_s == 0.0:
                return float( val_z_i )

            d_s = s_val - val_s_i

            return float( val_z_i + d_s * delta_val_z / delta_val_s )

        else:

            return None

    def s_subset(self,
                 s_start: numbers.Real,
                 s_end: Optional[numbers.Real] = None
                 ) -> Optional[np.ndarray]:
        """
        Return the s array values defined by the provided s range (extremes included).
        When the end value is not provided, a single-valued array is returned.

        :param s_start: the start s value (distance along the profile)
        :param s_end: the optional end s value (distance along the profile)
        :return: the s array subset, possibly with just a value

        Examples:
          >>> p = ZTrace(array('d', [ 0.0,  1.0,  2.0,  3.0, 3.14]), array('d', [10.0, 20.0, 0.0, 14.5, 17.9]))
          >>> p.s_subset(1.0)
          array([1.])
          >>> p.s_subset(0.0)
          array([0.])
          >>> p.s_subset(0.75)
          array([0.75])
          >>> p.s_subset(3.14)
          array([3.14])
          >>> p.s_subset(1.0, 2.0)
          array([1., 2.])
          >>> p.s_subset(0.75, 2.0)
          array([0.75, 1.  , 2.  ])
          >>> p.s_subset(0.75, 2.5)
          array([0.75, 1.  , 2.  , 2.5 ])
          >>> p.s_subset(0.75, 3.0)
          array([0.75, 1.  , 2.  , 3.  ])
          >>> p.s_subset(-1, 1)
          array([0., 1.])
          >>> p.s_subset(-1) is None
          True
          >>> p.s_subset(0.0, 10)
          array([0.  ,  1.  ,  2.  ,  3.  , 3.14])
          >>> p.s_subset(0.0, 3.14)
          array([0.  ,  1.  ,  2.  ,  3.  , 3.14])
        """

        # input data validity checks

        check_type(s_start, "Start s value", numbers.Real)
        if s_end is not None:
            check_type(s_end, "End s value", numbers.Real)

        if s_end is not None and s_end < s_start:
            raise Exception(f"Start is {s_start} while end is {s_end}")

        # result for single input value

        if s_end is None:
            if not self.s_min() <= s_start <= self.s_max():
                return None
            else:
                return np.array([s_start], dtype=np.float64)

        if s_start > self.s_max() or s_end < self.s_min():
            return None

        # fix ranges

        if s_start < self.s_min():
            s_start = self.s_min()
        if s_end > self.s_max():
            s_end = self.s_max()

        # particular case where start equal to end

        if s_end == s_start:
            return np.array([s_start], dtype=np.float64)

        # general case for x start < x end

        values = []

        s_start_upper_index_value = self.s_upper_ndx(s_start)

        if s_start < self.s(s_start_upper_index_value):
            values.append(s_start)

        s_end_upper_index_value = self.s_upper_ndx(s_end)

        for ndx in range(s_start_upper_index_value, s_end_upper_index_value):
            values.append(self.s(ndx))

        if s_end > self.s(s_end_upper_index_value - 1):
            values.append(s_end)

        return np.array(values, dtype=np.float64)

    def zs_from_s_range(self,
                        s_start: numbers.Real,
                        s_end: Optional[numbers.Real] = None
                        ) -> Optional[np.ndarray]:
        """
        Return the z array values defined by the provided s range (extremes included).
        When the end value is not provided, a single-valued array is returned.

        :param s_start: the start s value (distance along the profile)
        :param s_end: the optional end s value (distance along the profile)
        :return: the z array, possibly with just a value

        Examples:
          >>> p = ZTrace(array('d', [ 0.0,  1.0,  2.0,  3.0, 3.14]), array('d', [10.0, 20.0, 0.0, 14.5, 17.9]))
          >>> p.zs_from_s_range(1.0)
          array([20.])
          >>> p.zs_from_s_range(0.0)
          array([10.])
          >>> p.zs_from_s_range(0.75)
          array([17.5])
          >>> p.zs_from_s_range(3.14)
          array([17.9])
          >>> p.zs_from_s_range(1.0, 2.0)
          array([20., 0.])
          >>> p.zs_from_s_range(0.75, 2.0)
          array([17.5, 20. , 0. ])
          >>> p.zs_from_s_range(0.75, 2.5)
          array([17.5 , 20.  , 0.  , 7.25])
          >>> p.zs_from_s_range(0.75, 3.0)
          array([17.5, 20. , 0. , 14.5])
          >>> p.zs_from_s_range(-1, 1)
          array([10., 20.])
          >>> p.zs_from_s_range(-1) is None
          True
          >>> p.zs_from_s_range(0.0, 10)
          array([10. , 20. , 0. , 14.5, 17.9])
          >>> p.zs_from_s_range(0.0, 3.14)
          array([10. , 20. , 0. , 14.5, 17.9])
        """

        s_subset = self.s_subset(s_start, s_end)

        if s_subset is None:
            return None

        return np.array(list(map(self.z_linear_interpol, s_subset)), dtype=np.float64)

    def extend_in_place(self,
               another: 'ZTrace'
               ):
        """
        Extend an ZTrace instance in-place.
        Note that the last element of the first s-z array pair will be dropped
        and substituted by the first element of the second s-z array pair
        (no attempt to check values equality or to calculate a mean).
        Moverover, all the s values of the second s-z array pair will be incremented
        by the last s value of the first element.

        Examples:
          >>> f = ZTrace(np.array([ 0.0,  14.2,  20.0]), np.array([149.17, 132.4, 159.2]))
          >>> f.s_breaks()
          array([ 0., 20.])
          >>> f.s_max()
          20.0
          >>> s = ZTrace(np.array([ 0.0,  7.0,  11.0]), np.array([159.17, 180.1, 199.5]))
          >>> s.s_breaks()
          array([ 0., 11.])
          >>> s.s_breaks() + f.s_max()
          array([20., 31.])
          >>> f.extend_in_place(s)
          >>> f.s_arr()
          array([ 0. ,  14.2,  20. , 27. , 31. ])
          >>> f.z_arr()
          array([149.17, 132.4 , 159.17, 180.1 , 199.5 ])
          >>> f.s_breaks()
          array([ 0., 20., 31.])
        """

        check_type(another, "Second ZTrace instance", ZTrace)

        offset = self.s_max()
        self._s = np.append(self._s[:-1], another._s + offset)
        self._z = np.append(self._z[:-1], another._z)
        self._s_breaks = np.append(self._s_breaks[:-1], another._s_breaks + offset)

    def dir_slopes_deltas(self) -> Tuple[np.ndarray, np.ndarray]:

        delta_s = np.ediff1d(self._s, to_end=np.nan)
        delta_z = np.ediff1d(self._z, to_end=np.nan)

        return delta_s, delta_z

    def dir_slopes_ratios(self) -> np.ndarray:

        ds, dz = self.dir_slopes_deltas()

        return dz / ds

    def dir_slopes_percent(self) -> np.ndarray:

        return 100.0 * self.dir_slopes_ratios()

    def dir_slopes_radians(self) -> np.ndarray:

        ds, dz = self.dir_slopes_deltas()

        return np.arctan2(dz, ds)

    def dir_slopes_degrees(self) -> np.ndarray:

        return self.dir_slopes_radians() * 180.0 / np.pi

    def as_dir_slopes_degrees(self) -> 'ZTrace':

        return ZTrace(
            s_array=self._s,
            z_array=self.dir_slopes_degrees(),
            s_breaks=self._s_breaks
        )

    def as_abs_slopes_degrees(self) -> 'ZTrace':

        return ZTrace(
            s_array=self._s,
            z_array=np.fabs(self.dir_slopes_degrees()),
            s_breaks=self._s_breaks
        )


class ZTraces:
    """
    Class storing a set of z profiles.
    Important: it assumes and requires that all ZTraces have the same data length.
    """

    def __init__(self,
                 s_array: np.ndarray,
                 z_array: np.ndarray,
                 s_breaks: np.ndarray):
        """
        Instantiates a topographic profile set.

        :param s_array: the 1D-array of s values.
        :param z_array: the 2D-array of z values. Can store multiple z-series.
        :param s_breaks: the 1D-array of s break values (i.e., where the profile changes direction).
        """

        check_type(s_array, "S array", np.ndarray)
        check_type(z_array, "Z array", np.ndarray)
        check_type(s_breaks, "S breaks array", np.ndarray)

        if s_array.ndim != 1:
            raise Exception(f"S array must be 1D but {s_array.ndim} got")
        if z_array.ndim != 2:
            raise Exception(f"Z array must be 2D but {z_array.ndim} got")
        if s_breaks.ndim != 1:
            raise Exception(f"S breaks array must be 1D but {s_breaks.ndim} got")

        if not np.all(s_array[1:] > s_array[:-1]):
            raise Exception(f"S array must be strictly increasing")

        num_steps = len(s_array)
        if num_steps <= 1:
            raise Exception(f"At least a two-values S array is required but {num_steps} got")

        (num_profiles, num_profs_steps) = z_array.shape
        if num_profs_steps != num_steps:
            raise Exception(f"S array has {num_steps} steps while Z array has {num_profs_steps} steps")

        self._s_array = s_array.astype(np.float64)
        self._z_array = z_array.astype(np.float64)
        self._s_breaks_array = s_breaks.astype(np.float64)

    @classmethod
    def fromProfiles(cls,
                     z_profiles: List[ZTrace]):
        """
        Instantiates a topographic profile set.

        :param z_profiles: the topographic profile set.
        """

        check_type(z_profiles, "Topographic profiles set", List)
        for el in z_profiles:
            check_type(el, "Topographic profile", ZTrace)

        # ZTraces requires all profiles to share the same number of sampled
        # steps (see class docstring): fail here with a clear message instead
        # of letting np.array(...) below raise a cryptic "inhomogeneous shape"
        # ValueError when profiles have different lengths (e.g. parallel
        # profiles of differing length, or sampled with different distances)
        num_steps = [len(zprof.s_arr()) for zprof in z_profiles]
        if len(set(num_steps)) > 1:
            raise Exception(
                f"All profiles must have the same number of sampled steps to be "
                f"combined into a ZTraces instance, but got differing lengths: {num_steps}"
            )

        return cls(
            z_profiles[0].s_arr(),
            np.array([zprof.z_arr() for zprof in z_profiles]),
            z_profiles[0].s_breaks()
        )

    def __getitem__(self,
        ndx: numbers.Integral) -> ZTrace:
        """
        Returns the profile with the given index.
        """

        return self.make_ztrace_from_arr(self._z_array[ndx])

    def num_profiles(self) -> int:
        """
        Returns the number of profiles.
        """

        return self._z_array.shape[0]

    def s_breaks(self) -> np.ndarray:
        """
        Returns a copy of the S breaks array.
        """

        return np.copy(self._s_breaks_array)

    def s_min(self) -> float:
        """
        Returns the minimum s value in the topographic profiles.

        :return: the minimum s value in the profiles.
        """

        return float(self._s_array[0])

    def s_max(self) -> float:
        """
        Returns the maximum s value in the topographic profiles.

        :return: the maximum s value in the profiles.
        :rtype: optional numbers.Real.
        """

        return float(self._s_array[-1])

    def z_min(self) -> Optional[numbers.Real]:
        """
        Returns the minimum elevation value in the topographic profiles.

        :return: the minimum elevation value in the profiles.
        :rtype: optional numbers.Real.
        """

        return np.nanmin(self._z_array)

    def z_max(self) -> Optional[numbers.Real]:
        """
        Returns the maximum elevation value in the topographic profiles.

        :return: the maximum elevation value in the profiles.
        :rtype: optional numbers.Real.
        """

        return np.nanmax(self._z_array)

    def natural_elev_range(self) -> Tuple[numbers.Real, numbers.Real]:
        """
        Returns the elevation range of the profiles.

        :return: minimum and maximum values of the considered topographic profiles.
        :rtype: tuple of two floats.
        """

        return self.z_min(), self.z_max()

    def topoprofiles_params(self):
        """

        :return:
        """

        return self.s_min(), self.s_max(), self.z_min(), self.z_max()

    def max_arr(self) -> np.ndarray:
        """
        Return the array of the along-profiles max values.
        """

        return np.nanmax(self._z_array, axis=0)

    def min_arr(self) -> np.ndarray:
        """
        Return the array of the along-profiles min values.
        """

        return np.nanmin(self._z_array, axis=0)

    def mean_arr(self) -> np.ndarray:
        """
        Return the array of the along-profiles mean values.
        """

        return np.nanmean(self._z_array, axis=0)

    def median_arr(self) -> np.ndarray:
        """
        Return the array of the along-profiles median values.
        """

        return np.nanmedian(self._z_array, axis=0)

    def var_arr(self) -> np.ndarray:
        """
        Return the array of the along-profiles var values.
        """

        return np.nanvar(self._z_array, axis=0)

    def std_arr(self) -> np.ndarray:
        """
        Return the array of the along-profiles std values.
        """

        return np.nanstd(self._z_array, axis=0)

    def range_arr(self) -> np.ndarray:
        """
        Return the array of the along-profiles range values.
        """

        return np.ptp(self._z_array, axis=0)

    def middle_arr(self) -> np.ndarray:
        """
        Return the array of the along-profiles range values.
        """

        return self.min_arr() + 0.5 * self.range_arr()

    def make_ztrace_from_arr(self,
                             arr: np.ndarray) -> ZTrace:
        """
        Creates an ZTrace instance given a 1D z array.

        :param arr: the array representing the z values.
        :return: a new ZTrace
        """

        return ZTrace(
            self._s_array,
            arr,
            self._s_breaks_array
        )

    def profile_max(self) -> ZTrace:
        """
        Creates a profile of the max values.
        """

        return self.make_ztrace_from_arr(self.max_arr())

    def profile_min(self) -> ZTrace:
        """
        Creates a profile of the min values.
        """

        return self.make_ztrace_from_arr(self.min_arr())

    def profile_mean(self) -> ZTrace:
        """
        Creates a profile of the mean values.
        """

        return self.make_ztrace_from_arr(self.mean_arr())

    def profile_median(self) -> ZTrace:
        """
        Creates a profile of the median values.
        """

        return self.make_ztrace_from_arr(self.median_arr())

    def profile_std(self) -> ZTrace:
        """
        Creates a profile of the std values.
        """

        return self.make_ztrace_from_arr(self.std_arr())

    def profile_var(self) -> ZTrace:
        """
        Creates a profile of the var values.
        """

        return self.make_ztrace_from_arr(self.var_arr())

    def profile_range(self) -> ZTrace:
        """
        Creates a profile of the range values.
        """

        return self.make_ztrace_from_arr(self.range_arr())

    def profile_middle(self) -> ZTrace:
        """
        Creates a profile of the middle values.
        """

        return self.make_ztrace_from_arr(self.middle_arr())


def join_ztraces(
        *ztraces: ZTrace
) -> Tuple[Union[None, ZTrace], Error]:
    """
    Combines a tuple of ZTrace instances into an extended single ZTrace instance.

    Examples:
      >>> f = ZTrace(np.array([ 0.0,  14.2,  20.0]), np.array([149.17, 132.4, 159.2]))
      >>> s = ZTrace(np.array([ 0.0,  7.0,  11.0]), np.array([159.17, 180.1, 199.5]))
      >>> t = ZTrace(np.array([ 0.0,  22.0,  30.0]), np.array([199.5, 200.1, 179.1]))
      >>> combined, err = join_ztraces(f, s, t)
      >>> combined.s_arr()
      array([ 0. ,  14.2, 20. , 27. , 31. , 53. , 61. ])
      >>> combined.z_arr()
      array([149.17, 132.4 , 159.17, 180.1 , 199.5 , 200.1 , 179.1 ])
      >>> combined.s_breaks()
      array([ 0., 20., 31., 61.])
    """

    try:

        if len(ztraces) == 0:
            return None, Error(
                True,
                caller_name(),
                Exception("Length of ztraces is zero"),
                traceback.format_exc(),
            )

        combined_xy_arrays = ztraces[0].clone()

        for xy_array_pair in ztraces[1:]:
            combined_xy_arrays.extend_in_place(xy_array_pair)

        return combined_xy_arrays, Error()

    except Exception as e:

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


class ZArray:
    """
    Class storing a set of z profiles as arrays.
    Possible candidate substitute for ZTraces.
    """

    def __init__(self):

        self._s_array = None
        self._z_array = None
        self._s_breaks_array = None
        self._names = []

    def add_profile(self,
        z_trace: ZTrace,
        name: str = "undefined"):

        if self._s_array is None:
            self._s_array = z_trace.s_arr()

        if self._z_array is None:
            self._z_array = np.expand_dims(z_trace.z_arr(), axis=0)
        else:
            self._z_array = np.vstack((self._z_array, z_trace.z_arr()))

        if self._s_breaks_array is None:
            self._s_breaks_array = z_trace.s_breaks()

        self._names.append(name)

    def __getitem__(self,
        ndx: numbers.Integral) -> ZTrace:
        """
        Returns the profile with the given index.
        """

        return ZTrace(
            self._s_array,
            self._z_array[ndx, :],
            self._s_breaks_array,
            self._names[ndx]
        )

    def make_ztrace_from_arr(self,
        arr: np.ndarray,
        name: str = "undefined"
    ) -> ZTrace:
        """
        Creates an ZTrace instance given a 1D z array.

        :param arr: the array representing the z values.
        :return: a new ZTrace
        """

        return ZTrace(
            self._s_array,
            arr,
            self._s_breaks_array,
            name
        )

    def num_profiles(self) -> int:
        """
        Returns the number of profiles.
        """

        return self._z_array.shape[0]

    def s_breaks(self) -> np.ndarray:
        """
        Returns a copy of the S breaks array.
        """

        return np.copy(self._s_breaks_array)

    def s_min(self) -> float:
        """
        Returns the minimum s value in the topographic profiles.

        :return: the minimum s value in the profiles.
        """

        return float(self._s_array[0])

    def s_max(self) -> float:
        """
        Returns the maximum s value in the topographic profiles.

        :return: the maximum s value in the profiles.
        :rtype: optional numbers.Real.
        """

        return float(self._s_array[-1])

    def z_min(self) -> Optional[numbers.Real]:
        """
        Returns the minimum elevation value in the topographic profiles.

        :return: the minimum elevation value in the profiles.
        :rtype: optional numbers.Real.
        """

        return np.nanmin(self._z_array)

    def z_max(self) -> Optional[numbers.Real]:
        """
        Returns the maximum elevation value in the topographic profiles.

        :return: the maximum elevation value in the profiles.
        :rtype: optional numbers.Real.
        """

        return np.nanmax(self._z_array)

    def natural_elev_range(self) -> Tuple[numbers.Real, numbers.Real]:
        """
        Returns the elevation range of the profiles.

        :return: minimum and maximum values of the considered topographic profiles.
        :rtype: tuple of two floats.
        """

        return self.z_min(), self.z_max()

    def topoprofiles_params(self):
        """

        :return:
        """

        return self.s_min(), self.s_max(), self.z_min(), self.z_max()

    def max_arr(self) -> np.ndarray:
        """
        Return the array of the along-profiles max values.
        """

        return np.nanmax(self._z_array, axis=0)

    def min_arr(self) -> np.ndarray:
        """
        Return the array of the along-profiles min values.
        """

        return np.nanmin(self._z_array, axis=0)

    def mean_arr(self) -> np.ndarray:
        """
        Return the array of the along-profiles mean values.
        """

        return np.nanmean(self._z_array, axis=0)

    def median_arr(self) -> np.ndarray:
        """
        Return the array of the along-profiles median values.
        """

        return np.nanmedian(self._z_array, axis=0)

    def var_arr(self) -> np.ndarray:
        """
        Return the array of the along-profiles var values.
        """

        return np.nanvar(self._z_array, axis=0)

    def std_arr(self) -> np.ndarray:
        """
        Return the array of the along-profiles std values.
        """

        return np.nanstd(self._z_array, axis=0)

    def range_arr(self) -> np.ndarray:
        """
        Return the array of the along-profiles range values.
        """

        return np.ptp(self._z_array, axis=0)

    def middle_arr(self) -> np.ndarray:
        """
        Return the array of the along-profiles range values.
        """

        return self.min_arr() + 0.5 * self.range_arr()

    def profile_max(self) -> ZTrace:
        """
        Creates a profile of the max values.
        """

        return self.make_ztrace_from_arr(
            self.max_arr(),
            "maximum")

    def profile_min(self) -> ZTrace:
        """
        Creates a profile of the min values.
        """

        return self.make_ztrace_from_arr(
            self.min_arr(),
            "minimum")

    def profile_mean(self) -> ZTrace:
        """
        Creates a profile of the mean values.
        """

        return self.make_ztrace_from_arr(
            self.mean_arr(),
            "mean")

    def profile_median(self) -> ZTrace:
        """
        Creates a profile of the median values.
        """

        return self.make_ztrace_from_arr(
            self.median_arr(),
            "median")

    def profile_std(self) -> ZTrace:
        """
        Creates a profile of the std values.
        """

        return self.make_ztrace_from_arr(
            self.std_arr(),
            "standard deviation")

    def profile_var(self) -> ZTrace:
        """
        Creates a profile of the var values.
        """

        return self.make_ztrace_from_arr(
            self.var_arr(),
            "variance")

    def profile_range(self) -> ZTrace:
        """
        Creates a profile of the range values.
        """

        return self.make_ztrace_from_arr(
            self.range_arr(),
            "range")

    def profile_middle(self) -> ZTrace:
        """
        Creates a profile of the middle values.
        """

        return self.make_ztrace_from_arr(
            self.middle_arr(),
            "middle")


