import math
from typing import List

import numpy as np
from qgis.core import (
    QgsGeometry,
    QgsLineString,
    QgsProcessing,
    QgsFeature,
    QgsFeatureSink,
    QgsField,
    QgsProcessingException,
    QgsProcessingParameterFeatureSink,
    QgsProcessingParameterFeatureSource,
    QgsProcessingParameterField,
    QgsProcessingParameterNumber,
    QgsWkbTypes,
)
from qgis.PyQt.QtCore import QVariant

from PreCourlis.core.precourlis_file import PreCourlisFileLine
from PreCourlis.processing.precourlis_algorithm import PreCourlisAlgorithm


class InterpolateProfilesWithoutLateralResamplingAlgorithm(PreCourlisAlgorithm):
    SECTIONS = "SECTIONS"
    AXIS = "AXIS"
    LONG_STEP = "LONG_STEP"
    ATTR_CROSS_SECTION = "ATTR_CROSS_SECTION"
    AXIS_POINTS = "AXIS_POINTS"
    OUTPUT = "OUTPUT"

    def initAlgorithm(self, config=None):
        self.addParameter(
            QgsProcessingParameterFeatureSource(
                self.SECTIONS,
                self.tr("Sections"),
                types=[QgsProcessing.TypeVectorLine],
                defaultValue=None,
            )
        )
        self.addParameter(
            QgsProcessingParameterFeatureSource(
                self.AXIS,
                self.tr("Axis"),
                types=[QgsProcessing.TypeVectorLine],
                defaultValue=None,
            )
        )
        self.addParameter(
            QgsProcessingParameterNumber(
                self.LONG_STEP,
                self.tr("Longitudinal space step (in m)"),
                defaultValue=1.0,
                type=QgsProcessingParameterNumber.Double,
            )
        )
        self.addParameter(
            QgsProcessingParameterField(
                self.ATTR_CROSS_SECTION,
                self.tr("Attribute to identify cross-sections"),
                parentLayerParameterName=self.SECTIONS,
                defaultValue="sec_id",
            )
        )
        self.addParameter(
            QgsProcessingParameterFeatureSink(
                self.OUTPUT, self.tr("Interpolated profiles")
            )
        )

    def processAlgorithm(self, parameters, context, feedback):
        sections_source = self.parameterAsSource(parameters, self.SECTIONS, context)
        axis_source = self.parameterAsSource(parameters, self.AXIS, context)
        # axis_point_output = self.parameterAsOutputLayer(parameters, self.AXIS_POINTS, context)
        long_step = self.parameterAsDouble(parameters, self.LONG_STEP, context)

        file = PreCourlisFileLine(sections_source)
        layers = file.layers()

        axis: QgsFeature = next(axis_source.getFeatures())

        fields = PreCourlisFileLine.base_fields()
        for layer in file.layers():
            fields.append(QgsField(layer, QVariant.String, len=100000))
        sink, dest_id = self.parameterAsSink(
            parameters,
            self.OUTPUT,
            context,
            fields,
            QgsWkbTypes.LineString,
            sections_source.sourceCrs(),
        )
        if sink is None:
            raise QgsProcessingException(self.invalidSinkError(parameters, self.OUTPUT))

        # for each area between to sections
        total = (
            100.0 / sections_source.featureCount()
            if sections_source.featureCount()
            else 0
        )
        sec_id = 0
        sections = list(sections_source.getFeatures())
        for i_area in range(len(sections) - 1):
            first_section: QgsFeature = sections[i_area]
            last_section: QgsFeature = sections[i_area + 1]

            first_section_geom = first_section.geometry()
            first_section_geom.convertToSingleType()
            first_section_intersection = first_section_geom.intersection(
                axis.geometry()
            )
            first_section_abs_long = axis.geometry().lineLocatePoint(
                first_section_intersection
            )
            # abs_lat at intersection with axis
            first_section_abs_lat_offset = first_section.geometry().lineLocatePoint(
                first_section_intersection
            )

            last_section_geom = last_section.geometry()
            last_section_geom.convertToSingleType()
            last_section_intersection = last_section_geom.intersection(axis.geometry())
            last_section_abs_long = axis.geometry().lineLocatePoint(
                last_section_intersection
            )
            # abs_lat at intersection with axis
            last_section_abs_lat_offset = last_section.geometry().lineLocatePoint(
                last_section_intersection
            )

            if (
                first_section_geom.constGet().numPoints()
                != last_section_geom.constGet().numPoints()
            ):
                raise QgsProcessingException(
                    self.tr(
                        f"Profiles {first_section['sec_name']} and {last_section['sec_name']}"
                        " do not have the same number of points."
                    )
                )

            base_points: List[QgsGeometry] = []
            current_abs_long = first_section_abs_long + long_step
            while current_abs_long < last_section_abs_long:
                point = axis.geometry().interpolate(current_abs_long)
                base_points.append(point)
                current_abs_long += long_step

            # First area section
            sec_id += 1
            f = QgsFeature()
            f.setGeometry(first_section.geometry())
            f.setFields(fields)
            f.setAttribute("sec_id", sec_id)
            f.setAttribute("sec_name", first_section["sec_name"])
            f.setAttribute("abs_long", first_section["abs_long"])
            f.setAttribute("axis_x", first_section["axis_x"])
            f.setAttribute("axis_y", first_section["axis_y"])
            f.setAttribute("layers", first_section["layers"])
            f.setAttribute("p_id", first_section["p_id"])
            f.setAttribute("topo_bat", first_section["topo_bat"])
            f.setAttribute("abs_lat", first_section["abs_lat"])
            f.setAttribute("zfond", first_section["zfond"])
            for layer in layers:
                f.setAttribute(layer, first_section[layer])
            sink.addFeature(f, QgsFeatureSink.FastInsert)

            for intersection_point in base_points:
                sec_id += 1
                points = []

                current_section_abs_long = axis.geometry().lineLocatePoint(
                    intersection_point
                )

                def interpolate(first_value: float, last_value: float) -> float:
                    return np.interp(
                        [current_section_abs_long],
                        [first_section_abs_long, last_section_abs_long],
                        [first_value, last_value],
                    )[0]

                axis_angle = axis.geometry().interpolateAngle(current_section_abs_long)

                # Construct current section geometry
                for i_lat in range(first_section_geom.constGet().numPoints()):
                    first_point = QgsGeometry(
                        first_section_geom.constGet().pointN(i_lat).clone()
                    )
                    last_point = QgsGeometry(
                        last_section_geom.constGet().pointN(i_lat).clone()
                    )

                    first_relative_abs_lat = (
                        first_section.geometry().lineLocatePoint(first_point)
                        - first_section_abs_lat_offset
                    )
                    last_relative_abs_lat = (
                        last_section.geometry().lineLocatePoint(last_point)
                        - last_section_abs_lat_offset
                    )

                    relative_abs_lat = interpolate(
                        first_relative_abs_lat, last_relative_abs_lat
                    )

                    # TODO: proportional ?
                    angle = 90

                    new_point = intersection_point.constGet().project(
                        relative_abs_lat, angle + 180.0 / math.pi * axis_angle
                    )
                    points.append(new_point)

                line = QgsLineString([p.clone() for p in points])
                geom = QgsGeometry(line.clone())

                # Compute attributes
                abs_long = current_abs_long - first_section_abs_long
                abs_lats = [
                    geom.lineLocatePoint(QgsGeometry(p.clone())) for p in points
                ]
                first_section_zfond_list = [
                    float(s) for s in first_section.attribute("zfond").split(",")
                ]
                last_section_zfond_list = [
                    float(s) for s in last_section.attribute("zfond").split(",")
                ]
                zfond = [
                    interpolate(first_z, last_z)
                    for first_z, last_z in zip(
                        first_section_zfond_list,
                        last_section_zfond_list,
                    )
                ]
                layers_values = {}
                for layer in layers:
                    first_section_layer_values = [
                        float(s) for s in first_section.attribute(layer).split(",")
                    ]
                    last_section_layer_values = [
                        float(s) for s in last_section.attribute(layer).split(",")
                    ]
                    layers_values[layer] = [
                        interpolate(first_v, last_v)
                        for first_v, last_v in zip(
                            first_section_layer_values,
                            last_section_layer_values,
                        )
                    ]

                f = QgsFeature()
                f.setGeometry(geom)
                f.setFields(fields)
                f.setAttribute("sec_id", sec_id)
                f.setAttribute(
                    "sec_name",
                    "{}_{:04.3f}".format("P" if abs_long >= 0 else "N", abs_long),
                )
                f.setAttribute("abs_long", abs_long)
                f.setAttribute("axis_x", intersection_point.constGet().x())
                f.setAttribute("axis_y", intersection_point.constGet().y())
                f.setAttribute("layers", first_section.attribute("layers"))
                f.setAttribute("p_id", first_section.attribute("p_id"))
                f.setAttribute("topo_bat", first_section.attribute("topo_bat"))
                f.setAttribute("abs_lat", ",".join([str(v) for v in abs_lats]))
                f.setAttribute("zfond", ",".join([str(v) for v in zfond]))
                for layer in layers:
                    f.setAttribute(
                        layer, ",".join([str(v) for v in layers_values[layer]])
                    )
                sink.addFeature(f, QgsFeatureSink.FastInsert)

            feedback.setProgress(int(i_area * total))

        # Lastest section
        sec_id += 1
        f = QgsFeature()
        f.setGeometry(last_section.geometry())
        f.setFields(fields)
        f.setAttribute("sec_id", sec_id)
        f.setAttribute("sec_name", last_section["sec_name"])
        f.setAttribute("abs_long", last_section["abs_long"])
        f.setAttribute("axis_x", last_section["axis_x"])
        f.setAttribute("axis_y", last_section["axis_y"])
        f.setAttribute("layers", last_section["layers"])
        f.setAttribute("p_id", last_section["p_id"])
        f.setAttribute("topo_bat", last_section["topo_bat"])
        f.setAttribute("abs_lat", last_section["abs_lat"])
        f.setAttribute("zfond", last_section["zfond"])
        for layer in layers:
            f.setAttribute(layer, last_section[layer])
        sink.addFeature(f, QgsFeatureSink.FastInsert)

        return {"OUTPUT": dest_id}

    def name(self):
        return "interpolate_without_lateral_resampling"

    def displayName(self):
        return self.tr("Interpolate profiles without lateral resampling")

    def group(self):
        return self.tr("Interpolate")

    def groupId(self):
        return "Interpolate"

    def createInstance(self):
        return InterpolateProfilesWithoutLateralResamplingAlgorithm()

    def shortHelpString(self):
        return self.tr(
            "This algorithm interpolate a profiles layer by creating"
            " regular and perpendicular profiles along axis"
            " with fixed number of points per profile."
            "\n\nNote that the source profiles must have the same number of points"
            " and be perpendical to the axis."
        )
