In [19]:

    

from __future__ import division
import fiona
from shapely.geometry import shape, Point, Polygon
from descartes import PolygonPatch
from matplotlib import pylab as plt
import random as r

%matplotlib notebook


    

In [5]:

    

districts = {}
with fiona.open('data/point-in-polygon/Maryland_108_to_112.geojson', 'r') as source:
    for f in source:
        districts[f['properties']['district']] = shape(f['geometry'])


    

In [6]:

    

fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(111)
for d in districts:
    color = 'grey'
    alpha=0.5
    if d == '3':
        color='red'
        alpha=0.9
    for p in districts[d]:
        patch = PolygonPatch(p, facecolor=color, alpha=alpha)
        ax.add_patch(patch)

ax.autoscale_view(True, True, True)
ax.axis('off')
plt.show()


    

In [20]:

    

def generate_point_within_polygon(polygon):
    (minx, miny, maxx, maxy) = polygon.bounds
    i = 0
    while True:
        i+=1
        point = Point(r.uniform(minx, maxx), r.uniform(miny, maxy))
        if polygon.contains(point):
            return i, point


    

In [22]:

    

%time points = [generate_point_within_polygon(districts['3']) for _ in range(5000)]
print sum([p[0] for p in points])/len(points)


    

    




CPU times: user 20.6 s, sys: 0 ns, total: 20.6 s
Wall time: 20.6 s
3.5404

In [8]:

    

from meshpy.triangle import MeshInfo, build


    

In [9]:

    

points = [list(p.exterior.coords) for p in districts['3']]


    

In [10]:

    

triangles = []
for facet in points:
    range_points = range(len(facet))
    mesh_info = MeshInfo()
    mesh_info.set_points(facet)
    mesh_info.set_facets([(a,b) for a, b in zip(range_points[:-1], range_points[1:])])
    mesh = build(mesh_info)

    points_m = {}
    for i, p in enumerate(mesh.points):
        points_m[i] = p
    print "Point numbers in tetrahedra:"
    facets = []
    for i, t in enumerate(mesh.elements):
        facets.append([points_m[j] for j in t])

    ## Create triangles here
    triangles += [Polygon(f) for f in facets]

area_triangles = [t.area for t in triangles]


    

    




Point numbers in tetrahedra:
Point numbers in tetrahedra:

In [11]:

    

fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(111)
for d in triangles:
    patch = PolygonPatch(d, facecolor=color, alpha=alpha)
    ax.add_patch(patch)

ax.axis('off')
ax.autoscale_view(True, True, True)
plt.show()


    

In [12]:

    

import bisect
from math import sqrt

## Choose the triangle at random
cum = []
cumsum = 0
for m in area_triangles:
    cumsum += m
    cum.append(cumsum/sum(area_triangles))

points = []
for _ in range(100000):
    a = r.random()
    num =  bisect.bisect_left(cum, a)
    triangle = triangles[num]
    
    # Find random point in triangle
    a = r.random()
    b = r.random()
    
    coords = triangle.exterior.coords
    point = ((1-sqrt(a))*coords[0][0] + sqrt(a)*(1-b)*coords[1][0]+sqrt(a)*b*coords[2][0],
            (1-sqrt(a))*coords[0][1] + sqrt(a)*(1-b)*coords[1][1]+sqrt(a)*b*coords[2][1])
    points.append(point)


    

In [13]:

    

## Draw point at random  in a given triangle


    

In [14]:

    

fig  = plt.figure(figsize=(12,8))
ax = fig.add_subplot(111)
ax.scatter([p[0] for p in points], [p[1] for p in points], s=.01)
ax.axis('off')
plt.show()


    

In [ ]: