# coding=utf-8
#author:        4N
#createtime:    2022/3/8
#email:         nheweijun@sina.com


from osgeo import ogr
from osgeo.ogr import *
import math
import os
import copy
os.environ["CUDA_VISIBLE_DEVICES"] = "-1"
class util:
    @classmethod
    def envelop_2_polygon(cls,env):
        ring = ogr.Geometry(ogr.wkbLinearRing)
        ring.AddPoint(env[0], env[2])
        ring.AddPoint(env[0], env[3])
        ring.AddPoint(env[1], env[3])
        ring.AddPoint(env[1], env[2])
        ring.AddPoint(env[0], env[2])
        # Create polygon
        poly = ogr.Geometry(ogr.wkbPolygon)
        poly.AddGeometry(ring)
        return poly


class RepresentPoint(Geometry):
    user_data = None
    def SetUserData(self,object):
        self.user_data = object
    def GetUserData(self):
        return self.user_data

class DelauneyNode:

    point = None
    node = None
    children_nodes = []
    children_nodes_distance = []
    polygon = None

    def __init__(self, point):
        self.point = point

    def get_all_node(self):
        for node in self.children_nodes:
            yield node.get_all_node()
        yield self.node


    def cut_tree(self,threshold):

        trees = []

        for node in self.children_nodes:
            yield node.get_all_node()
        yield self.node


class ShapeUtil:

    driver = ogr.GetDriverByName("ESRI Shapefile")

    @classmethod
    def create_by_layer(cls,path,layer:Layer):
        data_source: DataSource = cls.driver.CreateDataSource(path)
        data_source.CopyLayer(layer,layer.GetName())
        data_source.Destroy()

    @classmethod
    def create_by_scheme(cls,path,name,sr,geo_type,scheme,features):
        data_source: DataSource = cls.driver.CreateDataSource(path)
        layer :Layer = data_source.CreateLayer(name, sr, geo_type)
        if scheme:
            layer.CreateFields(scheme)
        for feature in features:
            layer.CreateFeature(feature)
        data_source.Destroy()

    @classmethod
    def create_point(cls,path,name,point):
        data_source: DataSource = cls.driver.CreateDataSource(path)
        layer :Layer = data_source.CreateLayer(name, None, ogr.wkbPoint)

        feat_new = ogr.Feature(layer.GetLayerDefn())
        feat_new.SetGeometry(point)
        layer.CreateFeature(feat_new)
        data_source.Destroy()

    @classmethod
    def create_shp_fromwkts(cls,path,name,wkts):

        geo_type = None
        geoms = []
        for wkt in wkts:
            geom : Geometry = ogr.CreateGeometryFromWkt(wkt)
            if geo_type is None:
                geo_type = geom.GetGeometryType()
            geoms.append(geom)

        if os.path.exists(path):

            pre_name = ".".join(path.split(".")[0:-1])
            for bac in ["dbf","prj","cpg","shp","shx","sbn","sbx"]:
                try:
                    os.remove(pre_name+"."+bac)
                except Exception as e:
                    pass

        data_source: DataSource = cls.driver.CreateDataSource(path)
        layer :Layer = data_source.CreateLayer(name, None, geo_type)

        for geom in geoms:
            feat_new = ogr.Feature(layer.GetLayerDefn())
            feat_new.SetGeometry(geom)
            layer.CreateFeature(feat_new)
        data_source.Destroy()


class ShapeData:

    def __init__(self,path):
        driver: Driver = ogr.GetDriverByName("ESRI Shapefile")
        self.ds: DataSource = driver.Open(path, 0)
        if not self.ds:
            raise Exception("打开数据失败！")
        self.layer: Layer = self.ds.GetLayer(0)


    def get_polygons(self):

        polygons = []
        for feature in self.layer:
            f:Feature = feature
            geom : Geometry = copy.deepcopy(f.GetGeometryRef())
            if geom.GetGeometryType() == 3 or geom.GetGeometryType() == -2147483645:
                polygons.append(geom)
            if geom.GetGeometryType() == 6 or geom.GetGeometryType() == -2147483642:
                for i in range(geom.GetGeometryCount()):
                    polygons.append(geom.GetGeometryRef(i))
        return polygons

    def close(self):
        self.ds.Destroy()


