2022.03.22 备份

@@ -11,6 +11,8 @@ from ..models import Database,DES
 
 from app.util.component.ApiTemplate import ApiTemplate
 from app.util.component.ModelVisitor import ModelVisitor
+from app.util.component.UserAlias import UserAlias
+
 class Api(ApiTemplate):
     api_name = "获取数据库列表"
     def process(self):
@@ -32,6 +34,9 @@ class Api(ApiTemplate):
     
             database = database.limit(page_size).offset(page_index * page_size).all()
             res["data"]["list"] =[]
+
+            user_alias = UserAlias()
+
             for datab in database:
                 # 测试连接
                 try:
@@ -45,6 +50,8 @@ class Api(ApiTemplate):
                 datab_json = ModelVisitor.database_to_json(datab)
                 datab_json["table_count"] = table_count
                 datab_json["available"] = available
+                datab_json["display_name"] = user_alias.get_alias(datab_json["creator"])
+
                 del datab_json["connectstr"]
                 res["data"]["list"].append(datab_json)
 

@@ -4,12 +4,13 @@
 
 
 from ..models import Table,Catalog
-
+from app.modules.auth.models import User
 from sqlalchemy import or_,and_
-
+from sqlalchemy.orm import Query
 from app.util.component.ApiTemplate import ApiTemplate
-
+from app.util.component.UserAlias import UserAlias
 from app.util.component.ModelVisitor import ModelVisitor
+
 class Api(ApiTemplate):
     api_name = "表列表"
     def process(self):
@@ -37,10 +38,11 @@ class Api(ApiTemplate):
     
             sort = int(self.para.get("sort","1"))
             if sort.__eq__(1):
-                tables = Table.query.order_by(Table.update_time.desc(),Table.name)
+                tables:Query = Table.query.order_by(Table.update_time.desc(),Table.name)
             else:
                 tables = Table.query.order_by(Table.name)
 
+
             #过滤spatial_ref_sys
             tables = tables.filter(Table.name.notin_(["spatial_ref_sys"]))
 
@@ -75,8 +77,12 @@ class Api(ApiTemplate):
             tables = tables.limit(page_size).offset(page_index * page_size).all()
             res["data"]["list"]=[]
 
+
+            user_alias = UserAlias()
             for t in tables:
-                res["data"]["list"].append(ModelVisitor.table_to_json(t))
+                table_json =  ModelVisitor.table_to_json(t)
+                table_json["display_name"] = user_alias.get_alias(table_json["creator"])
+                res["data"]["list"].append(table_json)
 
             res["result"] = True
         except Exception as e:

@@ -10,7 +10,8 @@ from ..models import Task
 from app.util.component.ApiTemplate import ApiTemplate
 from app.util.component.ModelVisitor import ModelVisitor
 import json
-
+from app.modules.auth.models import User
+from app.util.component.UserAlias import UserAlias
 
 class Api(ApiTemplate):
     api_name = "任务列表"
@@ -33,7 +34,7 @@ class Api(ApiTemplate):
 
             state = parameter.get("state")
 
-            if creator:
+            if creator and creator !="admin":
                 tasks = tasks.filter_by(creator=creator)
             if start_time and end_time:
                 tasks = tasks.filter(Task.create_time.between(start_time, end_time))
@@ -52,6 +53,8 @@ class Api(ApiTemplate):
             res["data"]["list"] = [ModelVisitor.task_to_json(task) for task in tasks]
 
             try:
+                user_alias = UserAlias()
+
                 for task_info in res["data"]["list"]:
                     if task_info["task_type"]==1:
                         task_info["layers"] = []
@@ -60,6 +63,7 @@ class Api(ApiTemplate):
                             for key in m.get("layer").keys():
 
                                 task_info["layers"].append({"name":key,"new_name":m.get("layer").get(key)})
+                    task_info["display_name"] = user_alias.get_alias(task_info["creator"])
             except:
                 pass
 

@@ -39,9 +39,15 @@ class Api(ApiTemplate):
         try:
 
             this_time = datetime.datetime.now()
-            service_guid = uuid.uuid1().__str__ida()
+            service_guid = uuid.uuid1().__str__()
             map_service_guid = uuid.uuid1().__str__()
 
+            token = self.para.get("token")
+            creator = self.para.get("creator")
+            if creator is None:
+                raise Exception("用户检验失败！")
+
+
             # 逻辑是，先调用引擎接口，引擎说可以，才做持久化
             project_file = self.create_projectfile(self.para)
             para = {"name":self.para.get("name"),
@@ -72,7 +78,7 @@ class Api(ApiTemplate):
                 guid = service_guid,
                 name = self.para.get("name"),
                 title = self.para.get("title"),
-                creator=self.para.get("creator"),
+                creator=creator,
                 state = 1,
                 create_time = this_time,
                 update_time = this_time,
@@ -103,7 +109,9 @@ class Api(ApiTemplate):
             db.session.commit()
 
             res["data"] = service_guid
+            res["url"] = ""
             res["result"] = True
+
         except Exception as e:
             db.session.rollback()
             raise e
@@ -294,6 +302,12 @@ class Api(ApiTemplate):
              "type": "string",
              "description": "表名",
              "required": "true"},
+
+            {"name": "token",
+             "in": "formData",
+             "type": "string",
+             "description": "token",
+             "required": "true"},
         ],
         "responses": {
             200: {

@@ -11,6 +11,7 @@ from sqlalchemy import or_
 import requests
 import re
 from .util.ServiceType import ServiceType
+from app.util.component.UserAlias import UserAlias
 
 class Api(ApiTemplate):
     api_name = "服务列表"
@@ -99,6 +100,12 @@ class Api(ApiTemplate):
             services_json = sorted(services_json,key=lambda x:x["update_time"],reverse=True)
 
             res["data"]["list"] = services_json[page_index*page_size:(page_index+1)*page_size]
+
+            # 显示别名
+            user_alias = UserAlias()
+            for service_json in res["data"]["list"]:
+                service_json["display_name"] = user_alias.get_alias(service_json["creator"])
+
             res["result"] = True
 
         except Exception as e:

@@ -15,6 +15,7 @@ class ParameterUtil:
         parameter = dict()
         parameter.update(request.args.items())
         parameter.update(request.form.items())
+
         try:
             request_json = request.get_json()
             if json:

+# coding=utf-8
+#author:        4N
+#createtime:    2022/3/17
+#email:         nheweijun@sina.com
+
+from app.modules.auth.models import User
+
+class UserAlias:
+    def __init__(self):
+        self.creator_dict = {}
+
+    def get_alias(self,creator):
+        alias = self.creator_dict.get(creator)
+        if alias:
+            return alias
+        else:
+            user = User.query.filter_by(username=creator).one_or_none()
+            if user:
+                self.creator_dict[creator] = user.displayname
+            else:
+                self.creator_dict[creator] = ""
+            return self.creator_dict[creator]

+# coding=utf-8
+#author:        4N
+#createtime:    2022/3/14
+#email:         nheweijun@sina.com
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib.collections import PolyCollection,LineCollection
+import matplotlib as mpl
+from PIL import Image
+import os
+
+
+def draw(geom):
+
+
+
+    rectangle = [(0, 0), (0, 1), (1, 1), (1, 0),(0, 0)]
+    # hexagon = [(0, 0), (0, 1), (1, 2), (2, 1), (2, 0), (1, -1)]
+    l_shape = [(0, 0), (0, 3), (1, 3), (1, 1), (3, 1), (3, 0)]
+    # concave = [(0, 0), (0, 3), (1, 3), (1, 1), (2, 1), (2, 3), (3, 3), (3, 0)]
+
+    for points in [rectangle,l_shape]:
+        xs, ys = zip(*points)
+        # plt.plot(xs, ys, 'o')
+        plt.plot(xs, ys, '-')
+
+    plt.axis('off')
+    png_path = "F:/temp.png"
+    plt.savefig(png_path)
+    png = Image.open(png_path)
+    tif_path = "F:/temp.tif"
+    png.save(tif_path)
+    os.remove(png_path)
+
+
+if __name__ == '__main__':
+        draw()
\ No newline at end of file

+# coding=utf-8
+#author:        4N
+#createtime:    2022/3/15
+#email:         nheweijun@sina.com
+
+from osgeo import ogr
+from osgeo.ogr import *
+import math
+import os
+
+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
+
+
+p1= util.envelop_2_polygon([0,1,0,1])
+p2 = util.envelop_2_polygon([2,3,2,3])
+
+lll = p1.GetGeometryRef(0)
+print(lll.GetPoints())
+
+

+#!/bin/env python
+
+##
+## calculate the concave hull of a set of points
+## see: CONCAVE HULL: A K-NEAREST NEIGHBOURS APPROACH
+##      FOR THE COMPUTATION OF THE REGION OCCUPIED BY A
+##      SET OF POINTS
+##      Adriano Moreira and Maribel Yasmina Santos 2007
+##
+
+import numpy as np
+import scipy.spatial as spt
+import matplotlib.pyplot as plt
+from matplotlib.path import Path
+from test.zonghe import lineintersect as li
+
+def GetFirstPoint(dataset):
+    ''' Returns index of first point, which has the lowest y value '''
+    # todo: what if there is more than one point with lowest y?
+    imin = np.argmin(dataset[:,1])
+    return dataset[imin]
+
+def GetNearestNeighbors(dataset, point, k):
+    ''' Returns indices of k nearest neighbors of point in dataset'''
+    # todo: experiment with leafsize for performance
+    mytree = spt.cKDTree(dataset,leafsize=10)
+    distances, indices = mytree.query(point,k)
+    # todo: something strange here, we get more indices than points in dataset
+    #       so have to do this
+    return dataset[indices[:dataset.shape[0]]]
+
+def SortByAngle(kNearestPoints, currentPoint, prevPoint):
+    ''' Sorts the k nearest points given by angle '''
+    angles = np.zeros(kNearestPoints.shape[0])
+    i = 0
+    for NearestPoint in kNearestPoints:
+        # calculate the angle
+        angle = np.arctan2(NearestPoint[1]-currentPoint[1],
+                NearestPoint[0]-currentPoint[0]) - \
+                np.arctan2(prevPoint[1]-currentPoint[1],
+                prevPoint[0]-currentPoint[0])
+        angle = np.rad2deg(angle)
+        # only positive angles
+        angle = np.mod(angle+360,360)
+        #print NearestPoint[0], NearestPoint[1], angle
+        angles[i] = angle
+        i=i+1
+    return kNearestPoints[np.argsort(angles)]
+
+# def plotPoints(dataset):
+#     plt.plot(dataset[:,0],dataset[:,1],'o',markersize=10,markerfacecolor='0.75',
+#             markeredgewidth=1)
+#     plt.axis('equal')
+#     plt.axis([min(dataset[:,0])-0.5,max(dataset[:,0])+0.5,min(dataset[:,1])-0.5,
+#         max(dataset[:,1])+0.5])
+#     plt.show()
+#
+# def plotPath(dataset, path):
+#     plt.plot(dataset[:,0],dataset[:,1],'o',markersize=10,markerfacecolor='0.65',
+#             markeredgewidth=0)
+#     path = np.asarray(path)
+#     plt.plot(path[:,0],path[:,1],'o',markersize=10,markerfacecolor='0.55',
+#             markeredgewidth=0)
+#     plt.plot(path[:,0],path[:,1],'-',lw=1.4,color='k')
+#     plt.axis('equal')
+#     plt.axis([min(dataset[:,0])-0.5,max(dataset[:,0])+0.5,min(dataset[:,1])-0.5,
+#         max(dataset[:,1])+0.5])
+#     plt.axis('off')
+#     plt.savefig('./doc/figure_1.png', bbox_inches='tight')
+#     plt.show()
+
+def removePoint(dataset, point):
+    delmask = [np.logical_or(dataset[:,0]!=point[0],dataset[:,1]!=point[1])]
+    newdata = dataset[delmask]
+    return newdata
+
+
+def concaveHull(dataset, k):
+    assert k >= 3, 'k has to be greater or equal to 3.'
+    points = dataset
+    # todo: remove duplicate points from dataset
+    # todo: check if dataset consists of only 3 or less points
+    # todo: make sure that enough points for a given k can be found
+
+    firstpoint = GetFirstPoint(points)
+    # init hull as list to easily append stuff
+    hull = []
+    # add first point to hull
+    hull.append(firstpoint)
+    # and remove it from dataset
+    points = removePoint(points,firstpoint)
+    currentPoint = firstpoint
+    # set prevPoint to a Point righ of currentpoint (angle=0)
+    prevPoint = (currentPoint[0]+10, currentPoint[1])
+    step = 2
+
+    while ( (not np.array_equal(firstpoint, currentPoint) or (step==2)) and points.size > 0 ):
+        if ( step == 5 ): # we're far enough to close too early
+            points = np.append(points, [firstpoint], axis=0)
+        kNearestPoints = GetNearestNeighbors(points, currentPoint, k)
+        cPoints = SortByAngle(kNearestPoints, currentPoint, prevPoint)
+        # avoid intersections: select first candidate that does not intersect any
+        # polygon edge
+        its = True
+        i = 0
+        while ( (its==True) and (i<cPoints.shape[0]) ):
+                i=i+1
+                if ( np.array_equal(cPoints[i-1], firstpoint) ):
+                    lastPoint = 1
+                else:
+                    lastPoint = 0
+                j = 2
+                its = False
+                while ( (its==False) and (j<np.shape(hull)[0]-lastPoint) ):
+                    its = li.doLinesIntersect(hull[step-1-1], cPoints[i-1],
+                            hull[step-1-j-1],hull[step-j-1])
+                    j=j+1
+        if ( its==True ):
+            print("all candidates intersect -- restarting with k = ",k+1)
+            return concaveHull(dataset,k+1)
+        prevPoint = currentPoint
+        currentPoint = cPoints[i-1]
+        # add current point to hull
+        hull.append(currentPoint)
+        points = removePoint(points,currentPoint)
+        step = step+1
+    # check if all points are inside the hull
+    p = Path(hull)
+    pContained = p.contains_points(dataset, radius=0.0000000001)
+    if (not pContained.all()):
+        print("not all points of dataset contained in hull -- restarting with k = ",k+1)
+        return concaveHull(dataset, k+1)
+
+    print("finished with k = ",k)
+    return hull
+
+
+
+if __name__ == '__main__':
+
+
+
+    points_E = np.array([[1,1],[2,1],[3,1],[4,1],[5,1],[6,1],[6,2],[6,3],[5,3],[4,3],
+                    [3,3],[3,4],[3,5],[4,5],[5,5],[5,6],[5,7],[4,7],[3,7],[3,8],
+                    [3,9],[4,9],[5,9],[6,9],[6,10],[6,11],[5,11],[4,11],[3,11],[2,11],
+                    [1,11],[1,10],[1,9],[1,8],[1,7],[1,6],[1,5],[1,4],[1,3],[1,2],
+                    [5,2],[4,2],[3,2],[2,2],[2,3],[2,4],[2,5],[2,6],[2,7],[2,8],
+                    [2,9],[2,10],[3,10],[4,10],[5,10],[3,6],[4,6],[5,6],[4.5,7],[3,8.5],
+                    ])
+
+    print(concaveHull(points_E, 3))
+
+

+# coding=utf-8
+#author:        4N
+#createtime:    2022/3/15
+#email:         nheweijun@sina.com
+
+from osgeo import ogr
+from osgeo.ogr import *
+import math
+import os
+import copy
+from test.zonghe.RepresentPoint import RepresentPoint
+
+
+
+class DelauneyNode:
+
+    def __init__(self, rp):
+        self.rp = rp
+        self.children_nodes = list()
+        self.children_delauney_lines = list()
+
+    def get_all_node(self):
+        for node in self.children_nodes:
+            yield node.get_all_node()
+        yield self
+
+    def get_rps(self):
+        for node in self.children_nodes:
+            for rp in node.get_rps():
+                yield rp
+        yield self.rp
+
+    def get_polygons(self):
+        for node in self.children_nodes:
+            for rp in node.get_rps():
+                yield rp.polygon
+        yield self.rp.polygon
+
+    def get_delauney_lines(self):
+        for node in self.children_nodes:
+            for ll in node.get_delauney_lines():
+                yield ll
+        for l in self.children_delauney_lines:
+            yield l
+
+    def add_children(self,node,delauney_line):
+        self.children_nodes.append(node)
+        self.children_delauney_lines.append(delauney_line)
+
+class DelauneyLine:
+
+    def __init__(self,rp1,rp2):
+
+        self.rp1 = rp1
+        self.rp2 = rp2
+        self.distance = rp1.polygon.Distance(rp2.polygon)
+
+        self.line = ogr.Geometry(ogr.wkbLineString)
+        self.line.AddPoint(self.rp1.point.GetX(),self.rp1.point.GetY())
+        self.line.AddPoint(self.rp2.point.GetX(),self.rp2.point.GetY())
+
+
+class DelauneyMinTree:
+
+    def __init__(self):
+        self.delauney_lines_list = list()
+        self.rps = set()
+
+    def append(self,delauney_line):
+        # self.delauney_lines.add(delauney_line)
+
+        self.rps.add(delauney_line.rp1)
+        self.rps.add(delauney_line.rp2)
+
+        self.delauney_lines_list.append(delauney_line)
+
+    def get_rps(self):
+        return self.rps
+
+    def get_other_relate_delauney_lines(self,other_delauney_lines_kv):
+        other_relate_delauney_lines = []
+        for rp in self.rps:
+            other_relate_delauney_lines.extend(other_delauney_lines_kv.get(rp))
+        return other_relate_delauney_lines
+
+    def get_other_rp(self,delauney_line):
+        rp1 = delauney_line.rp1
+        rp2 = delauney_line.rp2
+
+        if not self.rps.__contains__(rp1):
+            return rp1
+        if not self.rps.__contains__(rp2):
+            return rp2
+        return None
+
+
+

+# coding=utf-8
+#author:        4N
+#createtime:    2022/3/15
+#email:         nheweijun@sina.com
+
+from osgeo import ogr
+from osgeo.ogr import *
+import math
+import os
+import copy
+class RepresentPoint:
+
+    def __init__(self,point,polygon):
+        self.point = point
+        self.polygon = polygon

+# coding=utf-8
+#author:        4N
+#createtime:    2022/3/15
+#email:         nheweijun@sina.com
+
+
+from osgeo import ogr
+from osgeo.ogr import *
+import math
+import os
+import copy
+
+class ShapeData:
+
+    driver = ogr.GetDriverByName("ESRI Shapefile")
+
+    def __init__(self,path):
+
+        self.ds: DataSource = self.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()
+
+
+
+
+    @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()
\ No newline at end of file

+# coding=utf-8
+#author:        4N
+#createtime:    2022/3/15
+#email:         nheweijun@sina.com
+
+from osgeo import ogr
+from osgeo.ogr import *
+import math
+import os
+import copy
+
+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
\ No newline at end of file

+# 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"
+
+from test.zonghe.ShapeData import ShapeData
+from test.zonghe.Delauney import DelauneyNode,DelauneyLine,DelauneyMinTree
+from test.zonghe.RepresentPoint import RepresentPoint
+from test.zonghe.Util import Util
+import cv2
+
+class Zonghe:
+
+    data:ShapeData
+    area_threshold = 100.0
+    buffer_threshold = 100.0
+    distance_buffer_threshold = 2.114691025217302e-04
+
+    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 = 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):
+            # print("\"" + center.ExportToWkt() + "\",")
+            return RepresentPoint(center, polygon)
+
+
+        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")
+
+        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)
+
+
+            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
+
+        line = self.division(polygon, center, in_hole, line_length, angle1, angle2, area_thre)
+
+        # 提前return
+        if in_hole:
+
+            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]
+            firt_short_length = short_length/20.0
+            radline,other_point = self.get_radline(center.GetX(),center.GetY(),nearest_point.GetX(),nearest_point.GetY(),firt_short_length)
+            # print("\"" + other_point.ExportToWkt() + "\",")
+
+
+            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)
+
+
+
+            return RepresentPoint(other_point,polygon)
+
+
+        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.ExportToWkt() + "\",")
+
+        rp = RepresentPoint(reprecent_point,polygon)
+
+        return rp
+
+
+
+
+
+    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):
+        a2 = math.pow((x1-x0),2)
+        b2 = math.pow((y1-y0),2)
+        len_size = len_size+ math.sqrt(a2+b2)
+
+        line: Geometry = ogr.Geometry(ogr.wkbLineString)
+        line.AddPoint(x0, y0)
+        tanx = abs(y1-y0) / abs(x1-x0)
+        cosx = 1 / math.sqrt(1 + tanx *tanx)
+        dx = cosx* len_size
+        dy = tanx * dx
+
+        if x1 == x0:
+            x2 = x0
+            if y1 > y0 :
+                y2 = y0 + abs(dy)
+            else:
+                y2 = y0 - abs(dy)
+        else:
+            if x1 < x0:
+                x2 = x0 - dx
+            else:
+                x2 = x0 + dx
+
+            if y1 < y0:
+                y2 = y0 - dy
+            else:
+                y2 = y0 + dy
+
+        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)
+
+
+            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)
+
+        one,other = self.jude(polygon,line,mid_angle,line_length)
+
+        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 create_delauney(self,rps):
+        multipoint: Geometry = ogr.Geometry(ogr.wkbMultiPoint)
+
+        rp_dict = {}
+        line_set = set()
+        delauney_lines = []
+
+        for rp in rps:
+            multipoint.AddGeometry(rp.point)
+
+            rp_dict[rp.point.GetPoint(0).__str__()] = rp
+
+        delauney:Geometry =  multipoint.DelaunayTriangulation()
+
+        for index in range(delauney.GetGeometryCount()):
+            triangle = delauney.GetGeometryRef(index)
+
+            triangle_line:Geometry = ogr.ForceToLineString(triangle)
+
+
+            for pi in range(triangle_line.GetPointCount()-1):
+                p = triangle_line.GetPoint(pi)
+                p_next = triangle_line.GetPoint(pi+1)
+
+                if p[0] < p_next[0]:
+                    key = "{}{}".format(p.__str__(),p_next.__str__())
+                elif p[0] > p_next[0]:
+                    key = "{}{}".format(p_next.__str__(), p.__str__())
+                else:
+                    if p[1] < p_next[1]:
+                        key = "{}{}".format(p.__str__(), p_next.__str__())
+                    else:
+                        key = "{}{}".format(p_next.__str__(), p.__str__())
+
+                if not line_set.__contains__(key):
+                    line_set.add(key)
+                    rp1 = rp_dict.get(p.__str__())
+                    rp2 = rp_dict.get(p_next.__str__())
+
+                    delauney_line = DelauneyLine(rp1,rp2)
+                    delauney_lines.append(delauney_line)
+
+        return delauney_lines
+
+
+
+    def create_min_delauney(self,delauney_lines,rps,polygons):
+        '''
+        使用prim算法计算最小生成树
+        :param delauney_lines:
+        :param polygons:
+        :return:
+        '''
+
+        kv = {}
+        rps_set = set(rps)
+
+        for delauney_line in delauney_lines:
+            if kv.get(delauney_line.rp1):
+
+                kv[delauney_line.rp1].add(delauney_line)
+            else:
+                kv[delauney_line.rp1] = set()
+                kv[delauney_line.rp1].add(delauney_line)
+
+            if kv.get(delauney_line.rp2):
+
+                kv[delauney_line.rp2].add(delauney_line)
+            else:
+
+                kv[delauney_line.rp2] = set()
+                kv[delauney_line.rp2].add(delauney_line)
+
+
+        #初始化树
+        delauney_lines = sorted(delauney_lines,key=lambda x:x.distance)
+
+
+        delauney_min_tree = DelauneyMinTree()
+        delauney_min_tree.append(delauney_lines[0])
+
+        rps_set.remove(delauney_lines[0].rp1)
+        rps_set.remove(delauney_lines[0].rp2)
+        kv[delauney_lines[0].rp1].remove(delauney_lines[0])
+        kv[delauney_lines[0].rp2].remove(delauney_lines[0])
+
+
+        while rps_set.__len__()>0:
+
+            relate_delauney_line = delauney_min_tree.get_other_relate_delauney_lines(kv)
+
+            relate_delauney_line_sort = sorted(relate_delauney_line,key=lambda x:x.distance)
+
+            min_distance_line = relate_delauney_line_sort[0]
+
+            #删除
+            other_rp = delauney_min_tree.get_other_rp(min_distance_line)
+            if other_rp is None:
+                kv[min_distance_line.rp1].remove(min_distance_line)
+                kv[min_distance_line.rp2].remove(min_distance_line)
+                continue
+
+            delauney_min_tree.append(min_distance_line)
+
+            rps_set.remove(other_rp)
+
+            kv[min_distance_line.rp1].remove(min_distance_line)
+            kv[min_distance_line.rp2].remove(min_distance_line)
+
+
+
+        for l in delauney_min_tree.delauney_lines_list:
+            line = l.line
+            # print("\"" + line.ExportToWkt() + "\",")
+
+        # print(delauney_min_tree.delauney_lines_list.__len__())
+
+        return delauney_min_tree
+
+    def cut_delauney_min_tree(self,delauney_min_tree:DelauneyMinTree,threshold):
+
+        trees = []
+        kv = {}
+
+        count = 0
+
+        for delauney_line in delauney_min_tree.delauney_lines_list:
+
+            dn1 = kv.get(delauney_line.rp1)
+            dn2 = kv.get(delauney_line.rp2)
+
+
+            if not dn1:
+                dn1 = DelauneyNode(delauney_line.rp1)
+                kv[delauney_line.rp1] = dn1
+
+                if not dn2:
+
+                    trees.append(dn1)
+
+                    dn2 = DelauneyNode(delauney_line.rp2)
+                    kv[delauney_line.rp2] = dn2
+
+                    if delauney_line.distance < threshold:
+                        dn1.add_children(dn2,delauney_line)
+                        count +=1
+
+                    else:
+                        #断开
+                        trees.append(dn2)
+                else:
+
+                    if delauney_line.distance < threshold:
+                        dn2.add_children(dn1,delauney_line)
+                        count += 1
+                    else:
+                        trees.append(dn1)
+            else:
+
+                dn2 = DelauneyNode(delauney_line.rp2)
+                kv[delauney_line.rp2] = dn2
+
+                if delauney_line.distance < threshold:
+                    dn1.add_children(dn2,delauney_line)
+                    count += 1
+                else:
+                    trees.append(dn2)
+
+
+        return trees
+
+    def merge(self,trees):
+
+        polygons = []
+        for tree in trees:
+
+            merge_polygon_raw = list()
+
+            dls = tree.get_delauney_lines()
+            for dl in dls:
+                p1 = dl.rp1.polygon
+                p2 = dl.rp2.polygon
+
+                merge_polygon_raw.append(p1)
+                merge_polygon_raw.append(p2)
+                line = dl.line
+
+                p1_lines = self.get_polygon_lines(p1)
+                p2_lines = self.get_polygon_lines(p2)
+
+                p1_intersect_lines = [l for l in p1_lines if l.Intersect(line)]
+                p2_intersect_lines = [l for l in p2_lines if l.Intersect(line)]
+
+                if p1_intersect_lines.__len__()>1:
+                    p1_intersect_lines = sorted(p1_intersect_lines,key=lambda x:p2.Distance(x))
+                    p1_intersect_line_rep = p1_intersect_lines[0]
+                elif p1_intersect_lines.__len__()==1:
+                    p1_intersect_line_rep = p1_intersect_lines[0]
+                else:
+                    continue
+
+
+                if p2_intersect_lines.__len__()>1:
+                    p2_intersect_lines = sorted(p2_intersect_lines,key=lambda x:p1.Distance(x))
+                    p2_intersect_line_rep = p2_intersect_lines[0]
+                elif p2_intersect_lines.__len__()==1:
+                    p2_intersect_line_rep = p2_intersect_lines[0]
+                else:
+                    continue
+
+
+                # #两条线没有距离跳过
+                # if p2_intersect_line_rep.Distance(p1_intersect_line_rep) == 0:
+                #     continue
+                # else:
+
+                mid_polygon = self.create_mid_polygon(p1_intersect_line_rep,p2_intersect_line_rep,p1,p2)
+
+                if mid_polygon.GetArea() > 0.0002116361000058077*0.0002116361000058077:
+                    continue
+
+                merge_polygon_raw.append(mid_polygon)
+
+            merge_polygon = None
+            for p in merge_polygon_raw:
+                if merge_polygon is None:
+                    merge_polygon = p
+                else:
+                    merge_polygon = merge_polygon.Union(p)
+
+            polygons.append(merge_polygon)
+
+
+        return polygons
+
+
+    def merge2(self,trees):
+
+        polygons = []
+
+
+
+
+        for tree in trees:
+
+            tree_polygons = tree.get_polygons()
+
+            points = []
+            for p in tree_polygons:
+
+                if p.GetGeometryCount() == 1:
+                    polygon_line: Geometry = ogr.ForceToLineString(p)
+
+                    polygon_line.Con
+                # 有孔Polygon
+                else:
+                    polygon_line: Geometry = ogr.ForceToLineString(p.GetGeometryRef(0))
+
+                points.extend(polygon_line.GetPoints()[:-1])
+
+            multipoint: Geometry = ogr.Geometry(ogr.wkbMultiPoint)
+
+            for p in  points:
+                poi = ogr.Geometry(ogr.wkbPoint)
+                poi.AddPoint(p[0],p[1])
+                multipoint.AddGeometry(poi)
+
+            delauney: Geometry = multipoint.DelaunayTriangulation()
+            merge_polygon = None
+
+            for p in delauney:
+                # if multi_polygon.Contains(p):
+                #     continue
+
+                pline:Geometry = ogr.ForceToLineString(p)
+                pline :list = pline.GetPoints()
+                pline.append(pline[1])
+
+
+                pass_or_not = False
+                for indx in range(3):
+                    p1 = pline[indx]
+                    p2 = pline[indx+1]
+
+                    ppline:Geometry = ogr.Geometry(ogr.wkbLineString)
+                    ppline.AddPoint(p1[0],p1[1])
+                    ppline.AddPoint(p2[0], p2[1])
+
+
+
+                    p3:Geometry = ogr.Geometry(ogr.wkbPoint)
+                    p3.AddPoint(pline[indx+2][0], pline[indx+2][1])
+
+                    if p3.Distance(ppline) > 0.0002116361000058077 or ppline.Length() > 0.0002116361000058077:
+                        pass_or_not = True
+                        break
+
+                if pass_or_not:
+                    continue
+
+                if merge_polygon is None:
+                    merge_polygon = p
+                else:
+                    merge_polygon = merge_polygon.Union(p)
+
+
+
+            polygons.append(merge_polygon)
+
+
+        return polygons
+
+    def merge_circle(self):
+        pass
+
+
+    def create_mid_polygon(self,p1_intersect_line_rep,p2_intersect_line_rep,p1,p2):
+        '''
+        创建衔接polygon
+        :param line:
+        :return:
+        '''
+
+        if p1_intersect_line_rep.Intersect(p2_intersect_line_rep):
+
+            p1_rep_p1 = p1_intersect_line_rep.GetPoints()[0]
+            p1_rep_p2 = p1_intersect_line_rep.GetPoints()[1]
+
+            p2_rep_p1 = p2_intersect_line_rep.GetPoints()[0]
+            p2_rep_p2 = p2_intersect_line_rep.GetPoints()[1]
+
+            polygon = self.create_polygon(p1_rep_p1, p2_rep_p2, p1_rep_p2, p2_rep_p1)
+
+        else:
+
+            p1_rep_p1 = p1_intersect_line_rep.GetPoints()[0]
+            p1_rep_p2 = p1_intersect_line_rep.GetPoints()[1]
+
+            p2_rep_p1 = p2_intersect_line_rep.GetPoints()[0]
+            p2_rep_p2 = p2_intersect_line_rep.GetPoints()[1]
+
+            line1 = ogr.Geometry(ogr.wkbLineString)
+
+
+            line1.AddPoint(p1_rep_p1[0], p1_rep_p1[1])
+            line1.AddPoint(p2_rep_p1[0], p2_rep_p1[1])
+
+            line2 = ogr.Geometry(ogr.wkbLineString)
+            line2.AddPoint(p1_rep_p2[0], p1_rep_p2[1])
+            line2.AddPoint(p2_rep_p2[0], p2_rep_p2[1])
+
+            if line1.Intersect(line2):
+                polygon = self.create_polygon(p1_rep_p1, p1_rep_p2, p2_rep_p1, p2_rep_p2)
+            else:
+
+                polygon = self.create_polygon(p1_rep_p1, p1_rep_p2, p2_rep_p2, p2_rep_p1)
+                if polygon.Intersect(p1) or polygon.Intersect(p2):
+                    polygon = self.create_polygon(p1_rep_p1, p1_rep_p2, p2_rep_p1, p2_rep_p2)
+
+
+        return polygon
+
+    def create_polygon(self,p1,p2,p3,p4):
+
+        ring = ogr.Geometry(ogr.wkbLinearRing)
+        ring.AddPoint(p1[0], p1[1])
+        ring.AddPoint(p2[0], p2[1])
+        ring.AddPoint(p3[0], p3[1])
+        ring.AddPoint(p4[0], p4[1])
+        ring.AddPoint(p1[0], p1[1])
+
+        poly = ogr.Geometry(ogr.wkbPolygon)
+        poly.AddGeometry(ring)
+        return poly
+
+
+    def get_polygon_lines(self,polygon):
+        #无孔Polygon
+        if polygon.GetGeometryCount() == 1:
+            polygon_line : Geometry = ogr.ForceToLineString(polygon)
+        # 有孔Polygon
+        else:
+            polygon_line : Geometry = ogr.ForceToLineString(polygon.GetGeometryRef(0))
+
+        points = polygon_line.GetPoints()
+
+        lines = []
+        for index in range(len(points)-1):
+
+            point_start = points[index]
+            point_end = points[index+1]
+            line = ogr.Geometry(ogr.wkbLineString)
+            line.AddPoint(point_start[0], point_start[1])
+            line.AddPoint(point_end[0], point_end[1])
+
+            lines.append(line)
+
+        return lines
+
+    def process(self):
+        polygons = self.data.get_polygons()
+        rps = []
+        for poly in polygons:
+            rp = self.get_polygon_reprecent_point(poly)
+            rps.append(rp)
+
+        delauney_lines = self.create_delauney(rps)
+
+        min_delauney = self.create_min_delauney(delauney_lines,rps,polygons)
+
+        trees = self.cut_delauney_min_tree(min_delauney,self.distance_buffer_threshold)
+
+
+        polygons = self.merge(trees)
+
+        for polygon in polygons:
+            print("\""+polygon.ExportToWkt() + "\",")
+
+if __name__ == '__main__':
+
+
+    # print(sd.get_polygons()[0])
+    # sd.close()
+
+    create = True
+    create = False
+
+    if create:
+
+        sd = ShapeData(r"J:\Data\制图综合\example.shp")
+
+        ext =  sd.layer.GetExtent()
+        threshold = (ext[1]-ext[0])/20
+        print(threshold)
+
+
+        zh = Zonghe(sd)
+        zh.process()
+    else:
+        wkts = [
+   ]
+
+        result = r"J:\Data\制图综合\zhonghe3.shp"
+        ShapeData.create_shp_fromwkts(result,"zh",wkts)
+

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

+#!/bin/env python
+
+##
+## determine if two line segments intersect
+## see: martin-thoma.com/how-to-check-if-two-line-segments-intersect/
+##
+
+import numpy as np
+
+def doBoundingBoxesIntersect(a, b, c, d):
+    '''
+    Check if bounding boxes do intersect. If one bounding box touches
+    the other, they do intersect.
+    First segment is of points a and b, second of c and d.
+    '''
+    ll1_x = min(a[0],b[0]); ll2_x = min(c[0],d[0])
+    ll1_y = min(a[1],b[1]); ll2_y = min(c[1],d[1])
+    ur1_x = max(a[0],b[0]); ur2_x = max(c[0],d[0])
+    ur1_y = max(a[1],b[1]); ur2_y = max(c[1],d[1])
+
+    return ll1_x <= ur2_x and \
+           ur1_x >= ll2_x and \
+           ll1_y <= ur2_y and \
+           ur1_y >= ll2_y
+
+def isPointOnLine(a,b,c):
+    '''
+    Check if a point is on a line.
+    '''
+    # move to origin
+    aTmp = (0,0)
+    bTmp = (b[0] - a[0], b[1] - a[1])
+    cTmp = (c[0] - a[0], c[1] - a[1])
+    r = np.cross(bTmp, cTmp)
+    return np.abs(r) < 0.0000000001
+
+def isPointRightOfLine(a,b,c):
+    '''
+    Check if a point (c) is right of a line (a-b).
+    If (c) is on the line, it is not right it.
+    '''
+    # move to origin
+    aTmp = (0,0)
+    bTmp = (b[0] - a[0], b[1] - a[1])
+    cTmp = (c[0] - a[0], c[1] - a[1])
+    return np.cross(bTmp, cTmp) < 0
+
+def lineSegmentTouchesOrCrossesLine(a,b,c,d):
+    '''
+    Check if line segment (a-b) touches or crosses
+    line segment (c-d).
+    '''
+    return isPointOnLine(a,b,c) or \
+           isPointOnLine(a,b,d) or \
+          (isPointRightOfLine(a,b,c) ^
+           isPointRightOfLine(a,b,d))
+
+def doLinesIntersect(a,b,c,d):
+    '''
+    Check if line segments (a-b) and (c-d) intersect.
+    '''
+    return doBoundingBoxesIntersect(a,b,c,d) and \
+           lineSegmentTouchesOrCrossesLine(a,b,c,d) and \
+           lineSegmentTouchesOrCrossesLine(c,d,a,b)
+
+
+##############################
+## Tests
+##############################
+
+def test_doBoundingBoxesIntersect():
+    A=(1,1); B=(2,2); C=(3,1); D=(4,2)
+    assert doBoundingBoxesIntersect(A,B,C,D) == False
+    A=(1,2); B=(3,4); C=(2,1); D=(4,3)
+    assert doBoundingBoxesIntersect(A,B,C,D) == True
+
+def test_isPointOnLine():
+    A=(1,1); B=(3,3); C=(2,2)
+    assert isPointOnLine(A,B,C) == True
+    A=(1,1); B=(3,3); C=(3,2)
+    assert isPointOnLine(A,B,C) == False
+
+
+def test_isPointRightOfLine():
+    A=(1,1); B=(3,3); C=(2,2)
+    assert isPointRightOfLine(A,B,C) == False
+    A=(1,1); B=(3,3); C=(3,2)
+    assert isPointRightOfLine(A,B,C) == True
+    A=(1,1); B=(3,3); C=(1,2)
+    assert isPointRightOfLine(A,B,C) == False
+
+# a lot of testcases to be tested with the final function
+
+def tcase(name):
+    if name == 'F1':
+        return (0,0), (7,7), (3,4), (4,5), False
+    elif name == 'F2':
+        return (-4,4), (-2,1), (-2,3), (0,0), False
+    elif name == 'F3':
+        return (0,0), (0,1), (2,2), (2,3), False
+    elif name == 'F4':
+        return (0,0), (0,1), (2,2), (3,2), False
+    elif name == 'F5':
+        return (-1,-1), (2,2), (3,3), (5,5), False
+    elif name == 'F6':
+        return (0,0), (1,1), (2,0), (0.5,2), False
+    elif name == 'F7':
+        return (1,1), (4,1), (2,2), (3,2), False
+    elif name == 'F8':
+        return (0,5), (6,0), (2,1), (2,2), False
+    elif name == 'T1':
+        return (0,-2), (0,2), (-2,0), (2,0), True
+    elif name == 'T2':
+        return (5,5), (0,0), (1,1), (8,2), True
+    elif name == 'T3':
+        return (-1,0), (0,0), (-1,-1), (-1,1), True
+    elif name == 'T4':
+        return (0,2), (2,2), (2,0), (2,4), True
+    elif name == 'T5':
+        return (0,0), (5,5), (1,1), (3,3), True
+    elif name == 'T6':
+        return (0,0), (3,3), (0,0), (3,3), True
+
+cases = ['F1', 'F2', 'F3', 'F4', 'F5', 'F6', 'F7', 'F8',
+         'T1', 'T2', 'T3', 'T4', 'T5', 'T6']
+
+def check_intersection(name):
+    A,B,C,D, result = tcase(name)
+    assert doLinesIntersect(A,B,C,D) == result
+
+def test_doLinesIntersect():
+    for case in cases:
+        yield check_intersection, case
+

+# 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"
+
+from test.zonghe.ShapeData import ShapeData
+from test.zonghe.ConcaveHull import concaveHull
+import numpy as np
+
+class Zonghe:
+
+    data:ShapeData
+    area_threshold = 100.0
+    buffer_threshold = 100.0
+    distance_buffer_threshold = 2.114691025217302e-04
+
+    def __init__(self,data):
+        self.data = data
+
+
+    def creat_point(self,tuple):
+        p: Geometry = ogr.Geometry(ogr.wkbPoint)
+        p.AddPoint(tuple[0], tuple[1])
+        return p
+
+
+
+    def merge2(self,p1,p2):
+
+        polygons = []
+        points = []
+
+        for p in [p1,p2]:
+
+            if p.GetGeometryCount() == 1:
+                polygon_line: Geometry = ogr.ForceToLineString(p)
+            # 有孔Polygon
+            else:
+                polygon_line: Geometry = ogr.ForceToLineString(p.GetGeometryRef(0))
+
+            points.extend(polygon_line.GetPoints()[:-1])
+
+        multipoint: Geometry = ogr.Geometry(ogr.wkbMultiPoint)
+
+
+        for p in  points:
+            poi = ogr.Geometry(ogr.wkbPoint)
+            poi.AddPoint(p[0],p[1])
+            multipoint.AddGeometry(poi)
+
+        delauney: Geometry = multipoint.DelaunayTriangulation()
+        merge_polygon = None
+
+        for p in delauney:
+
+            p_c = p.Centroid()
+
+            # print("\"" + p_c.ExportToWkt() + "\",")
+            if p1.Contains(p) or p2.Contains(p):
+                continue
+            print("\"" + p.ExportToWkt() + "\",")
+
+
+        return polygons
+
+
+    def merge(self,p1,p2):
+
+
+
+        if p1.GetGeometryCount() == 1:
+            polygon_line1: Geometry = ogr.ForceToLineString(p1)
+        # 有孔Polygon
+        else:
+            polygon_line1: Geometry = ogr.ForceToLineString(p1.GetGeometryRef(0))
+
+
+        p1_points =   polygon_line1.GetPoints()[:-1]
+
+        if p2.GetGeometryCount() == 1:
+            polygon_line2: Geometry = ogr.ForceToLineString(p2)
+        # 有孔Polygon
+        else:
+            polygon_line2: Geometry = ogr.ForceToLineString(p2.GetGeometryRef(0))
+
+        p2_points = polygon_line2.GetPoints()[:-1]
+
+
+        multipoint: Geometry = ogr.Geometry(ogr.wkbMultiPoint)
+
+
+
+
+        for p in p1_points:
+
+            p_t = False
+            for pp in p2_points:
+                line = ogr.Geometry(ogr.wkbLineString)
+                line.AddPoint(p[0], p[1])
+                line.AddPoint(pp[0], pp[1])
+
+                ggg:Geometry = line.Intersection(p1)
+
+
+                if ggg.GetGeometryType() in [1,-2147483647,3001]:
+                    p_t = True
+                    break
+            if p_t:
+                pppp = ogr.Geometry(ogr.wkbPoint)
+                pppp.AddPoint(p[0],p[1])
+
+                multipoint.AddGeometry(pppp)
+
+        for p in p2_points:
+
+            p_t = False
+            for pp in p1_points:
+                line = ogr.Geometry(ogr.wkbLineString)
+                line.AddPoint(p[0], p[1])
+                line.AddPoint(pp[0], pp[1])
+
+                ggg:Geometry = line.Intersection(p2)
+
+
+                if ggg.GetGeometryType() in [1,-2147483647,3001]:
+                    p_t = True
+                    break
+            if p_t:
+                pppp = ogr.Geometry(ogr.wkbPoint)
+                pppp.AddPoint(p[0],p[1])
+                multipoint.AddGeometry(pppp)
+
+
+        delauney: Geometry = multipoint.DelaunayTriangulation()
+        merge_polygon = None
+
+        for p in delauney:
+
+            print("\"" + p.ExportToWkt() + "\",")
+
+
+
+        return None
+
+    def merge3(self,polygons):
+
+        multi_poly = ogr.Geometry(ogr.wkbMultiPolygon)
+
+        points = []
+        for poly in polygons:
+            multi_poly.AddGeometry(poly)
+            if poly.GetGeometryCount() == 1:
+                poly_line: Geometry = ogr.ForceToLineString(poly)
+            # 有孔Polygon
+            else:
+                poly_line: Geometry = ogr.ForceToLineString(poly.GetGeometryRef(0))
+
+            points.extend(poly_line.GetPoints()[:-1])
+
+
+        points = np.array([[p[0],p[1]] for p in points])
+
+        merge_p_line = concaveHull(points, 3)
+
+        ring = ogr.Geometry(ogr.wkbLinearRing)
+
+        merge_p = ogr.Geometry(ogr.wkbPolygon)
+
+
+        for l in merge_p_line:
+            ring.AddPoint(l[0], l[1])
+
+        ring.AddPoint(merge_p_line[0][0], merge_p_line[0][1])
+
+        merge_p.AddGeometry(ring)
+
+
+        # for p in polygons:
+        #
+        #     merge_p = merge_p.Union(p)
+        #
+        #
+        #
+        # print(merge_p)
+
+        return [merge_p.ExportToWkt()]
+
+
+
+
+    def merge4(self,p1,p2):
+
+        polygons = []
+        points = []
+
+        mp: Geometry = ogr.Geometry(ogr.wkbMultiPolygon)
+        mp.AddGeometry(p1)
+        mp.AddGeometry(p2)
+
+
+        delauney: Geometry = mp.DelaunayTriangulation(bOnlyEdges=True)
+        merge_polygon = None
+
+        for p in delauney:
+
+            p_c = p.Centroid()
+
+            # print("\"" + p_c.ExportToWkt() + "\",")
+            # if p1.Contains(p) or p2.Contains(p):
+            #     continue
+            print("\"" + p.ExportToWkt() + "\",")
+
+
+        return polygons
+
+    def merge5(self,p1:Geometry,p2):
+
+
+
+        ext = p1.GetEnvelope()
+
+        tor = min(ext[1]-ext[0],ext[3]-ext[2])
+
+        pp = p1.Simplify(tor/20)
+
+        pp2 = p2.Simplify(tor / 20)
+
+        return [pp.ExportToWkt(),pp2.ExportToWkt() ]
+
+
+    def process(self):
+        polygons = self.data.get_polygons()
+
+        return self.merge3(polygons)
+
+
+if __name__ == '__main__':
+
+
+    # print(sd.get_polygons()[0])
+    # sd.close()
+
+    # create = True
+    # create = False
+    #
+    # if create:
+    #
+    #     sd = ShapeData(r"J:\Data\制图综合\me.shp")
+    #
+    #     ext =  sd.layer.GetExtent()
+    #     threshold = (ext[1]-ext[0])/20
+    #     print(threshold)
+    #
+    #
+    #     zh = Zonghe(sd)
+    #     zh.process()
+    # else:
+    #     wkts = ["POLYGON ((113.371086251 22.5101043861001,113.371311338 22.5099257505,113.371611445 22.5098560668001,113.371450698 22.5095291511,113.370775465 22.5096884259,113.371086251 22.5101043861001))"
+    #             ]
+    #     result = r"J:\Data\制图综合\mes.shp"
+    #     ShapeData.create_shp_fromwkts(result,"zh",wkts)
+    #
+
+
+    sd = ShapeData(r"J:\Data\制图综合\concave.shp")
+
+    ext =  sd.layer.GetExtent()
+    threshold = (ext[1]-ext[0])/20
+    print(threshold)
+
+
+    zh = Zonghe(sd)
+    wkts = zh.process()
+
+
+    result = r"J:\Data\制图综合\concave_merge.shp"
+    ShapeData.create_shp_fromwkts(result,"zh",wkts)
+

@@ -8,6 +8,7 @@ from osgeo import ogr
 from osgeo.ogr import *
 import math
 import os
+import copy
 os.environ["CUDA_VISIBLE_DEVICES"] = "-1"
 class util:
     @classmethod
@@ -24,6 +25,39 @@ class util:
         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")
@@ -88,7 +122,7 @@ class ShapeData:
 
     def __init__(self,path):
         driver: Driver = ogr.GetDriverByName("ESRI Shapefile")
-        self.ds: DataSource = driver.Open(path, 1)
+        self.ds: DataSource = driver.Open(path, 0)
         if not self.ds:
             raise Exception("打开数据失败！")
         self.layer: Layer = self.ds.GetLayer(0)
@@ -99,7 +133,7 @@ class ShapeData:
         polygons = []
         for feature in self.layer:
             f:Feature = feature
-            geom : Geometry = f.GetGeometryRef()
+            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:
@@ -114,7 +148,7 @@ class ShapeData:
 class Zonghe:
 
     data:ShapeData
-    area_threshold = 20.0
+    area_threshold = 100.0
     buffer_threshold = 100.0
 
     def __init__(self,data):
@@ -127,12 +161,20 @@ class Zonghe:
 
         polygon_env = polygon.GetEnvelope()
 
-        line_length = max(polygon_env[1]-polygon_env[0],polygon_env[3]-polygon_env[2])
+        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):
@@ -141,7 +183,10 @@ class Zonghe:
         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
 
@@ -160,9 +205,9 @@ class Zonghe:
                     in_hole = True
                     direction = range(0,190,10)
                     in_hole_line = ogr.ForceToLineString(hole_line)
-                    print("in_hole")
-
+                    # print("in_hole")
 
+        angle1,angle2=None,None
 
         for angle in direction:
 
@@ -171,45 +216,41 @@ class Zonghe:
             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, otherside_area = self.jude(polygon, line, angle, line_length)
+            # print(oneside_area, otherside_area)
 
-            oneside_area = 0
-            otherside_area = 0
+            # print("\"" + line.ExportToWkt() + "\",")
 
-            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:
+
+                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)
@@ -232,14 +273,14 @@ class Zonghe:
 
             radline,other_point = self.get_radline(center.GetX(),center.GetY(),nearest_point.GetX(),nearest_point.GetY(),firt_short_length)
 
-            print(other_point)
+            # 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")
+            # print("return hole_rep")
 
             return other_point
 
@@ -264,60 +305,33 @@ class Zonghe:
             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()+"\",")
-
+        print(reprecent_point)
 
-    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)
+        return reprecent_point
 
-        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)
+    def create_delauney(self,points):
+        multipoint: Geometry = ogr.Geometry(ogr.wkbMultiPoint)
+        for point in points:
+            multipoint.AddGeometry(point)
+        delauney =  multipoint.DelaunayTriangulation()
 
