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.