From 7e851045658b94a9376946245e9531c8f424d55d Mon Sep 17 00:00:00 2001 From: simonrp84 Date: Thu, 5 Dec 2019 16:32:22 +0000 Subject: [PATCH] Update to release version, as used in the go-around paper. --- GA_Detect.py | 139 ++++++++++++++++----- OS_Consts.py | 6 +- OS_Funcs.py | 286 ++++++++++++++++++++++++++++++-------------- OS_Output.py | 190 +++++++++++++++++++---------- OpenSky_Get_Data.py | 14 +-- metar_parse.py | 68 +++++++---- 6 files changed, 489 insertions(+), 214 deletions(-) diff --git a/GA_Detect.py b/GA_Detect.py index 8ac551f..0fe66d0 100644 --- a/GA_Detect.py +++ b/GA_Detect.py @@ -1,14 +1,21 @@ +"""A script to process OpenSky ADS-B data in an attempt to detect go-around events at an airport.""" from traffic.core import Traffic from datetime import timedelta -import metar_parse as MEP import multiprocessing as mp from OS_Airports import VABB import OS_Funcs as OSF import glob -def main(start_n, fidder): +def main(start_n, fidder, do_write): + """The main code for detecting go-arounds. + Arguments: + start_n -- The index of the first file to read + fidder -- the id of an open file to write log information into + do_write -- boolean flag specifying whether to output data to textfile + + """ # Total number of aircraft seen tot_n_ac = 0 @@ -32,16 +39,28 @@ def main(start_n, fidder): odir_da_ga = top_dir + 'OUT_DATA/PSGA/' odirs = [odir_pl_nm, odir_pl_ga, odir_da_nm, odir_da_ga] - - # Read METARs from disk - metars = MEP.get_metars('/home/proud/Desktop/GoAround_Paper/VABB_METAR') - - # File to save met info for g/a flights - metfid = open('/home/proud/Desktop/GoAround_Paper/GA_MET.csv','w') - metfid.write('ICAO24, Callsign, Time, Runway, Temp, Dewp, Wind_Spd,\ - Wind_Gust, Wind_Dir,CB, Pressure\n') - files = glob.glob(indir+'*.pkl') + # Output filenames for saving data about go-arounds + out_file_ga = 'GA_MET_NEW.csv' + # Output filenames for saving data about non-go-arounds + out_file_noga = 'GA_NOGA_NEW.csv' + + t_frmt = "%Y/%m/%d %H:%M:%S" + + # File to save met info for g/a flights + if (do_write): + metfid = open(out_file_ga, 'w') + nogfid = open(out_file_noga, 'w') + metfid.write('ICAO24, Callsign, GA_Time, L_Time, Runway, Heading, Alt, Lat, \ + Lon, gapt, rocvar, hdgvar, latvar, lonvar, gspvar, \ + Temp, Dewp, Wind_Spd, Wind_Gust, Wind_Dir,Cld_Base,\ + CB, Vis, Pressure\n') + nogfid.write('ICAO24, Callsign, GA_Time, L_Time, Runway, gapt, rocvar, \ + hdgvar, latvar, lonvar, gspvar, \ + Temp, Dewp, Wind_Spd, Wind_Gust, Wind_Dir,Cld_Base,\ + CB, Vis, Pressure\n') + files = [] + files = glob.glob(indir+'**.pkl') files.sort() fli_len = len(files) @@ -52,7 +71,7 @@ def main(start_n, fidder): # Number of files to open in one go n_files_proc = 55 - pool_proc = 56 + pool_proc = 100 f_data = [] pool = mp.Pool(processes=pool_proc) @@ -81,6 +100,8 @@ def main(start_n, fidder): traf_arr = Traffic.from_flights(f_data) p_list = [] f_data = [] + # Extend array end time if there's only one flight, else processing + # will fail as we check to ensure latest flight is > 5 mins from end if (len(traf_arr) == 1): end_time = traf_arr.end_time + timedelta(minutes=10) print("Extending timespan due to single aircraft") @@ -90,13 +111,13 @@ def main(start_n, fidder): # Now we process the results for flight in traf_arr: if (flight.stop + timedelta(minutes=5) < end_time): - p_list.append(pool.apply_async(OSF.proc_fl, args=(flight, - VABB.rwy_list, - odirs, - colormap, - metars, - True, - False,))) + p_list.append(pool.apply_async(OSF.proc_fl, + args=(flight, + VABB.rwy_list, + odirs, + colormap, + True, + False,))) else: f_data.append(flight) @@ -104,23 +125,75 @@ def main(start_n, fidder): t_res = p.get() if (t_res != -1): tot_n_ac += 1 + # If there's a go-around, this will be True if (t_res[0]): - metfid.write(t_res[1] + ',' + t_res[2] + ',') - metfid.write(t_res[3].strftime("%Y%m%d%H%M") + ',') - metfid.write(t_res[4] + ',') - metfid.write(str(t_res[5].temp) + ',') - metfid.write(str(t_res[5].dewp) + ',') - metfid.write(str(t_res[5].w_s) + ',') - metfid.write(str(t_res[5].w_g) + ',') - metfid.write(str(t_res[5].w_d) + ',') - metfid.write(str(t_res[5].cb) + ',') - metfid.write(str(t_res[5].vis) + ',') - metfid.write(str(t_res[5].pres) + '\n') + if (do_write): + metfid.write(t_res[1] + ',' + t_res[2] + ',') + metfid.write(t_res[3].strftime(t_frmt) + ',') + metfid.write(t_res[4].strftime(t_frmt) + ',') + metfid.write(t_res[5] + ',') + if (t_res[6] < 0): + t_res[6] = 360 + t_res[6] + metfid.write(str(t_res[6]) + ',') + metfid.write(str(t_res[7]) + ',') + metfid.write(str(t_res[8]) + ',') + metfid.write(str(t_res[9]) + ',') + metfid.write(str(t_res[10]) + ',') + metfid.write(str(t_res[11]) + ',') + metfid.write(str(t_res[12]) + ',') + metfid.write(str(t_res[13]) + ',') + metfid.write(str(t_res[14]) + ',') + metfid.write(str(t_res[15]) + ',') + metfid.write(str(t_res[16].temp) + ',') + metfid.write(str(t_res[16].dewp) + ',') + metfid.write(str(t_res[16].w_s) + ',') + metfid.write(str(t_res[16].w_g) + ',') + metfid.write(str(t_res[16].w_d) + ',') + metfid.write(str(t_res[16].cld) + ',') + if t_res[16].cb: + metfid.write('1,') + else: + metfid.write('0,') + metfid.write(str(t_res[16].vis) + ',') + metfid.write(str(t_res[16].pres) + '\n') tot_n_ga += 1 + # Otherwise, do the following (doesn't save g/a location etc). + else: + if (do_write): + nogfid.write(t_res[1] + ',' + t_res[2] + ',') + nogfid.write(t_res[3].strftime(t_frmt) + ',') + nogfid.write(t_res[4].strftime(t_frmt) + ',') + nogfid.write(t_res[5] + ',') + nogfid.write(str(t_res[10]) + ',') + nogfid.write(str(t_res[11]) + ',') + nogfid.write(str(t_res[12]) + ',') + nogfid.write(str(t_res[13]) + ',') + nogfid.write(str(t_res[14]) + ',') + nogfid.write(str(t_res[15]) + ',') + try: + outstr = str(t_res[16].temp) + ',' + outstr = outstr + str(t_res[16].dewp) + ',' + outstr = outstr + str(t_res[16].w_s) + ',' + outstr = outstr + str(t_res[16].w_g) + ',' + outstr = outstr + str(t_res[16].w_d) + ',' + outstr = outstr + str(t_res[16].cld) + ',' + if t_res[16].cb: + outstr = outstr + '1,' + else: + outstr = outstr + '0,' + outstr = outstr + str(t_res[16].vis) + ',' + outstr = outstr + str(t_res[16].pres) + except Exception as e: + print("No METAR data for this flight") + outstr = '' + nogfid.write(outstr + '\n') + print("\t-\tHave processed " + str(tot_n_ac) + " aircraft. Have seen " + str(tot_n_ga) + " go-arounds.") - - metfid.close() + + if (do_write): + metfid.close() + nogfid.close() # Use this to start processing from a given file number. @@ -129,6 +202,6 @@ def main(start_n, fidder): fid = open('/home/proud/Desktop/log.log', 'w') -main(init_num, fid) +main(init_num, fid, False) fid.close() diff --git a/OS_Consts.py b/OS_Consts.py index 30aaed2..908ef0e 100644 --- a/OS_Consts.py +++ b/OS_Consts.py @@ -1,3 +1,7 @@ +"""Some constants to use during go-around detection. + +These may need to be tuned to suit individual airports and/or situations. +""" # The threshold altitude used to determine if a go-around occurred, # check if points in a given time window (below) after a state change # are above this altitude @@ -35,7 +39,7 @@ # The threshold altitude for the state change, if change occurs above # this altitude then it's probably not a go-around -ga_st_alt_t = 1500. +ga_st_alt_t = 2500. # This is a list of icao24 addresses to exclude, for example general diff --git a/OS_Funcs.py b/OS_Funcs.py index 8de461b..956fa4f 100644 --- a/OS_Funcs.py +++ b/OS_Funcs.py @@ -1,6 +1,8 @@ +"""Core methods for processing ADS-B data and detecting go-arounds.""" from scipy.interpolate import UnivariateSpline as UniSpl from traffic.core import Traffic from datetime import timedelta +import metar_parse as MEP import pandas as pd import flightphase as flph @@ -9,16 +11,20 @@ import numpy as np +# Read METARs from disk +metars = MEP.get_metars('/home/proud/Desktop/GoAround_Paper/VABB_METAR') + + def estimate_rwy(df, rwy_list, verbose): - ''' - Guesses which runway a flight is attempting to land on + """Guess which runway a flight is attempting to land on. + Inputs: - df, a dict containing flight information - rwy_list, a list of runways to check, defined in OS_Airports - verbose, a bool specifying whether to verbosely print updates Returns: - A runway class from the list. - ''' + """ b_dist = 999. b_rwy = None b_pos = -1. @@ -77,13 +83,13 @@ def estimate_rwy(df, rwy_list, verbose): def get_flight(inf): - ''' - Load a series of flights from a file using Xavier's 'traffic' library + """Load a series of flights from a file using Xavier's 'traffic' library. + Input: - inf, the input filename Returns: - a list of flights - ''' + """ flist = [] # try: fdata = Traffic.from_file(inf).query("latitude == latitude") @@ -106,13 +112,13 @@ def get_flight(inf): def check_takeoff(df): - ''' - Checks if a flight is taking off. If so, we're not interested + """Check if a flight is taking off. If so, we're not interested. + Input: - df, a dict containing flight data Returns: - True if a takeoff, otherwise False - ''' + """ lf = len(df['gals']) # Check if there's enough data to process if (lf < 10): @@ -139,9 +145,8 @@ def check_takeoff(df): def get_future_time(times, c_pos, n_sec): - ''' - Find the array location that has the time closest to a - certain amount in the future. + """Find the array location that has the time closest to a certain amount in the future. + Inputs: - Times, an array of times in integer/float seconds - c_pos, the position in the array of the initial time @@ -149,7 +154,7 @@ def get_future_time(times, c_pos, n_sec): Returns: - An integer array position (>0) if time is found, or -1 if the time is not found - ''' + """ c_time = times[c_pos] f_time = c_time + n_sec diff = (np.abs(times - f_time)) @@ -161,8 +166,8 @@ def get_future_time(times, c_pos, n_sec): def check_ga(fd, verbose, first_pos=-1): - ''' - Check if a go-around occurred based on some simple tests + """Check if a go-around occurred based on some simple tests. + Inputs: - A dict containing flight data - A boolean for verbose mode. If True, a g/a warning is printed @@ -171,7 +176,7 @@ def check_ga(fd, verbose, first_pos=-1): Returns: - True if a go-around is likely to have occurred, false otherwise - An int specifying the array location where g/a occurred, or None - ''' + """ ga_flag = False bpt = None @@ -183,7 +188,7 @@ def check_ga(fd, verbose, first_pos=-1): if (labels[i] != labels[i-1]): cng[i] = True main_pts = (cng).nonzero() - if (np.all(cng == False)): + if np.all(cng is False): return ga_flag if (len(main_pts[0]) > 0): main_pts = main_pts[0] @@ -192,7 +197,7 @@ def check_ga(fd, verbose, first_pos=-1): if (pt < first_pos): continue # First check altitude of state change. G/A will be low alt - if (fd['gals'][pt] > CNS.ga_st_alt_t): + if (fd['alts'][pt] > CNS.ga_st_alt_t): continue # We only want points where initial state is Descent if (labels[pt-1] != "DE"): @@ -238,7 +243,7 @@ def check_ga(fd, verbose, first_pos=-1): vrts_p = (n_pts_vrt/n_pos)*100. if (n_pos > 10 and alts_p > 50 and vrts_p > 20): - if (verbose): + if verbose: ga_time = fd['strt'] + timedelta(seconds=int(fd['time'][pt])) print("\t-\tG/A warning:", fd['call'], @@ -250,10 +255,9 @@ def check_ga(fd, verbose, first_pos=-1): return ga_flag, bpt -def proc_fl(flight, check_rwys, odirs, colormap, metars, do_save, verbose): - ''' - The main processing routine, filters, assigns phases and determines - go-around status for a given flight. +def proc_fl(flight, check_rwys, odirs, colormap, do_save, verbose): + """Filter, assign phases and determine go-around status for a given flight. + Inputs: - A 'traffic' flight object - A list storing potential landing runways to check @@ -263,61 +267,79 @@ def proc_fl(flight, check_rwys, odirs, colormap, metars, do_save, verbose): - normal numpy data output - go-around numpy data output - A dict of colours used for flightpath labelling - - A dict of METARs used for correcting barometric altitude - A boolean specifying whether to save data or not - A boolean specifying whether to use verbose mode Returns: - Nothing - ''' - + """ + # First, check if a flight is not on exclusion list gd_fl = check_good_flight(flight) if (not gd_fl): if (verbose): print("\t-\tBad flight call:", flight.callsign) return -1 + + # Print some details if verbose if (verbose): print("\t-\tProcessing:", flight.callsign) + + # Resample trajectory to one second, this is used for runway estimation flight2 = flight.resample("1s") + + # Preprocess the data, sorting by time and putting into UNIX format fd = preproc_data(flight, verbose) fd2 = preproc_data(flight2, verbose) + + # If we don't have good data here, skip if (fd is None): if (verbose): print("\t-\tBad flight data:", flight.callsign, fd) return -1 if (fd2 is None): - if (verbose): + if verbose: print("\t-\tBad flight data:", flight.callsign, fd2) return -1 + + # We don't care about take-offs, so find and exclude takeoff = check_takeoff(fd) - if (takeoff == True): + if takeoff: return -1 + # Use Junzi's labelling method to get flight phases labels = do_labels(fd) if (np.all(labels == labels[0])): - if (verbose): + if verbose: print("\t-\tNo state change:", flight.callsign) return -1 fd['labl'] = labels + # Estimate which runway the flight is landing on (rwy), and at what + # point in the data arrays it does so (posser). rwy, posser = estimate_rwy(fd2, check_rwys, verbose) - if (rwy is None): - if (verbose): + rwy2, posser2 = estimate_rwy(fd, check_rwys, verbose) + + # If we can't estimate a runway, try using minimum altitude + # Either way, compute relative distances to the runway / minimum point + min_alt_pt = (np.nanmin(fd['alts']) == fd['alts']).nonzero() + if (len(min_alt_pt[0]) > 0): + min_alt_pt = min_alt_pt[0] + min_alt_pt = min_alt_pt[0] + if rwy is None: + if verbose: print('WARNING: Cannot find runway for flight ' + fd['call'] + ' ' + fd['ic24']) - pt = (np.nanmin(fd['alts']) == fd['alts']).nonzero() - if (len(pt[0]) > 0): - pt = pt[0] - pt = pt[0] fd['rwy'] = "None" - r_dis = np.sqrt((fd['lats'] - fd['lats'][pt]) * - (fd['lats'] - fd['lats'][pt]) + - (fd['lons'] - fd['lons'][pt]) * - (fd['lons'] - fd['lons'][pt])) + r_dis = np.sqrt((fd['lats'] - fd['lats'][min_alt_pt]) * + (fd['lats'] - fd['lats'][min_alt_pt]) + + (fd['lons'] - fd['lons'][min_alt_pt]) * + (fd['lons'] - fd['lons'][min_alt_pt])) else: fd['rwy'] = rwy.name r_dis = np.sqrt((fd['lats'] - rwy.rwy[0]) * (fd['lats'] - rwy.rwy[0]) + (fd['lons'] - rwy.rwy[1]) * (fd['lons'] - rwy.rwy[1])) + + # Convert degrees into km, not perfect but good enough r_dis = r_dis * 112. pt = (np.nanmin(r_dis) == r_dis).nonzero() if (len(pt[0]) > 0): @@ -325,7 +347,7 @@ def proc_fl(flight, check_rwys, odirs, colormap, metars, do_save, verbose): pt = pt[0] r_dis[0:pt] = r_dis[0:pt] * -1 fd['rdis'] = r_dis - + # Correct barometric altitudes t_alt = fd['alts'] l_time = fd['strt'] + (fd['dura'] / 2) @@ -338,40 +360,99 @@ def proc_fl(flight, check_rwys, odirs, colormap, metars, do_save, verbose): bmet, tdiff, l_time) fd['alts'] = t_alt + # Now the actual go-around check ga_flag, gapt = check_ga(fd, True) - - if (do_save): - spldict = create_spline(fd, bpos = gapt) + + # Make some plots if required, this needs a spline to smooth output + if do_save: + spldict = create_spline(fd, bpos=None) + # Choose output directory based upon go-around flag if (ga_flag): odir_pl = odirs[1] odir_np = odirs[3] else: odir_pl = odirs[0] odir_np = odirs[2] -# OSO.do_plots(fd, spldict, colormap, odir_pl, rwy=rwy) - if (ga_flag): - OSO.do_plots_dist(fd, spldict, colormap, odir_pl, rwy=rwy, bpos = gapt) - if (rwy != None): - garr = [ga_flag, fd['ic24'], fd['call'], l_time, rwy.name, bmet] + OSO.do_plots(fd, + spldict, + colormap, + odir_pl, + rwy=rwy, + bpos=None) + if (ga_flag): + ga_time = pd.Timestamp(fd['strt'] + + pd.Timedelta(seconds=fd['time'][gapt]), + tz='UTC') + else: + gapt = 0 + ga_time = fd['strt'] + + # Get the correct position of landing / go-around + gdpt = -1 + if (ga_flag): + if (gapt >= 10): + gdpt = gapt + elif (posser2 >= 10): + gdpt = posser2 + elif (min_alt_pt >= 10): + gdpt = min_alt_pt + else: + gdpt = -1 + + # Compute std dev for some variables immediately before g/a or landing + if (gdpt >= 10): + rocvar = np.std(fd['rocs'][gdpt-10:gdpt]) + hdgvar = np.std(fd['hdgs'][gdpt-10:gdpt]) + latvar = np.std(fd['lats'][gdpt-10:gdpt]) + lonvar = np.std(fd['lons'][gdpt-10:gdpt]) + gspvar = np.std(fd['spds'][gdpt-10:gdpt]) + timer = fd['time'][gdpt] - fd['time'][gdpt-10] + + rocvar = rocvar / timer + hdgvar = hdgvar / timer + latvar = latvar / timer + lonvar = lonvar / timer + gspvar = gspvar / timer + else: + rocvar = 0 + hdgvar = 0 + latvar = 0 + lonvar = 0 + gspvar = 0 + + # Return data to the calling function + if rwy is not None: + garr = [ga_flag, fd['ic24'], fd['call'], l_time, + ga_time, rwy.name, fd['hdgs'][gapt], + fd['alts'][gapt], fd['lats'][gapt], fd['lons'][gapt], + gapt, rocvar, hdgvar, latvar, lonvar, gspvar, + bmet] else: - garr = [ga_flag, fd['ic24'], fd['call'], l_time, 'None', bmet] -# OSO.to_numpy(fd, odir_np) + garr = [ga_flag, fd['ic24'], fd['call'], l_time, + ga_time, 'None', fd['hdgs'][gapt], + fd['alts'][gapt], fd['lats'][gapt], fd['lons'][gapt], + gapt, rocvar, hdgvar, latvar, lonvar, gspvar, + bmet] + fd['posser'] = posser2 + fd['gapt'] = gapt + fd['min_alt_pt'] = min_alt_pt + OSO.to_numpy(fd, odir_np) if (verbose): print("\t-\tDONE") return garr def check_good_flight(flight): - ''' - Checks if the flight callsign matches a series of - pre-defined 'bad' callsigns. Most of these are - ground vehicles. + """Check if the flight callsign matches a series of pre-defined 'bad' callsigns. + + Most of these are ground vehicles. Input: - A flight produced by the 'traffic' library. Returns: - False if callsign matches one of the 'bad' list - True if flight is not matched - ''' + + """ if (flight.icao24 in CNS.exclude_list): return False if (flight.callsign[0:7] == 'WILDLIF'): @@ -391,15 +472,16 @@ def check_good_flight(flight): def check_good_data(flight): - ''' - Checks that a flight has data suitable for inclusion in the study. + """Check that a flight has data suitable for inclusion in the study. + This discards ground-only flights, as well as those that do not appear to be attempting a landing. Input: - A flight produced by the 'traffic' library. Returns: - True if a flight is suitable, false otherwise. - ''' + + """ if (np.all(flight.data['geoaltitude'] > 3000)): return 'G_HIGH' if (np.all(flight.data['altitude'] > 3000)): @@ -415,8 +497,7 @@ def check_good_data(flight): def preproc_data(flight, verbose): - ''' - Preprocesses a flight into a format usable by Junzi's classifier + """Preprocesses a flight into a format usable by Junzi's classifier. Input: - A flight produced by the 'traffic' library. @@ -438,9 +519,10 @@ def preproc_data(flight, verbose): - strt: Time of first position in the flight datastream - stop: Time of last position in the flight datastream - dura: Reported duration of the flight - ''' + + """ isgd = check_good_data(flight) - if (isgd != True): + if (not isgd): if (verbose): print("Unsuitable flight:", flight.callsign, isgd) return None @@ -473,6 +555,22 @@ def preproc_data(flight, verbose): fdata['hdgs'] = hdgs fdata['rocs'] = f_data['vertical_rate'].values fdata['ongd'] = f_data['onground'].values + + # The next bit is needed in case a flight crosses two pkl files, which are + # usually one hour long. So a flight going from 00:59 -> 01:00 is in two + # files, and due to multiprocessing the two segments may be reversed. + df_tmp = pd.DataFrame(fdata) + df_new = df_tmp.sort_values(by=['time']) + + fdata['time'] = df_new['time'].values + fdata['lats'] = df_new['lats'].values + fdata['lons'] = df_new['lons'].values + fdata['alts'] = df_new['alts'].values + fdata['spds'] = df_new['spds'].values + fdata['gals'] = df_new['gals'].values + fdata['hdgs'] = df_new['hdgs'].values + fdata['rocs'] = df_new['rocs'].values + fdata['ongd'] = df_new['ongd'].values fdata['call'] = flight.callsign fdata['ic24'] = flight.icao24 fdata['strt'] = flight.start @@ -483,15 +581,15 @@ def preproc_data(flight, verbose): def find_closest_metar(l_time, metars): - ''' - Finds the best-fitting metar from a dict that matches a specified - time value + """Find the best-fitting metar from a dict that matches a specified time value. + Inputs: - The time to match (datetime) - A dict of METARS, each as a metobs class Returns: The best metar (as metobs) and the time difference in seconds - ''' + + """ tdiff = 1e8 bmet = None @@ -505,15 +603,16 @@ def find_closest_metar(l_time, metars): def correct_baro(balt, t0, p0): - ''' - A function to correct barometric altitude values from ISA to actual + """Correct barometric altitude values from ISA to actual. + Inputs: - balt: An array of baro altitudes - t0: A float with the surface temperature (C) - p0: A float with the surface pressure (hPa) Returns: - An array of corrected baro alts - ''' + + """ isa_t = 15.0 isa_p = 1013.25 # First, compute ISA pressure @@ -530,10 +629,9 @@ def correct_baro(balt, t0, p0): return alt -def create_spline(fd, bpos = None): - ''' - Creates the splines needed for plotting smoothed lines - on the output graphs +def create_spline(fd, bpos=None): + """Create the splines needed for plotting smoothed lines on the output graphs. + Input: - A dict of flight data, such as that returned by preproc_data() - An int speicfying the max array value to use @@ -544,34 +642,46 @@ def create_spline(fd, bpos = None): - rocspl - galspl - hdgspl - ''' + + """ spldict = {} - if (bpos == None): + if (bpos is None): bpos = len(fd['time']) - spldict['altspl'] = UniSpl(fd['time'][0: bpos], fd['alts'][0: bpos])(fd['time'][0: bpos]) - spldict['spdspl'] = UniSpl(fd['time'][0: bpos], fd['spds'][0: bpos])(fd['time'][0: bpos]) - spldict['rocspl'] = UniSpl(fd['time'][0: bpos], fd['rocs'][0: bpos])(fd['time'][0: bpos]) - spldict['galspl'] = UniSpl(fd['time'][0: bpos], fd['gals'][0: bpos])(fd['time'][0: bpos]) - spldict['hdgspl'] = UniSpl(fd['time'][0: bpos], fd['hdgs'][0: bpos])(fd['time'][0: bpos]) - spldict['latspl'] = UniSpl(fd['time'][0: bpos], fd['lats'][0: bpos])(fd['time'][0: bpos]) - spldict['lonspl'] = UniSpl(fd['time'][0: bpos], fd['lons'][0: bpos])(fd['time'][0: bpos]) + spldict['altspl'] = UniSpl(fd['time'][0: bpos], + fd['alts'][0: bpos])(fd['time'][0: bpos]) + spldict['spdspl'] = UniSpl(fd['time'][0: bpos], + fd['spds'][0: bpos])(fd['time'][0: bpos]) + spldict['rocspl'] = UniSpl(fd['time'][0: bpos], + fd['rocs'][0: bpos])(fd['time'][0: bpos]) + spldict['galspl'] = UniSpl(fd['time'][0: bpos], + fd['gals'][0: bpos])(fd['time'][0: bpos]) + spldict['hdgspl'] = UniSpl(fd['time'][0: bpos], + fd['hdgs'][0: bpos])(fd['time'][0: bpos]) + spldict['latspl'] = UniSpl(fd['time'][0: bpos], + fd['lats'][0: bpos])(fd['time'][0: bpos]) + spldict['lonspl'] = UniSpl(fd['time'][0: bpos], + fd['lons'][0: bpos])(fd['time'][0: bpos]) return spldict def do_labels(fd): - ''' - Perform the fuzzy labelling using Junzi's method. - Add an additional force label of aircraft with "onground=True" to + """Perform the fuzzy labelling using Junzi's method. + + Add an additional force label of aircraft with 'onground=True' to the 'GND' label category. Input: - A dict of flight data, such as that returned by preproc_data() Returns: - A numpy array containing categorised flight phases. - ''' - labels = flph.fuzzylabels(fd['time'], fd['gals'], - fd['spds'], fd['rocs'], twindow=15) + """ + try: + labels = flph.fuzzylabels(fd['time'], fd['alts'], + fd['spds'], fd['rocs'], twindow=15) + except Exception as e: + print("Error creating spline", e, fd['call']) + quit() pts = (fd['ongd']).nonzero() labels[pts] = 'GND' diff --git a/OS_Output.py b/OS_Output.py index 447cb9c..8dde7b0 100644 --- a/OS_Output.py +++ b/OS_Output.py @@ -1,14 +1,19 @@ -from matplotlib.lines import Line2D -import matplotlib.pyplot as plt +"""A set of functions to plot and/or save flight trajectory information.""" import numpy as np import os +# This line stops matplotlib messing up in terminal mode +os.environ['QT_QPA_PLATFORM'] = 'offscreen' + +from matplotlib.lines import Line2D +import matplotlib.pyplot as plt + + +def do_plots(fd, spld, cmap, outdir, app_ylim=True, + odpi=300, rwy=None, bpos=None): + """Creates and saves a series of plots showing relevant data for each flight that has been processed. -def do_plots(fd, spld, cmap, outdir, app_ylim=True, odpi=300, rwy=None): - ''' - Creates and saves a series of plots showing relevant data for each - flight that has been processed. Files are saved in a /YYYMMDD/ - subdirectory of the 'outdir' argument. + Files are saved in a /YYYMMDD/ subdirectory of the 'outdir' argument. This version saves data with time since first appearance in the datastream on the x-axis. Inputs: @@ -19,38 +24,83 @@ def do_plots(fd, spld, cmap, outdir, app_ylim=True, odpi=300, rwy=None): - (optional) A bool specifying to apply predefined axis limits. - (optional) An int specifying the desired output DPI - (optional) A runway class specifying the landing runway, or None + - (optional) An int specifying the max array position Returns: - Nothing - ''' - colors = [cmap[l] for l in fd['labl']] - fig, ax = plt.subplots(dpi=400) + """ + if bpos is None: + bpos = len(fd['time']) + + colors = [cmap[l] for l in fd['labl'][0:bpos]] fs = (20, 20) custom_lines = [] for color in cmap: lin = Line2D([0], [0], color=cmap[color], lw=0, marker='.') custom_lines.append(lin) - plt.subplot(311) - plt.plot(fd['time'], spld['galspl']/1000., '-', color='k', lw=0.1) - plt.scatter(fd['time'], fd['gals']/1000., marker='.', c=colors, lw=0) + fig, ax = plt.subplots(dpi=400) + fs = (20, 20) + + plt.subplot(411) + plt.plot(fd['time'][0:bpos], + spld['altspl']/1000., + '-', + color='k', + lw=0.1) + plt.scatter(fd['time'][0:bpos], + fd['alts'][0:bpos]/1000., + marker='.', + c=colors, + lw=0) plt.ylabel('altitude (kft)') if (app_ylim): plt.ylim(-0.5, 10.) - plt.subplot(312) - plt.plot(fd['time'], spld['rocspl']/1000., '-', color='k', lw=0.1) - plt.scatter(fd['time'], fd['rocs']/1000., marker='.', c=colors, lw=0) + plt.subplot(412) + plt.plot(fd['time'][0:bpos], + spld['rocspl']/1000., + '-', + color='k', + lw=0.1) + plt.scatter(fd['time'][0:bpos], + fd['rocs'][0:bpos]/1000., + marker='.', + c=colors, + lw=0) plt.ylabel('roc (kfpm)') if (app_ylim): - plt.ylim(-0.5, 10.) + plt.ylim(-1.5, 1.5) - plt.subplot(313) - plt.plot(fd['time'], spld['spdspl'], '-', color='k', lw=0.1) - plt.scatter(fd['time'], fd['spds'], marker='.', c=colors, lw=0) + plt.subplot(413) + plt.plot(fd['time'][0:bpos], + spld['spdspl'], + '-', + color='k', + lw=0.1) + plt.scatter(fd['time'][0:bpos], + fd['spds'][0:bpos], + marker='.', + c=colors, + lw=0) plt.ylabel('speed (kts)') if (app_ylim): plt.ylim(0., 400.) + plt.subplot(414) + plt.plot(fd['time'], + spld['hdgspl'], + '-', + color='k', + lw=0.1) + plt.scatter(fd['time'], + fd['hdgs'], + marker='.', + c=colors, + lw=0) + plt.ylabel('heading (deg)') + if (app_ylim): + plt.ylim(-180., 180.) + plt.legend(custom_lines, ['Ground', 'Climb', 'Cruise', 'Descent', 'Level', 'N/A'], bbox_to_anchor=(0., -0.33, 1., 0.102), @@ -64,8 +114,8 @@ def do_plots(fd, spld, cmap, outdir, app_ylim=True, odpi=300, rwy=None): if (not os.path.exists(odir)): try: os.mkdir(odir) - except: - None + except Exception: + pass timestr = fd['stop'].strftime("%Y%m%d%H%M") outf = odir + 'FLT_' + fd['ic24'] + '_' @@ -76,16 +126,17 @@ def do_plots(fd, spld, cmap, outdir, app_ylim=True, odpi=300, rwy=None): dpi=odpi, bbox_inches='tight', pad_inches=0) - plt.close() + plt.close(fig) def get_fig_outname(outdir, fd, figtype): + """Generates a name for the output figure based on type.""" odir = outdir + fd['stop'].strftime("%Y%m%d") + '/' if (not os.path.exists(odir)): try: os.mkdir(odir) - except: - None + except Exception: + pass timestr = fd['stop'].strftime("%Y%m%d%H%M") outf = odir + 'FLT_' + fd['ic24'] + '_' @@ -93,17 +144,17 @@ def get_fig_outname(outdir, fd, figtype): outf = outf + timestr + '_' + figtype + '.png' return outf - + + def make_yvals(dists, multis): - ''' - Return a list of y-values associated with a 6th order polynomial - and some x values. + """Return a list of y-values associated with a 6th order polynomial and some x values. + Inputs: - dists: List of x axis coords - multis: Polynomial coefficientsd Returns: - A list of y values - ''' + """ yvals = (np.power(dists, 6) * multis[0] + np.power(dists, 5) * multis[1] + np.power(dists, 4) * multis[2] + @@ -114,12 +165,13 @@ def make_yvals(dists, multis): return yvals + def do_plots_dist(fd, spld, cmap, outdir, - app_xlim=True, app_ylim=True, odpi=300, rwy=None, bpos = None): - ''' - Creates and saves a series of plots showing relevant data for each - flight that has been processed. Files are saved in a /YYYMMDD/ - subdirectory of the 'outdir' argument. + app_xlim=True, app_ylim=False, + odpi=300, rwy=None, bpos=None): + """Creates and saves a series of plots showing relevant data for each flight that has been processed. + + Files are saved in a /YYYMMDD/ subdirectory of the 'outdir' argument. This version saves data with distance to detected landing runway on the x-axis. Inputs: @@ -134,9 +186,8 @@ def do_plots_dist(fd, spld, cmap, outdir, - (optional) An int specifying the max array position Returns: - Nothing - ''' - - if (bpos == None): + """ + if bpos is None: bpos = len(fd['time']) xlims = [-10., 0.1] @@ -146,15 +197,14 @@ def do_plots_dist(fd, spld, cmap, outdir, else: hme = np.nanmean(fd['hdgs']) hdglim = [hme - 5, hme + 5] - + colors = [cmap[l] for l in fd['labl'][0:bpos]] fs = (20, 20) custom_lines = [] for color in cmap: lin = Line2D([0], [0], color=cmap[color], lw=0, marker='.') custom_lines.append(lin) - - + tmprd = np.copy(fd['rdis']) pts = (tmprd < xlims[0]).nonzero() tmprd[pts] = np.nan @@ -162,8 +212,11 @@ def do_plots_dist(fd, spld, cmap, outdir, tmprd[pts] = np.nan pts = (tmprd == tmprd).nonzero() -# plt.plot(fd['rdis'], spld['altspl']/1000., '-', color='k', lw=0.1) - plt.scatter(fd['rdis'][0:bpos], fd['alts'][0:bpos]/1000., marker='.', c=colors, lw=0) + plt.scatter(fd['rdis'][0:bpos], + fd['alts'][0:bpos]/1000., + marker='.', + c=colors, + lw=0) if (rwy is not None): yvals1 = make_yvals(distlist, rwy.alts1) yvalm = make_yvals(distlist, rwy.altm) @@ -183,7 +236,6 @@ def do_plots_dist(fd, spld, cmap, outdir, plt.legend(custom_lines, ['Ground', 'Climb', 'Cruise', 'Descent', 'Level', 'N/A'], loc='best', - #bbox_to_anchor=(0.5, -0.01), ncol=6, fancybox=False, mode="expand") @@ -192,8 +244,11 @@ def do_plots_dist(fd, spld, cmap, outdir, plt.savefig(outf, figsize=fs, dpi=odpi, bbox_inches='tight', pad_inches=0) plt.clf() -# plt.plot(fd['rdis'], spld['rocspl']/1000., '-', color='k', lw=0.1) - plt.scatter(fd['rdis'][0:bpos], fd['rocs'][0:bpos]/1000., marker='.', c=colors, lw=0) + plt.scatter(fd['rdis'][0:bpos], + fd['rocs'][0:bpos]/1000., + marker='.', + c=colors, + lw=0) if (rwy is not None): yvals1 = make_yvals(distlist, rwy.rocs1) yvalm = make_yvals(distlist, rwy.rocm) @@ -213,7 +268,6 @@ def do_plots_dist(fd, spld, cmap, outdir, plt.legend(custom_lines, ['Ground', 'Climb', 'Cruise', 'Descent', 'Level', 'N/A'], loc='best', - #bbox_to_anchor=(0.5, -0.01), ncol=6, fancybox=False, mode="expand") @@ -222,8 +276,11 @@ def do_plots_dist(fd, spld, cmap, outdir, plt.savefig(outf, figsize=fs, dpi=odpi, bbox_inches='tight', pad_inches=0) plt.clf() -# plt.plot(fd['rdis'], spld['spdspl'], '-', color='k', lw=0.1) - plt.scatter(fd['rdis'][0:bpos], fd['spds'][0:bpos], marker='.', c=colors, lw=0) + plt.scatter(fd['rdis'][0:bpos], + fd['spds'][0:bpos], + marker='.', + c=colors, + lw=0) plt.ylabel('Ground Speed (kts)') plt.xlabel('Distance to Runway (km)') if (app_xlim): @@ -236,7 +293,6 @@ def do_plots_dist(fd, spld, cmap, outdir, plt.legend(custom_lines, ['Ground', 'Climb', 'Cruise', 'Descent', 'Level', 'N/A'], loc='best', - #bbox_to_anchor=(0.5, -0.01), ncol=6, fancybox=False, mode="expand") @@ -245,8 +301,11 @@ def do_plots_dist(fd, spld, cmap, outdir, plt.savefig(outf, figsize=fs, dpi=odpi, bbox_inches='tight', pad_inches=0) plt.clf() -# plt.plot(fd['rdis'], spld['hdgspl'], '-', color='k', lw=0.1) - plt.scatter(fd['rdis'][0:bpos], fd['hdgs'][0:bpos], marker='.', c=colors, lw=0) + plt.scatter(fd['rdis'][0:bpos], + fd['hdgs'][0:bpos], + marker='.', + c=colors, + lw=0) if (rwy is not None): yvals1 = make_yvals(distlist, rwy.hdgs1) yvalm = make_yvals(distlist, rwy.hdgm) @@ -266,7 +325,6 @@ def do_plots_dist(fd, spld, cmap, outdir, plt.legend(custom_lines, ['Ground', 'Climb', 'Cruise', 'Descent', 'Level', 'N/A'], loc='best', - #bbox_to_anchor=(0.5, -0.01), ncol=6, fancybox=False, mode="expand") @@ -275,8 +333,11 @@ def do_plots_dist(fd, spld, cmap, outdir, plt.savefig(outf, figsize=fs, dpi=odpi, bbox_inches='tight', pad_inches=0) plt.clf() -# plt.plot(fd['rdis'], spld['lonspl'], '-', color='k', lw=0.1) - plt.scatter(fd['rdis'][0:bpos], fd['lons'][0:bpos], marker='.', c=colors, lw=0) + plt.scatter(fd['rdis'][0:bpos], + fd['lons'][0:bpos], + marker='.', + c=colors, + lw=0) if (rwy is not None): yvals1 = make_yvals(distlist, rwy.lons1) yvalm = make_yvals(distlist, rwy.lonm) @@ -295,7 +356,6 @@ def do_plots_dist(fd, spld, cmap, outdir, plt.legend(custom_lines, ['Ground', 'Climb', 'Cruise', 'Descent', 'Level', 'N/A'], loc='best', - #bbox_to_anchor=(0.5, -0.01), ncol=6, fancybox=False, mode="expand") @@ -304,8 +364,11 @@ def do_plots_dist(fd, spld, cmap, outdir, plt.savefig(outf, figsize=fs, dpi=odpi, bbox_inches='tight', pad_inches=0) plt.clf() -# plt.plot(fd['rdis'], spld['latspl'], '-', color='k', lw=0.1) - plt.scatter(fd['rdis'][0:bpos], fd['lats'][0:bpos], marker='.', c=colors, lw=0) + plt.scatter(fd['rdis'][0:bpos], + fd['lats'][0:bpos], + marker='.', + c=colors, + lw=0) if (rwy is not None): yvals1 = make_yvals(distlist, rwy.lats1) yvalm = make_yvals(distlist, rwy.latm) @@ -324,7 +387,6 @@ def do_plots_dist(fd, spld, cmap, outdir, plt.legend(custom_lines, ['Ground', 'Climb', 'Cruise', 'Descent', 'Level', 'N/A'], loc='best', - #bbox_to_anchor=(0.5, -0.01), ncol=6, fancybox=False, mode="expand") @@ -335,21 +397,21 @@ def do_plots_dist(fd, spld, cmap, outdir, def to_numpy(fd, outdir): - ''' - Save data for a single flight into a numpy pickle file. + """Save data for a single flight into a numpy pickle file. + Files are saved in a YYYYMMDD subdirectory. Inputs: - fd: Dict containing flight info - outdir: Location to store output Returns: - Nothing - ''' + """ odir = outdir + fd['stop'].strftime("%Y%m%d") + '/' if (not os.path.exists(odir)): try: os.mkdir(odir) - except: - None + except Exception: + pass outf = odir + 'FLT_' + fd['ic24'] + '_' outf = outf + fd['call'] + '_' outf = outf + fd['stop'].strftime("%Y%m%d%H%m") + '.pkl' diff --git a/OpenSky_Get_Data.py b/OpenSky_Get_Data.py index cc30425..33705f1 100644 --- a/OpenSky_Get_Data.py +++ b/OpenSky_Get_Data.py @@ -1,8 +1,7 @@ -''' +""" This script downloads data from the opensky library for a particular airport, a small perimeter is set up around the airport to catch the approach path - -''' +""" from datetime import datetime, timedelta from traffic.data import opensky @@ -29,12 +28,12 @@ def get_bounds(rwys): - ''' + """ This function computes the boundaries of the retrieval 'box' based upon the runways selected for processing. The box is approx. 0.25 degrees in each direction around the airport midpoint. - ''' + """ latlist = [] lonlist = [] for rwy in rwys: @@ -52,12 +51,11 @@ def get_bounds(rwys): def getter(init_time, bounds, timer, anam): - ''' + """ This function downloads the data, which is done in one hour segments. Each hour is downloaded separately using multiprocessing for efficiency. - ''' - + """ try: times = init_time + timedelta(hours=timer) dtst = times.strftime("%Y%m%d%H%M") diff --git a/metar_parse.py b/metar_parse.py index a93e69c..d556f91 100644 --- a/metar_parse.py +++ b/metar_parse.py @@ -1,11 +1,17 @@ +"""A script for processing METAR data. + +This uses the metar library available at: +https://github.com/python-metar/python-metar +""" + from datetime import datetime from metar import Metar import pytz class metobs: - ''' - A class to store METAR observations: + """A class to store METAR observations. + temp = temperature in C dewp = dewpoint in C w_s = wind speed in kts @@ -14,8 +20,11 @@ class metobs: cb = boolean specifying CBs near airport vis = visibility in km pres = pressure in hPa - ''' - def __init__(self, temp, dewp, w_s, w_d, w_g, cb, vis, pres): + cld = cloud base in ft + """ + + def __init__(self, temp, dewp, w_s, w_d, w_g, cb, vis, pres, cld): + """Setup the class.""" self.temp = temp self.dewp = dewp self.w_s = w_s @@ -24,17 +33,17 @@ def __init__(self, temp, dewp, w_s, w_d, w_g, cb, vis, pres): self.cb = cb self.vis = vis self.pres = pres + self.cld = cld -def get_metars(inf): - ''' - A function to parse metars from a file and convert into a - dict of metobs objets +def get_metars(inf, verbose): + """A function to parse metars from a file and convert into a dict of metobs objects. Input: - inf: The input file (as a string filename) Output: - a dict of metobs read from the file - ''' + """ + fid = open(inf, 'r') met_dict = {} @@ -48,31 +57,45 @@ def get_metars(inf): obs.time = metdate try: temp = obs.temp.value() - except: + except Exception as e: + if verbose: + print("ERROR: Bad temperature data", e) temp = 15 try: - dewp = obs.temp.value() - except: + dewp = obs.dewpt.value() + except Exception as e: + if verbose: + print("ERROR: Bad dewpoint data", e) dewp = 10 try: w_s = obs.wind_speed.value() - except: + except Exception as e: + if verbose: + print("ERROR: Bad wind speed data", e) w_s = 0 try: w_g = obs.wind_gust.value() - except: + except Exception as e: + if verbose: + print("ERROR: Bad wind gust data", e) w_g = 0 try: w_d = obs.wind_dir.value() - except: + except Exception as e: + if verbose: + print("ERROR: Bad wind direction data", e) w_d = 0 try: press = obs.press.value() - except: + except Exception as e: + if verbose: + print("ERROR: Bad pressure data", e) press = 1013.25 try: vis = obs.vis.value() - except: + except Exception as e: + if verbose: + print("ERROR: Bad visibility data", e) vis = 10000 cb = False @@ -80,11 +103,16 @@ def get_metars(inf): for wx in obs.sky: if (wx[2] == "CB"): cb = True - except: - None + except Exception as e: + if verbose: + print("ERROR: Bad CB data", e) - cur_obs = metobs(temp, dewp, w_s, w_d, w_g, cb, vis, press) + cld = 10000 + if (len(obs.sky) > 0): + if obs.sky[0][1] is not None: + cld = obs.sky[0][1].value() + cur_obs = metobs(temp, dewp, w_s, w_d, w_g, cb, vis, press, cld) met_dict[metdate] = cur_obs return met_dict