#!/usr/bin/env python
# coding: utf-8

# # Airport position Voronoi tesselation

# In[2]:


import numpy as np
import matplotlib.pyplot as plt
get_ipython().run_line_magic('pylab', 'inline')


# ### download data from openflights.org

# In[4]:


import urllib2
data = urllib2.urlopen("https://raw.githubusercontent.com/jpatokal/openflights/master/data/airports.dat")


# ### parse csv to dictionary

# In[5]:


import csv
data_str = data.read().split('\n')
reader = csv.DictReader(data_str,  fieldnames=('n', 'name1', 'name2', 'nation', 'id1', 'id2', 'lat', 'lon', 'x', 'y', 'z', 'q'), delimiter=',', quotechar='"')    


# ### select airports in a box region

# In[6]:


minlat = 33; maxlat = 65; minlon = -17; maxlon = 27;
airport_data = [x for x in reader if minlat-5 < float(x['lat']) < maxlat + 5 and minlon - 5 < float(x['lon']) < maxlon + 5]
lonlat = np.array([(float(x['lon']), float(x['lat'])) for x in airport_data])


# ### create map

# In[9]:


from mpl_toolkits.basemap import Basemap
m = Basemap(projection='merc', resolution = 'h', llcrnrlon=minlon, llcrnrlat=minlat, urcrnrlon=maxlon, urcrnrlat=maxlat )


# ### crate Voronoi tessellation

# In[12]:


x, y = m(lonlat[:,0], lonlat[:,1])
xy = np.dstack((x, y))[0]
from scipy.spatial import Voronoi
vor = Voronoi(xy)


# ### Draw

# In[13]:


fig = plt.figure(figsize=(15,15), dpi=100)
ax = fig.add_axes([0.1, 0.1, 0.8, 0.8]) 
m.drawcoastlines(color='white')
m.drawcountries(color='white')
m.fillcontinents(color = "0.8")
#for a in airport_data:
#    if minlon < float(a['lon']) < maxlon and minlat < float(a['lat']) < maxlat:
#        x, y = m(float(a['lon']), float(a['lat']))   
#        ax.text(x, y, a['id1'])
ax.plot(xy[:,0], xy[:,1], 'ko', ms=1.5)
for simplex in vor.ridge_vertices:
    simplex = np.asarray(simplex)
    if np.all(simplex >= 0):
        ax.plot(vor.vertices[simplex,0], vor.vertices[simplex,1], '-', color='darkred')
center = xy.mean(axis=0)
for pointidx, simplex in zip(vor.ridge_points, vor.ridge_vertices):
    simplex = np.asarray(simplex)
    if np.any(simplex < 0):
        i = simplex[simplex >= 0][0]
        t = xy[pointidx[1]] - xy[pointidx[0]]
        t /= np.linalg.norm(t)
        n = np.array([-t[1], t[0]])
        midpoint = xy[pointidx].mean(axis=0)
        far_point = vor.vertices[i] + np.sign(np.dot(midpoint - center, n)) * n * 10000000
        ax.plot([vor.vertices[i,0], far_point[0]], [vor.vertices[i,1], far_point[1]], '--', color='darkred')
ax.text(0.99, 0.01, 'data from OpenFlights.org',
        horizontalalignment='right',
        verticalalignment='bottom',
        transform=ax.transAxes)
plt.show()


# In[ ]: