# Prepare data & calculate convex hull¶

"The convex hull of a set of points is defined as the smallest convex polygon, that encloses all of the points in the set. Convex means, that the polygon has no corner that is bent inwards". From: https://medium.com/@pascal.sommer.ch/a-gentle-introduction-to-the-convex-hull-problem-62dfcabee90c

In [27]:
#Import the source data and libraries
import pandas as pd
import geopandas as gpd
import folium
import os
import folium
import json
import requests
from shapely.geometry import Polygon
from scipy.spatial import ConvexHull, convex_hull_plot_2d
import numpy
import numpy as np
import matplotlib.pyplot as plt
from functools import reduce
import operator
import math

dfx = df[df.groupby('tac')['latitude'].transform('nunique') > 2]

'''
For each TAC create a zip list of all lat, long pairs
group by tac, giving you a single row per TAC
Pass each row into the convex hull
Plot the original points
Plot the outline
'''

df1 = dfx[['tac','latitude', 'longitude']].copy()
df1.columns = ['TAC','Lat','Long']
df1['coordinates'] = list(zip(df1.Lat, df1.Long))
tac_grouped = df1.groupby('TAC').agg(list)[['coordinates']].reset_index()
shapes = []
polygons = []

tac = 0
for index, row in tac_grouped.iterrows():

#Set variables and create required empty lists
tac = tac +1
latz = []
longz = []

pointsx = np.array(tac_grouped.loc[index,'coordinates'])

hull = ConvexHull(pointsx)
plt.plot(pointsx[:,0], pointsx[:,1], 'o')

#Simplicies are indices of points forming the sides of the convex hull.
#Simplex is a point within the simplicies
for simplex in hull.simplices:
plt.plot(pointsx[simplex, 0], pointsx[simplex, 1], 'k-')

#Extract the null simplicies to get the lat longs of the co-ordinates around the perimeter
hull_indices = np.unique(hull.simplices.flat)
hull_pts = pointsx[hull_indices, :]

center = tuple(map(operator.truediv, reduce(lambda x, y: map(operator.add, x, y), hull_pts), [len(hull_pts)] * 2))
sorted_clockwise = (sorted(hull_pts, key=lambda coord: (-135 - math.degrees(math.atan2(*tuple(map(operator.sub, coord, center))[::-1]))) % 360))

#LOOP THROUGH THE SORTED POINTS & PUT THEM INTO A NUMPY ARRAY
points = []
for x in sorted_clockwise:
lat = x[0]
longit = x[1]
points.append([lat, longit])
points = numpy.array(points)

#Extract the lat and long values to pass into the mapping
lat_point_list = points[:,0]
lon_point_list = points[:,1]

#Create the polygomns
polygon_geom = Polygon(zip(lon_point_list, lat_point_list))
crs = {'init': 'epsg:4326'}
polygon = gpd.GeoDataFrame(index=[0], crs=crs, geometry=[polygon_geom])

#Each created polygon gets added to a list
polygons.append(polygon)

#output to files
polygon.to_file(filename=str(tac)+'polygon.geojson', driver='GeoJSON')
polygon.to_file(filename=str(tac)+'polygon.shp', driver="ESRI Shapefile")

In [68]:
    #plot on a map with central point being birmingham
m = folium.Map([51.509865, -0.118092], zoom_start=10, tiles='cartodbpositron')

#Loop through the list of polygons and plot them on a map
colours = ['Red', 'Black', 'Orange', 'Yellow', 'Blue', 'Purple', 'Green', 'Green', 'Green']

id = 1
for x in polygons:

#Set colour for this layer
try:
cul = colours[id]
except:
cul = 'Black'
id = id+1

#Plot the layer
folium.GeoJson(
x,

name='TAC Area',
style_function=lambda feature: {
'fillColor': cul,
'color': cul,
'weight': 1,
'dashArray': '5, 5',
'fillOpacity': 0.5,
}