script for generating a discrete number of random points within each of a set of polygons

points-within-polys.py

# assign a random location to every person in greater cinci
# inside their home census block

import psycopg2, random
from shapely import wkb
from shapely.geometry import Point

# connect to the census DB
conn_string = ("host='localhost' dbname='' user='' password=''")
census_conn = psycopg2.connect(conn_string)
census_cursor = census_conn.cursor()

# output
out = open('peeps.csv','w')
out.write('lon,lat\n')

peeps_done = 0

# get the population data
census_cursor.execute("""
	SELECT 
		total_population,
		ST_Transform(the_geom,4326)
	FROM blocks_2010
	WHERE total_population > 0;
""")
# for each block
for (tract_pop,geom) in census_cursor:
	geom = wkb.loads(geom,hex=True)
	minX,minY,maxX,maxY = geom.bounds
	rangeX,rangeY = (maxX-minX),(maxY-minY)
	# for each person in each block
	for person in range(0,tract_pop):
		x = minX + random.random() * rangeX
		y = minY + random.random() * rangeY
		p = Point(x,y)
		while not geom.contains(p):
			x = minX + random.random() * rangeX
			y = minY + random.random() * rangeY
			p = Point(x,y)
		# write the new point
		out.write(str(x)+','+str(y)+'\n')
		peeps_done += 1
		if peeps_done % 1000 == 0:
			print peeps_done