import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
import seaborn as sns
from IPython import display
from mpl_toolkits.basemap import Basemap # use for plotting lat lon & world maps
import pandas as pd
%matplotlib inline
# ^^ make plots appear in IPython notebook instead of separate window
#use below for HMM
from sklearn import linear_model
from sklearn.hmm import GaussianHMM, GMMHMM
f = open('runways' + '.txt', 'r') #renamed XPlane data file -- massive list of runway properties around the world
line_list = []
for line in f.readlines():
if '100 ' in line:
line.rstrip('\n')
line_list.append(line)
print "number of runways:", len(line_list)
number of runways: 34344
print line_list[0]
100 54.86 2 0 0.00 1 3 0 13 22.32526300 114.19222700 563.88 60.05 3 0 0 2 31 22.30395000 114.21587400 220.98 0.00 3 1 0 2
#example of 1 that won't work ... need to filter more
print line_list[76]
14 49.02711300 -122.37642100 90 0 Abbotsford Tower
new_line_list = []
for line in line_list:
if int(line[0]) == 1 and int(line[1]) == 0 and int(line[2]) == 0:
new_line_list.append(line)
print "number of runways:", len(new_line_list)
number of runways: 33527
coord1 = []
coord2 = []
coord3 = []
coord4 = []
for line in new_line_list:
temp = line.split()
try:
coord1.append(temp[9])
coord2.append(temp[10])
coord3.append(temp[18])
coord4.append(temp[19])
except:
break
#check to make sure all lengths equal
#print coord1[0:10]
print len(coord1)
print len(coord2)
print len(coord3)
print len(coord4)
33527 33527 33527 33527
print coord1[0:10]
print coord2[0:10]
print coord3[0:10]
print coord4[0:10]
['22.32526300', '29.88083567', '29.87181142', '-23.21480832', '-23.53246610', '-23.52883459', '29.65570843', '29.65007606', '32.77822536', '32.77741703'] ['114.19222700', '-103.70108085', '-103.70136737', '-051.19038734', '133.68587381', '133.66776261', '-096.57985673', '-096.58204544', '-104.38455756', '-104.37910499'] ['22.30395000', '29.86895821', '29.87798246', '-23.21352502', '-23.52755315', '-23.51384459', '29.64476323', '29.65039560', '32.78200130', '32.78280963'] ['114.21587400', '-103.69317249', '-103.69288597', '-051.18127931', '133.68954508', '133.69025847', '-096.57897605', '-096.57678734', '-104.37814356', '-104.38359613']
coordlist = np.zeros((len(new_line_list), 4))
for i in xrange(len(coord1)):
coordlist[i][0] = coord1[i]
coordlist[i][1] = coord2[i]
coordlist[i][2] = coord3[i]
coordlist[i][3] = coord4[i]
#print len(coordlist)
print coordlist[0:10]
[[ 22.325263 114.192227 22.30395 114.215874 ] [ 29.88083567 -103.70108085 29.86895821 -103.69317249] [ 29.87181142 -103.70136737 29.87798246 -103.69288597] [ -23.21480832 -51.19038734 -23.21352502 -51.18127931] [ -23.5324661 133.68587381 -23.52755315 133.68954508] [ -23.52883459 133.66776261 -23.51384459 133.69025847] [ 29.65570843 -96.57985673 29.64476323 -96.57897605] [ 29.65007606 -96.58204544 29.6503956 -96.57678734] [ 32.77822536 -104.38455756 32.7820013 -104.37814356] [ 32.77741703 -104.37910499 32.78280963 -104.38359613]]
print coordlist[0]
[ 22.325263 114.192227 22.30395 114.215874]
print coordlist[1][0]+coordlist[1][2], coordlist[1][1]+coordlist[1][3]
print coordlist[2][0]+coordlist[2][2], coordlist[2][1]+coordlist[2][3]
#probably are same runways, just in a criss cross format!
59.74979388 -207.39425334 59.74979388 -207.39425334
lat = []
lon = []
for i in xrange(len(coordlist)):
lat.append( (coordlist[i][0]+coordlist[i][2])/2.0 )
lon.append( (coordlist[i][1]+coordlist[i][3])/2.0 )
print lat[0:10]; print lon[0:10]
[22.3146065, 29.874896939999999, 29.874896939999999, -23.214166669999997, -23.530009624999998, -23.52133959, 29.65023583, 29.65023583, 32.780113329999999, 32.780113330000006] [114.20405049999999, -103.69712666999999, -103.69712667, -51.185833325000004, 133.687709445, 133.67901054000001, -96.579416390000006, -96.579416390000006, -104.38135056, -104.38135056]
coordlist_average = np.zeros((len(coordlist), 2))
for i in xrange(len(lon)):
coordlist_average[i][0] = lat[i]
coordlist_average[i][1] = lon[i]
print coordlist_average[0:10]
print len(coordlist_average) #check to make sure they're all still there!
[[ 22.3146065 114.2040505 ] [ 29.87489694 -103.69712667] [ 29.87489694 -103.69712667] [ -23.21416667 -51.18583333] [ -23.53000962 133.68770944] [ -23.52133959 133.67901054] [ 29.65023583 -96.57941639] [ 29.65023583 -96.57941639] [ 32.78011333 -104.38135056] [ 32.78011333 -104.38135056]] 33527
"""
Haversine equation.
Computes the great circle distance between two pairs of longitude/latitude.
Returns the distance in m or km depending on input (I use meters.)
"""
def haversine(r,lat1,lon1,lat2,lon2):
d = 2.0*r*np.arcsin(np.sqrt(np.sin(np.radians(lat2-lat1)/2.0)**2+np.cos(np.radians(lat1))*np.cos(np.radians(lat2))*np.sin(np.radians(lon2-lon1)/2.0)**2))
return d
earth_radius = 6371000 #in m, http://en.wikipedia.org/wiki/Earth
#testing first data set result
haversine(earth_radius,coordlist[0][0],coordlist[0][1],coordlist[0][2],coordlist[0][3])
3396.108094394885
coordlist_len = np.zeros((len(coordlist),1))
for i in xrange(len(coordlist_len)):
coordlist_len[i] = haversine(earth_radius,coordlist[i][0],coordlist[i][1],coordlist[i][2],coordlist[i][3])
print coordlist_len[0:10]
[[ 3396.10809439] [ 1525.02853835] [ 1067.51875355] [ 941.64587632] [ 662.21304719] [ 2835.27738788] [ 1220.0226416 ] [ 509.35849511] [ 732.01281279] [ 732.01324227]]
#convert ft to m
min_run_dist = 5000* 0.3048 # 5000 ft * 0.3048 ft/m
print "The minimum runway distance to land a Boeing 777 is", min_run_dist, "meters."
The minimum runway distance to land a Boeing 777 is 1524.0 meters.
runway777_indices = []
for i in xrange(len(coordlist_len)):
if coordlist_len[i] >= 1524.0 :
runway777_indices.append(i)
print len(runway777_indices)
print runway777_indices[0:50]
7357 [0, 1, 5, 11, 23, 27, 28, 29, 31, 35, 45, 46, 50, 51, 60, 61, 64, 65, 73, 74, 77, 79, 81, 83, 87, 89, 90, 93, 96, 98, 101, 102, 104, 117, 118, 119, 125, 127, 128, 130, 131, 132, 133, 134, 137, 138, 139, 155, 158, 159]
master_777_runway_coords = coordlist[runway777_indices]
master_777_runway_avg = coordlist_average[runway777_indices]
print len(master_777_runway_coords) == len(master_777_runway_avg)
print "-----------------------"
print master_777_runway_coords
print "-----------------------"
print master_777_runway_avg
True ----------------------- [[ 22.325263 114.192227 22.30395 114.215874 ] [ 29.88083567 -103.70108085 29.86895821 -103.69317249] [ -23.52883459 133.66776261 -23.51384459 133.69025847] ..., [ 37.09965032 -113.59303424 37.08153031 -113.593116 ] [ 55.771524 12.300378 55.767022 12.347614 ] [ 52.159333 4.404 52.173 4.432 ]] ----------------------- [[ 22.3146065 114.2040505 ] [ 29.87489694 -103.69712667] [ -23.52133959 133.67901054] ..., [ 37.09059032 -113.59307512] [ 55.769273 12.323996 ] [ 52.1661665 4.418 ]]
# The Inmarsat satellite is at 0,64.5 -- it's geostationary.
inmarsat = [0, 64.5]
#Now we plot the plane's known positions for this part of the model -- where it took off and where it disappeared
#Kuala Lumpur International Airport Coordinates: http://www.distancesfrom.com/my/Kuala-Lumpur-Airport-(KUL)-Malaysia-latitude-longitude-Kuala-Lumpur-Airport-(KUL)-Malaysia-latitude-/LatLongHistory/3308940.aspx
kualalumpur = [2.7544829, 101.7011363]
# Igari Waypoint. Source: # http://www.fallingrain.com/waypoint/SN/IGARI.html Given in hours,minutes,seconds.
igariwaypoint = [6. + 56./60. + 12./3600., 103. + 35./60. + 6./3600.]
print "inmarsat lat/lon:", inmarsat[0], inmarsat[1]
print "kuala lumpur int'l airport coord lat/lon:", kualalumpur[0],kualalumpur[1]
print "igari waypoint lat/lon:", igariwaypoint[0],igariwaypoint[1]
inmarsat lat/lon: 0 64.5 kuala lumpur int'l airport coord lat/lon: 2.7544829 101.7011363 igari waypoint lat/lon: 6.93666666667 103.585
eq_deg_km = 111.32 # number of km/degree at eq Source: http://en.wikipedia.org/wiki/Decimal_degrees
earth_radius = 6371 # overwrite previous earth radius in km, http://en.wikipedia.org/wiki/Earth
"""
ellipse function
derived from spherical trigonometry
(same from Part 1)
returns coords_array, which is in the form lat,lon
"""
def make_circle(radius,pts,lon_loc,lat_loc):
coords_array = np.zeros((pts, 2))
for i in xrange(pts):
#lat
coords_array[i][0] = radius*np.cos(np.radians(i)) + lat_loc
#lon
coords_array[i][1] = radius*np.sin(np.radians(i))/(np.cos(np.cos(np.radians(radius))*(lat_loc+radius*np.cos(np.radians(i)))*(np.pi/180.))) + lon_loc
return coords_array
cruise_speed = 905 # Boeing 777 cruising speed in km/hr
flight_time = 6.5 # hrs of flight time assumption
max_radius = (cruise_speed/eq_deg_km) * flight_time
#now make the circle, w/ 360 pts
max_circle = make_circle(max_radius,360,igariwaypoint[1],igariwaypoint[0])
runway_lats = []
runway_lons = []
for i in xrange(len(master_777_runway_avg)):
runway_lats.append(master_777_runway_avg[i][0])
runway_lons.append(master_777_runway_avg[i][1])
max_circle_lats = []
max_circle_lons = []
for i in xrange(len(max_circle)):
max_circle_lats.append(max_circle[i][0])
max_circle_lons.append(max_circle[i][1])
#load data and put in DataFrame
# http://en.wikipedia.org/wiki/Largest_cities_in_Asia
asiatoptencities = pd.read_csv('AsiaTopTenCities.csv')
#calculate city sum
cities_sum = np.sum(asiatoptencities['Population'])
#assemble fractions for cities
fractions = []
for i in xrange(len(asiatoptencities)):
fractions.append(np.float(asiatoptencities['Population'][i])/cities_sum)
#initialize new DataFrame column w/ the fractions
asiatoptencities['Weighted Population Fraction'] = fractions; asiatoptencities
City | Lat | Lon | Population | Country | Source | Notes | Weighted Population Fraction | |
---|---|---|---|---|---|---|---|---|
0 | Tokyo | 35.689506 | 139.691700 | 13185502 | Japan | http://en.wikipedia.org/wiki/Tokyo | NaN | 0.061747 |
1 | Seoul | 37.566536 | 126.977969 | 25620000 | South Korea | http://en.wikipedia.org/wiki/Seoul | NaN | 0.119977 |
2 | Shanghai | 31.200000 | 121.500000 | 21766000 | China | http://en.wikipedia.org/wiki/Shanghai | NaN | 0.101929 |
3 | Karachi | 24.860000 | 67.010000 | 23500000 | Pakistan | http://en.wikipedia.org/wiki/Karachi | NaN | 0.110050 |
4 | Beijing | 39.913889 | 116.391667 | 21150000 | China | http://en.wikipedia.org/wiki/Beijing | NaN | 0.099045 |
5 | Mumbai | 18.975000 | 72.825830 | 18414288 | India | http://en.wikipedia.org/wiki/Mumbai | NaN | 0.086233 |
6 | Jakarta | -6.200000 | 106.800000 | 28019545 | Indonesia | http://en.wikipedia.org/wiki/Jakarta | NaN | 0.131214 |
7 | Delhi | 28.610000 | 77.230000 | 22000000 | India | http://en.wikipedia.org/wiki/Delhi | NaN | 0.103025 |
8 | Osaka-Kobe-Kyoto | 35.011667 | 135.768333 | 18643915 | Japan | http://en.wikipedia.org/wiki/Osaka-Kobe-Kyoto | Used Kyoto lat/lon coords: http://en.wikipedia... | 0.087309 |
9 | Manila | 14.583333 | 120.966667 | 21241000 | Philippines | http://en.wikipedia.org/wiki/Manila | NaN | 0.099471 |
cities_lat = []
cities_lon = []
for i in xrange(len(asiatoptencities)):
cities_lat.append(asiatoptencities['Lat'][i])
cities_lon.append(asiatoptencities['Lon'][i])
#create lat/long grid of distance
lat_min = -45 #45 S
lat_max = 60 #60 N
lon_min = 50 #50 E
lon_max = 180 #180 E
lat_space = 5 #spacing for plotting latitudes and longitudes
lon_space = 5
"""
figure_function Bare bones basic map function for our HMM. We'll build on this in subsequent figures.
fig a Basemap figure that you initialize (i.e. assign the figure creation as a variable, fig)
"""
def figure_function(fig):
#Draw coasts
fig.drawcoastlines()
#Draw boundary
fig.drawmapboundary(fill_color='lightblue')
#Fill background
fig.fillcontinents(color='white',lake_color='lightblue')
#Draw parallels
parallels = np.arange(lat_min,lat_max,lat_space)
fig.drawparallels(np.arange(lat_min,lat_max,lat_space),labels=[1,1,0,1], fontsize=15)
#Draw meridians
meridians = np.arange(lon_min,lon_max,lon_space)
fig.drawmeridians(np.arange(lon_min,lon_max,lon_space),labels=[1,1,0,1], fontsize=15)
#Translate coords into map coord system to plot
#Known 777 Locs
x,y = fig(kualalumpur[1],kualalumpur[0]) #plotted as lon,lat NOT lat,lon -- watch out!!
x2,y2 = fig(igariwaypoint[1],igariwaypoint[0])
#Runway Locs
x3,y3 = fig(runway_lons,runway_lats)
#Max Circle
x4,y4 = fig(max_circle_lons,max_circle_lats)
#Cities Locs
x5,y5 = fig(cities_lon,cities_lat)
# plot coords w/ filled circles
fig.plot(x2,y2,'ko',markersize=15,label='Last Known MH370 Position')
fig.plot(x3,y3,'ro',markersize=5,label='Landable Runways')
fig.plot(x4,y4,'k--',markersize=10,label='Maximum Flight Reach')
fig.plot(x5,y5,'bo',markersize=8,label='Asian Cities, Top 10 Population')
#Draw arrows showing flight path
arrow1 = plt.arrow(x,y,x2-x,y2-y,linewidth=3,color='black',linestyle='dashed',label='flight path')
#Set figure size
fig = plt.figure(figsize=[30,20])
#Setup Basemap
fig = Basemap(width=13000000,height=15000000,projection='lcc',resolution='c',lat_0=5,lon_0=100,suppress_ticks=True)
#Call figure function
figure_function(fig)
#Put other plots in here -- they don't have to be in any strict order necessarily.
# (We don't have any right now)
#Make legend
legend = plt.legend(loc='lower left',fontsize=10,frameon=True,title='Legend',markerscale=1,prop={'size':15})
legend.get_title().set_fontsize('20')
#Add title
plt.title('Available Runways for MH370', fontsize=30)
#Show below
plt.show()
Electrical Fire or Depressurization (some type of movement for some fractional duration of chain followed by constant heading and speed)
Hijacked to land on a runway (heads towards a runway; allowed to switch)
Heading away from runway to ditch at sea (heads away from runways at all time steps)
Heads toward major cities (terrorism -- get coordinates of major cities)
In all scenarios add error to its actual trajectory so it isn't immediately obvious, and then run the Hidden Markov Model to detect most likely option, and see to what extent we're right for 10,000 trials (like in Lecture example)
master_777_runway_avg #coordinate list of 777 runways around the world that are landable
array([[ 22.3146065 , 114.2040505 ], [ 29.87489694, -103.69712667], [ -23.52133959, 133.67901054], ..., [ 37.09059032, -113.59307512], [ 55.769273 , 12.323996 ], [ 52.1661665 , 4.418 ]])
haversine_dists = np.zeros(len(master_777_runway_avg))
for i in xrange(len(master_777_runway_avg)):
haversine_dists[i] = haversine(earth_radius,igariwaypoint[0],igariwaypoint[1],master_777_runway_avg[i][0],master_777_runway_avg[i][1])
print "minimum runway distance in km is:", np.min(haversine_dists)
minimum runway distance in km is: 166.424481423
print "maximum runway distance in km is:", np.max(haversine_dists)
maximum runway distance in km is: 19971.897423
fig = plt.figure(figsize=[10,6])
plt.hist(haversine_dists,bins=25) #it's in km now instead of m -- earth_radius overwritten
plt.title('Histogram of Starting Distances to Runways',fontsize=15)
plt.xlabel('Distance in km')
plt.ylabel('# of Runways')
plt.show()
"""
get_heading gets heading relative to straight up which is 0 degrees (due north)
heading is direction plane has to go to, from its current location coords to its destination coords
lat_init starting latitude (current location)
lon_init starting longitude (current location)
dest_lat destination latitude
dest_lon destination longitude
"""
def get_heading(lat_init,lon_init,dest_lat,dest_lon):
#convert to x,y coordinate naming convention -- longitude is x axis and latitude is y axis
#(x0,y0) is center of circle where we start out
y0 = lat_init
x0 = lon_init
#(x,y) is along circumference of circle, where we want to end up
y = dest_lat
x = dest_lon
#test conditionally to see what half of the circle between starting location and destination, the angle will be on
#otherwise we don't get a unique answer for our angle -- there could be two possibilites for angles, one on either
#side.
if x-x0 >= 0:
#FYI -- be aware of difference between np.arctan and np.arctan2
angle_in_radians = np.pi/2 - np.arctan((y-y0)/(x-x0))
if x-x0 < 0:
angle_in_radians = 3*np.pi/2 - np.arctan((y-y0)/(x-x0))
angle_in_degrees = np.degrees(angle_in_radians)
return angle_in_degrees
max_hrs = 6.5
hrs = 2
sph = 60 # how many times we Sample Per Hour. once per minute is 60.
# (obviously there are 60 min/hr! the pt is our sampling rate might want to change)
#New note -- we use our make vector to advance the plane forward in the "ghost flight" scenario, for each time step.
#Write function that plots final destination from plane.
#Make a point at a distance on a heading to see where the plane would be if it continued on a straight line.
"""
make_vector
ang_radius how far the plane will advance in the time step -- the radius of the circle we construct
then we combine the distance (radius) with the heading to get the plane's new position
heading current heading of plane
lon_loc current longitude of plane
lat_loc current latitude of plane
"""
def make_vector(ang_radius,heading,lon_loc,lat_loc):
vec_lat = ang_radius*np.cos(np.radians(heading)) + lat_loc
vec_lon = ang_radius*np.sin(np.radians(heading))/(np.cos(np.cos(np.radians(ang_radius))*(lat_loc+ang_radius*np.cos(np.radians(heading)))*(np.pi/180.))) + lon_loc
return vec_lat,vec_lon
#Scenario 1
"""
scenario1 heads towards a runway it can get to -- right now, not allowed to switch
lat_init initial latitude of plane
lon_init initial longitude of plane
hrs hours to run simulation at
sph samples per hour. in other words, how often we update and plot the plane's position.
cruise_speed cruising speed of Boeing 777, 905 km/hr
eq_deg_km degrees per km at the equator
max_hrs amount of time plane can be in air
"""
def scenario1(lat_init,lon_init,hrs,sph,cruise_speed,eq_deg_km,max_hrs):
#number of samples we'll run simulation at
samples = hrs * sph
#max dist plane can fly
max_flying_dist = cruise_speed * max_hrs # km/hr * hrs = km
#dist that plane will fly during sampling period
sample_flying_dist = cruise_speed * hrs
runway_select = 0 #haven't chosen runway yet, so assign 0
while runway_select == 0:
#independent draw of index of runways - random runway draw
runway_draw = np.int(np.random.choice(len(master_777_runway_avg),size=1))
#now grab lats and lons of runway
dest_lat = master_777_runway_avg[runway_draw][0]
dest_lon = master_777_runway_avg[runway_draw][1]
#this runway will be the runway the plane tries to head to at all points in time throughout the simulation
#HOWEVER, we need to make sure this runway would be a realistic runway to get to given the TIME constraint, hrs
#so we compute the Haversine distance between the plane initial location and the runway we've chosen
#since earth_radius is in km, we get the dist_to_dist in km. *at first* in this notebook earth_radius was chosen
#to have units of meters since we were looking at runway lengths, but afterward we overwrote it with the earth
#radius in km.
dist_to_dest = haversine(earth_radius,lat_init,lon_init,dest_lat,dest_lon)
if dist_to_dest < max_flying_dist: #if we can fly there in the amount of time
if dist_to_dest > sample_flying_dist: #if we can't get there yet, in the amount of time to run simualtion
runway_select = 1 #we've chosen our runway
#we've chosen our runway.
#now we need to calculate the appropriate heading from 0-360 degrees of that destination
#then make_vector will take the heading, and the plane's current coords and advance it toward the runway
#we start by making the current plane location the center of our location, and the desired destination runway
#the edge of the circle, on the circumference. to get the circle's radius we use the Haversine distance formula.
circle_dest_radius = haversine(earth_radius,lat_init,lon_init,dest_lat,dest_lon)
circle_of_dest = make_circle(circle_dest_radius,360,lon_init,lat_init)
#we define 0 degrees being due north, i.e. straight up, and we know what the radius is, so using trigonometry
#we can calculate the angle of the plane's heading to advance the plane that way, and feed it into make_vector.
heading = np.around(get_heading(lat_init,lon_init,dest_lat,dest_lon))
#------------------------------------------------------------------------------#
#initialize where we'll store new plane position coords from simulation output
plane_lat = np.zeros(samples+1)
plane_lon = np.zeros(samples+1)
#initialize lat,lon,heading from input
plane_lat[0] = lat_init
plane_lon[0] = lon_init
#distance plane travels each time ... if we're doing 60 samples per hour, each time it only goes 1/60th
#distance it would travel if the units are km/hr
plane_dist = (cruise_speed/eq_deg_km)/sph
#main loop that calculates a new plane position each time
for i in xrange(samples): #for each sample -> we're going to MTP (move that plane)
#from the plane's location, make a circle with the 360 possible points it could head to
new_locations_circle = make_circle(plane_dist,360,plane_lon[0],plane_lat[0])
#we use the heading we'd need to take to reach that airport, from where the plane is currently located.
#and advance the plane in that direction
new_position = make_vector(plane_dist,heading, plane_lon[i], plane_lat[i])
plane_lat[i+1] = new_position[0] #i+1 because we had the initial coords be at 0th index, before 1st step
plane_lon[i+1] = new_position[1]
return plane_lat, plane_lon, dest_lat, dest_lon
#test function
testlats,testlons,testdestlat,testdestlon = scenario1(igariwaypoint[0],igariwaypoint[1],4,sph,cruise_speed,eq_deg_km, max_hrs)
#Set figure size
fig = plt.figure(figsize=[30,20])
#Setup Basemap
fig = Basemap(width=13000000,height=15000000,projection='lcc',resolution='c',lat_0=5,lon_0=100,suppress_ticks=True)
#Call figure function
figure_function(fig)
#Put other plots in here -- they don't have to be in any strict order necessarily.
x6,y6 = fig(testlons,testlats)
x7,y7 = fig(testdestlon,testdestlat)
fig.plot(x6,y6,'go',markersize=5,label='Plane Path')
fig.plot(x7,y7,'go',markersize=10,label='Target Runway')
#Make legend
legend = plt.legend(loc='lower left',fontsize=10,frameon=True,title='Legend',markerscale=1,prop={'size':15})
legend.get_title().set_fontsize('20')
#Add title
plt.title('Scenario 1', fontsize=30)
#Show below
plt.show()
#Scenario 2
"""
scenario2 examines scenario where plane heads to one of top 10 most populated cities in Asia.
lat_init initial latitude of plane
lon_init initial longitude of plane
hrs hours to run simulation at
sph samples per hour. in other words, how often we update and plot the plane's position.
cruise_speed cruising speed of Boeing 777, 905 km/hr
eq_deg_km degrees per km at the equator
max_hrs amount of time plane can be in air
"""
def scenario2(lat_init,lon_init,hrs,sph,cruise_speed,eq_deg_km,max_hrs):
samples = hrs * sph
max_flying_dist = cruise_speed * max_hrs # km/hr * hrs = km
sample_flying_dist = cruise_speed * hrs
#independent draw of index of Top Ten Asian cities - random city draw
city_draw = np.int(np.random.choice(len(asiatoptencities),size=None,p=fractions))
dest_lat = asiatoptencities['Lat'][city_draw]
dest_lon = asiatoptencities['Lon'][city_draw]
circle_dest_radius = haversine(earth_radius,lat_init,lon_init,dest_lat,dest_lon)
circle_of_dest = make_circle(circle_dest_radius,360,lon_init,lat_init)
heading = np.around(get_heading(lat_init,lon_init,dest_lat,dest_lon))
#------------------------------------------------------------------------------#
plane_lat = np.zeros(samples+1)
plane_lon = np.zeros(samples+1)
plane_lat[0] = lat_init
plane_lon[0] = lon_init
plane_dist = (cruise_speed/eq_deg_km)/sph
for i in xrange(samples):
new_locations_circle = make_circle(plane_dist,360,plane_lon[0],plane_lat[0])
new_position = make_vector(plane_dist,heading, plane_lon[i], plane_lat[i])
plane_lat[i+1] = new_position[0]
plane_lon[i+1] = new_position[1]
return plane_lat,plane_lon,dest_lat,dest_lon
#test function
testlats,testlons,testdestlat,testdestlon = scenario2(igariwaypoint[0],igariwaypoint[1],4,sph,cruise_speed,eq_deg_km, max_hrs)
#Set figure size
fig = plt.figure(figsize=[30,20])
#Setup Basemap
fig = Basemap(width=13000000,height=15000000,projection='lcc',resolution='c',lat_0=5,lon_0=100,suppress_ticks=True)
#Call figure function
figure_function(fig)
#Put other plots in here -- they don't have to be in any strict order necessarily.
x6,y6 = fig(testlons,testlats)
x7,y7 = fig(testdestlon,testdestlat)
fig.plot(x6,y6,'go',markersize=5,label='Plane Path')
fig.plot(x7,y7,'go',markersize=10,label='Target City')
#Make legend
legend = plt.legend(loc='lower left',fontsize=10,frameon=True,title='Legend',markerscale=1,prop={'size':15})
legend.get_title().set_fontsize('20')
#Add title
plt.title('Scenario 2', fontsize=30)
#Show below
plt.show()
#Scenario 3
"""
scenario3 Heading *toward* Top 10 cities, to ditch at sea (at all time steps)
lat_init initial latitude of plane
lon_init initial longitude of plane
hrs hours to run simulation at
sph samples per hour. in other words, how often we update and plot the plane's position.
cruise_speed cruising speed of Boeing 777, 905 km/hr
eq_deg_km degrees per km at the equator
max_hrs amount of time plane can be in air
"""
def scenario3(lat_init,lon_init,hrs,sph,cruise_speed,eq_deg_km,max_hrs):
samples = hrs * sph
plane_dist = (cruise_speed/eq_deg_km)/sph
#------------------------------------------------------------------------------#
plane_lat = np.zeros(samples+1)
plane_lon = np.zeros(samples+1)
plane_lat[0] = lat_init
plane_lon[0] = lon_init
for i in xrange(samples):
#make a circle with radius equal to the distance traveled to evaluate the distances over
#I believe we're doing it so that a 0 degree heading is first, and then around clockwise
#so we can just grab the index of the minimum, and then choose that for our next heading.
new_locs = make_circle(plane_dist,360,plane_lon[i],plane_lat[i])
cumsum_citydist = 10e5 #pick insanely high number
#now iterate through the possibilities, evaluating the distance between the Top 10 Cities
for j in xrange(len(new_locs)):
temp_citydist = 0
for k in xrange(len(asiatoptencities)): #iterate through cities
temp_citydist += haversine(earth_radius,new_locs[j][0],new_locs[j][1],
asiatoptencities['Lat'][k],asiatoptencities['Lon'][k])
if temp_citydist < cumsum_citydist:
cumsum_citydist = temp_citydist
dest_lat = new_locs[j][0]
dest_lon = new_locs[j][1]
#print cumsum_citydist #testing to make sure distance is going down
plane_lat[i+1] = dest_lat
plane_lon[i+1] = dest_lon
return plane_lat, plane_lon
#test function
testlats,testlons = scenario3(igariwaypoint[0],igariwaypoint[1],4,sph,cruise_speed,eq_deg_km,max_hrs)
#Set figure size
fig = plt.figure(figsize=[30,20])
#Setup Basemap
fig = Basemap(width=13000000,height=15000000,projection='lcc',resolution='c',lat_0=5,lon_0=100,suppress_ticks=True)
#Call figure function
figure_function(fig)
#Put other plots in here -- they don't have to be in any strict order necessarily.
x6,y6 = fig(testlons,testlats)
fig.plot(x6,y6,'go',markersize=5,label='Plane Path')
#Make legend
legend = plt.legend(loc='lower left',fontsize=10,frameon=True,title='Legend',markerscale=1,prop={'size':15})
legend.get_title().set_fontsize('20')
#Add title
plt.title('Scenario 3', fontsize=30)
#Show below
plt.show()
#Scenario 4
"""
scenario4 Heading *away* from Top 10 cities (at all time steps)
lat_init initial latitude of plane
lon_init initial longitude of plane
hrs hours to run simulation at
sph samples per hour. in other words, how often we update and plot the plane's position.
cruise_speed cruising speed of Boeing 777, 905 km/hr
eq_deg_km degrees per km at the equator
max_hrs amount of time plane can be in air
"""
def scenario4(lat_init,lon_init,hrs,sph,cruise_speed,eq_deg_km,max_hrs):
samples = hrs * sph
plane_dist = (cruise_speed/eq_deg_km)/sph
#------------------------------------------------------------------------------#
plane_lat = np.zeros(samples+1)
plane_lon = np.zeros(samples+1)
plane_lat[0] = lat_init
plane_lon[0] = lon_init
for i in xrange(samples):
new_locs = make_circle(plane_dist,360,plane_lon[i],plane_lat[i])
cumsum_citydist = 0 #pick insanely low number
for j in xrange(len(new_locs)):
temp_citydist = 0
for k in xrange(len(asiatoptencities)):
temp_citydist += haversine(earth_radius,new_locs[j][0],new_locs[j][1],
asiatoptencities['Lat'][k],asiatoptencities['Lon'][k])
if temp_citydist > cumsum_citydist:
cumsum_citydist = temp_citydist
dest_lat = new_locs[j][0]
dest_lon = new_locs[j][1]
plane_lat[i+1] = dest_lat
plane_lon[i+1] = dest_lon
return plane_lat, plane_lon
#test function
testlats,testlons = scenario4(igariwaypoint[0],igariwaypoint[1],4,sph,cruise_speed,eq_deg_km,max_hrs)
#Set figure size
fig = plt.figure(figsize=[30,20])
#Setup Basemap
fig = Basemap(width=13000000,height=15000000,projection='lcc',resolution='c',lat_0=5,lon_0=100,suppress_ticks=True)
#Call figure function
figure_function(fig)
#Put other plots in here -- they don't have to be in any strict order necessarily.
x6,y6 = fig(testlons,testlats)
fig.plot(x6,y6,'go',markersize=5,label='Plane Path')
#Make legend
legend = plt.legend(loc='lower left',fontsize=10,frameon=True,title='Legend',markerscale=1,prop={'size':15})
legend.get_title().set_fontsize('20')
#Add title
plt.title('Scenario 4', fontsize=30)
#Show below
plt.show()
#Scenario 5
"""
scenario5 Heading *away* from runways in [hrs] radius (at all time steps)
lat_init initial latitude of plane
lon_init initial longitude of plane
hrs hours to run simulation at
sph samples per hour. in other words, how often we update and plot the plane's position.
cruise_speed cruising speed of Boeing 777, 905 km/hr
eq_deg_km degrees per km at the equator
max_hrs amount of time plane can be in air
"""
def scenario5(lat_init,lon_init,hrs,sph,cruise_speed,eq_deg_km,max_hrs):
samples = hrs * sph
sample_flying_dist = cruise_speed * hrs
plane_dist = (cruise_speed/eq_deg_km)/sph
#------------------------------------------------------------------------------#
plane_lat = np.zeros(samples+1)
plane_lon = np.zeros(samples+1)
plane_lat[0] = lat_init
plane_lon[0] = lon_init
#speed up calculation -- grab all runways that are within sample flying distance from initial plane loc
close_runways = []
for i in xrange(len(master_777_runway_avg)):
if haversine(earth_radius,plane_lat[0],plane_lon[0],master_777_runway_avg[i][0],
master_777_runway_avg[i][1]) < sample_flying_dist:
close_runways.append(master_777_runway_avg[i])
for i in xrange(samples):
new_locs = make_circle(plane_dist,360,plane_lon[i],plane_lat[i])
cumsum_runwaydist = 0 #pick insanely low number
for j in xrange(len(new_locs)):
temp_runwaydist = 0
for k in xrange(len(close_runways)):
temp_runwaydist += haversine(earth_radius,new_locs[j][0],new_locs[j][1],
close_runways[k][0],close_runways[k][1])
if temp_runwaydist > cumsum_runwaydist:
cumsum_runwaydist = temp_runwaydist
dest_lat = new_locs[j][0]
dest_lon = new_locs[j][1]
plane_lat[i+1] = dest_lat
plane_lon[i+1] = dest_lon
return plane_lat, plane_lon
#test function
testlats,testlons = scenario5(igariwaypoint[0],igariwaypoint[1],4,sph,cruise_speed,eq_deg_km,max_hrs)
#Set figure size
fig = plt.figure(figsize=[30,20])
#Setup Basemap
fig = Basemap(width=13000000,height=15000000,projection='lcc',resolution='c',lat_0=5,lon_0=100,suppress_ticks=True)
#Call figure function
figure_function(fig)
#Put other plots in here -- they don't have to be in any strict order necessarily.
x6,y6 = fig(testlons,testlats)
fig.plot(x6,y6,'go',markersize=5,label='Plane Path')
#Make legend
legend = plt.legend(loc='lower left',fontsize=10,frameon=True,title='Legend',markerscale=1,prop={'size':15})
legend.get_title().set_fontsize('20')
#Add title
plt.title('Scenario 5', fontsize=30)
#Show below
plt.show()
#Scenario 6
"""
scenario6 Heading *toward* runways in [hrs] radius (at all time steps)
lat_init initial latitude of plane
lon_init initial longitude of plane
hrs hours to run simulation at
sph samples per hour. in other words, how often we update and plot the plane's position.
cruise_speed cruising speed of Boeing 777, 905 km/hr
eq_deg_km degrees per km at the equator
max_hrs amount of time plane can be in air
"""
def scenario6(lat_init,lon_init,hrs,sph,cruise_speed,eq_deg_km,max_hrs):
samples = hrs * sph
sample_flying_dist = cruise_speed * hrs
plane_dist = (cruise_speed/eq_deg_km)/sph
#------------------------------------------------------------------------------#
plane_lat = np.zeros(samples+1)
plane_lon = np.zeros(samples+1)
plane_lat[0] = lat_init
plane_lon[0] = lon_init
#speed up calculation -- grab all runways that are within sample flying distance from initial plane loc
close_runways = []
for i in xrange(len(master_777_runway_avg)):
if haversine(earth_radius,plane_lat[0],plane_lon[0],master_777_runway_avg[i][0],
master_777_runway_avg[i][1]) < sample_flying_dist:
close_runways.append(master_777_runway_avg[i])
for i in xrange(samples):
new_locs = make_circle(plane_dist,360,plane_lon[i],plane_lat[i])
cumsum_runwaydist = 1e20 #pick insanely high number
for j in xrange(len(new_locs)):
temp_runwaydist = 0
for k in xrange(len(close_runways)):
temp_runwaydist += haversine(earth_radius,new_locs[j][0],new_locs[j][1],
close_runways[k][0],close_runways[k][1])
if temp_runwaydist < cumsum_runwaydist:
cumsum_runwaydist = temp_runwaydist
dest_lat = new_locs[j][0]
dest_lon = new_locs[j][1]
plane_lat[i+1] = dest_lat
plane_lon[i+1] = dest_lon
return plane_lat, plane_lon
#test function
testlats,testlons = scenario6(igariwaypoint[0],igariwaypoint[1],4,sph,cruise_speed,eq_deg_km,max_hrs)
#Set figure size
fig = plt.figure(figsize=[30,20])
#Setup Basemap
fig = Basemap(width=13000000,height=15000000,projection='lcc',resolution='c',lat_0=5,lon_0=100,suppress_ticks=True)
#Call figure function
figure_function(fig)
#Put other plots in here -- they don't have to be in any strict order necessarily.
x6,y6 = fig(testlons,testlats)
fig.plot(x6,y6,'go',markersize=5,label='Plane Path')
#Make legend
legend = plt.legend(loc='lower left',fontsize=10,frameon=True,title='Legend',markerscale=1,prop={'size':15})
legend.get_title().set_fontsize('20')
#Add title
plt.title('Scenario 6', fontsize=30)
#Show below
plt.show()
#uncomment/comment the below to see the local maximum achieved
print "latitude length: " , len(testlats)
print "longitude length: ", len(testlons)
print "------------------"
print "last 20 latitudes: ", testlats[len(testlats)-20:]
print "last 20 longitudes: ", testlons[len(testlons)-20:]
latitude length: 241 longitude length: 241 ------------------ last 20 latitudes: [ 14.34143266 14.21965032 14.34143266 14.21965032 14.34143266 14.21965032 14.34143266 14.21965032 14.34143266 14.21965032 14.34143266 14.21965032 14.34143266 14.21965032 14.34143266 14.21965032 14.34143266 14.21965032 14.34143266 14.21965032] last 20 longitudes: [ 103.51758231 103.5788569 103.51754915 103.57882374 103.51751599 103.57879058 103.51748283 103.57875742 103.51744967 103.57872426 103.51741651 103.5786911 103.51738335 103.57865794 103.5173502 103.57862478 103.51731704 103.57859162 103.51728388 103.57855847]