geospatialpython/DensityProfiler_pure.py at master · GeospatialPython/geospatialpython · GitHub

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
importrandom
importshapefile
importpngcanvas

defpip(x,y,poly):
n=len(poly)
inside=False
p1x,p1y=poly[0]
foriinrange(n+1):
p2x,p2y=poly[i%n]
ify>min(p1y,p2y):
ify<=max(p1y,p2y):
ifx<=max(p1x,p2x):
ifp1y!=p2y:
xinters= (y-p1y)*(p2x-p1x)/(p2y-p1y)+p1x
ifp1x==p2xorx<=xinters:
inside=notinside
p1x,p1y=p2x,p2y
returninside

# Source shapefile
r=shapefile.Reader("GIS_CensusTract_poly.shp")

# pixel to coordinate info
xdist=r.bbox[2] -r.bbox[0]
ydist=r.bbox[3] -r.bbox[1]
iwidth=600
iheight=500
xratio=iwidth/xdist
yratio=iheight/ydist

c=pngcanvas.PNGCanvas(iwidth,iheight,color=[255,255,255,0xff])

# background color
c.filledRectangle(0,0,iwidth,iheight)

# Pen color
c.color= [139,137,137,0xff]

# Draw the polygons
forshapeinr.shapes():
pixels= []
forx,yinshape.points:
px=int(iwidth- ((r.bbox[2] -x) *xratio))
py=int((r.bbox[3] -y) *yratio)
pixels.append([px,py])
c.polyline(pixels)

rnum=0
trnum=len(r.shapeRecords())
forsrinr.shapeRecords():
rnum+=1
#print rnum, " of ", trnum
density=sr.record[20]
total=int(density/50)
count=0
minx, miny, maxx, maxy=sr.shape.bbox
whilecount<total:
x=random.uniform(minx,maxx)
y=random.uniform(miny,maxy)
ifpip(x,y,sr.shape.points):
count+=1
#print " ", count, " of ", total
px=int(iwidth- ((r.bbox[2] -x) *xratio))
py=int((r.bbox[3] -y) *yratio)
c.point(px,py,color=[255,0,0,0xff])

f=file("density_pure.png", "wb")
f.write(c.dump())
f.close()