from scipy import interpolate
import dem
import ezdxf
from data_types import WGS84Point
import transformer
import numpy as np
from functools import singledispatch
from typing import Tuple


class Plane:
    # dem文件坐标系统
    def __init__(self, dem_system, plane_system, dem_path):
        self._dem_path = dem_path
        self._dem_system = dem_system
        self._start_84_point = None
        self._end_84_point = None
        self._plane_system = plane_system
        self._transformer_84_to_dem_ = transformer.Projector(
            transformer.EPSG_WGS_84, dem_system
        )
        self._transformer_dem_to_84 = transformer.Projector(
            dem_system, transformer.EPSG_WGS_84
        )

    @singledispatch
    def set_line(self, start: WGS84Point, end: WGS84Point):
        self._start_84_point = start
        self._end_84_point = end
        pass

    # def set_line_xy(self, start: Tuple[float], end: Tuple[float]):
    @set_line.register(tuple)
    def set_line_yx(self, start: Tuple[float, float], end: Tuple[float, float]):
        _transformer = self._transformer_dem_to_84
        start_wgs_84 = self._transformer_dem_to_84.transform(*start)
        end_wgs_84 = self._transformer_dem_to_84.transform(*end)
        self._start_84_point = WGS84Point(start_wgs_84[0], start_wgs_84[1], 0)
        self._end_84_point = WGS84Point(end_wgs_84[0], end_wgs_84[1], 0)
        pass

    def get_elevation(self, n):
        point_start = self._start_84_point
        point_end = self._end_84_point
        all_wgs_84_point = self._generate_wgs_84_point(point_start, point_end, n)
        # 转换为DEM的坐标
        transformer_dem_to_84 = self._transformer_84_to_dem_
        all_dem_point = np.array(
            [
                transformer_dem_to_84.transform(y_lat, x_lon)
                for y_lat, x_lon in all_wgs_84_point[:, [0, 1]]
            ]
        )
        dem_gcs = dem.Dem(self._dem_path)
        elevation_in_line = dem_gcs.get_elevation(all_dem_point)
        return np.array(elevation_in_line)

        # 计算档距

    def cal_span(self):
        start = self._start_84_point
        end = self._end_84_point
        _transformer_84_yx = transformer.Projector(
            transformer.EPSG_WGS_84, transformer.EPSG_Xian_1980_114E
        )
        start_point = np.array(
            _transformer_84_yx.transform(start.latitude, start.longitude)
        )
        end_point = np.array(_transformer_84_yx.transform(end.latitude, end.longitude))
        return int(np.round(np.sum((start_point - end_point) ** 2) ** 0.5))

    def _generate_wgs_84_point(
        self, point_start: WGS84Point, point_end: WGS84Point, n_point
    ):
        plane_system = self._plane_system
        _transformer_84_yx = transformer.Projector(
            transformer.EPSG_WGS_84, plane_system
        )
        start_y, start_x = _transformer_84_yx.transform(
            point_start.latitude, point_start.longitude
        )
        end_y, end_x = _transformer_84_yx.transform(
            point_end.latitude, point_end.longitude
        )
        # 插入点
        x_linspace = np.linspace(start_x, end_x, n_point)
        y_linspace = np.linspace(start_y, end_y, n_point)
        elevation_linspace = np.linspace(
            point_start.elevation, point_end.elevation, n_point
        )
        yxz_point = np.array([y_linspace, x_linspace, elevation_linspace]).transpose()
        _transformer_yx_84 = transformer.Projector(
            plane_system, transformer.EPSG_WGS_84
        )
        wgs_84_point = np.array(
            [
                np.array(_transformer_yx_84.transform(y, x))
                for y, x, _elevation in yxz_point
            ]
        )
        return wgs_84_point


def draw_dxf(elevation_in_line, interval, dxf_path):
    doc = ezdxf.new(dxfversion="R2004")
    msp = doc.modelspace()
    x_axis = (
        np.linspace(0, len(elevation_in_line) * interval, len(elevation_in_line)) / 5
    )
    points = np.array([x_axis, (elevation_in_line - elevation_in_line[0]) * 2])
    new_x_axis = (
        np.linspace(
            0, len(elevation_in_line) * interval, int(len(elevation_in_line) / 12.5)
        )
        / 5
    )  # 平滑用
    func = interpolate.interp1d(points[0, :], points[1, :], kind="cubic")
    # msp.add_polyline2d([(x + 10, func(x)) for x in new_x_axis])
    msp.add_polyline2d([(x + 10, y) for x, y in points.transpose()])
    doc.saveas(dxf_path)


def main():
    _plane = Plane(
        transformer.EPSG_Xian_1980_114E, transformer.EPSG_Xian_1980_114E, r"d:\307.tif"
    )
    # start_point = WGS84Point(33.24643289, 114.58574497, 28)
    # end_point = WGS84Point(33.22873308, 114.57951634, 28)
    _plane.set_line_yx((3680231.27, 554470.91), (3678264.93, 553901.32))
    span = _plane.cal_span()
    print("档距 {span}".format(span=span))
    elevation_in_line = _plane.get_elevation(int(span / 5))
    # dem_gcs = dem.Dem()
    # dem_gcs.get_dem_info(if_print=False)
    #
    # span = int(cal_span(start_point, end_point))
    # print("档距 {span}".format(span=span))
    # wgs_84_point_in_line = generate_wgs_84_point(start_point, end_point, span)
    # elevation_in_line = dem_gcs.get_elevation(wgs_84_point_in_line)
    draw_dxf(elevation_in_line, 5, "a.dxf")