Great circle flights from Beijing, arcpy version

In the previous post, we found a pure Python solution to generating great circles from a common point, specifically from Beijing airports to all other airports in the world that maintain commercial flight routes from it. This situation is modeled by finding geodesic arcs, great circles from point to point in the spheroid.

In this post we present a solution implemented only with arcpy (ESRI's ArcGIS 10.0 Python geoprocessing package):

	# Great circle flights from Beijing - arcpy version
	# @author rnz0
	# See http://bit.ly/HeLThl
	import arcpy
	import os
	arcpy.env.workspace = "C:/pythonGIS/great_circles"
	inTables = ["data/PEK.csv", "data/airports_new.csv"]
	# Create out geodatabase in output/gc.db
	# This file may already exist so wrap it in a try/except block
	geodb = ("C:/pythonGIS/great_circles/output", "gc.gdb")
	try:
	arcpy.CreateFileGDB_management(geodb[0], geodb[1])
	except:
	pass
	geodb = '/'.join(geodb)
	prj = "Coordinate Systems/Geographic Coordinate Systems/World/WGS 1984.prj"
	# We first need to convert csv files to dBase tables for being
	# able to obtain group by, join, etc.
	arcpy.TableToDBASE_conversion(inTables, arcpy.env.workspace)
	PEK_table = "C:/pythonGIS/great_circles/PEK.dbf"
	AIRPORTS_table = "C:/pythonGIS/great_circles/airports_new.dbf"
	# Join flights from Beijing table (PEK.dbf) with all airports (airports_new.dbf)
	# to add latitude and longitude columns information
	# http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#/Join_Field/001700000065000000/
	arcpy.JoinField_management(
	PEK_table, # The table or feature class to which the join table will be joined.
	"To_OID", # The field in the input table on which the join will be based.
	AIRPORTS_table, # The table to be joined to the input table.
	"ID", # The field in the input table on which the join will be based.
	"#" # The fields in the input table on which the join will be based. (optional)
	)
	# Join again with PEK to repeat Beijing's coordinates on each row.
	# Longitude and latitude columns will be not be replaced, new latitud_1,
	# and longitude_1 columns will be created containing the value 40.08, 111.48 repeated.
	# Notice here latitude_1 is not created but rather latitud_1 to keep column name
	# to 10 characters wide
	arcpy.JoinField_management(
	PEK_table,
	"From_OID",
	AIRPORTS_table,
	"ID",
	"latitude;longitude"
	)
	# Create Geodesic arcs from PEK to any other airport
	great_arcs = os.path.join(geodb, "great_arcs")
	# http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//0017000000tv000000
	arcpy.XYToLine_management(
	PEK_table, # Table with fields for starting and ending X and Y coordinates.
	great_arcs, # Feature class that will contain the output lines.
	"longitud_1", # Starting point's X coordinate field.
	"latitude_1", # Starting point's Y coordinate field.
	"longitude", # Ending point's X coordinate field.
	"latitude", # Ending point's Y coordinate field.
	"GEODESIC",# Type of two-point lines to construct. Straight line based on a spheroid.
	"#", # ID field from the input table. T
	prj # Spatial Reference of the input coordinates
	)

view raw great_circles_arc.py hosted with ❤ by GitHub

Here you can download the complete script with data. Change the root of your workspace from "C:\pythonGIS" accordingly and unzip the files in the new root defined.

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