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


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 = 10.0
    buffer_threshold = 100.0

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



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

        polygon_area = polygon.GetArea()

        polygon_env = polygon.GetEnvelope()

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

        center:Geometry = polygon.Centroid()
        if polygon.Contains(center):
            return center


        left_area = 0
        right_area = polygon_area

        left_or_right = None
        fen = False

        if polygon.GetGeometryCount() == 1:
            polygon_line = ogr.ForceToLineString(polygon)
        # 有孔Polygon
        else:
            polygon_line = ogr.ForceToLineString(polygon.GetGeometryRef(0))

        # print(len(polygon_line.GetPoints()))
        point_i = 0

        line = None
        next_point = polygon_line.GetPoint(point_i)
        last_point = None
        last_last_point = None

        angle = []

        calculate_times = 0

        while abs(left_area - right_area) > (polygon_area / self.area_threshold):
            # print(calculate_times)
            calculate_times += 1

            line,other_point = self.get_line(center.GetX(),center.GetY(),next_point[0],next_point[1],line_length)

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


            if clip_geom.GetGeometryType() == 3 or clip_geom.GetGeometryType() == -2147483645:

                cl = ogr.ForceToLineString(clip_geom)
                p1 = cl.GetPoint(0)
                p2 = cl.GetPoint(1)
                tp = [(p1[0] + p2[0])/2.0, (p1[1] + p2[1])/2.0]
                # clip_geom_center = clip_geom.Centroid()


                side_value = self.side(center,other_point,self.creat_point(tp))
                if side_value > 0:
                    left_area = clip_geom.GetArea()
                    right_area = 0
                else:
                    right_area = clip_geom.GetArea()
                    left_area = 0



            else:

                left_area = 0
                right_area = 0

                for gi in range(clip_geom.GetGeometryCount()):
                    clip_geom_i :Geometry = clip_geom.GetGeometryRef(gi)
                    clip_geom_i_center = clip_geom_i.Centroid()
                    side_value = self.side(center, other_point, clip_geom_i_center)

                    if side_value > 0:
                        left_area += clip_geom_i.GetArea()
                    else:
                        right_area  += clip_geom_i.GetArea()


            print(left_area,right_area)

            last_last_point = last_point
            last_point = next_point
            left_or_right = 1 if left_area > right_area else -1


            # 进入二分法步骤
            if fen:

                alternate:Geometry = ogr.Geometry(ogr.wkbPoint)
                alternate.AddPoint(angle[0][0],angle[0][1])
                # print(angle)
                alternate_side = self.side(center, other_point, alternate)
                print(alternate_side)
                if left_or_right * alternate_side > 0:
                    angle = [angle[0],next_point]
                    next_point = ((next_point[0]+angle[0][0])/2.0,(next_point[1]+angle[0][1])/2.0)

                    # print(abs(next_point[0]-angle[0][0]))
                else:
                    angle = [angle[1],next_point]
                    next_point = ((next_point[0]+angle[1][0])/2.0,(next_point[1]+angle[1][1])/2.0)

                    # print(abs(next_point[0] - angle[1][0]))

            else:
                # print(point_i)
                point_i += 1
                next_point = polygon_line.GetPoint(point_i)

                next_point_side = self.side(center, other_point, self.creat_point(next_point))

                if left_or_right * next_point_side > 0 :
                    pass
                else:
                    next_point = ((last_point[0]+last_last_point[0])/2.0,(last_point[1]+last_last_point[1])/2.0)
                    angle = [last_last_point,last_point]
                    fen = True

        print(calculate_times)

        point_intersection:Geometry = line.Intersection(polygon_line)

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

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

        return self.creat_point(reprecent_point)



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

        line: Geometry = ogr.Geometry(ogr.wkbLineString)
        line.AddPoint(x0, y0)
        x2 = x0 + len_size

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

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

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



    def side(self,A,B,P):


        side = (B.GetX() - A.GetX()) * (P.GetY() - A.GetY()) - (B.GetY() - A.GetY()) * (P.GetX() - A.GetX())
        return side




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

    driver: Driver = ogr.GetDriverByName("ESRI Shapefile")
    ds: DataSource = driver.Open(r"J:\Data\矢量数据\佛山\制图综合\t2.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])
    print(pre)

    result="J:\Data\矢量数据\佛山\制图综合\zhresult.shp"
    print(pre)

    ShapeUtil.create_point(result,"zh",pre)