CreateData.py
1.4 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
# 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
layer_name = "random"
driver:Driver = ogr.GetDriverByName("ESRI Shapefile")
data_source: DataSource = driver.CreateDataSource(r"E:/Data/bigdata/{}.shp".format(layer_name))
outLayer = data_source.CreateLayer(layer_name, geom_type=ogr.wkbPoint )
t_point = time.time()
#打开广东边界
ds: DataSource = driver.Open("E:\Data\广东边界\gdsample.shp", 1)
if not ds:
raise Exception("打开数据失败！")
layer: Layer = ds.GetLayer(0)
feat :Feature = layer.GetFeature(0)
geo:Geometry=feat.GetGeometryRef()
extent = geo.GetEnvelope()
i=0
while i<100000:
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):
i+=1
# print(i)
featureDefn = outLayer.GetLayerDefn()
outFeature = ogr.Feature(featureDefn)
outFeature.SetGeometry(point)
outLayer.CreateFeature(outFeature)
if i%10000==0:
print(i)
print(time.time() - t_point)
t_point=time.time()
outFeature=None