In [2]:

    

import numpy as np
from scipy.spatial import Voronoi, voronoi_plot_2d
import shapely
import geopandas as gpd
import sys
import matplotlib.pyplot as plt
import matplotlib as mpl
%matplotlib inline

mpl.rcParams['figure.figsize'] = 15,15  #set the default map size
mpl.rcParams['patch.linewidth'] = 0.5  #set default polygon line width

np.random.seed(352)

# Local path on user's machine
path = '/Users/jgaboardi/Algorithms_etc./Data/'


    

    




/Users/jgaboardi/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/matplotlib/font_manager.py:273: UserWarning: Matplotlib is building the font cache using fc-list. This may take a moment.
  warnings.warn('Matplotlib is building the font cache using fc-list. This may take a moment.')

In [3]:

    

def in_box(households, bounding_box):
    return np.logical_and(np.logical_and(bounding_box[0] <= households[:, 0],
                                         households[:, 0] <= bounding_box[1]),
                          np.logical_and(bounding_box[2] <= households[:, 1],
                                         households[:, 1] <= bounding_box[3]))


    

In [16]:

    

def voronoi(households, bounding_box):
    # Select households inside the bounding box
    i = in_box(households, bounding_box)
    print i
    # Mirror points
    points_center = households[i, :]
    print points_center
    points_left = np.copy(points_center)
    points_left[:, 0] = bounding_box[0] - (points_left[:, 0] - bounding_box[0])
    points_right = np.copy(points_center)
    points_right[:, 0] = bounding_box[1] + (bounding_box[1] - points_right[:, 0])
    points_down = np.copy(points_center)
    points_down[:, 1] = bounding_box[2] - (points_down[:, 1] - bounding_box[2])
    points_up = np.copy(points_center)
    points_up[:, 1] = bounding_box[3] + (bounding_box[3] - points_up[:, 1])
    points = np.append(points_center,
                       np.append(np.append(points_left,
                                           points_right,
                                           axis=0),
                                 np.append(points_down,
                                           points_up,
                                           axis=0),
                                 axis=0),
                       axis=0)
    # Compute Voronoi
    vor = Voronoi(points)
    # Filter regions
    regions = []
    for region in vor.regions:
        flag = True
        for index in region:
            if index == -1:
                flag = False
                break
            else:
                x = vor.vertices[index, 0]
                y = vor.vertices[index, 1]
                if not(bounding_box[0] - eps <= x and x <= bounding_box[1] + eps and
                       bounding_box[2] - eps <= y and y <= bounding_box[3] + eps):
                    flag = False
                    break
        if region != [] and flag:
            regions.append(region)
    vor.filtered_points = points_center
    vor.filtered_regions = regions
    return vor


    

In [17]:

    

def centroid_region(vertices):
    # Polygon's signed area
    A = 0
    # Centroid's x
    C_x = 0
    # Centroid's y
    C_y = 0
    for i in range(0, len(vertices) - 1):
        s = (vertices[i, 0] * vertices[i + 1, 1] - vertices[i + 1, 0] * vertices[i, 1])
        A = A + s
        C_x = C_x + (vertices[i, 0] + vertices[i + 1, 0]) * s
        C_y = C_y + (vertices[i, 1] + vertices[i + 1, 1]) * s
    A = 0.5 * A
    C_x = (1.0 / (6.0 * A)) * C_x
    C_y = (1.0 / (6.0 * A)) * C_y
    return np.array([[C_x, C_y]])


    

In [18]:

    

eps = sys.float_info.epsilon

n_households = 5
households = []
for i in range(n_households):
    households.append([float(np.random.uniform(0,10,1)), float(np.random.uniform(0,10,1))])
households = np.array(households)
#households = np.random.rand(n_households, 2)
#bounding_box = np.array([0., 1., 0., 1.]) # [x_min, x_max, y_min, y_max]

bounding_box = np.array([0., 10., 0., 10.]) # [x_min, x_max, y_min, y_max]
vor = voronoi(households, bounding_box)


    

    




[ True  True  True  True  True]
[[ 6.17038772  6.8657563 ]
 [ 7.6635099   9.04430825]
 [ 1.97552373  8.3309021 ]
 [ 8.95717121  9.92134941]
 [ 7.29467212  9.36690738]]

In [15]:

    

fig = plt.figure()
ax = fig.gca()

# Plot initial points
ax.plot(vor.filtered_points[:, 0], vor.filtered_points[:, 1], 'b.')

# Plot ridges points
for region in vor.filtered_regions:
    vertices = vor.vertices[region, :]
    ax.plot(vertices[:, 0], vertices[:, 1], 'go')

# Plot ridges
for region in vor.filtered_regions:
    vertices = vor.vertices[region + [region[0]], :]
    ax.plot(vertices[:, 0], vertices[:, 1], 'k-')

# Compute and plot centroids
centroids = []
for region in vor.filtered_regions:
    vertices = vor.vertices[region + [region[0]], :]
    centroid = centroid_region(vertices)
    centroids.append(list(centroid[0, :]))
    ax.plot(centroid[:, 0], centroid[:, 1], 'r.')

ax.set_xlim([-0.1, 10.1])
ax.set_ylim([-0.1, 10.1])
plt.savefig("bounded_voronoi.png")

voronoi_plot_2d(vor)
plt.savefig("voronoi.png")


    

In [ ]:

    

len(vor.filtered_points) == len(vor.filtered_regions)


    

In [ ]:

    

simulated_lines = [shapely.geometry.LineString(vor.filtered_points[line])
                    for line in vor.filtered_regions
                    if -1 not in line]


    

In [ ]:

    

simulated_gdf = gpd.GeoSeries([p for p in shapely.ops.polygonize(vor.filtered_regions)])


    

In [ ]:

    

simulated_gdf.plot()


    

In [ ]:

    

vor.filtered_regions


    

In [ ]:

    

vor.filtered_points


    

In [ ]: