Airport position Voronoi tesselation

In [1]:
import numpy as np
import matplotlib.pyplot as plt
%pylab inline

Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.zmq.pylab.backend_inline].
For more information, type 'help(pylab)'.

download data from openflights.org

In [2]:
import urllib2
data = urllib2.urlopen("http://sourceforge.net/p/openflights/code/HEAD/tree/openflights/data/airports.dat?format=raw")

parse csv to dictionary

In [3]:
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 [4]:
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 [5]:
from mpl_toolkits.basemap import Basemap
m = Basemap(projection='merc', resolution = 'h', llcrnrlon=minlon, llcrnrlat=minlat, urcrnrlon=maxlon, urcrnrlat=maxlat )

crate Voronoi tessellation

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

Draw

In [7]:
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()
Back to top