
from math import atan

from gst.qgis.projections import *


def calculate_local_geographic_north_azimuth_from_crs(
    current_qgis_crs: QgsCoordinateReferenceSystem,
    current_crs_local_pt: Point,
) -> Tuple[Union[None, numbers.Real], Error]:
    """
    Returns the azimuth of the North direction, clockwise in degrees,
    calculated from the top (Y) map direction.
    """

    try:

        # Calculates dip direction correction with respect to project CRS y-axis orientation

        current_crs_local_pt_x = current_crs_local_pt.x
        current_crs_local_pt_y = current_crs_local_pt.y

        epsg_4326_local_pt_lon, epsg_4326_local_pt_lat = qgs_project_xy(
            x=current_crs_local_pt_x,
            y=current_crs_local_pt_y,
            src_crs=current_qgis_crs
        )

        epsg_4326_north_dummy_pt_lon = epsg_4326_local_pt_lon  # no change
        epsg_4326_north_dummy_pt_lat = epsg_4326_local_pt_lat + (1.0 / 1200.0)  # add 3 minute-seconds (approximately 90 meters)

        current_crs_north_dummy_pt_x, current_crs_north_dummy_pt_y = qgs_project_xy(
            x=epsg_4326_north_dummy_pt_lon,
            y=epsg_4326_north_dummy_pt_lat,
            dest_crs=current_qgis_crs
        )

        current_crs_north_dummy_pt = Point(
            current_crs_north_dummy_pt_x,
            current_crs_north_dummy_pt_y
        )

        current_crs_north_directed_segment = Segment(
            start_pt=current_crs_local_pt.as_point2d(),
            end_pt=current_crs_north_dummy_pt)

        current_crs_north_directed_vector = current_crs_north_directed_segment.as_vector2d()

        geographic_north_azimuth_cw = current_crs_north_directed_vector.azimuth_degr()

        return geographic_north_azimuth_cw, Error()

    except Exception as e:

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


def calculate_local_geographic_north_azimuth_from_epsg_4326(
    longitude: numbers.Real,
    latitude: numbers.Real,
    planar_qgis_crs: QgsCoordinateReferenceSystem,
) -> Tuple[Union[None, numbers.Real], Error]:
    """
    Returns the azimuth of the North direction, clockwise in degrees,
    calculated from the top (Y) map direction.
    """

    try:

        # Calculates dip direction correction with respect to project CRS y-axis orientation

        north_dummy_pt_longitude = longitude  # no change
        north_dummy_pt_latitude = latitude + (1.0 / 1200.0)  # add 3 minute-seconds (approximately 90 meters)

        current_crs_north_dummy_pt_x, current_crs_north_dummy_pt_y = qgs_project_xy(
            x=north_dummy_pt_longitude,
            y=north_dummy_pt_latitude,
            dest_crs=planar_qgis_crs
        )

        current_crs_north_dummy_pt = Point(
            current_crs_north_dummy_pt_x,
            current_crs_north_dummy_pt_y
        )

        current_crs_local_pt_x, current_crs_local_pt_y = qgs_project_xy(
            x=longitude,
            y=latitude,
            dest_crs=planar_qgis_crs
        )

        current_crs_local_pt = Point(
            current_crs_local_pt_x,
            current_crs_local_pt_y,
        )

        current_crs_north_directed_segment = Segment(
            start_pt=current_crs_local_pt.as_point2d(),
            end_pt=current_crs_north_dummy_pt)

        current_crs_north_directed_vector = current_crs_north_directed_segment.as_vector2d()

        geographic_north_azimuth_cw = current_crs_north_directed_vector.azimuth_degr()

        return geographic_north_azimuth_cw, Error()

    except Exception as e:

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


