-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathwork_stab.py
107 lines (90 loc) · 3.05 KB
/
work_stab.py
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
import pysal
from pandas import *
import fiona
from shapely.geometry import Point
import gdal, ogr, os, osr
from point_fre import point_fre
from polyraster import polyraster
from union import union
from polysplit import polysplit
#import rpy2.robjects as robjects
#from rpy2.robjects.packages import importr
#import rpy2.robjects.numpy2ri
#rpy2.robjects.numpy2ri.activate()
# Set up our R namespaces
#R = rpy2.robjects.r #run R in python namespaces
#spatstat = importr('spatstat') #import R packages
iteration=range(1001)
setwd="D:/crime geocoding/SOM/stability/"
db= pysal.open('D:/crime geocoding/SOM/tulsa_crime0911.dbf', 'r')
d = {col: db.by_col(col) for col in db.header}
table1=DataFrame(d)
polycount=[]
for it in iteration[1:21]:
#point frequency in raster
fre=setwd+'quadrat'+str(it)+".tif"
cenpoint=setwd+'pointcen'+str(it)+".shp"
aggpoint=setwd+'aggpoint'+str(it-1)+".shp"
point_fre(repeat=1, gridsize=20, table1=table1,outfile=fre, pointfile=cenpoint)
"""
#aggregate points
if it == 1:
mypoints=[]
for gf in fiona.open(cenpoint):
mypoints.append(Point(gf['geometry']['coordinates']))
crime=[]
for gf in fiona.open('D:/crime geocoding/SOM/tulsa_crime0911.shp'):
crime.append(Point(gf['geometry']['coordinates']))
mypoints=crime+mypoints
elif it > 1:
points=[]
for gf in fiona.open(cenpoint):
points.append(Point(gf['geometry']['coordinates']))
mypoints=points+mypoints #combine
spatialReference = osr.SpatialReference()
spatialReference.ImportFromProj4('+proj=utm +zone=15 +ellps=WGS84 +datum=WGS84 +units=m +no_defs ')
driver = ogr.GetDriverByName('Esri Shapefile')
ds = driver.CreateDataSource(aggpoint)
layer = ds.CreateLayer('', spatialReference, geom_type=ogr.wkbPoint)
# Add one attribute
layer.CreateField(ogr.FieldDefn('id', ogr.OFTInteger))
defn = layer.GetLayerDefn()
for i in mypoints:
feat = ogr.Feature(defn)
feat.SetField('id', 0)
geom = ogr.CreateGeometryFromWkb(i.wkb)
feat.SetGeometry(geom)
layer.CreateFeature(feat)
feat = geom = None # destroy these
"""
#polygonize
#inin=setwd+'quadrat'+str(it)+".tif"
outout=setwd+'quadrat'+str(it)+".shp"
polyraster(infile=fre, outfile=outout)
#polygon union and polygon split
if it == 1:
in_cade=outout
mypoly=[]
for count in fiona.open(outout):
mypoly.append(count['geometry'])
polycount.append(len(mypoly))
elif it > 1 and it == 2:
in_cade=[in_cade]+[outout]
out_cade=setwd+'cascade_union'+str(it)+".shp"
out_split=setwd+'cascade_split'+str(it)+".shp"
union(infile=in_cade, outfile=out_cade)
polysplit(infile=out_cade, outfile=out_split)
mypoly=[]
for count in fiona.open(out_split):
mypoly.append(count['geometry'])
polycount.append(len(mypoly))
elif it > 2:
in_cade=[out_split]+[outout]
union(infile=in_cade, outfile=setwd+'cascade_union'+str(it)+".shp")
out_cade=setwd+'cascade_union'+str(it)+".shp"
out_split=setwd+'cascade_split'+str(it)+".shp"
polysplit(infile=out_cade, outfile=out_split)
mypoly=[]
for count in fiona.open(out_split):
mypoly.append(count['geometry'])
polycount.append(len(mypoly))