2022.03.31 修复注册空表问题

@@ -205,7 +205,7 @@ class EntryDataVacuate:
 
 
             # 复制原图层的属性
             # 复制原图层的属性
             # 去掉fid的属性
             # 去掉fid的属性
-            schema = [sche for sche in layer.schema if not sche.name.__eq__(fid)]
+            schema = [sche for sche in layer.schema if not sche.name.lower().__eq__(fid)]
             pg_layer.CreateFields(schema)
             pg_layer.CreateFields(schema)
 
 
             #创建抽稀过程
             #创建抽稀过程

@@ -38,14 +38,7 @@ class Api(ApiTemplate):
                     append_dict["sr_wkt"] = sr.ExportToWkt()
                     append_dict["sr_wkt"] = sr.ExportToWkt()
                     append_dict["sr_proj4"] = sr.ExportToProj4()
                     append_dict["sr_proj4"] = sr.ExportToProj4()
 
 
-                srid_sql = '''select st_srid("{}") from public."{}" limit 1'''.format(layer.GetGeometryColumn(),layer.GetName())
-                srid_layer:Layer = pg_ds.ExecuteSQL(srid_sql)
-                if srid_layer:
-                    srid_feature : Feature = srid_layer.GetNextFeature()
-                    if srid_feature:
-                        if srid_feature.GetField(0):
-                            append_dict["srid"] = str(srid_feature.GetField(0))
-
+                append_dict["srid"] = str(PGUtil.get_srid(pg_ds,layer.GetName()))
             else:
             else:
                 append_dict["exist"]=0
                 append_dict["exist"]=0
             columns = table.relate_columns.order_by(Columns.name).all()
             columns = table.relate_columns.order_by(Columns.name).all()

@@ -5,7 +5,6 @@
 
 
 from app.models import db
 from app.models import db
 from sqlalchemy import Column, Integer, String, ForeignKey, Text, DateTime, Time,Float,Binary,Sequence
 from sqlalchemy import Column, Integer, String, ForeignKey, Text, DateTime, Time,Float,Binary,Sequence
-from sqlalchemy import Sequence
 from sqlalchemy.orm import relationship
 from sqlalchemy.orm import relationship
 import base64
 import base64
 from pyDes import *
 from pyDes import *

@@ -38,7 +38,7 @@ class Api(ApiTemplate):
             map_service_update = {}
             map_service_update = {}
 
 
             for key in self.para.keys():
             for key in self.para.keys():
-                if key in ["name","title","state","description","overview","catalog_guid"]:
+                if key in ["name","title","description","overview","catalog_guid"]:
                     service_update[key] = self.para.get(key)
                     service_update[key] = self.para.get(key)
                 if key in ["name","tile","project","capabilities"]:
                 if key in ["name","tile","project","capabilities"]:
                     map_service_update[key] = self.para.get(key)
                     map_service_update[key] = self.para.get(key)
@@ -47,7 +47,9 @@ class Api(ApiTemplate):
                     "title":self.para.get("title"),
                     "title":self.para.get("title"),
                     "type":"mapserver",
                     "type":"mapserver",
                     "capabilities":2,
                     "capabilities":2,
-                    "project":self.para.get("project")}
+                    "project":self.para.get("project"),
+                    "state":service.state
+                    }
 
 
             map_service_register_url = "{}/dmap/api/manager/regservice".format(configure.dmap_engine)
             map_service_register_url = "{}/dmap/api/manager/regservice".format(configure.dmap_engine)
             resp: requests.Response  = requests.post(map_service_register_url,data=json.dumps(para),
             resp: requests.Response  = requests.post(map_service_register_url,data=json.dumps(para),

@@ -199,11 +199,16 @@ class Api(ApiTemplate):
 
 
         extent = [float(x) for x in table.extent.split(",")]
         extent = [float(x) for x in table.extent.split(",")]
 
 
+        minx = para.get("minx")
+        miny = para.get("miny")
+        maxx = para.get("maxx")
+        maxy = para.get("maxy")
+
         xml = xml_format.format(name = para.get("name"),
         xml = xml_format.format(name = para.get("name"),
-                                xmin = extent[0],
-                                xmax = extent[1],
-                                ymin = extent[2],
-                                ymax = extent[3],
+                                xmin = minx if minx else extent[0],
+                                xmax = maxx if maxx else extent[1],
+                                ymin = miny if miny else extent[2],
+                                ymax = maxy if maxy else extent[3],
                                 feature_count = table.feature_count,
                                 feature_count = table.feature_count,
                                 database_guid = table.database_guid,
                                 database_guid = table.database_guid,
                                 gemetry_type = gemetry_type,
                                 gemetry_type = gemetry_type,
@@ -302,12 +307,23 @@ class Api(ApiTemplate):
              "type": "string",
              "type": "string",
              "description": "表名",
              "description": "表名",
              "required": "true"},
              "required": "true"},
-
             {"name": "token",
             {"name": "token",
              "in": "formData",
              "in": "formData",
              "type": "string",
              "type": "string",
              "description": "token",
              "description": "token",
              "required": "true"},
              "required": "true"},
+            {"name": "minx",
+             "in": "formData",
+             "type": "string"},
+            {"name": "miny",
+             "in": "formData",
+             "type": "string"},
+            {"name": "maxx",
+             "in": "formData",
+             "type": "string"},
+            {"name": "maxy",
+             "in": "formData",
+             "type": "string"},
         ],
         ],
         "responses": {
         "responses": {
             200: {
             200: {

@@ -3,7 +3,7 @@
 #createtime:    2021/10/26
 #createtime:    2021/10/26
 #email:         nheweijun@sina.com
 #email:         nheweijun@sina.com
 
 
-from sqlalchemy import Column, Integer, String, ForeignKey, Text, DateTime, Time,Float,Binary
+from sqlalchemy import Column, Integer, String, ForeignKey, Text, DateTime
 from sqlalchemy.orm import relationship
 from sqlalchemy.orm import relationship
 from app.models import db
 from app.models import db
 
 

@@ -104,6 +104,7 @@ class Api(ApiTemplate):
             # 显示别名
             # 显示别名
             user_alias = UserAlias()
             user_alias = UserAlias()
             for service_json in res["data"]["list"]:
             for service_json in res["data"]["list"]:
+
                 service_json["display_name"] = user_alias.get_alias(service_json["creator"])
                 service_json["display_name"] = user_alias.get_alias(service_json["creator"])
 
 
             res["result"] = True
             res["result"] = True

@@ -6,6 +6,8 @@
 from app.util.component.ApiTemplate import ApiTemplate
 from app.util.component.ApiTemplate import ApiTemplate
 from .models import Service
 from .models import Service
 from .util.ServiceType import ServiceType
 from .util.ServiceType import ServiceType
+
+
 class Api(ApiTemplate):
 class Api(ApiTemplate):
     api_name = "注册服务"
     api_name = "注册服务"
     def process(self):
     def process(self):

@@ -27,7 +27,6 @@ class Api(ApiTemplate):
                 if service:
                 if service:
                     service.state = state
                     service.state = state
 
 
-
                     operate = "startService" if state == 1 else "stopService"
                     operate = "startService" if state == 1 else "stopService"
                     server_type = "tileserver" if s_type == "电子地图" else "mapserver"
                     server_type = "tileserver" if s_type == "电子地图" else "mapserver"
                     dmap_engine_url = "{}/dmap/api/manager/{}?servicename={}&servername={}".format(configure.dmap_engine,
                     dmap_engine_url = "{}/dmap/api/manager/{}?servicename={}&servername={}".format(configure.dmap_engine,

@@ -35,7 +35,7 @@ class Api(ApiTemplate):
             service_update = {"update_time":this_time}
             service_update = {"update_time":this_time}
             tile_update = {}
             tile_update = {}
             for key in self.para.keys():
             for key in self.para.keys():
-                if key in ["name","title","state","description","overview","catalog_guid"]:
+                if key in ["name","title","description","overview","catalog_guid"]:
                     service_update[key] = self.para.get(key)
                     service_update[key] = self.para.get(key)
                 if key in ["name","tile_type","vendor","crs","datasource","description",
                 if key in ["name","tile_type","vendor","crs","datasource","description",
                            "layer_name","layer_alias","layer_title","layer_style","layer_format",
                            "layer_name","layer_alias","layer_title","layer_style","layer_format",
@@ -61,7 +61,8 @@ class Api(ApiTemplate):
 
 
             para = {"name":self.para.get("name") if self.para.get("name") else se.name,
             para = {"name":self.para.get("name") if self.para.get("name") else se.name,
                     "title":self.para.get("title") if self.para.get("title") else se.title,
                     "title":self.para.get("title") if self.para.get("title") else se.title,
-                    "type":"tileserver","capabilities":1 if tile_type == "WMTS" else 16,"project":project_file}
+                    "type":"tileserver","capabilities":1 if tile_type == "WMTS" else 16,"project":project_file,
+                    "state": service.state}
 
 
 
 
             tile_service_edit_url = "{}/dmap/api/manager/RegService".format(configure.dmap_engine)
             tile_service_edit_url = "{}/dmap/api/manager/RegService".format(configure.dmap_engine)

@@ -7,7 +7,6 @@ from flask import current_app,request
 import traceback
 import traceback
 import json
 import json
 from app.models import db
 from app.models import db
-from app.modules.auth.models import OAuth2Token,User
 import app
 import app
 
 
 class ApiTemplate:
 class ApiTemplate:

@@ -3,6 +3,7 @@
 #createtime:    2021/5/24
 #createtime:    2021/5/24
 #email:         nheweijun@sina.com
 #email:         nheweijun@sina.com
 from osgeo import ogr
 from osgeo import ogr
+from osgeo.ogr import Layer
 from sqlalchemy import create_engine
 from sqlalchemy import create_engine
 from sqlalchemy.orm import sessionmaker,Session
 from sqlalchemy.orm import sessionmaker,Session
 from app.util.component.StructurePrint import StructurePrint
 from app.util.component.StructurePrint import StructurePrint
@@ -126,16 +127,19 @@ class PGUtil:
 
 
     @classmethod
     @classmethod
     def get_srid(cls,pg_ds,table_name):
     def get_srid(cls,pg_ds,table_name):
-        layer = pg_ds.GetLayerByName(table_name)
+        layer:Layer = pg_ds.GetLayerByName(table_name)
         if not layer:
         if not layer:
             return None
             return None
-        srid_sql = '''select st_srid("{}") from public."{}" limit 1'''.format(layer.GetGeometryColumn(), layer.GetName())
+
+        srid_sql ="SELECT srid FROM public.geometry_columns where f_table_name='{}'".format(table_name)
+
         srid_layer = pg_ds.ExecuteSQL(srid_sql)
         srid_layer = pg_ds.ExecuteSQL(srid_sql)
         srid_feature = srid_layer.GetNextFeature()
         srid_feature = srid_layer.GetNextFeature()
         if srid_feature:
         if srid_feature:
             if srid_feature.GetField(0):
             if srid_feature.GetField(0):
                 return int(srid_feature.GetField(0))
                 return int(srid_feature.GetField(0))
         else:
         else:
+
             return None
             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
+from test.zonghe.RepresentPoint import RepresentPoint
+import time
+
+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()
+        self.touch_delauney_lines = set()
+
+    def append(self,delauney_line,other_delauney_lines_kv):
+
+
+        candicates = []
+        if not self.rps.__contains__(delauney_line.rp1):
+            self.rps.add(delauney_line.rp1)
+
+            candicates.extend(other_delauney_lines_kv.get(delauney_line.rp1))
+
+        if not self.rps.__contains__(delauney_line.rp2):
+            self.rps.add(delauney_line.rp2)
+
+            candicates.extend(other_delauney_lines_kv.get(delauney_line.rp2))
+
+        for candicate in candicates:
+            if (candicate.rp1 in self.rps) and (candicate.rp2 in self.rps):
+                #初始化
+                try:
+                    self.touch_delauney_lines.remove(candicate)
+                except:
+                    print(111)
+
+            else:
+                self.touch_delauney_lines.add(candicate)
+
+        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,rps_set):
+
+        other_relate_delauney_lines = []
+
+        for rp in self.rps:
+            candicates = other_delauney_lines_kv.get(rp)
+            candicates = [c for c in candicates if c.rp1 in rps_set or c.rp2 in rps_set]
+            other_relate_delauney_lines.extend(candicates)
+
+        # other_relate_delauney_lines.append()
+
+        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/31
+#email:         nheweijun@sina.com
+
+# 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.MapSynthesize.ShapeData import ShapeData
+from test.MapSynthesize.Delauney import DelauneyNode,DelauneyLine,DelauneyMinTree
+from test.MapSynthesize.RepresentPoint import RepresentPoint
+from test.MapSynthesize.PolygonMerge import PolygonMerge
+import cv2
+import time
+import numpy as np
+
+class MapSynthesize:
+
+    data:ShapeData
+    area_threshold = 10.0
+    buffer_threshold = 100.0
+    distance_buffer_threshold = 1.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
+        line = None
+        oneside_area, otherside_area = 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 oneside_area > 0 and otherside_area >0:
+                if abs(oneside_area-otherside_area) < area_thre:
+                    break
+
+            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
+
+        #第一条线就是均分线
+        if angle1 is None or angle2 is None:
+            pass
+        else:
+            #如果面积差大于阈值，二分法查找
+            if abs(oneside_area - otherside_area) > area_thre:
+                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):
+        '''
+        二分法查找
+        :param polygon:
+        :param center:
+        :param in_hole:
+        :param line_length:
+        :param angle1:
+        :param angle2:
+        :param area_thre:
+        :return:
+        '''
+
+        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 = {}
+        delauney_lines = []
+
+        rps_sub = set()
+        for rp in rps:
+            multipoint.AddGeometry(rp.point)
+            # 多边形代表点相同，需要融合多边形
+            if rp_dict.get(rp.point.GetPoint(0).__str__()):
+                rp_exist = rp_dict.get(rp.point.GetPoint(0).__str__())
+                rp_exist.polygon = rp_exist.polygon.Union(rp.polygon)
+            else:
+                rp_dict[rp.point.GetPoint(0).__str__()] = rp
+                rps_sub.add(rp)
+
+        delauney:Geometry =  multipoint.DelaunayTriangulation(bOnlyEdges=True)
+
+        for index in range(delauney.GetGeometryCount()):
+            triangle_line:Geometry = delauney.GetGeometryRef(index)
+            p = triangle_line.GetPoint(0)
+            p_next = triangle_line.GetPoint(1)
+
+            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,rps_sub
+
+    def create_min_delauney(self,delauney_lines,rps):
+        '''
+        使用prim算法计算最小生成树
+        :param delauney_lines:
+        :param polygons:
+        :return:
+        '''
+
+        kv = dict()
+        print(len(rps))
+
+
+        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)
+
+        print(len(kv))
+
+        print("完成三角边倒排索引")
+
+        #初始化树
+        delauney_lines = sorted(delauney_lines,key=lambda x:x.distance)
+        delauney_min_tree = DelauneyMinTree()
+
+        delauney_min_tree.append(delauney_lines[0],kv)
+
+        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])
+
+        print("完成最小生成树初始化")
+
+        while rps_set.__len__()>0:
+
+            print(rps_set.__len__())
+            t1 = time.time()
+            # relate_delauney_line = delauney_min_tree.get_other_relate_delauney_lines(kv,rps_set)
+
+            relate_delauney_line = delauney_min_tree.touch_delauney_lines
+
+
+            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,kv)
+
+            rps_set.remove(other_rp)
+
+            kv[min_distance_line.rp1].remove(min_distance_line)
+            kv[min_distance_line.rp2].remove(min_distance_line)
+
+            # print(time.time() - t2)
+        print("完成最小生成树")
+
+        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 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() in [0,1]:
+            polygon_line : Geometry = ogr.ForceToLineString(polygon)
+        # 有孔Polygon
+        else:
+            polygon_line : Geometry = ogr.ForceToLineString(polygon.GetGeometryRef(0))
+
+        return polygon_line
+
+    def merge_concave(self, trees, distance_threshold):
+        merge_polygons = []
+        count = 1
+
+        for tree in trees:
+            print(count)
+            count+=1
+            polygons = list(tree.get_polygons())
+
+            merge_raw: Geometry = polygons[0]
+            for p in polygons:
+                merge_raw = merge_raw.Union(p)
+
+            points = []
+
+
+            for i in range(merge_raw.GetGeometryCount()):
+
+                poly = merge_raw.GetGeometryRef(i)
+                poly_line :Geometry = self.get_polygon_lines(poly)
+                if poly_line:
+                    for indx in range(poly_line.GetPointCount() - 1):
+                        poi = poly_line.GetPoint(indx)
+
+                        points.append(poi)
+
+                        poi_next = poly_line.GetPoint(indx + 1)
+
+                        dis = math.sqrt(math.pow(poi[0] - poi_next[0], 2) + math.pow(poi[1] - poi_next[1], 2))
+                        # print(dis)
+                        if dis > distance_threshold:
+                            times = int(dis / distance_threshold)
+                            # print(times)
+                            for time in range(1, times + 1):
+                                x, y = self.get_subpoint(poi[0], poi[1], poi_next[0], poi_next[1], (time * distance_threshold))
+                                points.append([x, y])
+
+            points = np.array([[p[0], p[1]] for p in points])
+
+
+            cp = []
+            for pp in points:
+                poi = ogr.Geometry(ogr.wkbPoint)
+                poi.AddPoint(pp[0], pp[1])
+                cp.append(poi)
+
+            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)
+            merge_p = merge_p.Simplify(self.distance_buffer_threshold/10)
+            merge_polygons.append(merge_p)
+
+        return merge_polygons
+
+    def get_subpoint(self, x0, y0, x1, y1, len_size):
+
+        len_size = len_size
+        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
+
+        return x2, y2
+
+    def process(self):
+        polygons = self.data.get_polygons()
+        rps = []
+
+        for poly in polygons:
+            rp = self.get_polygon_reprecent_point(poly)
+            rps.append(rp)
+        print("完成代表点计算")
+
+        #代表点要融合
+        delauney_lines,rps_sub = self.create_delauney(rps)
+
+        print("完成树生成计算")
+        min_delauney = self.create_min_delauney(delauney_lines,rps_sub)
+        print("完成最小生成树计算")
+        trees = self.cut_delauney_min_tree(min_delauney,self.distance_buffer_threshold)
+        print("完成树裁剪")
+
+        polygon_result = []
+        for index,tree in enumerate(trees):
+
+            polygons = tree.get_polygons()
+            print(index)
+            polygons = list(polygons)
+            print(len(polygons))
+            polygon_result.append(PolygonMerge.merge(polygons,self.distance_buffer_threshold))
+
+        # polygons = self.merge_concave(trees,self.distance_buffer_threshold)
+
+        return [p.ExportToWkt() for p in polygon_result]
+
+if __name__ == '__main__':
+
+
+
+    t1 = time.time()
+    sd = ShapeData(r"J:\Data\制图综合result\zhongshansub.shp")
+
+    mapsynthesize = MapSynthesize(sd)
+    wkts = mapsynthesize.process()
+    result = r"J:\Data\制图综合result\zhongshansub_merge22.shp"
+    ShapeData.create_shp_fromwkts(result,"zh",wkts)
+
+    print(time.time()-t1)

+# coding=utf-8
+#author:        4N
+#createtime:    2022/3/30
+#email:         nheweijun@sina.com
+
+
+from osgeo import ogr,gdal
+from osgeo.ogr import *
+
+from test.zonghe.ShapeData import ShapeData
+import copy
+
+from rtreelib import RTree,Rect
+
+
+
+class RTreeData:
+
+    def __init__(self,triangle,state):
+        self.state = state
+        self.triangle = triangle
+
+    def set_state(self,state):
+        self.state = state
+
+
+def env2rect(env):
+    return Rect(env[0],env[2],env[1],env[3])
+
+
+class PolygonMerge:
+
+    @classmethod
+    def get_delauney(cls,polygon):
+
+        polygon_line = cls.get_polygon_lines(polygon)
+        polygon_line_points:list = polygon_line.GetPoints()
+
+        if len(polygon_line_points)==4:
+            tris = [polygon]
+        else:
+            d = 0
+            for ind in range(len(polygon_line_points) - 1):
+                p = polygon_line_points[ind]
+                p_n = polygon_line_points[ind + 1]
+
+                d += -0.5 * (p_n[1] + p[1]) * (p_n[0] - p[0])
+
+            if d < 0:
+                polygon_line_points.reverse()
+
+            tris = list(cls.get_tri(polygon_line_points))
+
+        return tris
+
+    @classmethod
+    def get_tri(cls,polygon_line_points:list):
+
+        if len(polygon_line_points) == 4 or len(polygon_line_points) == 3:
+            yield cls.create_polygon(polygon_line_points)
+        else:
+            for index,point in enumerate(polygon_line_points):
+                if index == 0:
+                   continue
+                else:
+                    try:
+                        point_front = polygon_line_points[index-1]
+                        point_next = polygon_line_points[index+1]
+                    except Exception as e:
+                        print(e)
+                        kk=1
+                        d=1
+
+                check = (point[0]-point_front[0]) * (point_next[1]-point[1]) - (point_next[0]-point[0])*(point[1]-point_front[1])
+
+                if check > 0:
+                    triangle = cls.create_polygon([point_front,point,point_next,point_front])
+                    yield triangle
+
+                    polygon_line_points.pop(index)
+                    for tri in cls.get_tri(polygon_line_points):
+                        yield tri
+                    break
+
+    @classmethod
+    def merge(cls,polygons:list,distance_buffer_threshold):
+
+        rtree = RTree()
+
+        polygon:Geometry = polygons[0]
+        # multipolygon = ogr.Geometry(ogr.wkbMultiPolygon)
+        # for poly in polygons:
+        #     multipolygon.AddGeometry(poly)
+
+        for poly in polygons:
+            polygon = polygon.Union(poly)
+
+
+
+        delauney: Geometry = polygon.DelaunayTriangulation()
+
+        for index in range(delauney.GetGeometryCount()):
+            de:Geometry = copy.deepcopy(delauney.GetGeometryRef(index))
+            de_extent = de.GetEnvelope()
+            rtree.insert(RTreeData(de,1), env2rect(de_extent))
+
+        bians: list = cls.get_bian(polygon)
+
+
+        for bian in bians:
+
+            bian_ext = bian.GetEnvelope()
+
+            query_result = rtree.query(env2rect(bian_ext))
+
+            intersection_data = []
+
+
+            for entry in query_result:
+
+                # 已被抛弃的不要
+                if entry.data.state == 0:
+                    continue
+                intersection_geom: Geometry = bian.Intersection(entry.data.triangle)
+
+                if not intersection_geom.IsEmpty():
+                    if not intersection_geom.Equals(bian):
+                        if intersection_geom.GetGeometryType() in [2, 5, -2147483643, -2147483646, 3002, 3005]:
+                            intersection_data.append(entry.data)
+
+
+            for data in intersection_data:
+                data.set_state(0)
+
+            if intersection_data:
+
+                inter_tri_merge:Geometry = intersection_data[0].triangle
+
+                for ind in range(1,len(intersection_data)):
+                    inter_tri_merge = inter_tri_merge.Union(intersection_data[ind].triangle)
+
+
+                inter_tri_merge_line:Geometry = cls.get_polygon_lines(inter_tri_merge)
+
+                inter_tri_merge_line_points:list = inter_tri_merge_line.GetPoints()
+
+                bian_ps = bian.GetPoints()
+
+                if len(inter_tri_merge_line_points[0]) == 3:
+                    try:
+                        first_index = inter_tri_merge_line_points.index(bian_ps[0])
+                        second_index = inter_tri_merge_line_points.index(bian_ps[1])
+                    except:
+                        rtree.insert(RTreeData(inter_tri_merge, 1), env2rect(inter_tri_merge.GetEnvelope()))
+                        continue
+                else:
+                    try:
+                        first_index = inter_tri_merge_line_points.index((bian_ps[0][0],bian_ps[0][1]))
+                        second_index = inter_tri_merge_line_points.index((bian_ps[1][0],bian_ps[1][1]))
+                    except:
+                        rtree.insert(RTreeData(inter_tri_merge, 1), env2rect(inter_tri_merge.GetEnvelope()))
+                        continue
+
+
+
+                small_index = min(first_index,second_index)
+                big_index = max(first_index,second_index)
+
+                small_points:list = inter_tri_merge_line_points[:small_index+1]
+                small_points.extend(inter_tri_merge_line_points[big_index:])
+                small_polygon:Geometry = cls.create_polygon(small_points)
+
+                delauney = cls.get_delauney(small_polygon)
+                for de in delauney:
+                    de_extent = de.GetEnvelope()
+                    rtree.insert(RTreeData(de, 1), env2rect(de_extent))
+
+                big_points = inter_tri_merge_line_points[small_index:big_index+1]
+                big_points.append(inter_tri_merge_line_points[small_index])
+
+                big_polygon:Geometry = cls.create_polygon(big_points)
+
+                delauney = cls.get_delauney(big_polygon)
+                for de in delauney:
+                    de_extent = de.GetEnvelope()
+                    rtree.insert(RTreeData(de, 1), env2rect(de_extent))
+
+        for entry in rtree.get_leaf_entries():
+            triangle = entry.data.triangle
+            if entry.data.state==0:
+                continue
+            if polygon.Contains(triangle):
+                pass
+            else:
+                tri_center,tri_center2 = cls.get_center(triangle)
+                distance = tri_center.Distance(polygon)
+                distance2 = tri_center2.Distance(polygon)
+
+                if distance > distance_buffer_threshold :
+                    continue
+                # if distance2 > distance_buffer_threshold:
+                #     continue
+
+                polygon = polygon.Union(triangle)
+                # polygon.AddGeometry(triangle)
+
+        result :Geometry = polygon
+
+        return result
+
+
+    @classmethod
+    def get_center(cls,tri:Geometry):
+
+        tri_line: Geometry = cls.get_polygon_lines(tri)
+        tri_line_points = tri_line.GetPoints()
+
+        lines = []
+
+        for index in range(len(tri_line_points)-1):
+            first = tri_line_points[index]
+            second = tri_line_points[index+1]
+            line: Geometry = ogr.Geometry(ogr.wkbLineString)
+            line.AddPoint(first[0], first[1])
+            line.AddPoint(second[0], second[1])
+            lines.append(line)
+
+        lines = sorted(lines,key=lambda l:l.Length())
+
+        return lines[-1].Centroid(),lines[-2].Centroid()
+
+    @classmethod
+    def get_polygon_lines(cls,polygon):
+        #无孔Polygon
+        if polygon.GetGeometryCount() in [0,1]:
+            polygon_line : Geometry = ogr.ForceToLineString(polygon)
+        # 有孔Polygon
+        else:
+            polygon_line : Geometry = ogr.ForceToLineString(polygon.GetGeometryRef(0))
+
+        return polygon_line
+
+    @classmethod
+    def get_bian(cls,polygon):
+
+        bians = []
+
+        polygon_lines = []
+
+        if polygon.GetGeometryCount() in [0, 1]:
+            polygon_line: Geometry = ogr.ForceToLineString(polygon)
+            polygon_lines.append(polygon_line)
+        # 有孔Polygon
+        else:
+            for index in range(polygon.GetGeometryCount()):
+                polygon_line: Geometry = ogr.ForceToLineString(polygon.GetGeometryRef(index))
+                polygon_lines.append(polygon_line)
+
+        for polygon_line in polygon_lines:
+            for index in range(polygon_line.GetPointCount()-1):
+                p = polygon_line.GetPoint(index)
+                p_next = polygon_line.GetPoint(index+1)
+
+                line: Geometry = ogr.Geometry(ogr.wkbLineString)
+
+                if len(p) ==2:
+                    line.AddPoint(p[0], p[1],0)
+                    line.AddPoint(p_next[0], p_next[1],0)
+                else:
+                    line.AddPoint(p[0], p[1])
+                    line.AddPoint(p_next[0], p_next[1])
+                bians.append(line)
+
+        return bians
+
+    @classmethod
+    def create_polygon(cls,ps):
+
+        ring = ogr.Geometry(ogr.wkbLinearRing)
+        for p in ps:
+            ring.AddPoint(p[0], p[1])
+        poly = ogr.Geometry(ogr.wkbPolygon)
+        poly.AddGeometry(ring)
+        return poly
+
+
+if __name__ == '__main__':
+    import time
+    t1 = time.time()
+
+
+    # sd = ShapeData(r"J:\Data\制图综合result\t4.shp")
+    # mp = ogr.Geometry(ogr.wkbMultiPolygon)
+    # polygons = sd.get_polygons()
+    #
+    # wkts = PolygonMerge.merge(polygons)
+    #
+    # result = r"J:\Data\制图综合result\t4_merge.shp"
+    #
+    # ShapeData.create_shp_fromwkts(result,"zh",wkts)
+    #
+    # print(time.time()-t1)
+
+
+    p = PolygonMerge.create_polygon([(113.37029204400005, 22.526576828100076, 0.0), (113.36972403800007, 22.526369320699985, 0.0), (113.36974547200009, 22.526379207800005, 0.0), (113.37025989200004, 22.52656695200011, 0.0), (113.37029204400005, 22.526576828100076, 0.0)])
+
+    pl:Geometry = ogr.ForceToLineString(p)
+
+    ps:list = pl.GetPoints()
+    ps.reverse()
+
+    PolygonMerge.get_tri(ps)
+
+    # wkts =[p.ExportToWkt()]
+    #
+    # result = r"J:\Data\制图综合result\error.shp"
+    #
+    # ShapeData.create_shp_fromwkts(result,"zh",wkts)
+    #
+    # print(time.time()-t1)

+# coding=utf-8
+#author:        4N
+#createtime:    2022/3/15
+#email:         nheweijun@sina.com
+
+import uuid
+class RepresentPoint:
+
+    def __init__(self,point,polygon):
+        self.point = point
+        self.polygon = polygon
+        self.uuid = uuid.uuid1().__str__()

+# 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, 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 = copy.deepcopy(f.GetGeometryRef())
+            if geom is None:
+                continue
+            if geom.GetGeometryType() in [3,-2147483645,3001]:
+                polygons.append(ogr.ForceToPolygon(geom))
+            if geom.GetGeometryType() in [6 ,-2147483642 ,3006]:
+                for i in range(geom.GetGeometryCount()):
+                    polygons.append(ogr.ForceToPolygon(geom.GetGeometryRef(i)))
+
+        # p_cs = []
+        # for p in polygons:
+        #     polygon_line: Geometry = self.get_polygon_lines(p)
+        #     print(polygon_line)
+        #     p_cs.append(self.create_polygon(polygon_line.GetPoints()))
+        #
+        # return p_cs
+        return polygons
+
+
+    def get_polygon_lines(self,polygon):
+        #无孔Polygon
+        if polygon.GetGeometryCount() in [0,1]:
+            polygon_line : Geometry = ogr.ForceToLineString(polygon)
+        # 有孔Polygon
+        else:
+            polygon_line : Geometry = ogr.ForceToLineString(polygon.GetGeometryRef(0))
+
+        return polygon_line
+
+
+    def create_polygon(self,ps):
+
+        ring = ogr.Geometry(ogr.wkbLinearRing)
+        for p in ps:
+            ring.AddPoint(p[0], p[1],0)
+        poly = ogr.Geometry(ogr.wkbPolygon)
+        poly.AddGeometry(ring)
+        return poly
+
+
+
+    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()
+
+
+if __name__ == '__main__':
+    sd = ShapeData(r"J:\Data\制图综合result\test.shp")
+    # polygons = sd.get_polygons()
+    # wkts = []
+    # for p in polygons:
+    #     # print(p.ExportToWkt())
+    #     wkts.append(p.ExportToWkt())
+    # result = r"J:\Data\制图综合result\zhongshan_ronghe.shp"
+    # ShapeData.create_shp_fromwkts(result,"zh",wkts)
+
+    layer = sd.layer
+
+    fn = "PG: user=postgres password=chinadci host=172.26.60.101 port=5432 dbname=ceshi "
+    driver = ogr.GetDriverByName("PostgreSQL")
+    if driver is None:
+        raise Exception("打开PostgreSQL驱动失败，可能是当前GDAL未支持PostgreSQL驱动！")
+    ds:DataSource = driver.Open(fn, 1)
+
+
+    pg_layer = ds.CreateLayer("test", layer.GetSpatialRef(), layer.GetGeomType(), ["overwrite=yes"])
+    schema =  layer.schema
+    pg_layer.CreateFields(schema)
+    for feature in layer:
+        out_feature: Feature = copy.copy(feature)
+        out_feature.SetFID(out_feature.GetFID() + 1)
+        pg_layer.CreateFeature(out_feature)
+
+
+    ds.Destroy()

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

@@ -10,7 +10,6 @@
 
 
 import numpy as np
 import numpy as np
 import scipy.spatial as spt
 import scipy.spatial as spt
-import matplotlib.pyplot as plt
 from matplotlib.path import Path
 from matplotlib.path import Path
 from test.zonghe import lineintersect as li
 from test.zonghe import lineintersect as li
 
 
@@ -47,27 +46,6 @@ def SortByAngle(kNearestPoints, currentPoint, prevPoint):
         i=i+1
         i=i+1
     return kNearestPoints[np.argsort(angles)]
     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):
 def removePoint(dataset, point):
     delmask = [np.logical_or(dataset[:,0]!=point[0],dataset[:,1]!=point[1])]
     delmask = [np.logical_or(dataset[:,0]!=point[0],dataset[:,1]!=point[1])]
@@ -89,9 +67,13 @@ def concaveHull(dataset, k):
     hull.append(firstpoint)
     hull.append(firstpoint)
     # and remove it from dataset
     # and remove it from dataset
     points = removePoint(points,firstpoint)
     points = removePoint(points,firstpoint)
+
     currentPoint = firstpoint
     currentPoint = firstpoint
     # set prevPoint to a Point righ of currentpoint (angle=0)
     # set prevPoint to a Point righ of currentpoint (angle=0)
+
     prevPoint = (currentPoint[0]+10, currentPoint[1])
     prevPoint = (currentPoint[0]+10, currentPoint[1])
+
+
     step = 2
     step = 2
 
 
     while ( (not np.array_equal(firstpoint, currentPoint) or (step==2)) and points.size > 0 ):
     while ( (not np.array_equal(firstpoint, currentPoint) or (step==2)) and points.size > 0 ):
@@ -116,7 +98,7 @@ def concaveHull(dataset, k):
                             hull[step-1-j-1],hull[step-j-1])
                             hull[step-1-j-1],hull[step-j-1])
                     j=j+1
                     j=j+1
         if ( its==True ):
         if ( its==True ):
-            print("all candidates intersect -- restarting with k = ",k+1)
+            # print("all candidates intersect -- restarting with k = ",k+1)
             return concaveHull(dataset,k+1)
             return concaveHull(dataset,k+1)
         prevPoint = currentPoint
         prevPoint = currentPoint
         currentPoint = cPoints[i-1]
         currentPoint = cPoints[i-1]
@@ -127,11 +109,12 @@ def concaveHull(dataset, k):
     # check if all points are inside the hull
     # check if all points are inside the hull
     p = Path(hull)
     p = Path(hull)
     pContained = p.contains_points(dataset, radius=0.0000000001)
     pContained = p.contains_points(dataset, radius=0.0000000001)
+
     if (not pContained.all()):
     if (not pContained.all()):
-        print("not all points of dataset contained in hull -- restarting with k = ",k+1)
+        # print("not all points of dataset contained in hull -- restarting with k = ",k+1)
         return concaveHull(dataset, k+1)
         return concaveHull(dataset, k+1)
 
 
-    print("finished with k = ",k)
+    # print("finished with k = ",k)
     return hull
     return hull
 
 
 
 

@@ -9,8 +9,7 @@ import math
 import os
 import os
 import copy
 import copy
 from test.zonghe.RepresentPoint import RepresentPoint
 from test.zonghe.RepresentPoint import RepresentPoint
-
-
+import time
 
 
 class DelauneyNode:
 class DelauneyNode:
 
 
@@ -65,22 +64,52 @@ class DelauneyMinTree:
     def __init__(self):
     def __init__(self):
         self.delauney_lines_list = list()
         self.delauney_lines_list = list()
         self.rps = set()
         self.rps = set()
+        self.touch_delauney_lines = set()
+
+    def append(self,delauney_line,other_delauney_lines_kv):
+
+
+        candicates = []
+        if not self.rps.__contains__(delauney_line.rp1):
+            self.rps.add(delauney_line.rp1)
+
+            candicates.extend(other_delauney_lines_kv.get(delauney_line.rp1))
+
+        if not self.rps.__contains__(delauney_line.rp2):
+            self.rps.add(delauney_line.rp2)
+
+            candicates.extend(other_delauney_lines_kv.get(delauney_line.rp2))
 
 
-    def append(self,delauney_line):
-        # self.delauney_lines.add(delauney_line)
+        for candicate in candicates:
+            if (candicate.rp1 in self.rps) and (candicate.rp2 in self.rps):
+                #初始化
+                try:
+                    self.touch_delauney_lines.remove(candicate)
+                except:
+                    print(111)
 
 
-        self.rps.add(delauney_line.rp1)
-        self.rps.add(delauney_line.rp2)
+            else:
+                self.touch_delauney_lines.add(candicate)
 
 
         self.delauney_lines_list.append(delauney_line)
         self.delauney_lines_list.append(delauney_line)
 
 
+
+
     def get_rps(self):
     def get_rps(self):
         return self.rps
         return self.rps
 
 
-    def get_other_relate_delauney_lines(self,other_delauney_lines_kv):
+
+    def get_other_relate_delauney_lines(self,other_delauney_lines_kv,rps_set):
+
         other_relate_delauney_lines = []
         other_relate_delauney_lines = []
+
         for rp in self.rps:
         for rp in self.rps:
-            other_relate_delauney_lines.extend(other_delauney_lines_kv.get(rp))
+            candicates = other_delauney_lines_kv.get(rp)
+            candicates = [c for c in candicates if c.rp1 in rps_set or c.rp2 in rps_set]
+            other_relate_delauney_lines.extend(candicates)
+
+        # other_relate_delauney_lines.append()
+
         return other_relate_delauney_lines
         return other_relate_delauney_lines
 
 
     def get_other_rp(self,delauney_line):
     def get_other_rp(self,delauney_line):

+# coding=utf-8
+#author:        4N
+#createtime:    2022/3/30
+#email:         nheweijun@sina.com
+
+import math
+
+from osgeo import ogr,gdal
+from osgeo.ogr import *
+
+from test.zonghe.ShapeData import ShapeData
+import copy
+
+from scipy import spatial
+from rtreelib import RTree,Rect
+
+
+
+class RTreeData:
+
+    def __init__(self,triangle,state):
+        self.state = state
+        self.triangle = triangle
+
+    def set_state(self,state):
+        self.state = state
+
+
+def env2rect(env):
+    return Rect(env[0],env[2],env[1],env[3])
+
+
+class PolygonDivision:
+
+    @classmethod
+    def get_delauney(cls,polygon):
+
+        polygon_line = cls.get_polygon_lines(polygon)
+        polygon_line_points:list = polygon_line.GetPoints()
+
+        if len(polygon_line_points)==4:
+            tris = [polygon]
+        else:
+            d = 0
+            for ind in range(len(polygon_line_points) - 1):
+                p = polygon_line_points[ind]
+                p_n = polygon_line_points[ind + 1]
+
+                d += -0.5 * (p_n[1] + p[1]) * (p_n[0] - p[0])
+
+            if d < 0:
+                polygon_line_points.reverse()
+
+            tris = list(cls.get_tri(polygon_line_points))
+
+        return tris
+
+    @classmethod
+    def get_tri(cls,polygon_line_points:list):
+
+        if len(polygon_line_points) == 4:
+            yield cls.create_polygon(polygon_line_points)
+        else:
+            for index,point in enumerate(polygon_line_points):
+                if index == 0:
+                   continue
+                # elif index == len(polygon_line_points)-1 :
+                #     pass
+                else:
+                    point_front = polygon_line_points[index-1]
+                    point_next = polygon_line_points[index+1]
+
+                check = (point[0]-point_front[0]) * (point_next[1]-point[1]) - (point_next[0]-point[0])*(point[1]-point_front[1])
+
+                if check > 0:
+                    triangle = cls.create_polygon([point_front,point,point_next,point_front])
+                    yield triangle
+                    polygon_line_points.pop(index)
+                    for tri in cls.get_tri(polygon_line_points):
+                        yield tri
+                    break
+
+    @classmethod
+    def triangle_division(cls,polygon:Geometry):
+
+        rtree = RTree()
+        polygon = polygon.UnionCascaded()
+        delauney: Geometry = polygon.DelaunayTriangulation()
+
+        for index in range(delauney.GetGeometryCount()):
+            de:Geometry = copy.deepcopy(delauney.GetGeometryRef(index))
+            de_extent = de.GetEnvelope()
+            rtree.insert(RTreeData(de,1), env2rect(de_extent))
+
+
+        # return [entry.data.triangle.ExportToWkt() for entry in rtree.get_leaf_entries()]
+        bians: list = cls.get_bian(polygon)
+
+
+        for bian in bians:
+
+            bian_ext = bian.GetEnvelope()
+
+            query_result = rtree.query(env2rect(bian_ext))
+
+            intersection_data = []
+
+
+            for entry in query_result:
+
+                # 已被抛弃的不要
+                if entry.data.state == 0:
+                    continue
+                intersection_geom: Geometry = bian.Intersection(entry.data.triangle)
+
+                if not intersection_geom.IsEmpty():
+                    if not intersection_geom.Equals(bian):
+                        if intersection_geom.GetGeometryType() in [2, 5, -2147483643, -2147483646, 3002, 3005]:
+                            intersection_data.append(entry.data)
+
+
+            for data in intersection_data:
+                data.set_state(0)
+
+            if intersection_data:
+
+                inter_tri_merge:Geometry = intersection_data[0].triangle
+
+                for ind in range(1,len(intersection_data)):
+                    inter_tri_merge = inter_tri_merge.Union(intersection_data[ind].triangle)
+
+
+
+                inter_tri_merge_line:Geometry = cls.get_polygon_lines(inter_tri_merge)
+
+                inter_tri_merge_line_points:list = inter_tri_merge_line.GetPoints()
+
+                bian_ps = bian.GetPoints()
+
+                if len(inter_tri_merge_line_points[0]) == 3:
+                    try:
+                        first_index = inter_tri_merge_line_points.index(bian_ps[0])
+                        second_index = inter_tri_merge_line_points.index(bian_ps[1])
+                    except:
+                        print(11)
+                        print(bian)
+                        print(inter_tri_merge)
+
+                        rtree.insert(RTreeData(inter_tri_merge, 1), env2rect(inter_tri_merge.GetEnvelope()))
+
+                        continue
+                else:
+                    try:
+                        first_index = inter_tri_merge_line_points.index((bian_ps[0][0],bian_ps[0][1]))
+                        second_index = inter_tri_merge_line_points.index((bian_ps[1][0],bian_ps[1][1]))
+                    except:
+                        print(11)
+                        print(bian)
+                        print(inter_tri_merge)
+
+                        rtree.insert(RTreeData(inter_tri_merge, 1), env2rect(inter_tri_merge.GetEnvelope()))
+                        continue
+
+
+
+                small_index = min(first_index,second_index)
+                big_index = max(first_index,second_index)
+
+                small_points:list = inter_tri_merge_line_points[:small_index+1]
+                small_points.extend(inter_tri_merge_line_points[big_index:])
+                small_polygon:Geometry = cls.create_polygon(small_points)
+
+
+
+                # delauney: Geometry = small_polygon.DelaunayTriangulation()
+                #
+                # for index in range(delauney.GetGeometryCount()):
+                #     de: Geometry = copy.deepcopy(delauney.GetGeometryRef(index))
+                #
+                #     if small_polygon.Contains(de):
+                #         de_extent = de.GetEnvelope()
+                #         rtree.insert(RTreeData(de, 1), env2rect(de_extent))
+
+                delauney = cls.get_delauney(small_polygon)
+                for de in delauney:
+                    de_extent = de.GetEnvelope()
+                    rtree.insert(RTreeData(de, 1), env2rect(de_extent))
+
+
+                big_points = inter_tri_merge_line_points[small_index:big_index+1]
+                big_points.append(inter_tri_merge_line_points[small_index])
+
+                big_polygon:Geometry = cls.create_polygon(big_points)
+
+                # delauney: Geometry = big_polygon.DelaunayTriangulation()
+                #
+                # for index in range(delauney.GetGeometryCount()):
+                #     de: Geometry = copy.deepcopy(delauney.GetGeometryRef(index))
+                #
+                #     if big_polygon.Contains(de):
+                #
+                #         de_extent = de.GetEnvelope()
+                #         rtree.insert(RTreeData(de, 1), env2rect(de_extent))
+
+                delauney = cls.get_delauney(big_polygon)
+                for de in delauney:
+                    de_extent = de.GetEnvelope()
+                    rtree.insert(RTreeData(de, 1), env2rect(de_extent))
+
+        tris = []
+
+        for entry in rtree.get_leaf_entries():
+            if entry.data.state==0:
+                continue
+            if polygon.Contains(entry.data.triangle):
+                pass
+            else:
+                tris.append(entry.data.triangle)
+
+        for tri in tris:
+
+            tri_center,long_line = cls.get_center(tri)
+            dd = tri_center.Distance(polygon)
+            if dd > 1.114691025217302e-04:
+                continue
+
+            polygon.AddGeometry(tri)
+        result :Geometry = polygon.UnionCascaded()
+        return [result.ExportToWkt()]
+
+
+    @classmethod
+    def get_center(cls,tri:Geometry):
+        tri_line: Geometry = ogr.ForceToLineString(tri)
+        tri_line_points = tri_line.GetPoints()
+
+        lines = []
+
+        for index in range(len(tri_line_points)-1):
+            first = tri_line_points[index]
+            second = tri_line_points[index+1]
+            line: Geometry = ogr.Geometry(ogr.wkbLineString)
+            line.AddPoint(first[0], first[1])
+            line.AddPoint(second[0], second[1])
+            lines.append(line)
+
+        lines = sorted(lines,key=lambda l:l.Length())
+
+        return lines[-1].Centroid(),lines[-1]
+
+    @classmethod
+    def get_polygon_lines(cls,polygon):
+        #无孔Polygon
+        if polygon.GetGeometryCount() in [0,1]:
+            polygon_line : Geometry = ogr.ForceToLineString(polygon)
+        # 有孔Polygon
+        else:
+            polygon_line : Geometry = ogr.ForceToLineString(polygon.GetGeometryRef(0))
+
+        return polygon_line
+
+    @classmethod
+    def get_bian(cls,polygon):
+
+        bians = []
+
+        polygon_lines = []
+
+        if polygon.GetGeometryCount() in [0, 1]:
+            polygon_line: Geometry = ogr.ForceToLineString(polygon)
+            polygon_lines.append(polygon_line)
+        # 有孔Polygon
+        else:
+            for index in range(polygon.GetGeometryCount()):
+                polygon_line: Geometry = ogr.ForceToLineString(polygon.GetGeometryRef(index))
+                polygon_lines.append(polygon_line)
+
+        for polygon_line in polygon_lines:
+            for index in range(polygon_line.GetPointCount()-1):
+                p = polygon_line.GetPoint(index)
+                p_next = polygon_line.GetPoint(index+1)
+
+                line: Geometry = ogr.Geometry(ogr.wkbLineString)
+
+                if len(p) ==2:
+                    line.AddPoint(p[0], p[1],0)
+                    line.AddPoint(p_next[0], p_next[1],0)
+                else:
+                    line.AddPoint(p[0], p[1])
+                    line.AddPoint(p_next[0], p_next[1])
+                bians.append(line)
+
+        return bians
+
+    @classmethod
+    def create_polygon(cls,ps):
+
+        ring = ogr.Geometry(ogr.wkbLinearRing)
+        for p in ps:
+            ring.AddPoint(p[0], p[1])
+        poly = ogr.Geometry(ogr.wkbPolygon)
+        poly.AddGeometry(ring)
+        return poly
+
+
+if __name__ == '__main__':
+    import time
+    t1 = time.time()
+    gdal.UseExceptions()
+
+    sd = ShapeData(r"J:\Data\制图综合result\t4.shp")
+    mp = ogr.Geometry(ogr.wkbMultiPolygon)
+    polygons = sd.get_polygons()
+    for p in polygons:
+        mp.AddGeometry(p)
+
+    wkts = PolygonDivision.triangle_division(mp)
+
+    result = r"J:\Data\制图综合result\t4_merge.shp"
+
+    ShapeData.create_shp_fromwkts(result,"zh",wkts)
+
+    print(time.time()-t1)
+
+

@@ -3,13 +3,10 @@
 #createtime:    2022/3/15
 #createtime:    2022/3/15
 #email:         nheweijun@sina.com
 #email:         nheweijun@sina.com
 
 
-from osgeo import ogr
-from osgeo.ogr import *
-import math
-import os
-import copy
+import uuid
 class RepresentPoint:
 class RepresentPoint:
 
 
     def __init__(self,point,polygon):
     def __init__(self,point,polygon):
         self.point = point
         self.point = point
         self.polygon = polygon
         self.polygon = polygon
+        self.uuid = uuid.uuid1().__str__()

@@ -16,7 +16,7 @@ class ShapeData:
 
 
     def __init__(self,path):
     def __init__(self,path):
 
 
-        self.ds: DataSource = self.driver.Open(path, 0)
+        self.ds: DataSource = self.driver.Open(path, 1)
         if not self.ds:
         if not self.ds:
             raise Exception("打开数据失败！")
             raise Exception("打开数据失败！")
         self.layer: Layer = self.ds.GetLayer(0)
         self.layer: Layer = self.ds.GetLayer(0)
@@ -28,13 +28,46 @@ class ShapeData:
         for feature in self.layer:
         for feature in self.layer:
             f:Feature = feature
             f:Feature = feature
             geom : Geometry = copy.deepcopy(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:
+            if geom is None:
+                continue
+            if geom.GetGeometryType() in [3,-2147483645,3001]:
+                polygons.append(ogr.ForceToPolygon(geom))
+            if geom.GetGeometryType() in [6 ,-2147483642 ,3006]:
                 for i in range(geom.GetGeometryCount()):
                 for i in range(geom.GetGeometryCount()):
-                    polygons.append(geom.GetGeometryRef(i))
+                    polygons.append(ogr.ForceToPolygon(geom.GetGeometryRef(i)))
+
+        # p_cs = []
+        # for p in polygons:
+        #     polygon_line: Geometry = self.get_polygon_lines(p)
+        #     print(polygon_line)
+        #     p_cs.append(self.create_polygon(polygon_line.GetPoints()))
+        #
+        # return p_cs
         return polygons
         return polygons
 
 
+
+    def get_polygon_lines(self,polygon):
+        #无孔Polygon
+        if polygon.GetGeometryCount() in [0,1]:
+            polygon_line : Geometry = ogr.ForceToLineString(polygon)
+        # 有孔Polygon
+        else:
+            polygon_line : Geometry = ogr.ForceToLineString(polygon.GetGeometryRef(0))
+
+        return polygon_line
+
+
+    def create_polygon(self,ps):
+
+        ring = ogr.Geometry(ogr.wkbLinearRing)
+        for p in ps:
+            ring.AddPoint(p[0], p[1],0)
+        poly = ogr.Geometry(ogr.wkbPolygon)
+        poly.AddGeometry(ring)
+        return poly
+
+
+
     def close(self):
     def close(self):
         self.ds.Destroy()
         self.ds.Destroy()
 
 
@@ -94,4 +127,35 @@ class ShapeData:
             feat_new = ogr.Feature(layer.GetLayerDefn())
             feat_new = ogr.Feature(layer.GetLayerDefn())
             feat_new.SetGeometry(geom)
             feat_new.SetGeometry(geom)
             layer.CreateFeature(feat_new)
             layer.CreateFeature(feat_new)
-        data_source.Destroy()
+        data_source.Destroy()
+
+
+if __name__ == '__main__':
+    sd = ShapeData(r"J:\Data\制图综合result\test.shp")
+    # polygons = sd.get_polygons()
+    # wkts = []
+    # for p in polygons:
+    #     # print(p.ExportToWkt())
+    #     wkts.append(p.ExportToWkt())
+    # result = r"J:\Data\制图综合result\zhongshan_ronghe.shp"
+    # ShapeData.create_shp_fromwkts(result,"zh",wkts)
+
+    layer = sd.layer
+
+    fn = "PG: user=postgres password=chinadci host=172.26.60.101 port=5432 dbname=ceshi "
+    driver = ogr.GetDriverByName("PostgreSQL")
+    if driver is None:
+        raise Exception("打开PostgreSQL驱动失败，可能是当前GDAL未支持PostgreSQL驱动！")
+    ds:DataSource = driver.Open(fn, 1)
+
+
+    pg_layer = ds.CreateLayer("test", layer.GetSpatialRef(), layer.GetGeomType(), ["overwrite=yes"])
+    schema =  layer.schema
+    pg_layer.CreateFields(schema)
+    for feature in layer:
+        out_feature: Feature = copy.copy(feature)
+        out_feature.SetFID(out_feature.GetFID() + 1)
+        pg_layer.CreateFeature(out_feature)
+
+
+    ds.Destroy()

@@ -10,19 +10,21 @@ import math
 import os
 import os
 import copy
 import copy
 os.environ["CUDA_VISIBLE_DEVICES"] = "-1"
 os.environ["CUDA_VISIBLE_DEVICES"] = "-1"
-
+from test.zonghe.ConcaveHull import concaveHull
 from test.zonghe.ShapeData import ShapeData
 from test.zonghe.ShapeData import ShapeData
 from test.zonghe.Delauney import DelauneyNode,DelauneyLine,DelauneyMinTree
 from test.zonghe.Delauney import DelauneyNode,DelauneyLine,DelauneyMinTree
 from test.zonghe.RepresentPoint import RepresentPoint
 from test.zonghe.RepresentPoint import RepresentPoint
 from test.zonghe.Util import Util
 from test.zonghe.Util import Util
 import cv2
 import cv2
+import time
+import numpy as np
 
 
 class Zonghe:
 class Zonghe:
 
 
     data:ShapeData
     data:ShapeData
-    area_threshold = 100.0
+    area_threshold = 10.0
     buffer_threshold = 100.0
     buffer_threshold = 100.0
-    distance_buffer_threshold = 2.114691025217302e-04
+    distance_buffer_threshold = 1.114691025217302e-04
 
 
     def __init__(self,data):
     def __init__(self,data):
         self.data = data
         self.data = data
@@ -37,11 +39,8 @@ class Zonghe:
 
 
         line_length = math.sqrt(math.pow(polygon_env[1] - polygon_env[0],2) + math.pow(polygon_env[3] - polygon_env[2],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])
         short_length = min(polygon_env[1]-polygon_env[0],polygon_env[3]-polygon_env[2])
 
 
-
-
         center:Geometry = polygon.Centroid()
         center:Geometry = polygon.Centroid()
         # print(center)
         # print(center)
 
 
@@ -80,7 +79,8 @@ class Zonghe:
                     # print("in_hole")
                     # print("in_hole")
 
 
         angle1,angle2=None,None
         angle1,angle2=None,None
-
+        line = None
+        oneside_area, otherside_area = None,None
         for angle in direction:
         for angle in direction:
 
 
             if in_hole:
             if in_hole:
@@ -91,6 +91,10 @@ class Zonghe:
             oneside_area, otherside_area = self.jude(polygon, line, angle, line_length)
             oneside_area, otherside_area = self.jude(polygon, line, angle, line_length)
 
 
 
 
+            if oneside_area > 0 and otherside_area >0:
+                if abs(oneside_area-otherside_area) < area_thre:
+                    break
+
             if last_state is None:
             if last_state is None:
 
 
                 angle2 = angle
                 angle2 = angle
@@ -115,7 +119,14 @@ class Zonghe:
                 else:
                 else:
                     last_state = now_state
                     last_state = now_state
 
 
-        line = self.division(polygon, center, in_hole, line_length, angle1, angle2, area_thre)
+        #第一条线就是均分线
+        if angle1 is None or angle2 is None:
+            pass
+        else:
+            #如果面积差大于阈值，二分法查找
+            if abs(oneside_area - otherside_area) > area_thre:
+                line = self.division(polygon, center, in_hole, line_length, angle1, angle2, area_thre)
+
 
 
         # 提前return
         # 提前return
         if in_hole:
         if in_hole:
@@ -156,6 +167,7 @@ class Zonghe:
             ylist.append(pi.GetY())
             ylist.append(pi.GetY())
             plist.append(pi)
             plist.append(pi)
 
 
+
         reprecent_point = [(min(xlist)+max(xlist))/2.0,(min(ylist)+max(ylist))/2.0]
         reprecent_point = [(min(xlist)+max(xlist))/2.0,(min(ylist)+max(ylist))/2.0]
 
 
         reprecent_point = self.creat_point(reprecent_point)
         reprecent_point = self.creat_point(reprecent_point)
@@ -298,6 +310,8 @@ class Zonghe:
     def division(self,polygon,center,in_hole,line_length,angle1,angle2,area_thre):
     def division(self,polygon,center,in_hole,line_length,angle1,angle2,area_thre):
 
 
         mid_angle = (angle1 + angle2)/2.0
         mid_angle = (angle1 + angle2)/2.0
+
+
         if in_hole:
         if in_hole:
             line = self.get_doubleline(center.GetX(), center.GetY(), mid_angle, line_length)
             line = self.get_doubleline(center.GetX(), center.GetY(), mid_angle, line_length)
         else:
         else:
@@ -373,55 +387,41 @@ class Zonghe:
         return p
         return p
 
 
 
 
-
-
     def create_delauney(self,rps):
     def create_delauney(self,rps):
         multipoint: Geometry = ogr.Geometry(ogr.wkbMultiPoint)
         multipoint: Geometry = ogr.Geometry(ogr.wkbMultiPoint)
 
 
         rp_dict = {}
         rp_dict = {}
-        line_set = set()
         delauney_lines = []
         delauney_lines = []
 
 
+        rps_sub = set()
         for rp in rps:
         for rp in rps:
             multipoint.AddGeometry(rp.point)
             multipoint.AddGeometry(rp.point)
+            # 多边形代表点相同，需要融合多边形
+            if rp_dict.get(rp.point.GetPoint(0).__str__()):
+                rp_exist = rp_dict.get(rp.point.GetPoint(0).__str__())
+                rp_exist.polygon = rp_exist.polygon.Union(rp.polygon)
+            else:
+                rp_dict[rp.point.GetPoint(0).__str__()] = rp
+                rps_sub.add(rp)
 
 
-            rp_dict[rp.point.GetPoint(0).__str__()] = rp
-
-        delauney:Geometry =  multipoint.DelaunayTriangulation()
+        delauney:Geometry =  multipoint.DelaunayTriangulation(bOnlyEdges=True)
 
 
         for index in range(delauney.GetGeometryCount()):
         for index in range(delauney.GetGeometryCount()):
-            triangle = delauney.GetGeometryRef(index)
-
-            triangle_line:Geometry = ogr.ForceToLineString(triangle)
-
+            triangle_line:Geometry = delauney.GetGeometryRef(index)
+            p = triangle_line.GetPoint(0)
+            p_next = triangle_line.GetPoint(1)
 
 
-            for pi in range(triangle_line.GetPointCount()-1):
-                p = triangle_line.GetPoint(pi)
-                p_next = triangle_line.GetPoint(pi+1)
+            rp1 = rp_dict.get(p.__str__())
+            rp2 = rp_dict.get(p_next.__str__())
 
 
-                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)
 
 
-                    delauney_line = DelauneyLine(rp1,rp2)
-                    delauney_lines.append(delauney_line)
+        return delauney_lines,rps_sub
 
 
-        return delauney_lines
 
 
 
 
-
-    def create_min_delauney(self,delauney_lines,rps,polygons):
+    def create_min_delauney(self,delauney_lines,rps):
         '''
         '''
         使用prim算法计算最小生成树
         使用prim算法计算最小生成树
         :param delauney_lines:
         :param delauney_lines:
@@ -429,68 +429,72 @@ class Zonghe:
         :return:
         :return:
         '''
         '''
 
 
-        kv = {}
+        kv = dict()
+        print(len(rps))
+
+
         rps_set = set(rps)
         rps_set = set(rps)
 
 
         for delauney_line in delauney_lines:
         for delauney_line in delauney_lines:
             if kv.get(delauney_line.rp1):
             if kv.get(delauney_line.rp1):
-
                 kv[delauney_line.rp1].add(delauney_line)
                 kv[delauney_line.rp1].add(delauney_line)
             else:
             else:
                 kv[delauney_line.rp1] = set()
                 kv[delauney_line.rp1] = set()
                 kv[delauney_line.rp1].add(delauney_line)
                 kv[delauney_line.rp1].add(delauney_line)
 
 
             if kv.get(delauney_line.rp2):
             if kv.get(delauney_line.rp2):
-
                 kv[delauney_line.rp2].add(delauney_line)
                 kv[delauney_line.rp2].add(delauney_line)
             else:
             else:
-
                 kv[delauney_line.rp2] = set()
                 kv[delauney_line.rp2] = set()
                 kv[delauney_line.rp2].add(delauney_line)
                 kv[delauney_line.rp2].add(delauney_line)
 
 
+        print(len(kv))
+
+        print("完成三角边倒排索引")
 
 
         #初始化树
         #初始化树
         delauney_lines = sorted(delauney_lines,key=lambda x:x.distance)
         delauney_lines = sorted(delauney_lines,key=lambda x:x.distance)
-
-
         delauney_min_tree = DelauneyMinTree()
         delauney_min_tree = DelauneyMinTree()
-        delauney_min_tree.append(delauney_lines[0])
+
+        delauney_min_tree.append(delauney_lines[0],kv)
 
 
         rps_set.remove(delauney_lines[0].rp1)
         rps_set.remove(delauney_lines[0].rp1)
         rps_set.remove(delauney_lines[0].rp2)
         rps_set.remove(delauney_lines[0].rp2)
+
         kv[delauney_lines[0].rp1].remove(delauney_lines[0])
         kv[delauney_lines[0].rp1].remove(delauney_lines[0])
         kv[delauney_lines[0].rp2].remove(delauney_lines[0])
         kv[delauney_lines[0].rp2].remove(delauney_lines[0])
 
 
+        print("完成最小生成树初始化")
 
 
         while rps_set.__len__()>0:
         while rps_set.__len__()>0:
 
 
-            relate_delauney_line = delauney_min_tree.get_other_relate_delauney_lines(kv)
+            print(rps_set.__len__())
+            t1 = time.time()
+            # relate_delauney_line = delauney_min_tree.get_other_relate_delauney_lines(kv,rps_set)
+
+            relate_delauney_line = delauney_min_tree.touch_delauney_lines
+
 
 
             relate_delauney_line_sort = sorted(relate_delauney_line,key=lambda x:x.distance)
             relate_delauney_line_sort = sorted(relate_delauney_line,key=lambda x:x.distance)
 
 
             min_distance_line = relate_delauney_line_sort[0]
             min_distance_line = relate_delauney_line_sort[0]
-
             #删除
             #删除
             other_rp = delauney_min_tree.get_other_rp(min_distance_line)
             other_rp = delauney_min_tree.get_other_rp(min_distance_line)
+
             if other_rp is None:
             if other_rp is None:
                 kv[min_distance_line.rp1].remove(min_distance_line)
                 kv[min_distance_line.rp1].remove(min_distance_line)
                 kv[min_distance_line.rp2].remove(min_distance_line)
                 kv[min_distance_line.rp2].remove(min_distance_line)
                 continue
                 continue
 
 
-            delauney_min_tree.append(min_distance_line)
+            delauney_min_tree.append(min_distance_line,kv)
 
 
             rps_set.remove(other_rp)
             rps_set.remove(other_rp)
 
 
             kv[min_distance_line.rp1].remove(min_distance_line)
             kv[min_distance_line.rp1].remove(min_distance_line)
             kv[min_distance_line.rp2].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__())
+            # print(time.time() - t2)
+        print("完成最小生成树")
 
 
         return delauney_min_tree
         return delauney_min_tree
 
 
@@ -546,277 +550,154 @@ class Zonghe:
 
 
         return trees
         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 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])
 
 
-    def merge2(self,trees):
+        poly = ogr.Geometry(ogr.wkbPolygon)
+        poly.AddGeometry(ring)
+        return poly
 
 
-        polygons = []
 
 
+    def get_polygon_lines(self,polygon):
+        #无孔Polygon
+        if polygon.GetGeometryCount() in [0,1]:
+            polygon_line : Geometry = ogr.ForceToLineString(polygon)
+        # 有孔Polygon
+        else:
+            polygon_line : Geometry = ogr.ForceToLineString(polygon.GetGeometryRef(0))
 
 
+        return polygon_line
 
 
+    def merge_concave(self, trees, distance_threshold):
+        merge_polygons = []
+        count = 1
 
 
         for tree in trees:
         for tree in trees:
+            print(count)
+            count+=1
+            polygons = list(tree.get_polygons())
 
 
-            tree_polygons = tree.get_polygons()
+            merge_raw: Geometry = polygons[0]
+            for p in polygons:
+                merge_raw = merge_raw.Union(p)
 
 
             points = []
             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)
 
 
 
 
+            for i in range(merge_raw.GetGeometryCount()):
 
 
-            polygons.append(merge_polygon)
+                poly = merge_raw.GetGeometryRef(i)
+                poly_line :Geometry = self.get_polygon_lines(poly)
+                if poly_line:
+                    for indx in range(poly_line.GetPointCount() - 1):
+                        poi = poly_line.GetPoint(indx)
 
 
+                        points.append(poi)
 
 
-        return polygons
+                        poi_next = poly_line.GetPoint(indx + 1)
 
 
-    def merge_circle(self):
-        pass
+                        dis = math.sqrt(math.pow(poi[0] - poi_next[0], 2) + math.pow(poi[1] - poi_next[1], 2))
+                        # print(dis)
+                        if dis > distance_threshold:
+                            times = int(dis / distance_threshold)
+                            # print(times)
+                            for time in range(1, times + 1):
+                                x, y = self.get_subpoint(poi[0], poi[1], poi_next[0], poi_next[1], (time * distance_threshold))
+                                points.append([x, y])
 
 
+            points = np.array([[p[0], p[1]] for p in points])
 
 
-    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):
+            cp = []
+            for pp in points:
+                poi = ogr.Geometry(ogr.wkbPoint)
+                poi.AddPoint(pp[0], pp[1])
+                cp.append(poi)
 
 
-            p1_rep_p1 = p1_intersect_line_rep.GetPoints()[0]
-            p1_rep_p2 = p1_intersect_line_rep.GetPoints()[1]
+            merge_p_line = concaveHull(points, 3)
 
 
-            p2_rep_p1 = p2_intersect_line_rep.GetPoints()[0]
-            p2_rep_p2 = p2_intersect_line_rep.GetPoints()[1]
+            ring = ogr.Geometry(ogr.wkbLinearRing)
 
 
-            polygon = self.create_polygon(p1_rep_p1, p2_rep_p2, p1_rep_p2, p2_rep_p1)
+            merge_p = ogr.Geometry(ogr.wkbPolygon)
 
 
-        else:
+            for l in merge_p_line:
+                ring.AddPoint(l[0], l[1])
 
 
-            p1_rep_p1 = p1_intersect_line_rep.GetPoints()[0]
-            p1_rep_p2 = p1_intersect_line_rep.GetPoints()[1]
+            ring.AddPoint(merge_p_line[0][0], merge_p_line[0][1])
 
 
-            p2_rep_p1 = p2_intersect_line_rep.GetPoints()[0]
-            p2_rep_p2 = p2_intersect_line_rep.GetPoints()[1]
+            merge_p.AddGeometry(ring)
+            merge_p = merge_p.Simplify(self.distance_buffer_threshold/10)
+            merge_polygons.append(merge_p)
 
 
-            line1 = ogr.Geometry(ogr.wkbLineString)
+        return merge_polygons
 
 
 
 
-            line1.AddPoint(p1_rep_p1[0], p1_rep_p1[1])
-            line1.AddPoint(p2_rep_p1[0], p2_rep_p1[1])
+    def get_subpoint(self, x0, y0, x1, y1, len_size):
 
 
-            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])
+        len_size = len_size
+        tanx = abs(y1 - y0) / abs(x1 - x0)
+        cosx = 1 / math.sqrt(1 + tanx * tanx)
+        dx = cosx * len_size
+        dy = tanx * dx
 
 
-            if line1.Intersect(line2):
-                polygon = self.create_polygon(p1_rep_p1, p1_rep_p2, p2_rep_p1, p2_rep_p2)
+        if x1 == x0:
+            x2 = x0
+            if y1 > y0:
+                y2 = y0 + abs(dy)
             else:
             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
+                y2 = y0 - abs(dy)
         else:
         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])
+            if x1 < x0:
+                x2 = x0 - dx
+            else:
+                x2 = x0 + dx
 
 
-            lines.append(line)
+            if y1 < y0:
+                y2 = y0 - dy
+            else:
+                y2 = y0 + dy
 
 
-        return lines
+        return x2, y2
 
 
     def process(self):
     def process(self):
         polygons = self.data.get_polygons()
         polygons = self.data.get_polygons()
         rps = []
         rps = []
+
         for poly in polygons:
         for poly in polygons:
             rp = self.get_polygon_reprecent_point(poly)
             rp = self.get_polygon_reprecent_point(poly)
             rps.append(rp)
             rps.append(rp)
+        print("完成代表点计算")
 
 
-        delauney_lines = self.create_delauney(rps)
-
-        min_delauney = self.create_min_delauney(delauney_lines,rps,polygons)
+        #代表点要融合
+        delauney_lines,rps_sub = self.create_delauney(rps)
 
 
+        print("完成树生成计算")
+        min_delauney = self.create_min_delauney(delauney_lines,rps_sub)
+        print("完成最小生成树计算")
         trees = self.cut_delauney_min_tree(min_delauney,self.distance_buffer_threshold)
         trees = self.cut_delauney_min_tree(min_delauney,self.distance_buffer_threshold)
+        print("完成树裁剪")
 
 
+        polygons = self.merge_concave(trees,self.distance_buffer_threshold)
 
 
-        polygons = self.merge(trees)
-
-        for polygon in polygons:
-            print("\""+polygon.ExportToWkt() + "\",")
+        return [p.ExportToWkt() for p in polygons]
 
 
 if __name__ == '__main__':
 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 = [
-   ]
+    t1 = time.time()
+    sd = ShapeData(r"J:\Data\制图综合result\zhongshansub.shp")
 
 
-        result = r"J:\Data\制图综合\zhonghe3.shp"
-        ShapeData.create_shp_fromwkts(result,"zh",wkts)
+    zh = Zonghe(sd)
+    wkts = zh.process()
+    result = r"J:\Data\制图综合result\zhongshan_single.shp"
+    ShapeData.create_shp_fromwkts(result,"zh",wkts)
 
 
+    print(time.time()-t1)

+# 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.ConcaveHull import concaveHull
+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
+import time
+import numpy as np
+
+class Zonghe:
+
+    data:ShapeData
+    area_threshold = 10.0
+    buffer_threshold = 100.0
+    distance_buffer_threshold = 1.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
+        line = None
+        oneside_area, otherside_area = 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 oneside_area > 0 and otherside_area >0:
+                if abs(oneside_area-otherside_area) < area_thre:
+                    break
+
+            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
+
+        #第一条线就是均分线
+        if angle1 is None or angle2 is None:
+            pass
+        else:
+            #如果面积差大于阈值，二分法查找
+            if abs(oneside_area - otherside_area) > area_thre:
+                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 = {}
+        delauney_lines = []
+
+        rps_sub = set()
+        for rp in rps:
+            multipoint.AddGeometry(rp.point)
+            # 多边形代表点相同，需要融合多边形
+            if rp_dict.get(rp.point.GetPoint(0).__str__()):
+                rp_exist = rp_dict.get(rp.point.GetPoint(0).__str__())
+                rp_exist.polygon = rp_exist.polygon.Union(rp.polygon)
+            else:
+                rp_dict[rp.point.GetPoint(0).__str__()] = rp
+                rps_sub.add(rp)
+
+        delauney:Geometry =  multipoint.DelaunayTriangulation(bOnlyEdges=True)
+
+        for index in range(delauney.GetGeometryCount()):
+            triangle_line:Geometry = delauney.GetGeometryRef(index)
+            p = triangle_line.GetPoint(0)
+            p_next = triangle_line.GetPoint(1)
+
+            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,rps_sub
+
+    def create_min_delauney(self,delauney_lines,rps):
+        '''
+        使用prim算法计算最小生成树
+        :param delauney_lines:
+        :param polygons:
+        :return:
+        '''
+
+        kv = dict()
+        # print(len(rps))
+
+
+        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)
+
+
+
+        print("完成三角边倒排索引")
+
+        #初始化树
+        delauney_lines = sorted(delauney_lines,key=lambda x:x.distance)
+        delauney_min_tree = DelauneyMinTree()
+
+        delauney_min_tree.append(delauney_lines[0],kv)
+
+        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])
+
+        print("完成最小生成树初始化")
+
+        while rps_set.__len__()>0:
+
+            relate_delauney_line = delauney_min_tree.touch_delauney_lines
+
+
+            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,kv)
+
+            rps_set.remove(other_rp)
+
+            kv[min_distance_line.rp1].remove(min_distance_line)
+            kv[min_distance_line.rp2].remove(min_distance_line)
+
+
+        print("完成最小生成树")
+
+        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 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() in [0,1]:
+            polygon_line : Geometry = ogr.ForceToLineString(polygon)
+        # 有孔Polygon
+        else:
+            polygon_line : Geometry = ogr.ForceToLineString(polygon.GetGeometryRef(0))
+
+        return polygon_line
+
+    def merge_concave(self, trees, distance_threshold):
+        merge_polygons = []
+        count = 1
+
+        for tree in trees:
+            print(count)
+            count+=1
+            polygons = list(tree.get_polygons())
+
+            merge_raw: Geometry = polygons[0]
+            for p in polygons:
+                merge_raw = merge_raw.Union(p)
+
+            points = []
+
+
+            for i in range(merge_raw.GetGeometryCount()):
+
+                poly = merge_raw.GetGeometryRef(i)
+                poly_line :Geometry = self.get_polygon_lines(poly)
+                if poly_line:
+                    for indx in range(poly_line.GetPointCount() - 1):
+                        poi = poly_line.GetPoint(indx)
+
+                        points.append(poi)
+
+                        poi_next = poly_line.GetPoint(indx + 1)
+
+                        dis = math.sqrt(math.pow(poi[0] - poi_next[0], 2) + math.pow(poi[1] - poi_next[1], 2))
+                        # print(dis)
+                        if dis > distance_threshold:
+                            times = int(dis / distance_threshold)
+                            # print(times)
+                            for ts in range(1, times + 1):
+                                x, y = self.get_subpoint(poi[0], poi[1], poi_next[0], poi_next[1], (ts * distance_threshold))
+                                points.append([x, y])
+
+            points = np.array([[p[0], p[1]] for p in points])
+
+
+            cp = []
+            for pp in points:
+                poi = ogr.Geometry(ogr.wkbPoint)
+                poi.AddPoint(pp[0], pp[1])
+                cp.append(poi)
+
+
+            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)
+            merge_p = merge_p.Simplify(self.distance_buffer_threshold/10)
+            merge_polygons.append(merge_p)
+
+        return merge_polygons
+        # return cp
+
+    def get_subpoint(self, x0, y0, x1, y1, len_size):
+
+        len_size = len_size
+        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
+
+        return x2, y2
+
+    def merge_delauney(self,trees):
+
+        merge_polygons = []
+        count = 1
+
+        for tree in trees:
+
+            polygons = list(tree.get_polygons())
+
+            mp: Geometry = ogr.Geometry(ogr.wkbMultiPolygon)
+            for p in polygons:
+                mp.AddGeometry(p)
+
+
+
+            delauney: Geometry = mp.DelaunayTriangulation()
+
+            intp = []
+            for p in delauney:
+
+                res = p
+                ispass = False
+                for pp in polygons:
+                    check:Geometry = res.Difference(pp)
+                    if check.IsEmpty():
+                        ispass = True
+                        break
+                    else:
+                        res = check
+                if not ispass:
+                    res_extend = self.extend_tri(res)
+                    intp.append(res_extend)
+
+            # polygons.extend(intp)
+            #
+            # mp: Geometry = ogr.Geometry(ogr.wkbMultiPolygon)
+            #
+            # for p in polygons:
+            #     mp.AddGeometry(p)
+            # merge:Geometry = mp.UnionCascaded()
+            # merge_polygons.append(merge)
+
+        return intp
+
+
+
+    def extend_tri(self,triangle:Geometry):
+
+        center:Geometry = triangle.Centroid()
+        tri_line: Geometry = self.get_polygon_lines(triangle)
+
+        result_geo = ogr.Geometry(ogr.wkbPolygon)
+        result_geo_ring:Geometry = ogr.Geometry(ogr.wkbLinearRing)
+        for point in tri_line.GetPoints():
+            radpoint = self.get_radpoint(center.GetX(),center.GetY(),point[0],point[1],0.00000001)
+
+            result_geo_ring.AddPoint(radpoint.GetX(),radpoint.GetY())
+        result_geo.AddGeometry(result_geo_ring)
+
+        return result_geo
+
+    def get_radpoint(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)
+
+        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
+
+        other_point = ogr.Geometry(ogr.wkbPoint)
+        other_point.AddPoint(x2, y2)
+
+        return  other_point
+
+
+    def process(self):
+        polygons = self.data.get_polygons()
+        rps = []
+
+        for poly in polygons:
+            rp = self.get_polygon_reprecent_point(poly)
+            rps.append(rp)
+        print("完成代表点计算")
+
+        #代表点要融合
+        delauney_lines,rps_sub = self.create_delauney(rps)
+
+        print("完成树生成计算")
+        min_delauney = self.create_min_delauney(delauney_lines,rps_sub)
+        print("完成最小生成树计算")
+        trees = self.cut_delauney_min_tree(min_delauney,self.distance_buffer_threshold)
+        print("完成树裁剪")
+
+        polygons = self.merge_concave(trees,self.distance_buffer_threshold)
+        # polygons = self.merge_delauney(trees)
+
+        return [p.ExportToWkt() for p in polygons]
+
+if __name__ == '__main__':
+
+
+    # print(sd.get_polygons()[0])
+    # sd.close()
+
+    t1 = time.time()
+    sd = ShapeData(r"J:\Data\制图综合result\zhongshansub2.shp")
+
+    zh = Zonghe(sd)
+    wkts = zh.process()
+    result = r"J:\Data\制图综合result\zhongshansub2_m2.shp"
+    ShapeData.create_shp_fromwkts(result,"zh",wkts)
+
+    print(time.time()-t1)

+# 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.ConcaveHull import concaveHull
+from test.zonghe.ShapeData import ShapeData
+from test.zonghe.Delauney import DelauneyNode,DelauneyLine,DelauneyMinTree
+from test.zonghe.RepresentPoint import RepresentPoint
+
+import time
+
+import numpy as np
+from multiprocessing import Process
+from threading import Thread
+class MyThread(Thread):
+    def __init__(self,func,args=()):
+        super(MyThread,self).__init__()
+        self.func = func
+        self.args = args
+    def run(self):
+        self.result = self.func(*self.args)
+    def get_result(self):
+        try:
+            return self.result
+        except Exception:
+            return None
+
+class Zonghe:
+
+    data:ShapeData
+    area_threshold = 10.0
+    buffer_threshold = 100.0
+    distance_buffer_threshold = 1.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
+        line = None
+        oneside_area, otherside_area = 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 oneside_area > 0 and otherside_area >0:
+                if abs(oneside_area-otherside_area) < area_thre:
+                    break
+
+            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
+
+        #第一条线就是均分线
+        if angle1 is None or angle2 is None:
+            pass
+        else:
+            #如果面积差大于阈值，二分法查找
+            if abs(oneside_area - otherside_area) > area_thre:
+                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 = {}
+        delauney_lines = []
+
+        rps_sub = set()
+        for rp in rps:
+            multipoint.AddGeometry(rp.point)
+            # 多边形代表点相同，需要融合多边形
+            if rp_dict.get(rp.point.GetPoint(0).__str__()):
+                rp_exist = rp_dict.get(rp.point.GetPoint(0).__str__())
+                rp_exist.polygon = rp_exist.polygon.Union(rp.polygon)
+            else:
+                rp_dict[rp.point.GetPoint(0).__str__()] = rp
+                rps_sub.add(rp)
+
+        delauney:Geometry =  multipoint.DelaunayTriangulation(bOnlyEdges=True)
+
+        for index in range(delauney.GetGeometryCount()):
+            triangle_line:Geometry = delauney.GetGeometryRef(index)
+            p = triangle_line.GetPoint(0)
+            p_next = triangle_line.GetPoint(1)
+
+            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,rps_sub
+
+    def create_min_delauney(self,delauney_lines,rps):
+        '''
+        使用prim算法计算最小生成树
+        :param delauney_lines:
+        :param polygons:
+        :return:
+        '''
+
+        kv = dict()
+        print(len(rps))
+
+
+        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)
+
+        print(len(kv))
+
+        print("完成三角边倒排索引")
+
+        #初始化树
+        delauney_lines = sorted(delauney_lines,key=lambda x:x.distance)
+        delauney_min_tree = DelauneyMinTree()
+
+        delauney_min_tree.append(delauney_lines[0],kv)
+
+        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])
+
+        print("完成最小生成树初始化")
+
+        while rps_set.__len__()>0:
+
+            print(rps_set.__len__())
+            t1 = time.time()
+            # relate_delauney_line = delauney_min_tree.get_other_relate_delauney_lines(kv,rps_set)
+
+            relate_delauney_line = delauney_min_tree.touch_delauney_lines
+
+
+            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,kv)
+
+            rps_set.remove(other_rp)
+
+            kv[min_distance_line.rp1].remove(min_distance_line)
+            kv[min_distance_line.rp2].remove(min_distance_line)
+
+            # print(time.time() - t2)
+        print("完成最小生成树")
+
+        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 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() in [0,1]:
+            polygon_line : Geometry = ogr.ForceToLineString(polygon)
+        # 有孔Polygon
+        else:
+            polygon_line : Geometry = ogr.ForceToLineString(polygon.GetGeometryRef(0))
+
+        return polygon_line
+
+    def get_subpoint(self, x0, y0, x1, y1, len_size):
+
+        len_size = len_size
+        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
+
+        return x2, y2
+
+    def merge_concave_parallel(self, trees, distance_threshold):
+        merge_polygons = []
+
+        tree_set = set(trees)
+        thread_set = set()
+
+        while len(tree_set) > 0 :
+            while len(thread_set) < 16 and len(tree_set) > 0:
+
+                print(len(tree_set))
+                tree = tree_set.pop()
+                polygons = list(tree.get_polygons())
+                merge_threading = MyThread(self.merge,args=(polygons,distance_threshold))
+                merge_threading.start()
+                thread_set.add(merge_threading)
+
+
+            remove_list = []
+            for thread in thread_set:
+
+
+                if not thread.is_alive():
+                    merge_polygons.append(thread.get_result())
+                    remove_list.append(thread)
+            for re in remove_list:
+                thread_set.remove(re)
+                del re
+
+        for thread in thread_set:
+            thread.join()
+            merge_polygons.append(thread.get_result())
+
+        return merge_polygons
+
+    def merge(self,polygons,distance_threshold):
+        merge_raw: Geometry = polygons[0]
+        for p in polygons:
+            merge_raw = merge_raw.Union(p)
+
+        points = []
+
+        for i in range(merge_raw.GetGeometryCount()):
+
+            poly = merge_raw.GetGeometryRef(i)
+            poly_line: Geometry = self.get_polygon_lines(poly)
+            if poly_line:
+                for indx in range(poly_line.GetPointCount() - 1):
+                    poi = poly_line.GetPoint(indx)
+
+                    points.append(poi)
+
+                    poi_next = poly_line.GetPoint(indx + 1)
+
+                    dis = math.sqrt(math.pow(poi[0] - poi_next[0], 2) + math.pow(poi[1] - poi_next[1], 2))
+                    # print(dis)
+                    if dis > distance_threshold:
+                        times = int(dis / distance_threshold)
+                        # print(times)
+                        for time in range(1, times + 1):
+                            x, y = self.get_subpoint(poi[0], poi[1], poi_next[0], poi_next[1],
+                                                     (time * distance_threshold))
+                            points.append([x, y])
+
+        points = np.array([[p[0], p[1]] for p in points])
+
+        cp = []
+        for pp in points:
+            poi = ogr.Geometry(ogr.wkbPoint)
+            poi.AddPoint(pp[0], pp[1])
+            cp.append(poi)
+
+        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)
+        merge_p = merge_p.Simplify(self.distance_buffer_threshold / 10)
+
+        return merge_p
+
+    def process(self):
+        polygons = self.data.get_polygons()
+        rps = []
+
+        for poly in polygons:
+            rp = self.get_polygon_reprecent_point(poly)
+            rps.append(rp)
+        print("完成代表点计算")
+
+        #代表点要融合
+        delauney_lines,rps_sub = self.create_delauney(rps)
+
+        print("完成树生成计算")
+        min_delauney = self.create_min_delauney(delauney_lines,rps_sub)
+        print("完成最小生成树计算")
+        trees = self.cut_delauney_min_tree(min_delauney,self.distance_buffer_threshold)
+        print("完成树裁剪")
+
+        polygons = self.merge_concave_parallel(trees,self.distance_buffer_threshold)
+
+        return [p.ExportToWkt() for p in polygons]
+
+if __name__ == '__main__':
+
+
+    # print(sd.get_polygons()[0])
+    # sd.close()
+
+    t1 = time.time()
+    sd = ShapeData(r"J:\Data\制图综合result\zhongshansub.shp")
+
+    zh = Zonghe(sd)
+    wkts = zh.process()
+    result = r"J:\Data\制图综合result\zhongshansub_pm.shp"
+    ShapeData.create_shp_fromwkts(result,"zh",wkts)
+
+    print(time.time()-t1)

+# coding=utf-8
+#author:        4N
+#createtime:    2022/3/8
+#email:         nheweijun@sina.com
+
+import math
+
+from osgeo import ogr
+from osgeo.ogr import *
+
+from test.zonghe.ShapeData import ShapeData
+
+
+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)
+
+        p1_int = []
+        for indx,p in enumerate(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])
+                p1_int.append(p)
+                multipoint.AddGeometry(pppp)
+
+
+        for indx,p in enumerate(p2_points):
+
+            p_t = False
+            for pp in p1_int:
+                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()
+
+        for indx in range(delauney.GetGeometryCount()):
+            tri:Geometry = delauney.GetGeometryRef(indx)
+            dif1:Geometry = tri.Difference(p1)
+            dif2:Geometry = tri.Difference(p2)
+
+        # return [multipoint.ExportToWkt()]
+
+        return [p.ExportToWkt() for p in delauney]
+
+
+
+    #这个是融合面
+    def mergebac(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)
+
+
+        p1_point_start = None
+        p1_point_end = None
+
+        for indx,p in enumerate(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])
+
+                if p1_point_start is None:
+                    p1_point_start = indx
+
+                multipoint.AddGeometry(pppp)
+                p1_point_end = indx
+
+        p1_int = p1_points[p1_point_start:p1_point_end+1]
+
+
+        p2_point_start = None
+        p2_point_end = None
+
+
+        for indx,p in enumerate(p2_points):
+
+            p_t = False
+            for pp in p1_int:
+                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])
+
+                if p2_point_start is None:
+                    p2_point_start = indx
+
+                multipoint.AddGeometry(pppp)
+
+                p2_point_end = indx
+
+        delauney: Geometry = multipoint.DelaunayTriangulation()
+
+        p2_int = p2_points[p2_point_start:p2_point_end + 1]
+        if p2_point_end == len(p2_points)-1:
+            p2_int = p2_points[p2_point_start:p2_point_end + 1]
+
+        p1_int.extend(p2_int)
+
+        print(p2_point_start,p2_point_end)
+
+        multipoint: Geometry = ogr.Geometry(ogr.wkbMultiPoint)
+
+        for p in p1_int:
+            pppp = ogr.Geometry(ogr.wkbPoint)
+            pppp.AddPoint(p[0], p[1])
+            multipoint.AddGeometry(pppp)
+
+        return [multipoint.ExportToWkt()]
+
+        # return [p.ExportToWkt() for p in delauney]
+
+
+    def merge4(self,polygons):
+
+
+        mp: Geometry = ogr.Geometry(ogr.wkbMultiPolygon)
+        for p in polygons:
+            mp.AddGeometry(p)
+
+
+
+        delauney: Geometry = mp.DelaunayTriangulation()
+
+        intp = []
+        for p in delauney:
+
+            res = p
+            ispass = False
+            for pp in polygons:
+                check:Geometry = res.Difference(pp)
+                if check.IsEmpty():
+                    ispass = True
+                    break
+                else:
+                    res = check
+            if not ispass:
+                res_buffer = self.extend_tri(res)
+
+                intp.append(res_buffer)
+
+
+        polygons.extend(intp)
+
+        mp: Geometry = ogr.Geometry(ogr.wkbMultiPolygon)
+
+        for p in polygons:
+            mp.AddGeometry(p)
+
+
+
+        merge:Geometry = mp.UnionCascaded()
+
+        return [merge.ExportToWkt()]
+
+        # return [p.ExportToWkt() for p in intp]
+
+
+    def extend_tri(self,triangle:Geometry):
+
+        center:Geometry = triangle.Centroid()
+        tri_line: Geometry = ogr.ForceToLineString(triangle)
+
+        result_geo = ogr.Geometry(ogr.wkbPolygon)
+        result_geo_ring:Geometry = ogr.Geometry(ogr.wkbLinearRing)
+        for point in tri_line.GetPoints():
+            radpoint = self.get_radpoint(center.GetX(),center.GetY(),point[0],point[1],0.00000001)
+
+            result_geo_ring.AddPoint(radpoint.GetX(),radpoint.GetY())
+        result_geo.AddGeometry(result_geo_ring)
+
+        return result_geo
+
+    def get_radpoint(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)
+
+        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
+
+        other_point = ogr.Geometry(ogr.wkbPoint)
+        other_point.AddPoint(x2, y2)
+
+        return  other_point
+
+
+
+    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.merge(polygons[0],polygons[1])
+
+        return self.merge4(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\制图综合result\ex1.shp")
+
+    ext =  sd.layer.GetExtent()
+    threshold = (ext[1]-ext[0])/20
+
+    zh = Zonghe(sd)
+    wkts = zh.process()
+
+    result = r"J:\Data\制图综合result\ex1_1.shp"
+    ShapeData.create_shp_fromwkts(result,"zh",wkts)
+

@@ -150,26 +150,53 @@ class Zonghe:
 
 
         return None
         return None
 
 
-    def merge3(self,polygons):
+    def merge3(self,polygons,distance_threshold):
+
+        merge_raw :Geometry = polygons[0]
+        for p in polygons:
+            merge_raw = merge_raw.Union(p)
 
 
-        multi_poly = ogr.Geometry(ogr.wkbMultiPolygon)
 
 
         points = []
         points = []
-        for poly in polygons:
-            multi_poly.AddGeometry(poly)
+
+        for i in range(merge_raw.GetGeometryCount()):
+
+            poly = merge_raw.GetGeometryRef(i)
             if poly.GetGeometryCount() == 1:
             if poly.GetGeometryCount() == 1:
                 poly_line: Geometry = ogr.ForceToLineString(poly)
                 poly_line: Geometry = ogr.ForceToLineString(poly)
             # 有孔Polygon
             # 有孔Polygon
             else:
             else:
                 poly_line: Geometry = ogr.ForceToLineString(poly.GetGeometryRef(0))
                 poly_line: Geometry = ogr.ForceToLineString(poly.GetGeometryRef(0))
 
 
-            points.extend(poly_line.GetPoints()[:-1])
 
 
+            for indx in range(poly_line.GetPointCount() -1) :
+                poi = poly_line.GetPoint(indx)
+
+                points.append(poi)
+
+                poi_next = poly_line.GetPoint(indx+1)
+
+                dis = math.sqrt(math.pow(poi[0] - poi_next[0],2) + math.pow(poi[1] - poi_next[1],2))
+                # print(dis)
+                if dis > distance_threshold:
+                    times = int(dis/distance_threshold)
+                    # print(times)
+                    for time in range(1,times+1):
+                        x,y = self.get_radline(poi[0],poi[1],poi_next[0],poi_next[1],(time*distance_threshold))
+                        points.append([x,y])
 
 
         points = np.array([[p[0],p[1]] for p in points])
         points = np.array([[p[0],p[1]] for p in points])
 
 
+
+        for pp in points:
+            poi = ogr.Geometry(ogr.wkbPoint)
+            poi.AddPoint(pp[0],pp[1])
+            print("\"" + poi.ExportToWkt() + "\",")
+
+
         merge_p_line = concaveHull(points, 3)
         merge_p_line = concaveHull(points, 3)
 
 
+
         ring = ogr.Geometry(ogr.wkbLinearRing)
         ring = ogr.Geometry(ogr.wkbLinearRing)
 
 
         merge_p = ogr.Geometry(ogr.wkbPolygon)
         merge_p = ogr.Geometry(ogr.wkbPolygon)
@@ -194,32 +221,46 @@ class Zonghe:
         return [merge_p.ExportToWkt()]
         return [merge_p.ExportToWkt()]
 
 
 
 
+    def get_radline(self,x0,y0,x1,y1,len_size):
 
 
+        len_size = len_size
+        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
 
 
-    def merge4(self,p1,p2):
+            if y1 < y0:
+                y2 = y0 - dy
+            else:
+                y2 = y0 + dy
 
 
-        polygons = []
-        points = []
+        return x2,y2
 
 
-        mp: Geometry = ogr.Geometry(ogr.wkbMultiPolygon)
-        mp.AddGeometry(p1)
-        mp.AddGeometry(p2)
+    def merge4(self,polygons):
 
 
 
 
-        delauney: Geometry = mp.DelaunayTriangulation(bOnlyEdges=True)
-        merge_polygon = None
+        mp: Geometry = ogr.Geometry(ogr.wkbMultiPolygon)
+        for p in polygons:
+            mp.AddGeometry(p)
 
 
-        for p in delauney:
 
 
-            p_c = p.Centroid()
 
 
-            # print("\"" + p_c.ExportToWkt() + "\",")
-            # if p1.Contains(p) or p2.Contains(p):
-            #     continue
-            print("\"" + p.ExportToWkt() + "\",")
+        delauney: Geometry = mp.DelaunayTriangulation()
 
 
 
 
-        return polygons
+        return [p.ExportToWkt() for p in delauney]
 
 
     def merge5(self,p1:Geometry,p2):
     def merge5(self,p1:Geometry,p2):
 
 
@@ -239,7 +280,9 @@ class Zonghe:
     def process(self):
     def process(self):
         polygons = self.data.get_polygons()
         polygons = self.data.get_polygons()
 
 
-        return self.merge3(polygons)
+        # return self.merge3(polygons,0.00007)
+
+        return self.merge4(polygons)
 
 
 
 
 if __name__ == '__main__':
 if __name__ == '__main__':
@@ -270,7 +313,7 @@ if __name__ == '__main__':
     #
     #
 
 
 
 
-    sd = ShapeData(r"J:\Data\制图综合\concave.shp")
+    sd = ShapeData(r"J:\Data\制图综合\concave5.shp")
 
 
     ext =  sd.layer.GetExtent()
     ext =  sd.layer.GetExtent()
     threshold = (ext[1]-ext[0])/20
     threshold = (ext[1]-ext[0])/20
@@ -280,7 +323,6 @@ if __name__ == '__main__':
     zh = Zonghe(sd)
     zh = Zonghe(sd)
     wkts = zh.process()
     wkts = zh.process()
 
 
-
-    result = r"J:\Data\制图综合\concave_merge.shp"
+    result = r"J:\Data\制图综合\concave5_merge.shp"
     ShapeData.create_shp_fromwkts(result,"zh",wkts)
     ShapeData.create_shp_fromwkts(result,"zh",wkts)