2022.03.22 修改模块回调问题

@@ -6,8 +6,10 @@ from app.util.component.ParameterUtil import ParameterUtil
 from flask import current_app,request
 import traceback
 import json
-from app.modules.auth.models import OAuth2Token,User,db
+from app.models import db
+from app.modules.auth.models import OAuth2Token,User
 import app
+
 class ApiTemplate:
     #模板方法
     para = dict()
@@ -28,10 +30,11 @@ class ApiTemplate:
             if app.GLOBAL_DIC.get(token):
                 self.para["creator"] = app.GLOBAL_DIC.get(token)
             else:
-                user = User.query.join(OAuth2Token).filter(OAuth2Token.access_token == token).one_or_none()
-                if user:
-                    self.para["creator"] = user.username
-                    app.GLOBAL_DIC[token] = user.username
+                sql = "select username from dmap_user join dmap_oauth2_token on dmap_user.id=dmap_oauth2_token.user_id where dmap_oauth2_token.access_token='{}'".format(token)
+                result = db.session.execute(sql).fetchone()
+                if result:
+                    self.para["creator"] = result[0]
+                    app.GLOBAL_DIC[token] = result[0]
 
     def para_check(self):
         pass

+# 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 numpy as np
+
+class Zonghe:
+
+    data:ShapeData
+    area_threshold = 100.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
+
+        for angle in direction:
+
+            if in_hole:
+                line = self.get_doubleline(center.GetX(),center.GetY(),angle,line_length)
+            else:
+                line = self.get_line(center.GetX(),center.GetY(),angle,line_length)
+
+            oneside_area, otherside_area = self.jude(polygon, line, angle, line_length)
+
+
+            if last_state is None:
+
+                angle2 = angle
+
+                if oneside_area == 0 and otherside_area==0:
+                    last_state = 0
+                else:
+                    last_state = 1 if oneside_area < otherside_area else -1
+            else:
+
+                angle1 = angle2
+                angle2 = angle
+
+                if oneside_area == 0 and otherside_area==0:
+                    now_state = 0
+                else:
+                    now_state = 1 if oneside_area < otherside_area else -1
+
+                if last_state* now_state < 0:
+                    # print(angle)
+                    break
+                else:
+                    last_state = now_state
+
+        line = self.division(polygon, center, in_hole, line_length, angle1, angle2, area_thre)
+
+        # 提前return
+        if in_hole:
+
+            point_intersection: Geometry = line.Intersection(in_hole_line)
+            plist = []
+            for i in range(point_intersection.GetGeometryCount()):
+                pi = point_intersection.GetGeometryRef(i)
+                plist.append(pi)
+            plist_d = [(pi,center.Distance(pi)) for pi in plist]
+            plist_d.sort(key=lambda x:x[1])
+
+            nearest_point= plist_d[0][0]
+            firt_short_length = short_length/20.0
+            radline,other_point = self.get_radline(center.GetX(),center.GetY(),nearest_point.GetX(),nearest_point.GetY(),firt_short_length)
+            # print("\"" + other_point.ExportToWkt() + "\",")
+
+
+            while not polygon.Contains(other_point):
+
+                firt_short_length = firt_short_length / 2
+                radline, other_point = self.get_radline(center.GetX(), center.GetY(), nearest_point.GetX(),
+                                                        nearest_point.GetY(), firt_short_length)
+
+
+
+            return RepresentPoint(other_point,polygon)
+
+
+        point_intersection:Geometry = line.Intersection(polygon_line)
+        xlist = []
+        ylist = []
+        plist = []
+
+        for i in range(point_intersection.GetGeometryCount()):
+            pi = point_intersection.GetGeometryRef(i)
+            xlist.append(pi.GetX())
+            ylist.append(pi.GetY())
+            plist.append(pi)
+
+        reprecent_point = [(min(xlist)+max(xlist))/2.0,(min(ylist)+max(ylist))/2.0]
+
+        reprecent_point = self.creat_point(reprecent_point)
+
+        if not polygon.Contains(reprecent_point):
+            plist_d = [(pi,center.Distance(pi)) for pi in plist]
+            plist_d.sort(key=lambda x:x[1])
+            reprecent_point = self.creat_point(((plist_d[0][0].GetX()+ plist_d[1][0].GetX())/2.0, (plist_d[0][0].GetY()+plist_d[1][0].GetY())/2.0))
+
+        # print("\""+reprecent_point.ExportToWkt() + "\",")
+
+        rp = RepresentPoint(reprecent_point,polygon)
+
+        return rp
+
+
+
+
+
+    def get_line(self, x0, y0 , angle, r):
+
+        line: Geometry = ogr.Geometry(ogr.wkbLineString)
+        line.AddPoint(x0, y0)
+        hudu = angle/360.0 * 2 * math.pi
+        dx = math.sin(hudu) * r
+        dy = math.cos(hudu) * r
+        line.AddPoint(x0+dx, y0+dy)
+
+        return line
+
+    def get_doubleline(self, x0, y0 , angle, r):
+
+        line: Geometry = ogr.Geometry(ogr.wkbLineString)
+
+        hudu = angle/360.0 * 2 * math.pi
+        dx = math.sin(hudu) * r
+        dy = math.cos(hudu) * r
+
+        hudu = (angle+180)/360.0 * 2 * math.pi
+        dx2 = math.sin(hudu) * r
+        dy2 = math.cos(hudu) * r
+
+        line.AddPoint(x0 + dx2, y0 + dy2)
+        line.AddPoint(x0+dx, y0+dy)
+
+        return line
+
+
+    def get_radline(self,x0,y0,x1,y1,len_size):
+        a2 = math.pow((x1-x0),2)
+        b2 = math.pow((y1-y0),2)
+        len_size = len_size+ math.sqrt(a2+b2)
+
+        line: Geometry = ogr.Geometry(ogr.wkbLineString)
+        line.AddPoint(x0, y0)
+        tanx = abs(y1-y0) / abs(x1-x0)
+        cosx = 1 / math.sqrt(1 + tanx *tanx)
+        dx = cosx* len_size
+        dy = tanx * dx
+
+        if x1 == x0:
+            x2 = x0
+            if y1 > y0 :
+                y2 = y0 + abs(dy)
+            else:
+                y2 = y0 - abs(dy)
+        else:
+            if x1 < x0:
+                x2 = x0 - dx
+            else:
+                x2 = x0 + dx
+
+            if y1 < y0:
+                y2 = y0 - dy
+            else:
+                y2 = y0 + dy
+
+        line.AddPoint(x2, y2)
+
+        other_point = ogr.Geometry(ogr.wkbPoint)
+        other_point.AddPoint(x2, y2)
+
+        return line , other_point
+
+    def get_direction(self,polygon,center,line_length):
+        '''
+        判断旋转方向
+        :param polygon:
+        :param center:
+        :param line_length:
+        :return:
+        '''
+
+        line = self.get_line(center.GetX(),center.GetY(),0,line_length)
+
+        oneside_area, otherside_area = self.jude(polygon,line,0,line_length)
+
+        if oneside_area == 0 and otherside_area == 0 :
+            return range(0, 390, 30)
+        else:
+            line = self.get_line(center.GetX(), center.GetY(), 10, line_length)
+            oneside_area_contrast, otherside_area_contrast = self.jude(polygon, line, 30, line_length)
+
+            if oneside_area_contrast == 0 and otherside_area_contrast == 0:
+                return range(360, -30, -30)
+            else:
+                if abs(oneside_area - otherside_area) > abs(oneside_area_contrast - otherside_area_contrast):
+                    return range(0, 390, 30)
+                else:
+                    return range(360, -30, -30)
+
+
+    def jude(self,polygon,line,angle,line_length):
+
+        line_buffer: Geometry = line.Buffer(line_length / self.buffer_threshold)
+        clip_geom: Geometry = polygon.Difference(line_buffer)
+
+        oneside_area = 0
+        otherside_area = 0
+
+        gc = clip_geom.GetGeometryCount()
+        if gc <= 1:
+            oneside_area = 0
+            otherside_area = 0
+        else:
+            triangle1, triangle2 = self.get_side_triangle(line, angle)
+
+
+            for gi in range(clip_geom.GetGeometryCount()):
+                clip_geom_i: Geometry = clip_geom.GetGeometryRef(gi)
+                it = clip_geom_i.Intersection(triangle1)
+
+                if clip_geom_i.Intersect(triangle1) and it.Buffer(line_length / self.buffer_threshold).Intersect(
+                        line_buffer):
+                    oneside_area += clip_geom_i.GetArea()
+                else:
+                    otherside_area += clip_geom_i.GetArea()
+        return oneside_area, otherside_area
+
+    def division(self,polygon,center,in_hole,line_length,angle1,angle2,area_thre):
+
+        mid_angle = (angle1 + angle2)/2.0
+        if in_hole:
+            line = self.get_doubleline(center.GetX(), center.GetY(), mid_angle, line_length)
+        else:
+            line = self.get_line(center.GetX(), center.GetY(), mid_angle, line_length)
+
+        one,other = self.jude(polygon,line,mid_angle,line_length)
+
+        while abs(one - other) > area_thre:
+            if one > other:
+                if angle1 > angle2:
+                    return self.division(polygon, center, in_hole, line_length, angle2, mid_angle, area_thre)
+                else:
+                    return self.division(polygon, center, in_hole, line_length, mid_angle, angle1, area_thre)
+            else:
+                if angle1 > angle2:
+                    return self.division(polygon, center, in_hole, line_length, mid_angle, angle1, area_thre)
+                else:
+                    return self.division(polygon, center, in_hole, line_length, angle2, mid_angle, area_thre)
+        return line
+
+
+    def get_side_triangle(self,line:Geometry,angle):
+        '''
+        获取方向三角形
+        :param line:
+        :param angle:
+        :return:
+        '''
+
+        start = line.GetPoint(0)
+        end = line.GetPoint(1)
+
+
+        angle_t = angle-30
+        hudu = angle_t/360.0 * 2 * math.pi
+
+        r = line.Length() / (2*math.cos( math.pi/6))
+        dx = math.sin(hudu) * r
+        dy = math.cos(hudu) * r
+
+
+        ring1 = ogr.Geometry(ogr.wkbLinearRing)
+        ring1.AddPoint(start[0], start[1])
+        ring1.AddPoint(start[0] + dx, start[1] + dy)
+        ring1.AddPoint(end[0], end[1])
+        ring1.AddPoint(start[0], start[1])
+        triangle1 = ogr.Geometry(ogr.wkbPolygon)
+        triangle1.AddGeometry(ring1)
+
+
+        angle_t = angle+30
+        hudu = angle_t/360.0 * 2 * math.pi
+
+        r = line.Length() / (2*math.cos( math.pi/6))
+
+        dx = math.sin(hudu) * r
+        dy = math.cos(hudu) * r
+
+        ring2 = ogr.Geometry(ogr.wkbLinearRing)
+        ring2.AddPoint(start[0], start[1])
+        ring2.AddPoint(start[0]+dx, start[1]+dy)
+        ring2.AddPoint(end[0], end[1])
+        ring2.AddPoint(start[0], start[1])
+        triangle2 = ogr.Geometry(ogr.wkbPolygon)
+        triangle2.AddGeometry(ring2)
+
+        return triangle1,triangle2
+
+
+    def creat_point(self,tuple):
+        p: Geometry = ogr.Geometry(ogr.wkbPoint)
+        p.AddPoint(tuple[0], tuple[1])
+        return p
+
+
+
+
+    def create_delauney(self,rps):
+        multipoint: Geometry = ogr.Geometry(ogr.wkbMultiPoint)
+
+        rp_dict = {}
+        line_set = set()
+        delauney_lines = []
+
+        for rp in rps:
+            multipoint.AddGeometry(rp.point)
+
+            rp_dict[rp.point.GetPoint(0).__str__()] = rp
+
+        delauney:Geometry =  multipoint.DelaunayTriangulation()
+
+        for index in range(delauney.GetGeometryCount()):
+            triangle = delauney.GetGeometryRef(index)
+
+            triangle_line:Geometry = ogr.ForceToLineString(triangle)
+
+
+            for pi in range(triangle_line.GetPointCount()-1):
+                p = triangle_line.GetPoint(pi)
+                p_next = triangle_line.GetPoint(pi+1)
+
+                if p[0] < p_next[0]:
+                    key = "{}{}".format(p.__str__(),p_next.__str__())
+                elif p[0] > p_next[0]:
+                    key = "{}{}".format(p_next.__str__(), p.__str__())
+                else:
+                    if p[1] < p_next[1]:
+                        key = "{}{}".format(p.__str__(), p_next.__str__())
+                    else:
+                        key = "{}{}".format(p_next.__str__(), p.__str__())
+
+                if not line_set.__contains__(key):
+                    line_set.add(key)
+                    rp1 = rp_dict.get(p.__str__())
+                    rp2 = rp_dict.get(p_next.__str__())
+
+                    delauney_line = DelauneyLine(rp1,rp2)
+                    delauney_lines.append(delauney_line)
+
+        return delauney_lines
+
+
+
+    def create_min_delauney(self,delauney_lines,rps,polygons):
+        '''
+        使用prim算法计算最小生成树
+        :param delauney_lines:
+        :param polygons:
+        :return:
+        '''
+
+        kv = {}
+        rps_set = set(rps)
+
+        for delauney_line in delauney_lines:
+            if kv.get(delauney_line.rp1):
+
+                kv[delauney_line.rp1].add(delauney_line)
+            else:
+                kv[delauney_line.rp1] = set()
+                kv[delauney_line.rp1].add(delauney_line)
+
+            if kv.get(delauney_line.rp2):
+
+                kv[delauney_line.rp2].add(delauney_line)
+            else:
+
+                kv[delauney_line.rp2] = set()
+                kv[delauney_line.rp2].add(delauney_line)
+
+
+        #初始化树
+        delauney_lines = sorted(delauney_lines,key=lambda x:x.distance)
+
+
+        delauney_min_tree = DelauneyMinTree()
+        delauney_min_tree.append(delauney_lines[0])
+
+        rps_set.remove(delauney_lines[0].rp1)
+        rps_set.remove(delauney_lines[0].rp2)
+        kv[delauney_lines[0].rp1].remove(delauney_lines[0])
+        kv[delauney_lines[0].rp2].remove(delauney_lines[0])
+
+
+        while rps_set.__len__()>0:
+
+            relate_delauney_line = delauney_min_tree.get_other_relate_delauney_lines(kv)
+
+            relate_delauney_line_sort = sorted(relate_delauney_line,key=lambda x:x.distance)
+
+            min_distance_line = relate_delauney_line_sort[0]
+
+            #删除
+            other_rp = delauney_min_tree.get_other_rp(min_distance_line)
+            if other_rp is None:
+                kv[min_distance_line.rp1].remove(min_distance_line)
+                kv[min_distance_line.rp2].remove(min_distance_line)
+                continue
+
+            delauney_min_tree.append(min_distance_line)
+
+            rps_set.remove(other_rp)
+
+            kv[min_distance_line.rp1].remove(min_distance_line)
+            kv[min_distance_line.rp2].remove(min_distance_line)
+
+
+
+        for l in delauney_min_tree.delauney_lines_list:
+            line = l.line
+            # print("\"" + line.ExportToWkt() + "\",")
+
+        # print(delauney_min_tree.delauney_lines_list.__len__())
+
+        return delauney_min_tree
+
+    def cut_delauney_min_tree(self,delauney_min_tree:DelauneyMinTree,threshold):
+
+        trees = []
+        kv = {}
+
+        count = 0
+
+        for delauney_line in delauney_min_tree.delauney_lines_list:
+
+            dn1 = kv.get(delauney_line.rp1)
+            dn2 = kv.get(delauney_line.rp2)
+
+
+            if not dn1:
+                dn1 = DelauneyNode(delauney_line.rp1)
+                kv[delauney_line.rp1] = dn1
+
+                if not dn2:
+
+                    trees.append(dn1)
+
+                    dn2 = DelauneyNode(delauney_line.rp2)
+                    kv[delauney_line.rp2] = dn2
+
+                    if delauney_line.distance < threshold:
+                        dn1.add_children(dn2,delauney_line)
+                        count +=1
+
+                    else:
+                        #断开
+                        trees.append(dn2)
+                else:
+
+                    if delauney_line.distance < threshold:
+                        dn2.add_children(dn1,delauney_line)
+                        count += 1
+                    else:
+                        trees.append(dn1)
+            else:
+
+                dn2 = DelauneyNode(delauney_line.rp2)
+                kv[delauney_line.rp2] = dn2
+
+                if delauney_line.distance < threshold:
+                    dn1.add_children(dn2,delauney_line)
+                    count += 1
+                else:
+                    trees.append(dn2)
+
+
+        return trees
+
+    def merge(self,trees):
+
+        polygons = []
+        for tree in trees:
+
+            merge_polygon_raw = list()
+
+            dls = tree.get_delauney_lines()
+            for dl in dls:
+                p1 = dl.rp1.polygon
+                p2 = dl.rp2.polygon
+
+                merge_polygon_raw.append(p1)
+                merge_polygon_raw.append(p2)
+                line = dl.line
+
+                p1_lines = self.get_polygon_lines(p1)
+                p2_lines = self.get_polygon_lines(p2)
+
+                p1_intersect_lines = [l for l in p1_lines if l.Intersect(line)]
+                p2_intersect_lines = [l for l in p2_lines if l.Intersect(line)]
+
+                if p1_intersect_lines.__len__()>1:
+                    p1_intersect_lines = sorted(p1_intersect_lines,key=lambda x:p2.Distance(x))
+                    p1_intersect_line_rep = p1_intersect_lines[0]
+                elif p1_intersect_lines.__len__()==1:
+                    p1_intersect_line_rep = p1_intersect_lines[0]
+                else:
+                    continue
+
+
+                if p2_intersect_lines.__len__()>1:
+                    p2_intersect_lines = sorted(p2_intersect_lines,key=lambda x:p1.Distance(x))
+                    p2_intersect_line_rep = p2_intersect_lines[0]
+                elif p2_intersect_lines.__len__()==1:
+                    p2_intersect_line_rep = p2_intersect_lines[0]
+                else:
+                    continue
+
+
+                # #两条线没有距离跳过
+                # if p2_intersect_line_rep.Distance(p1_intersect_line_rep) == 0:
+                #     continue
+                # else:
+
+                mid_polygon = self.create_mid_polygon(p1_intersect_line_rep,p2_intersect_line_rep,p1,p2)
+
+                if mid_polygon.GetArea() > 0.0002116361000058077*0.0002116361000058077:
+                    continue
+
+                merge_polygon_raw.append(mid_polygon)
+
+            merge_polygon = None
+            for p in merge_polygon_raw:
+                if merge_polygon is None:
+                    merge_polygon = p
+                else:
+                    merge_polygon = merge_polygon.Union(p)
+
+            polygons.append(merge_polygon)
+
+
+        return polygons
+
+
+    def merge2(self,trees):
+
+        polygons = []
+
+
+
+
+        for tree in trees:
+
+            tree_polygons = tree.get_polygons()
+
+            points = []
+            for p in tree_polygons:
+
+                if p.GetGeometryCount() == 1:
+                    polygon_line: Geometry = ogr.ForceToLineString(p)
+
+                    polygon_line.Con
+                # 有孔Polygon
+                else:
+                    polygon_line: Geometry = ogr.ForceToLineString(p.GetGeometryRef(0))
+
+                points.extend(polygon_line.GetPoints()[:-1])
+
+            multipoint: Geometry = ogr.Geometry(ogr.wkbMultiPoint)
+
+            for p in  points:
+                poi = ogr.Geometry(ogr.wkbPoint)
+                poi.AddPoint(p[0],p[1])
+                multipoint.AddGeometry(poi)
+
+            delauney: Geometry = multipoint.DelaunayTriangulation()
+            merge_polygon = None
+
+            for p in delauney:
+                # if multi_polygon.Contains(p):
+                #     continue
+
+                pline:Geometry = ogr.ForceToLineString(p)
+                pline :list = pline.GetPoints()
+                pline.append(pline[1])
+
+
+                pass_or_not = False
+                for indx in range(3):
+                    p1 = pline[indx]
+                    p2 = pline[indx+1]
+
+                    ppline:Geometry = ogr.Geometry(ogr.wkbLineString)
+                    ppline.AddPoint(p1[0],p1[1])
+                    ppline.AddPoint(p2[0], p2[1])
+
+
+
+                    p3:Geometry = ogr.Geometry(ogr.wkbPoint)
+                    p3.AddPoint(pline[indx+2][0], pline[indx+2][1])
+
+                    if p3.Distance(ppline) > 0.0002116361000058077 or ppline.Length() > 0.0002116361000058077:
+                        pass_or_not = True
+                        break
+
+                if pass_or_not:
+                    continue
+
+                if merge_polygon is None:
+                    merge_polygon = p
+                else:
+                    merge_polygon = merge_polygon.Union(p)
+
+
+
+            polygons.append(merge_polygon)
+
+
+        return polygons
+
+    def merge_circle(self):
+        pass
+
+
+    def create_mid_polygon(self,p1_intersect_line_rep,p2_intersect_line_rep,p1,p2):
+        '''
+        创建衔接polygon
+        :param line:
+        :return:
+        '''
+
+        if p1_intersect_line_rep.Intersect(p2_intersect_line_rep):
+
+            p1_rep_p1 = p1_intersect_line_rep.GetPoints()[0]
+            p1_rep_p2 = p1_intersect_line_rep.GetPoints()[1]
+
+            p2_rep_p1 = p2_intersect_line_rep.GetPoints()[0]
+            p2_rep_p2 = p2_intersect_line_rep.GetPoints()[1]
+
+            polygon = self.create_polygon(p1_rep_p1, p2_rep_p2, p1_rep_p2, p2_rep_p1)
+
+        else:
+
+            p1_rep_p1 = p1_intersect_line_rep.GetPoints()[0]
+            p1_rep_p2 = p1_intersect_line_rep.GetPoints()[1]
+
+            p2_rep_p1 = p2_intersect_line_rep.GetPoints()[0]
+            p2_rep_p2 = p2_intersect_line_rep.GetPoints()[1]
+
+            line1 = ogr.Geometry(ogr.wkbLineString)
+
+
+            line1.AddPoint(p1_rep_p1[0], p1_rep_p1[1])
+            line1.AddPoint(p2_rep_p1[0], p2_rep_p1[1])
+
+            line2 = ogr.Geometry(ogr.wkbLineString)
+            line2.AddPoint(p1_rep_p2[0], p1_rep_p2[1])
+            line2.AddPoint(p2_rep_p2[0], p2_rep_p2[1])
+
+            if line1.Intersect(line2):
+                polygon = self.create_polygon(p1_rep_p1, p1_rep_p2, p2_rep_p1, p2_rep_p2)
+            else:
+
+                polygon = self.create_polygon(p1_rep_p1, p1_rep_p2, p2_rep_p2, p2_rep_p1)
+                if polygon.Intersect(p1) or polygon.Intersect(p2):
+                    polygon = self.create_polygon(p1_rep_p1, p1_rep_p2, p2_rep_p1, p2_rep_p2)
+
+
+        return polygon
+
+    def create_polygon(self,p1,p2,p3,p4):
+
+        ring = ogr.Geometry(ogr.wkbLinearRing)
+        ring.AddPoint(p1[0], p1[1])
+        ring.AddPoint(p2[0], p2[1])
+        ring.AddPoint(p3[0], p3[1])
+        ring.AddPoint(p4[0], p4[1])
+        ring.AddPoint(p1[0], p1[1])
+
+        poly = ogr.Geometry(ogr.wkbPolygon)
+        poly.AddGeometry(ring)
+        return poly
+
+
+    def get_polygon_lines(self,polygon):
+        #无孔Polygon
+        if polygon.GetGeometryCount() == 1:
+            polygon_line : Geometry = ogr.ForceToLineString(polygon)
+        # 有孔Polygon
+        else:
+            polygon_line : Geometry = ogr.ForceToLineString(polygon.GetGeometryRef(0))
+
+        points = polygon_line.GetPoints()
+
+        lines = []
+        for index in range(len(points)-1):
+
+            point_start = points[index]
+            point_end = points[index+1]
+            line = ogr.Geometry(ogr.wkbLineString)
+            line.AddPoint(point_start[0], point_start[1])
+            line.AddPoint(point_end[0], point_end[1])
+
+            lines.append(line)
+
+        return lines
+
+    def merge_concave(self, trees, distance_threshold):
+        merge_polygons = []
+        for tree in trees:
+            polygons = list(tree.get_polygons())
+
+            merge_raw: Geometry = polygons[0]
+            for p in polygons:
+                merge_raw = merge_raw.Union(p)
+
+            points = []
+            dd= merge_raw.GetGeometryCount()
+
+            for i in range(merge_raw.GetGeometryCount()):
+
+                poly = merge_raw.GetGeometryRef(i)
+                if poly.GetGeometryCount() == 1:
+                    poly_line: Geometry = ogr.ForceToLineString(poly)
+                # 有孔Polygon
+                else:
+                    poly_line: Geometry = ogr.ForceToLineString(poly.GetGeometryRef(0))
+                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])
+
+            for pp in points:
+                poi = ogr.Geometry(ogr.wkbPoint)
+                poi.AddPoint(pp[0], pp[1])
+
+            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_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)
+
+        delauney_lines = self.create_delauney(rps)
+
+        min_delauney = self.create_min_delauney(delauney_lines,rps,polygons)
+
+        trees = self.cut_delauney_min_tree(min_delauney,self.distance_buffer_threshold)
+
+        print(len(trees))
+        polygons = self.merge_concave(trees,0.00007)
+
+        return [p.ExportToWkt() for p in polygons]
+
+if __name__ == '__main__':
+
+
+    # print(sd.get_polygons()[0])
+    # sd.close()
+
+
+    sd = ShapeData(r"J:\Data\制图综合\example.shp")
+
+    ext =  sd.layer.GetExtent()
+    threshold = (ext[1]-ext[0])/20
+    print(threshold)
+
+
+    zh = Zonghe(sd)
+    wkts = zh.process()
+
+    result = r"J:\Data\制图综合\zhonghe_concave2.shp"
+    ShapeData.create_shp_fromwkts(result,"zh",wkts)
+