Voronoi.py
2.8 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
# coding=utf-8
#author: 4N
#createtime: 2021/6/9
#email: nheweijun@sina.com
from osgeo.ogr import *
from osgeo.gdal import *
#
from osgeo import ogr
import random
import time
from scipy.spatial import Voronoi,voronoi_plot_2d
import numpy as np
import random
import os
import time
print(os.path.abspath(__file__))
begin_time = time.time()
create_feature=100000
layer_name="voronoixian"
work_dir = r"E:/Data/bigdata"
# work_dir =os.path.dirname(os.path.abspath(__file__))
driver:Driver = ogr.GetDriverByName("ESRI Shapefile")
#打开广东边界
ds: DataSource = driver.Open("E:\Data\广东边界\gdxian.shp", 1)
if not ds:
raise Exception("打开数据失败！")
layer: Layer = ds.GetLayer(0)
each_piece_features = create_feature*1.0/layer.GetFeatureCount()
data_source=None
outLayer=None
for feat in layer:
geo:Geometry=feat.GetGeometryRef()
extent = geo.GetEnvelope()
points = []
i=0
while i<each_piece_features:
point = ogr.Geometry(ogr.wkbPoint)
x = (extent[1]-extent[0])*random.random()+extent[0]
y = (extent[3]-extent[2])*random.random()+extent[2]
point.AddPoint_2D(x,y)
# print(geo.Intersect(point))
if geo.Intersect(point):
points.append([x,y])
i+=1
array = np.array(points)
# 沃罗诺伊图
vor = Voronoi(array,furthest_site=False, incremental=True, qhull_options=None)
# 泰森多边形的顶点
vertices = vor.vertices
regions = vor.regions
if not data_source:
data_source: DataSource = driver.CreateDataSource(os.path.join(work_dir,"{}.shp".format(layer_name)))
if not outLayer:
outLayer = data_source.CreateLayer(layer_name, geom_type=ogr.wkbPolygon )
print(time.time()-begin_time)
print("here")
vonoi_time = time.time()
for index,r in enumerate(regions):
if len(r) == 0:
continue
if -1 in r:
continue
angulars = []
for id in r:
angulars.append(vertices[id])
angulars.append(vertices[r[0]])
ring = ogr.Geometry(ogr.wkbLinearRing)
for point in angulars:
ring.AddPoint(point[0],point[1])
# Create polygon
poly = ogr.Geometry(ogr.wkbPolygon)
poly.AddGeometry(ring)
poly:Geometry= geo.Intersection(poly)
if poly:
if poly.IsEmpty():
continue
featureDefn = outLayer.GetLayerDefn()
outFeature = ogr.Feature(featureDefn)
outFeature.SetGeometry(poly)
outLayer.CreateFeature(outFeature)
outFeature = None
print(time.time()-vonoi_time)
data_source.Destroy()
ds.Destroy()