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


df = pd.read_csv('gldn.csv')
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')
    #folium.GeoJson(polygon).add_to(m)
    
    #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,
            }
        ).add_to(m)
        

    folium.LayerControl().add_to(m)
    folium.LatLngPopup().add_to(m)
    m
Out[68]:
In [ ]: