demographics-point-builder.py

import numpy as np
import json
import os

def bbox(coords):
    sep = zip(*coords)
    Xmax = max(sep[0])
    Xmin = min(sep[0])
    Ymax = max(sep[1])
    Ymin = min(sep[1])
    return [Xmax,Ymax,Xmin,Ymin]

def randPoint(coords):
    bb = bbox(coords)
    randX = (np.random.rand() * (bb[0] - bb[2])) + bb[2]
    randY = (np.random.rand() * (bb[1] - bb[3])) + bb[3]
    return [randX, randY]

def point_in_poly(x,y,poly):
    n = len(poly)
    inside = False
    p1x,p1y = poly[0]
    for i in range(n+1):
        p2x,p2y = poly[i % n]
        if y > min(p1y,p2y):
            if y <= max(p1y,p2y):
                if x <= max(p1x,p2x):
                    if p1y != p2y:
                        xints = (y-p1y)*(p2x-p1x)/(p2y-p1y)+p1x
                    if p1x == p2x or x <= xints:
                        inside = not inside
        p1x,p1y = p2x,p2y
    return inside

def fillPoly(coords,num,name):
    return_points = []
    success = 0
    tries = 0
    while success < num:
        test_point = randPoint(coords)
        if point_in_poly(test_point[0],test_point[1],coords):
            test_point.append(name)
            return_points.append(test_point)
            success += 1
        tries += 1
    return return_points

def makePoints(polygons, features):
    # assuming features is a multi-dim array in the form [[feature_vector_one],[feature_vector_two]]
    accepted_points = []
    for i in range(0,len(polygons)):
        poly = polygons[i]
        for j in range(0,len(features)):
            feature = features[j]
            points = fillPoly(poly,int(feature['value'][i]),feature['name'])
            accepted_points.append(points)
    return accepted_points
    
for filename in os.listdir('geojson'):
    with open('./geojson/' + filename) as f:
        data = json.load(f)
        
        poly_coords = []
        population = []
        renters = []
        owners = []
        for feature in data['features']:
            poly_coords.append(feature['geometry']['coordinates'][0][0])
            population.append(feature['properties']['POPULATION'])
            owners.append(feature['properties']['O'])
            renters.append(feature['properties']['R'])
            
        feature_matrix = [{'name':'O','value':owners},{'name':'R','value':renters}]
        people_points = makePoints(poly_coords,feature_matrix)
        
        geodict = {'type':'FeatureCollection','features':[]}
        for block in people_points:
            for point in block:
                geopoint = {'type':'Feature','geometry':{'type':'Point','coordinates':[point[0],point[1]]},'properties':{'ro':point[2]}}
                geodict['features'].append(geopoint)
        with open('./output/' + filename,'w') as geofile:
            json.dump(geodict, geofile)
        print 'Finished: ' + filename