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


from osgeo import ogr,gdal
from osgeo.ogr import *
from test.MapSynthesize.ShapeData import ShapeData
import copy
import math
from rtreelib import RTree,Rect



class RTreeData:

    def __init__(self,triangle,state):
        self.state = state
        self.triangle = triangle

    def set_state(self,state):
        self.state = state


def env2rect(env):
    return Rect(env[0],env[2],env[1],env[3])


class PolygonMerge:

    @classmethod
    def get_delauney(cls,polygon):

        polygon_line = cls.get_polygon_lines(polygon)
        polygon_line_points:list = polygon_line.GetPoints()

        if len(polygon_line_points)==4:
            tris = [polygon]
        else:
            d = 0
            for ind in range(len(polygon_line_points) - 1):
                p = polygon_line_points[ind]
                p_n = polygon_line_points[ind + 1]

                d += -0.5 * (p_n[1] + p[1]) * (p_n[0] - p[0])

            if d < 0:
                polygon_line_points.reverse()

            tris = list(cls.simple_polygon_division(polygon_line_points))

        return tris

    @classmethod
    def simple_polygon_division(cls,polygon_line_points:list):
        '''
        简单多边形三角剖分算法
        :param polygon_line_points: 多边形逆时针点序列
        :return:
        '''

        if len(polygon_line_points) == 4:
            yield cls.create_polygon(polygon_line_points)
        else:
            for index,point in enumerate(polygon_line_points):
                try:
                    if index == 0:
                        point_front = polygon_line_points[-1]
                        point_next = polygon_line_points[1]
                    elif index == len(polygon_line_points) - 1:
                        point_front = polygon_line_points[index - 1]
                        point_next = polygon_line_points[0]
                    else:
                        point_front = polygon_line_points[index - 1]
                        point_next = polygon_line_points[index + 1]

                except:
                    dd=1
                check = (point[0] - point_front[0]) * (point_next[1] -point[1]) - (point_next[0] - point[0]) * (point[1] - point_front[1])

                #如果是凸点
                if check > 0:
                    triangle = cls.create_polygon([point_front,point,point_next,point_front])

                    #判断是否包含其他点
                    contains = False
                    for i,poi_xy in enumerate(polygon_line_points):

                        if index == 0:
                            if i in [0, 1, len(polygon_line_points) - 1]:
                                continue

                        elif index == len(polygon_line_points) - 1:
                            if i in [index, 0 ,index - 1]:
                                continue
                        else:
                            if i in [index, index-1, index+1]:
                                continue

                        poi:Geometry = ogr.Geometry(ogr.wkbPoint)
                        poi.AddPoint(poi_xy[0],poi_xy[1])
                        if triangle.Contains(poi):
                            contains = True

                    #如果不包含其他点，为有效剖分,结束循环
                    if not contains:
                        yield triangle
                        polygon_line_points.pop(index)
                        for tri in cls.simple_polygon_division(polygon_line_points):
                            yield tri
                        break


    @classmethod
    def merge(cls,polygons:list,distance_buffer_threshold):

        rtree = RTree()

        polygon:Geometry = polygons[0]
        # multipolygon = ogr.Geometry(ogr.wkbMultiPolygon)
        # for poly in polygons:
        #     multipolygon.AddGeometry(poly)

        tris = []
        for poly in polygons:
            polygon = polygon.Union(poly)

        delauney: Geometry = polygon.DelaunayTriangulation()

        for index in range(delauney.GetGeometryCount()):
            de:Geometry = copy.deepcopy(delauney.GetGeometryRef(index))
            de_extent = de.GetEnvelope()
            rtree.insert(RTreeData(de,1), env2rect(de_extent))
            tris.append(de)

        bians: list = cls.get_bian(polygon)

        for bian in bians:

            bian_ext = bian.GetEnvelope()
            query_result = list(rtree.query(env2rect(bian_ext)))

            intersection_data = []

            for entry in query_result:

                # 已被抛弃的不要
                if entry.data.state == 0:
                    continue
                intersection_geom: Geometry = bian.Intersection(entry.data.triangle)


                if not intersection_geom.IsEmpty():
                    if not intersection_geom.Equals(bian):
                        if intersection_geom.GetGeometryType() in [2, 5, -2147483643, -2147483646, 3002, 3005]:
                            intersection_data.append(entry.data)



            for data in intersection_data:
                data.set_state(0)

            if intersection_data:


                inter_tri_merge:Geometry = intersection_data[0].triangle

                for ind in range(1,len(intersection_data)):
                    inter_tri_merge = inter_tri_merge.Union(intersection_data[ind].triangle)

                try:
                    inter_tri_merge_line:Geometry = cls.get_polygon_lines(inter_tri_merge)

                    inter_tri_merge_line_points:list = inter_tri_merge_line.GetPoints()

                    bian_ps = bian.GetPoints()

                    if len(inter_tri_merge_line_points[0]) == 3:
                        first_index = inter_tri_merge_line_points.index(bian_ps[0])
                        second_index = inter_tri_merge_line_points.index(bian_ps[1])
                    else:

                        first_index = inter_tri_merge_line_points.index((bian_ps[0][0],bian_ps[0][1]))
                        second_index = inter_tri_merge_line_points.index((bian_ps[1][0],bian_ps[1][1]))

                    small_index = min(first_index,second_index)
                    big_index = max(first_index,second_index)

                    small_points:list = inter_tri_merge_line_points[:small_index+1]
                    small_points.extend(inter_tri_merge_line_points[big_index:])
                    small_polygon:Geometry = cls.create_polygon(small_points)

                    delauney = cls.get_delauney(small_polygon)
                    for de in delauney:
                        de_extent = de.GetEnvelope()
                        rtree.insert(RTreeData(de, 1), env2rect(de_extent))

                    big_points = inter_tri_merge_line_points[small_index:big_index+1]
                    big_points.append(inter_tri_merge_line_points[small_index])

                    big_polygon:Geometry = cls.create_polygon(big_points)

                    delauney = cls.get_delauney(big_polygon)
                    for de in delauney:
                        de_extent = de.GetEnvelope()
                        rtree.insert(RTreeData(de, 1), env2rect(de_extent))
                except:
                    print(inter_tri_merge)
                    print(bian)


        for entry in rtree.get_leaf_entries():
            triangle = entry.data.triangle
            if entry.data.state==0:
                continue
            if polygon.Contains(triangle):
                pass
            else:
                centers = cls.get_center(triangle)
                # querycenter :Geometry= tri_center.Buffer(distance_buffer_threshold)

                # query_result = list(rtree.query(env2rect(querycenter.GetEnvelope())))

                pass_ = False
                for center in centers:
                    distance = center.Distance(polygon)
                    if distance > distance_buffer_threshold :
                        pass_ = True
                        break

                if pass_:
                    continue

                # if len(query_result)==0:
                #     continue
                polygon:Geometry = polygon.Union(triangle)

                # polygon.AddGeometry(triangle)

        # result: Geometry = polygon

        # 去岛
        result :Geometry = cls.remove_hole(polygon,distance_buffer_threshold)

        return result

    @classmethod
    def remove_hole(cls,polygon:Geometry,distance_buffer_threshold):

        polygons = []

        if polygon.GetGeometryType() in [6, -2147483642, 3006]:
            for index in range(polygon.GetGeometryCount()):
                each = polygon.GetGeometryRef(index)
                polygons.append(each)
        else:
            polygons.append(polygon)

        polygons = [cls.remove_hole_simple(p,distance_buffer_threshold) for p in polygons]

        if len(polygons) == 1:
            return polygons[0]
        else:
            mp:Geometry = ogr.Geometry(ogr.wkbMultiPolygon)
            for p in polygons:
                mp.AddGeometry(p)
            return mp

    @classmethod
    def remove_hole_simple(cls,polygon:Geometry,distance_buffer_threshold):

        if polygon.GetGeometryCount() >1:
            for index in range(1,polygon.GetGeometryCount()):

                each:Geometry = copy.deepcopy(polygon.GetGeometryRef(index))

                if each.GetArea() < distance_buffer_threshold*distance_buffer_threshold*5:
                    polygon.RemoveGeometry(index)
                    cls.remove_hole_simple(polygon,distance_buffer_threshold)
                    break

        return polygon

    @classmethod
    def remove_hole_simple2(cls,polygon:Geometry,distance_buffer_threshold):

        if polygon.GetGeometryCount() >1:
            for index in range(1,polygon.GetGeometryCount()):

                each:Geometry = polygon.GetGeometryRef(index)
                print(polygon.GetGeometryCount())
                if each:
                    if each.GetArea() < distance_buffer_threshold*distance_buffer_threshold*18:
                        polygon.RemoveGeometry(index)

        return polygon

    @classmethod
    def get_center2(cls,tri:Geometry):

        tri_line: Geometry = cls.get_polygon_lines(tri)
        tri_line_points = tri_line.GetPoints()



        max = 0
        poi_tuple = []


        center1 = ogr.Geometry(ogr.wkbPoint)
        center2 = ogr.Geometry(ogr.wkbPoint)
        center3 = ogr.Geometry(ogr.wkbPoint)


        for index in range(len(tri_line_points)-1):

            first = tri_line_points[index]
            second = tri_line_points[index+1]

            dis = math.sqrt(math.pow(first[0]-second[0],2) + math.pow(first[1]-second[1],2))

            if dis > max:
                max = dis
                poi_tuple = [first, second]



        center1.AddPoint((poi_tuple[0][0]+poi_tuple[1][0])/2,(poi_tuple[0][1]+poi_tuple[1][1])/2)
        center2.AddPoint((poi_tuple[0][0] + center1.GetX()) / 2, (poi_tuple[0][1] + center1.GetY()) / 2)
        center3.AddPoint((poi_tuple[1][0] + center1.GetX()) / 2, (poi_tuple[1][1] + center1.GetY()) / 2)



        return [center1,center2,center3]

    @classmethod
    def get_center(cls, tri: Geometry):

        center_set = set()

        tri_line: Geometry = cls.get_polygon_lines(tri)
        tri_line_points = tri_line.GetPoints()



        max = 0
        poi_tuple = []


        center1 = ogr.Geometry(ogr.wkbPoint)
        center2 = ogr.Geometry(ogr.wkbPoint)
        center3 = ogr.Geometry(ogr.wkbPoint)


        for index in range(len(tri_line_points)-1):

            first = tri_line_points[index]
            second = tri_line_points[index+1]

            dis = math.sqrt(math.pow(first[0]-second[0],2) + math.pow(first[1]-second[1],2))

            if dis > max:
                max = dis
                poi_tuple = [first, second]



        center1.AddPoint((poi_tuple[0][0]+poi_tuple[1][0])/2,(poi_tuple[0][1]+poi_tuple[1][1])/2)
        center2.AddPoint((poi_tuple[0][0] + center1.GetX()) / 2, (poi_tuple[0][1] + center1.GetY()) / 2)
        center3.AddPoint((poi_tuple[1][0] + center1.GetX()) / 2, (poi_tuple[1][1] + center1.GetY()) / 2)

        center_set.add(center1)
        center_set.add(center2)
        center_set.add(center3)

        center4 = ogr.Geometry(ogr.wkbPoint)
        center5 = ogr.Geometry(ogr.wkbPoint)
        center6 = ogr.Geometry(ogr.wkbPoint)

        center4.AddPoint((tri_line_points[0][0] + tri_line_points[1][0]) / 2, (tri_line_points[0][1] + tri_line_points[1][1]) / 2)
        center5.AddPoint((tri_line_points[1][0] + tri_line_points[2][0]) / 2, (tri_line_points[1][1] + tri_line_points[2][1]) / 2)
        center6.AddPoint((tri_line_points[2][0] + tri_line_points[3][0]) / 2, (tri_line_points[2][1] + tri_line_points[3][1]) / 2)

        center_set.add(center4)
        center_set.add(center5)
        center_set.add(center6)

        return center_set

    @classmethod
    def get_polygon_lines(cls,polygon):
        #无孔Polygon
        if polygon.GetGeometryCount() in [0,1]:
            polygon_line : Geometry = ogr.ForceToLineString(polygon)
        # 有孔Polygon
        else:
            polygon_line : Geometry = ogr.ForceToLineString(polygon.GetGeometryRef(0))

        return polygon_line

    @classmethod
    def get_bian(cls,polygon):

        bians = []

        polygon_lines = []

        if polygon.GetGeometryCount() in [0, 1]:
            polygon_line: Geometry = ogr.ForceToLineString(polygon)
            polygon_lines.append(polygon_line)
        # 有孔Polygon
        else:
            for index in range(polygon.GetGeometryCount()):
                polygon_line: Geometry = ogr.ForceToLineString(polygon.GetGeometryRef(index))
                polygon_lines.append(polygon_line)

        for polygon_line in polygon_lines:
            for index in range(polygon_line.GetPointCount()-1):
                p = polygon_line.GetPoint(index)
                p_next = polygon_line.GetPoint(index+1)

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

                if len(p) ==2:
                    line.AddPoint(p[0], p[1],0)
                    line.AddPoint(p_next[0], p_next[1],0)
                else:
                    line.AddPoint(p[0], p[1])
                    line.AddPoint(p_next[0], p_next[1])
                bians.append(line)

        return bians

    @classmethod
    def create_polygon(cls,ps):

        ring = ogr.Geometry(ogr.wkbLinearRing)
        for p in ps:
            ring.AddPoint(p[0], p[1])
        poly = ogr.Geometry(ogr.wkbPolygon)
        poly.AddGeometry(ring)
        return poly


if __name__ == '__main__':
    import time
    t1 = time.time()


    sd = ShapeData(r"J:\Data\论文\rh.shp")

    polygons = sd.get_polygons()

    p_rh = [PolygonMerge.remove_hole(p,1.114691025217302e-04) for p in polygons]

    wkts = [p.ExportToWkt() for p in p_rh]

    result = r"J:\Data\论文\rhresult.shp"

    ShapeData.create_shp_fromwkts(result,"zh",wkts)