# 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
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 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, 1)
        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 = 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 = 20.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()

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


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

        center:Geometry = polygon.Centroid()

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

        last_state = None
        line = None
        calculate_time = 0
        direction = self.get_direction(polygon, center, line_length)
        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")



        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)

            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:
                calculate_time += 1
                triangle1, triangle2 = self.get_side_triangle(line, angle)

                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()
            print(oneside_area,otherside_area)
            if last_state is None:
                if oneside_area == 0 and otherside_area==0:
                    last_state = 0
                else:
                    last_state = 1 if oneside_area < otherside_area else -1
            else:
                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:
                    break
                else:
                    last_state = now_state

        # 提前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))

        return reprecent_point


    def create_delauney(self,polygon:Geometry):
        de =  polygon.DelaunayTriangulation()

        for i in range(de.GetGeometryCount()):
            d:Geometry = de.GetGeometryRef(i)

            print("\""+d.ExportToWkt()+"\",")


    def create_delauney_test(self):

        mp:Geometry = ogr.Geometry(ogr.wkbMultiPoint)

        p1: Geometry = ogr.Geometry(ogr.wkbPoint)
        p1.AddPoint(1,1)

        mp.AddGeometry(p1)

        p2: Geometry = ogr.Geometry(ogr.wkbPoint)
        p2.AddPoint(1.5, 1.7)

        mp.AddGeometry(p2)


        p3: Geometry = ogr.Geometry(ogr.wkbPoint)
        p3.AddPoint(2, 10)

        mp.AddGeometry(p3)

        p4: Geometry = ogr.Geometry(ogr.wkbPoint)
        p4.AddPoint(3, 5)
        mp.AddGeometry(p4)

        p5: Geometry = ogr.Geometry(ogr.wkbPoint)
        p5.AddPoint(4, 2)
        mp.AddGeometry(p5)

        p6: Geometry = ogr.Geometry(ogr.wkbPoint)

        p6.AddPoint(3, 1)

        mp.AddGeometry(p6)

        de =  mp.DelaunayTriangulation()

        print(mp)

        for i in range(de.GetGeometryCount()):
            d:Geometry = de.GetGeometryRef(i)

            print("\""+d.ExportToWkt()+"\",")

    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)
        line_buffer:Geometry = line.Buffer(line_length / self.buffer_threshold)
        clip_geom: Geometry = polygon.Difference(line_buffer)

        gc = clip_geom.GetGeometryCount()
        oneside_area = 0
        otherside_area = 0

        if gc <= 1:
            return range(0,370,10)
        else:
            triangle1, triangle2 = self.get_side_triangle(line, 0)

            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()

            oneside_area_contrast = 0
            otherside_area_contrast = 0

            line = self.get_line(center.GetX(), center.GetY(), 10, line_length)
            line_buffer: Geometry = line.Buffer(line_length / self.buffer_threshold)
            clip_geom: Geometry = polygon.Difference(line_buffer)

            gc = clip_geom.GetGeometryCount()
            if gc <= 1:
                return range(360, -10, -10)
            else:
                triangle1, triangle2 = self.get_side_triangle(line, 10)

                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_contrast += clip_geom_i.GetArea()
                    else:
                        otherside_area_contrast += clip_geom_i.GetArea()

                if abs(oneside_area - otherside_area) > abs(oneside_area_contrast - otherside_area_contrast):
                    return range(0, 370, 10)
                else:
                    return range(360, -10, -10)

    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 process(self):
        polygons = self.data.get_polygons()
        for poly in polygons:
            reprecent_point = self.get_polygon_reprecent_point(poly)
            print(reprecent_point)



if __name__ == '__main__':

    # sd = ShapeData(r"J:\Data\矢量数据\佛山\制图综合.shp")
    # 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\制图综合\t1.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
        zh = Zonghe("d")
        pre = zh.get_polygon_reprecent_point(polygons[0])



        zh.create_delauney_test()


    else:


        wkts = ["POLYGON ((2 10 0,1 1 0,1.5 1.7 0,2 10 0))",
"POLYGON ((2 10 0,1.5 1.7 0,3 5 0,2 10 0))",
"POLYGON ((2 10 0,3 5 0,4 2 0,2 10 0))",
"POLYGON ((3 1 0,4 2 0,1.5 1.7 0,3 1 0))",
"POLYGON ((3 1 0,1.5 1.7 0,1 1 0,3 1 0))",
"POLYGON ((1.5 1.7 0,4 2 0,3 5 0,1.5 1.7 0))",]
        result = r"J:\Data\制图综合\mpde.shp"
        ShapeUtil.create_shp_fromwkts(result,"zh",wkts)