class Zonghe:

    data:ShapeData
    area_threshold = 100.0
    buffer_threshold = 100.0

    def __init__(self,data):
        self.data = data



    #只处理Polygon，MultiPolygon要拆开
    def get_polygon_reprecent_point(self,polygon:Geometry):

        polygon_env = polygon.GetEnvelope()

        area_thre = polygon.GetArea()/self.area_threshold

        # line_length = max(polygon_env[1]-polygon_env[0],polygon_env[3]-polygon_env[2])

        line_length = math.sqrt(math.pow(polygon_env[1] - polygon_env[0],2) + math.pow(polygon_env[3] - polygon_env[2],2))


        short_length = min(polygon_env[1]-polygon_env[0],polygon_env[3]-polygon_env[2])



        center:Geometry = polygon.Centroid()
        # print(center)


        # 提前return
        if polygon.Contains(center):
            return center

        last_state = None
        line = None
        calculate_time = 0

        direction = self.get_direction(polygon, center, line_length)

        # print(direction)
        in_hole = False
        in_hole_line = None


        #无孔Polygon
        if polygon.GetGeometryCount() == 1:
            polygon_line = ogr.ForceToLineString(polygon)

        # 有孔Polygon,要分2种情况
        else:
            polygon_line = ogr.ForceToLineString(polygon.GetGeometryRef(0))
            for i in range(1,polygon.GetGeometryCount()):
                hole_line = polygon.GetGeometryRef(i)
                hole_polygon : Geometry = ogr.ForceToPolygon(hole_line)
                if hole_polygon.Contains(center):
                    in_hole = True
                    direction = range(0,190,10)
                    in_hole_line = ogr.ForceToLineString(hole_line)
                    # print("in_hole")

        angle1,angle2=None,None

        for angle in direction:

            if in_hole:
                line = self.get_doubleline(center.GetX(),center.GetY(),angle,line_length)
            else:
                line = self.get_line(center.GetX(),center.GetY(),angle,line_length)

            oneside_area, otherside_area = self.jude(polygon, line, angle, line_length)
            # print(oneside_area, otherside_area)

            # print("\"" + line.ExportToWkt() + "\",")

            if last_state is None:

                angle2 = angle

                if oneside_area == 0 and otherside_area==0:
                    last_state = 0
                else:
                    last_state = 1 if oneside_area < otherside_area else -1
            else:

                angle1 = angle2
                angle2 = angle

                if oneside_area == 0 and otherside_area==0:
                    now_state = 0
                else:
                    now_state = 1 if oneside_area < otherside_area else -1

                if last_state* now_state < 0:
                    # print(angle)
                    break
                else:
                    last_state = now_state


        # print("start division")
        # print(angle1, angle2)

        line = self.division(polygon, center, in_hole, line_length, angle1, angle2, area_thre)

        # 提前return
        if in_hole:
            # print(line)
            # print(in_hole_line)
            point_intersection: Geometry = line.Intersection(in_hole_line)
            plist = []
            for i in range(point_intersection.GetGeometryCount()):
                pi = point_intersection.GetGeometryRef(i)
                plist.append(pi)
            plist_d = [(pi,center.Distance(pi)) for pi in plist]
            plist_d.sort(key=lambda x:x[1])

            nearest_point= plist_d[0][0]
            # print(line_length)
            # print(short_length)
            # print(nearest_point)

            firt_short_length = short_length/20.0
            # print(firt_short_length)

            radline,other_point = self.get_radline(center.GetX(),center.GetY(),nearest_point.GetX(),nearest_point.GetY(),firt_short_length)

            # print(other_point)

            while not polygon.Contains(other_point):
                firt_short_length = firt_short_length / 2
                radline, other_point = self.get_radline(center.GetX(), center.GetY(), nearest_point.GetX(),
                                                        nearest_point.GetY(), firt_short_length)

            # print("return hole_rep")

            return other_point


        point_intersection:Geometry = line.Intersection(polygon_line)
        xlist = []
        ylist = []
        plist = []

        for i in range(point_intersection.GetGeometryCount()):
            pi = point_intersection.GetGeometryRef(i)
            xlist.append(pi.GetX())
            ylist.append(pi.GetY())
            plist.append(pi)

        reprecent_point = [(min(xlist)+max(xlist))/2.0,(min(ylist)+max(ylist))/2.0]

        reprecent_point = self.creat_point(reprecent_point)

        if not polygon.Contains(reprecent_point):
            plist_d = [(pi,center.Distance(pi)) for pi in plist]
            plist_d.sort(key=lambda x:x[1])
            reprecent_point = self.creat_point(((plist_d[0][0].GetX()+ plist_d[1][0].GetX())/2.0, (plist_d[0][0].GetY()+plist_d[1][0].GetY())/2.0))

        print(reprecent_point)

        return reprecent_point


    def create_delauney(self,points):
        multipoint: Geometry = ogr.Geometry(ogr.wkbMultiPoint)
        for point in points:
            multipoint.AddGeometry(point)
        delauney =  multipoint.DelaunayTriangulation()





    # def delauney_min_tree(self,points,delauney):
    #     lines = []
    #     for i in range(delauney.GetGeometryCount()):
    #         triangle:Geometry = delauney.GetGeometryRef(i)
    #
    #         triangle_line :Geometry = ogr.ForceToLineString(triangle)
    #         triangle_line_points = triangle_line.GetPoints()
    #
    #         distance =
    #         bian1 = (triangle_line_points[0],triangle_line_points[1],)



    def get_line(self, x0, y0 , angle, r):

        line: Geometry = ogr.Geometry(ogr.wkbLineString)
        line.AddPoint(x0, y0)
        hudu = angle/360.0 * 2 * math.pi
        dx = math.sin(hudu) * r
        dy = math.cos(hudu) * r
        line.AddPoint(x0+dx, y0+dy)

        return line

    def get_doubleline(self, x0, y0 , angle, r):

        line: Geometry = ogr.Geometry(ogr.wkbLineString)

        hudu = angle/360.0 * 2 * math.pi
        dx = math.sin(hudu) * r
        dy = math.cos(hudu) * r

        hudu = (angle+180)/360.0 * 2 * math.pi
        dx2 = math.sin(hudu) * r
        dy2 = math.cos(hudu) * r

        line.AddPoint(x0 + dx2, y0 + dy2)
        line.AddPoint(x0+dx, y0+dy)

        return line


    def get_radline(self,x0,y0,x1,y1,len_size):

        line: Geometry = ogr.Geometry(ogr.wkbLineString)
        line.AddPoint(x0, y0)
        tanx = (y1-y0) / (x1-x0)
        cosx = 1 / math.sqrt(1 + tanx *tanx)
        dx = cosx* len_size



        x2 = x1 + dx

        if x1 != x0:
            y2 = ((y1-y0)*(x2-x0))/(x1-x0) + y0
        else:
            x2 = x0
            y2 = y0 + dx
        line.AddPoint(x2, y2)

        other_point = ogr.Geometry(ogr.wkbPoint)
        other_point.AddPoint(x2, y2)

        return line , other_point

    def get_direction(self,polygon,center,line_length):
        '''
        判断旋转方向
        :param polygon:
        :param center:
        :param line_length:
        :return:
        '''

        line = self.get_line(center.GetX(),center.GetY(),0,line_length)

        oneside_area, otherside_area = self.jude(polygon,line,0,line_length)

        if oneside_area == 0 and otherside_area == 0 :
            return range(0, 390, 30)
        else:
            line = self.get_line(center.GetX(), center.GetY(), 10, line_length)
            oneside_area_contrast, otherside_area_contrast = self.jude(polygon, line, 30, line_length)

            if oneside_area_contrast == 0 and otherside_area_contrast == 0:
                return range(360, -30, -30)
            else:
                if abs(oneside_area - otherside_area) > abs(oneside_area_contrast - otherside_area_contrast):
                    return range(0, 390, 30)
                else:
                    return range(360, -30, -30)


    def jude(self,polygon,line,angle,line_length):

        line_buffer: Geometry = line.Buffer(line_length / self.buffer_threshold)
        clip_geom: Geometry = polygon.Difference(line_buffer)

        oneside_area = 0
        otherside_area = 0

        gc = clip_geom.GetGeometryCount()
        if gc <= 1:
            oneside_area = 0
            otherside_area = 0
        else:
            triangle1, triangle2 = self.get_side_triangle(line, angle)
            # print(angle)
            # print(line)
            # print(triangle1)
            # print(triangle2)

            for gi in range(clip_geom.GetGeometryCount()):
                clip_geom_i: Geometry = clip_geom.GetGeometryRef(gi)
                it = clip_geom_i.Intersection(triangle1)

                if clip_geom_i.Intersect(triangle1) and it.Buffer(line_length / self.buffer_threshold).Intersect(
                        line_buffer):
                    oneside_area += clip_geom_i.GetArea()
                else:
                    otherside_area += clip_geom_i.GetArea()
        return oneside_area, otherside_area

    def division(self,polygon,center,in_hole,line_length,angle1,angle2,area_thre):

        mid_angle = (angle1 + angle2)/2.0
        if in_hole:
            line = self.get_doubleline(center.GetX(), center.GetY(), mid_angle, line_length)
        else:
            line = self.get_line(center.GetX(), center.GetY(), mid_angle, line_length)

        # print("\""+line.ExportToWkt()+"\",")

        one,other = self.jude(polygon,line,mid_angle,line_length)
        # print(mid_angle)
        # print(one,other)

        while abs(one - other) > area_thre:
            if one > other:
                if angle1 > angle2:
                    return self.division(polygon, center, in_hole, line_length, angle2, mid_angle, area_thre)
                else:
                    return self.division(polygon, center, in_hole, line_length, mid_angle, angle1, area_thre)
            else:
                if angle1 > angle2:
                    return self.division(polygon, center, in_hole, line_length, mid_angle, angle1, area_thre)
                else:
                    return self.division(polygon, center, in_hole, line_length, angle2, mid_angle, area_thre)
        return line


    def get_side_triangle(self,line:Geometry,angle):
        '''
        获取方向三角形
        :param line:
        :param angle:
        :return:
        '''

        start = line.GetPoint(0)
        end = line.GetPoint(1)


        angle_t = angle-30
        hudu = angle_t/360.0 * 2 * math.pi

        r = line.Length() / (2*math.cos( math.pi/6))
        dx = math.sin(hudu) * r
        dy = math.cos(hudu) * r


        ring1 = ogr.Geometry(ogr.wkbLinearRing)
        ring1.AddPoint(start[0], start[1])
        ring1.AddPoint(start[0] + dx, start[1] + dy)
        ring1.AddPoint(end[0], end[1])
        ring1.AddPoint(start[0], start[1])
        triangle1 = ogr.Geometry(ogr.wkbPolygon)
        triangle1.AddGeometry(ring1)


        angle_t = angle+30
        hudu = angle_t/360.0 * 2 * math.pi

        r = line.Length() / (2*math.cos( math.pi/6))

        dx = math.sin(hudu) * r
        dy = math.cos(hudu) * r

        ring2 = ogr.Geometry(ogr.wkbLinearRing)
        ring2.AddPoint(start[0], start[1])
        ring2.AddPoint(start[0]+dx, start[1]+dy)
        ring2.AddPoint(end[0], end[1])
        ring2.AddPoint(start[0], start[1])
        triangle2 = ogr.Geometry(ogr.wkbPolygon)
        triangle2.AddGeometry(ring2)

        return triangle1,triangle2


    def creat_point(self,tuple):
        p: Geometry = ogr.Geometry(ogr.wkbPoint)
        p.AddPoint(tuple[0], tuple[1])
        return p


    def merge(self,delauney_node:DelauneyNode):
        nodes =  delauney_node.get_all_node()
        points = [node.point for node in nodes]

        return None


    def process(self):
        polygons = self.data.get_polygons()

        Rps = []
        for poly in polygons:
            reprecent_point = self.get_polygon_reprecent_point(poly)




if __name__ == '__main__':


    # print(sd.get_polygons()[0])
    # sd.close()

    create = True
    # create = False

    if create:

        # driver: Driver = ogr.GetDriverByName("ESRI Shapefile")
        # ds: DataSource = driver.Open(r"J:\Data\制图综合\left.shp",0)
        #
        # layer: Layer = ds.GetLayer(0)
        #
        # polygons = []
        # for feature in layer:
        #     f:Feature = feature
        #     geom : Geometry = f.GetGeometryRef()
        #     if geom.GetGeometryType() == 3 or geom.GetGeometryType() == -2147483645:
        #         polygons.append(geom)
        #     if geom.GetGeometryType() == 6 or geom.GetGeometryType() == -2147483642:
        #         for i in range(geom.GetGeometryCount()):
        #             polygons.append(geom.GetGeometryRef(i))
        # del ds

        sd = ShapeData(r"J:\Data\制图综合\t5.shp")
        zh = Zonghe(sd)
        zh.process()



    else:


        wkts = ["POINT (113.211177054031 22.8681300466504)"
                ]

        result = r"J:\Data\制图综合\left_center.shp"
        ShapeUtil.create_shp_fromwkts(result,"zh",wkts)