-        p6: Geometry = ogr.Geometry(ogr.wkbPoint)
 
-        p6.AddPoint(3, 1)
 
-        mp.AddGeometry(p6)
 
-        de =  mp.DelaunayTriangulation()
 
-        print(mp)
+    # 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],)
 
-        for i in range(de.GetGeometryCount()):
-            d:Geometry = de.GetGeometryRef(i)
 
-            print("\""+d.ExportToWkt()+"\",")
 
     def get_line(self, x0, y0 , angle, r):
 
@@ -356,6 +370,8 @@ class Zonghe:
         cosx = 1 / math.sqrt(1 + tanx *tanx)
         dx = cosx* len_size
 
+
+
         x2 = x1 + dx
 
         if x1 != x0:
@@ -380,17 +396,42 @@ class Zonghe:
         '''
 
         line = self.get_line(center.GetX(),center.GetY(),0,line_length)
-        line_buffer:Geometry = line.Buffer(line_length / self.buffer_threshold)
+
+        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)
 
-        gc = clip_geom.GetGeometryCount()
         oneside_area = 0
         otherside_area = 0
 
+        gc = clip_geom.GetGeometryCount()
         if gc <= 1:
-            return range(0,370,10)
+            oneside_area = 0
+            otherside_area = 0
         else:
-            triangle1, triangle2 = self.get_side_triangle(line, 0)
+            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)
@@ -401,34 +442,35 @@ class Zonghe:
                     oneside_area += clip_geom_i.GetArea()
                 else:
                     otherside_area += clip_geom_i.GetArea()
+        return oneside_area, otherside_area
 
-            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)
+    def division(self,polygon,center,in_hole,line_length,angle1,angle2,area_thre):
 
-            gc = clip_geom.GetGeometryCount()
-            if gc <= 1:
-                return range(360, -10, -10)
-            else:
-                triangle1, triangle2 = self.get_side_triangle(line, 10)
+        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)
 
-                for gi in range(clip_geom.GetGeometryCount()):
-                    clip_geom_i: Geometry = clip_geom.GetGeometryRef(gi)
-                    it = clip_geom_i.Intersection(triangle1)
+        # print("\""+line.ExportToWkt()+"\",")
 
-                    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()
+        one,other = self.jude(polygon,line,mid_angle,line_length)
+        # print(mid_angle)
+        # print(one,other)
 
-                if abs(oneside_area - otherside_area) > abs(oneside_area_contrast - otherside_area_contrast):
-                    return range(0, 370, 10)
+        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 range(360, -10, -10)
+                    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):
         '''
@@ -483,57 +525,63 @@ class Zonghe:
         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)
-            print(reprecent_point)
+
 
 
 
 if __name__ == '__main__':
 
-    # sd = ShapeData(r"J:\Data\矢量数据\佛山\制图综合.shp")
+
     # print(sd.get_polygons()[0])
     # sd.close()
 
     create = True
-    create = False
+    # 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])
+        # 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()
 
 
-        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"
+        wkts = ["POINT (113.211177054031 22.8681300466504)"
+                ]
+
+        result = r"J:\Data\制图综合\left_center.shp"
         ShapeUtil.create_shp_fromwkts(result,"zh",wkts)