Generating great circles from a common origin with Python

Following Claudia Engel's blog post: Great circles on a recentered worldmap, in ggplot, in which a solution in R is provided, I was asked to build an alternative solution in Python. I used the following components:

Python 2.7
pandas 0.7.0.dev
pyproj 1.9.0

	# See: http://bit.ly/GKgHU2
	import pandas
	import pyproj
	from string import Template
	# Dataframe from table of airports
	airports = pandas.read_csv("data/airports_new.csv")
	# Dataframe from table of flights departing from Beijing airport (PEK)
	flights = pandas.read_csv("data/PEK.csv")
	# Count number of flights from PEK
	counts = flights.groupby(['From','To'])['Distance'].count().reset_index()
	vector_columns = ['latitude','longitude','altitude','IATA','name']
	# Dataframe merging counts with additional destination airport information
	flights_ll = pandas.merge(counts, airports, left_on='To', right_on='IATA')[[0]+vector_columns]
	# lat, lon, z and name simulated vector of PEK airport
	PEK = airports.ix[ airports['IATA']=='PEK', vector_columns].as_matrix()[0]
	all_airports_from_pek = flights_ll[vector_columns].as_matrix()
	# Geodetic converter
	geod = pyproj.Geod(ellps='WGS84')
	# Number of points of linestring that approximates each great circle
	STEPS = 100
	# KML output.
	# NOTE: I could have used pyKML but preferred manual string templating assembly.
	def toKML(great_circles):
	linestyle_tpl = Template( """
	<Style id="linestyle$count">
	<LineStyle>
	<color>7f0000ff</color>
	<width>$count</width>
	<gx:labelVisibility>1</gx:labelVisibility>
	</LineStyle>
	</Style>"""
	)
	airports_tpl = Template("""
	<Placemark>
	<name>$name ($IATA)</name>
	<styleUrl>#sn_shaded_dot</styleUrl>
	<Point>
	<extrude>true</extrude>
	<altitudeMode>relativeToGround</altitudeMode>
	<coordinates>$lon,$lat,2</coordinates>
	</Point>
	</Placemark>"""
	)
	linestring_tpl = Template("""
	<Placemark>
	<name>PEK to $IATA</name>
	<styleUrl>#linestyle$count</styleUrl>
	<LineString>
	<extrude>1</extrude>
	<tessellate>1</tessellate>
	<coordinates>
	$coordinates
	</coordinates>
	</LineString>
	</Placemark>"""
	)
	KML_tpl = Template("""<?xml version="1.0" encoding="UTF-8"?>
	<kml xmlns="http://www.opengis.net/kml/2.2">
	<Document>
	<name>Flights from PEK</name>
	<open>1</open>
	<Style id="sn_shaded_dot">
	<IconStyle>
	<scale>1.2</scale>
	<Icon>
	<href>http://maps.google.com/mapfiles/kml/shapes/shaded_dot.png</href>
	</Icon>
	</IconStyle>
	</Style>
	$linestyles
	<Folder>
	<name>Airports</name>
	$airports
	</Folder>
	<Folder>
	<name>Flights</name>
	$linestrings
	</Folder>
	</Document>
	</kml>"""
	)
	return KML_tpl.substitute(
	dict(
	linestyles = ''.join([ linestyle_tpl.substitute(dict(count = i)) for i in xrange(1,11)]),
	airports = ''.join([ airports_tpl.substitute(dict(name=a[-1],IATA=a[-2],lat=a[0],lon=a[1])) for a in all_airports_from_pek]),
	linestrings =''.join([ linestring_tpl.substitute(gc) for gc in great_circles])
	)
	)
	def great_arc(FROM, TO):
	# geod.npts takes points as (lon, lat), not as (lat, lon)
	great_circle = geod.npts(FROM[1], FROM[0], TO[1], TO[0], STEPS)
	#az1, az2, D = geod.inv(FROM[1], FROM[0], TO[1], TO[0])
	for i, (lat, lon) in enumerate(great_circle):
	# Altitude is 0, irrelevant at this scale
	yield (lat, lon, 0)
	def main():
	great_circles = [
	dict(
	count = F[0],
	name = F[-1],
	IATA = F[-2],
	coordinates = ' '.join(['%s,%s,%s' % p for p in great_arc(PEK, F[1:])])
	) for F in flights_ll.as_matrix() # F[0] is the flight count and F[:1] destination vector
	]
	return toKML(great_circles)
	if __name__ == '__main__':
	print main()

view raw great_circles.py hosted with ❤ by GitHub

The CSV data files with airport locations and flight information from Beijing (PEK) airport are airports_new.csv and PEK.csv, respectively. These data were derived from openflights.org.

The KML generated is easily visualized in Google Earth:

As expected, the vector layer is displayed in other GIS desktop clients unpleasantly. See for instance how it looks in QGIS:

I decided not to alter the coordinates and recenter the world map around Beijing and fix the split polygons problem. Google Earth does the trick.

Produce the output KML file from command line (UNIX) like so:

$ python great_circles.py > out.kml

Here is the compressed KML output file.

Alternative partial implementations

Shown next is how I use pyKML (0.1.0) to generate 'almost' the same output above:

	# Programmatic KML output generation using
	# pyKML http://pypi.python.org/pypi/pykml
	def toKML_pykml(great_circles):
	from pykml.factory import KML_ElementMaker as KML
	from pykml.factory import GX_ElementMaker as GX
	from lxml import etree
	doc = KML.kml(
	KML.Document(
	KML.name('Flights from PEK'),
	KML.Style(
	KML.IconStyle(
	KML.scale('1.2'),
	KML.Icon(
	KML.href('http://maps.google.com/mapfiles/kml/shapes/shaded_dot.png')
	)
	),
	id = 'sn_shaded_dot'
	),
	KML.Folder(
	KML.name('Airports'),
	*[
	KML.Placemark(
	KML.name("%s (%s)"% (a[-1], a[-2])),
	KML.styleUrl("#sn_shaded_dot"),
	KML.Point(
	KML.extrude('true'),
	KML.altitudeMode('relativeToGround'),
	KML.coordinates('%s,%s,2' % (a[1],a[0]))
	)
	) for a in all_airports_from_pek
	],
	id = 'a'
	),
	KML.Folder(
	KML.name('Flights'),
	id = 'f',
	*[
	KML.Placemark(
	KML.name('PEK to %s' % gc['IATA']),
	KML.styleUrl('#linestyle%s' % gc['count']),
	KML.LineString(
	KML.extrude('1'),
	KML.tessellate('1'),
	KML.coordinates(gc['coordinates'])
	)
	) for gc in great_circles
	]
	),
	*[
	KML.Style(
	KML.LineStyle(
	KML.color('7f0000ff'),
	KML.width('%d' % i),
	GX.labelVisibility('1'),
	),
	id='linestyle%d' % i
	) for i in xrange(1,11)
	])
	)
	return etree.tostring(doc, pretty_print=True)

view raw toKML_pykml.py hosted with ❤ by GitHub