def convert_plane_orientation_to_crs(
    dip_direction_as_north_azimuth: numbers.Real,
    source_dip_angle: numbers.Real,
    local_pt_longitude: numbers.Real,
    local_pt_latitude: numbers.Real,
    dest_crs: QgsCoordinateReferenceSystem,
    #dummy_distance: numbers.Real = 100  # meters
) -> Tuple[Union[None, Plane], Error]:

    try:

        geographic_north_direction_azimuth_cw, err = calculate_local_geographic_north_azimuth_from_epsg_4326(
            longitude=local_pt_longitude,
            latitude=local_pt_latitude,
            planar_qgis_crs=dest_crs,
        )

        if err:
            return None, err

        print(f"Azimuth correction (degrees, clockwise): {geographic_north_direction_azimuth_cw}")

        dest_crs_dip_direction = (dip_direction_as_north_azimuth + geographic_north_direction_azimuth_cw) % 360.0
        dest_crs_dip_angle = source_dip_angle

        """ 20230827: not sure whether the following processing, modified as needed, as necessary
        
        src_crs_direction = Direct(
            az=src_crs_plane.dipazim,
            pl=0.0)

        src_crs_direction_versor = src_crs_direction.as_versor()

        src_crs_direction_dummy_vector = src_crs_direction_versor.scale(dummy_distance)

        # NOTE: it assumes that the source CRS is not in lon-lat !!

        src_crs_shifted_pt = src_crs_local_pt.shift2d_by_vect(
            v=src_crs_direction_dummy_vector
        )

        src_crs_dummy_direction_end_pt_x, src_crs_dummy_direction_end_pt_y = src_crs_shifted_pt.xy()

        src_crs_dummy_z_range = dummy_distance * tan(radians(src_crs_plane.dipang))

        dest_crs_local_pt_x, dest_crs_local_pt_y = qgs_project_xy(
            x=src_crs_local_pt.x,
            y=src_crs_local_pt.y,
            src_crs=src_crs,
            dest_crs=dest_crs,
        )

        dest_crs_local_pt = Point(
            dest_crs_local_pt_x,
            dest_crs_local_pt_y
        )

        dest_crs_dummy_direction_end_pt_x, dest_crs_dummy_direction_end_pt_y = qgs_project_xy(
            x=src_crs_dummy_direction_end_pt_x,
            y=src_crs_dummy_direction_end_pt_y,
            src_crs=src_crs,
            dest_crs=dest_crs,
        )

        dest_crs_dummy_direction_end_pt = Point(
            dest_crs_dummy_direction_end_pt_x,
            dest_crs_dummy_direction_end_pt_y
        )

        # geoplane attitude in destination CRS

        dest_crs_dummy_direction_segment = Segment(
            start_pt=dest_crs_local_pt,
            end_pt=dest_crs_dummy_direction_end_pt,
        )

        dest_crs_dummy_direction_vector = dest_crs_dummy_direction_segment.as_vector2d()

        dest_crs_dip_direction = dest_crs_dummy_direction_vector.azimuth_degr()

        dest_crs_dummy_direction_vector_length = dest_crs_dummy_direction_vector.length
        dest_crs_dummy_z_range = src_crs_dummy_z_range

        dest_crs_dip_angle = degrees(atan(dest_crs_dummy_z_range / dest_crs_dummy_direction_vector_length))
        """

        return Plane(
            dest_crs_dip_direction,
            dest_crs_dip_angle
        ), Error()

    except Exception as e:

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


def convert_plane_orientation_between_projections(
    src_crs_plane: Plane,
    src_crs_local_pt: Point,
    src_crs: QgsCoordinateReferenceSystem,
    dest_crs: QgsCoordinateReferenceSystem,
    dummy_distance: numbers.Real = 100  # meters
) -> Tuple[Union[None, Plane], Error]:

    try:

        src_crs_direction = Direct(
            az=src_crs_plane.dipazim,
            pl=0.0)

        src_crs_direction_versor = src_crs_direction.as_versor()

        src_crs_direction_dummy_vector = src_crs_direction_versor.scale(dummy_distance)

        # NOTE: it assumes that the source CRS is not in lon-lat !!

        src_crs_shifted_pt = src_crs_local_pt.shift2d_by_vect(
            v=src_crs_direction_dummy_vector
        )

        src_crs_dummy_direction_end_pt_x, src_crs_dummy_direction_end_pt_y = src_crs_shifted_pt.xy()

        src_crs_dummy_z_range = dummy_distance * tan(radians(src_crs_plane.dipang))

        dest_crs_local_pt_x, dest_crs_local_pt_y = qgs_project_xy(
            x=src_crs_local_pt.x,
            y=src_crs_local_pt.y,
            src_crs=src_crs,
            dest_crs=dest_crs,
        )

        dest_crs_local_pt = Point(
            dest_crs_local_pt_x,
            dest_crs_local_pt_y
        )

        dest_crs_dummy_direction_end_pt_x, dest_crs_dummy_direction_end_pt_y = qgs_project_xy(
            x=src_crs_dummy_direction_end_pt_x,
            y=src_crs_dummy_direction_end_pt_y,
            src_crs=src_crs,
            dest_crs=dest_crs,
        )

        dest_crs_dummy_direction_end_pt = Point(
            dest_crs_dummy_direction_end_pt_x,
            dest_crs_dummy_direction_end_pt_y
        )

        # geoplane attitude in destination CRS

        dest_crs_dummy_direction_segment = Segment(
            start_pt=dest_crs_local_pt,
            end_pt=dest_crs_dummy_direction_end_pt,
        )

        dest_crs_dummy_direction_vector = dest_crs_dummy_direction_segment.as_vector2d()

        dest_crs_dip_direction = dest_crs_dummy_direction_vector.azimuth_degr()

        dest_crs_dummy_direction_vector_length = dest_crs_dummy_direction_vector.length
        dest_crs_dummy_z_range = src_crs_dummy_z_range

        dest_crs_dip_angle = degrees(atan(dest_crs_dummy_z_range / dest_crs_dummy_direction_vector_length))

        return Plane(
            dest_crs_dip_direction,
            dest_crs_dip_angle
        ), Error()

    except Exception as e:

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


def correct_plane_azimuth_to_geographic_north(
    src_crs_plane: Plane,
    src_crs_local_pt: Point,
    src_crs: QgsCoordinateReferenceSystem,
) -> Tuple[Union[None, Plane], Error]:

    try:

        north_orientation, err = calculate_local_geographic_north_azimuth_from_crs(
            current_qgis_crs=src_crs,
            current_crs_local_pt=src_crs_local_pt,
        )

        if err:
            return None, err

        corrected_north_plane_azimuth = (src_crs_plane.dipazim - north_orientation) % 360.0

        corrected_plane = Plane(
            azim=corrected_north_plane_azimuth,
            dip_ang=src_crs_plane.dipang
        )

        return corrected_plane, Error()

    except Exception as e:

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



