2022.03.11 合并备份

-
-from osgeo import ogr
-
-def envelop_2_polygon(env):
-    ring = ogr.Geometry(ogr.wkbLinearRing)
-    ring.AddPoint(env[0], env[2])
-    ring.AddPoint(env[0], env[3])
-    ring.AddPoint(env[1], env[3])
-    ring.AddPoint(env[1], env[2])
-    ring.AddPoint(env[0], env[2])
-    # Create polygon
-    poly = ogr.Geometry(ogr.wkbPolygon)
-    poly.AddGeometry(ring)
-    return poly
-
-
-env = [1,3,1,3]
-
-poly = envelop_2_polygon(env)
-line =  ogr.Geometry(ogr.wkbLineString)
-line.AddPoint(0, 2)
-line.AddPoint(4, 2)
-line.AddPoint(5, 2)
-
-print(poly)
-print(line)
-print(line.Simplify(1))

+# coding=utf-8
+#author:        4N
+#createtime:    2022/3/8
+#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
+os.environ["CUDA_VISIBLE_DEVICES"] = "-1"
+class util:
+    @classmethod
+    def envelop_2_polygon(cls,env):
+        ring = ogr.Geometry(ogr.wkbLinearRing)
+        ring.AddPoint(env[0], env[2])
+        ring.AddPoint(env[0], env[3])
+        ring.AddPoint(env[1], env[3])
+        ring.AddPoint(env[1], env[2])
+        ring.AddPoint(env[0], env[2])
+        # Create polygon
+        poly = ogr.Geometry(ogr.wkbPolygon)
+        poly.AddGeometry(ring)
+        return poly
+
+
+class ShapeUtil:
+
+    driver = ogr.GetDriverByName("ESRI Shapefile")
+
+    @classmethod
+    def create_by_layer(cls,path,layer:Layer):
+        data_source: DataSource = cls.driver.CreateDataSource(path)
+        data_source.CopyLayer(layer,layer.GetName())
+        data_source.Destroy()
+
+    @classmethod
+    def create_by_scheme(cls,path,name,sr,geo_type,scheme,features):
+        data_source: DataSource = cls.driver.CreateDataSource(path)
+        layer :Layer = data_source.CreateLayer(name, sr, geo_type)
+        if scheme:
+            layer.CreateFields(scheme)
+        for feature in features:
+            layer.CreateFeature(feature)
+        data_source.Destroy()
+
+    @classmethod
+    def create_point(cls,path,name,point):
+        data_source: DataSource = cls.driver.CreateDataSource(path)
+        layer :Layer = data_source.CreateLayer(name, None, ogr.wkbPoint)
+
+        feat_new = ogr.Feature(layer.GetLayerDefn())
+        feat_new.SetGeometry(point)
+        layer.CreateFeature(feat_new)
+        data_source.Destroy()
+
+
+class ShapeData:
+
+    def __init__(self,path):
+        driver: Driver = ogr.GetDriverByName("ESRI Shapefile")
+        self.ds: DataSource = driver.Open(path, 1)
+        if not self.ds:
+            raise Exception("打开数据失败！")
+        self.layer: Layer = self.ds.GetLayer(0)
+
+
+    def get_polygons(self):
+
+        polygons = []
+        for feature in self.layer:
+            f:Feature = feature
+            geom : Geometry = f.GetGeometryRef()
+            if geom.GetGeometryType() == 3 or geom.GetGeometryType() == -2147483645:
+                polygons.append(geom)
+            if geom.GetGeometryType() == 6 or geom.GetGeometryType() == -2147483642:
+                for i in range(geom.GetGeometryCount()):
+                    polygons.append(geom.GetGeometryRef(i))
+        return polygons
+
+    def close(self):
+        self.ds.Destroy()
+
+
+class Zonghe:
+
+    data:ShapeData
+    area_threshold = 10.0
+    buffer_threshold = 100.0
+
+    def __init__(self,data):
+        self.data = data
+
+
+
+    #只处理Polygon，MultiPolygon要拆开
+    def get_polygon_reprecent_point(self,polygon:Geometry):
+
+        polygon_area = polygon.GetArea()
+
+        polygon_env = polygon.GetEnvelope()
+
+        line_length = max(polygon_env[1]-polygon_env[0],polygon_env[3]-polygon_env[2])
+
+        center:Geometry = polygon.Centroid()
+        if polygon.Contains(center):
+            return center
+
+
+        left_area = 0
+        right_area = polygon_area
+
+        left_or_right = None
+        fen = False
+
+        if polygon.GetGeometryCount() == 1:
+            polygon_line = ogr.ForceToLineString(polygon)
+        # 有孔Polygon
+        else:
+            polygon_line = ogr.ForceToLineString(polygon.GetGeometryRef(0))
+
+        # print(len(polygon_line.GetPoints()))
+        point_i = 0
+
+        line = None
+        next_point = polygon_line.GetPoint(point_i)
+        last_point = None
+        last_last_point = None
+
+        angle = []
+
+        calculate_times = 0
+
+        while abs(left_area - right_area) > (polygon_area / self.area_threshold):
+            # print(calculate_times)
+            calculate_times += 1
+
+            line,other_point = self.get_line(center.GetX(),center.GetY(),next_point[0],next_point[1],line_length)
+
+            clip_geom:Geometry = polygon.Difference(line.Buffer(line_length/self.buffer_threshold))
+
+
+            if clip_geom.GetGeometryType() == 3 or clip_geom.GetGeometryType() == -2147483645:
+
+                cl = ogr.ForceToLineString(clip_geom)
+                p1 = cl.GetPoint(0)
+                p2 = cl.GetPoint(1)
+                tp = [(p1[0] + p2[0])/2.0, (p1[1] + p2[1])/2.0]
+                # clip_geom_center = clip_geom.Centroid()
+
+
+                side_value = self.side(center,other_point,self.creat_point(tp))
+                if side_value > 0:
+                    left_area = clip_geom.GetArea()
+                    right_area = 0
+                else:
+                    right_area = clip_geom.GetArea()
+                    left_area = 0
+
+
+
+            else:
+
+                left_area = 0
+                right_area = 0
+
+                for gi in range(clip_geom.GetGeometryCount()):
+                    clip_geom_i :Geometry = clip_geom.GetGeometryRef(gi)
+                    clip_geom_i_center = clip_geom_i.Centroid()
+                    side_value = self.side(center, other_point, clip_geom_i_center)
+
+                    if side_value > 0:
+                        left_area += clip_geom_i.GetArea()
+                    else:
+                        right_area  += clip_geom_i.GetArea()
+
+
+            print(left_area,right_area)
+
+            last_last_point = last_point
+            last_point = next_point
+            left_or_right = 1 if left_area > right_area else -1
+
+
+            # 进入二分法步骤
+            if fen:
+
+                alternate:Geometry = ogr.Geometry(ogr.wkbPoint)
+                alternate.AddPoint(angle[0][0],angle[0][1])
+                # print(angle)
+                alternate_side = self.side(center, other_point, alternate)
+                print(alternate_side)
+                if left_or_right * alternate_side > 0:
+                    angle = [angle[0],next_point]
+                    next_point = ((next_point[0]+angle[0][0])/2.0,(next_point[1]+angle[0][1])/2.0)
+
+                    # print(abs(next_point[0]-angle[0][0]))
+                else:
+                    angle = [angle[1],next_point]
+                    next_point = ((next_point[0]+angle[1][0])/2.0,(next_point[1]+angle[1][1])/2.0)
+
+                    # print(abs(next_point[0] - angle[1][0]))
+
+            else:
+                # print(point_i)
+                point_i += 1
+                next_point = polygon_line.GetPoint(point_i)
+
+                next_point_side = self.side(center, other_point, self.creat_point(next_point))
+
+                if left_or_right * next_point_side > 0 :
+                    pass
+                else:
+                    next_point = ((last_point[0]+last_last_point[0])/2.0,(last_point[1]+last_last_point[1])/2.0)
+                    angle = [last_last_point,last_point]
+                    fen = True
+
+        print(calculate_times)
+
+        point_intersection:Geometry = line.Intersection(polygon_line)
+
+        xlist = []
+        ylist = []
+        for i in range(point_intersection.GetGeometryCount()):
+            pi = point_intersection.GetGeometryRef(i)
+            xlist.append(pi.GetX())
+            ylist.append(pi.GetY())
+
+        reprecent_point = [(min(xlist)+max(xlist))/2.0,(min(ylist)+max(ylist))/2.0]
+
+        return self.creat_point(reprecent_point)
+
+
+
+    def get_line(self,x0,y0,x1,y1,len_size):
+
+        line: Geometry = ogr.Geometry(ogr.wkbLineString)
+        line.AddPoint(x0, y0)
+        x2 = x0 + len_size
+
+        if x1 != x0:
+            y2 = ((y1-y0)*(x2-x0))/(x1-x0) + y0
+        else:
+            x2 = x0
+            y2 = y0 + len_size
+        line.AddPoint(x2, y2)
+
+        other_point = ogr.Geometry(ogr.wkbPoint)
+        other_point.AddPoint(x2, y2)
+
+        return line, other_point
+    def creat_point(self,tuple):
+        p: Geometry = ogr.Geometry(ogr.wkbPoint)
+        p.AddPoint(tuple[0], tuple[1])
+        return p
+
+
+
+    def side(self,A,B,P):
+
+
+        side = (B.GetX() - A.GetX()) * (P.GetY() - A.GetY()) - (B.GetY() - A.GetY()) * (P.GetX() - A.GetX())
+        return side
+
+
+
+
+    def process(self):
+        polygons = self.data.get_polygons()
+        for poly in polygons:
+            reprecent_point = self.get_polygon_reprecent_point(poly)
+            print(reprecent_point)
+
+
+
+if __name__ == '__main__':
+
+    # sd = ShapeData(r"J:\Data\矢量数据\佛山\制图综合.shp")
+    # print(sd.get_polygons()[0])
+    # sd.close()
+
+    driver: Driver = ogr.GetDriverByName("ESRI Shapefile")
+    ds: DataSource = driver.Open(r"J:\Data\矢量数据\佛山\制图综合\t2.shp",0)
+
+    layer: Layer = ds.GetLayer(0)
+
+    polygons = []
+    for feature in layer:
+        f:Feature = feature
+        geom : Geometry = f.GetGeometryRef()
+        if geom.GetGeometryType() == 3 or geom.GetGeometryType() == -2147483645:
+            polygons.append(geom)
+        if geom.GetGeometryType() == 6 or geom.GetGeometryType() == -2147483642:
+            for i in range(geom.GetGeometryCount()):
+                polygons.append(geom.GetGeometryRef(i))
+    del ds
+    zh = Zonghe("d")
+    pre = zh.get_polygon_reprecent_point(polygons[0])
+    print(pre)
+
+    result="J:\Data\矢量数据\佛山\制图综合\zhresult.shp"
+    print(pre)
+
+    ShapeUtil.create_point(result,"zh",pre)
+