New retrieve_runway_list dervices runway data from OurAirports db

GA_Detect.py

@@ -2,17 +2,30 @@
 from traffic.core import Traffic
 from datetime import timedelta
 import multiprocessing as mp
-from OS_Airports import VABB
+import metar_parse as MEP
+from OS_Airports.RWY import get_runway_list
 import OS_Funcs as OSF
+import os
 import glob
-
-
-def main(start_n, fidder, do_write):
+import click
+
+
+@click.command()
+@click.option('--top-dir', default='./')
+@click.option('--start-n', default=0)
+@click.option('--do-write/--no-write', default=True)
+@click.option('--do-plot/--no-plot', default=True)
+@click.option('--metars-file', default='VABB_METAR')
+@click.option('--airport', default='VABB')
+@click.option('--n-files-proc', default=55)
+@click.option('--pool-proc', default=16)
+@click.option('--verbose', default=False)
+def main(top_dir, start_n, do_write, do_plot, metars_file, airport,
+         n_files_proc, pool_proc, verbose):
     """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
 
     """
@@ -22,7 +35,6 @@ def main(start_n, fidder, do_write):
     # Of which go-arounds
     tot_n_ga = 0
 
-    top_dir = '/gf2/eodg/SRP002_PROUD_ADSBREP/GO_AROUNDS/VABB/'
     # indir stores the opensky data
     indir = top_dir + 'INDATA/'
 
@@ -39,11 +51,13 @@ def main(start_n, fidder, do_write):
     odir_da_ga = top_dir + 'OUT_DATA/PSGA/'
 
     odirs = [odir_pl_nm, odir_pl_ga, odir_da_nm, odir_da_ga]
+    for odir in odirs:
+        os.makedirs(odir, exist_ok=True)
 
     # Output filenames for saving data about go-arounds
-    out_file_ga = 'GA_MET_NEW.csv'
+    out_file_ga = top_dir + 'GA_MET_NEW.csv'
     # Output filenames for saving data about non-go-arounds
-    out_file_noga = 'GA_NOGA_NEW.csv'
+    out_file_noga = top_dir + 'GA_NOGA_NEW.csv'
 
     t_frmt = "%Y/%m/%d %H:%M:%S"
 
@@ -60,18 +74,16 @@ def main(start_n, fidder, do_write):
                       Temp, Dewp, Wind_Spd, Wind_Gust, Wind_Dir,Cld_Base,\
                       CB, Vis, Pressure\n')
     files = []
-    files = glob.glob(indir+'**.pkl')
+    files = glob.glob(indir+'*.pkl') + glob.glob(indir+'*/*.pkl')
     files.sort()
 
     fli_len = len(files)
 
-    colormap = {'GND': 'black', 'CL': 'green', 'CR': 'blue',
+    colormap = {'GND': 'black', 'GN': 'black', 'CL': 'green', 'CR': 'blue',
                 'DE': 'orange', 'LVL': 'purple', 'NA': 'red'}
 
-    # Number of files to open in one go
-    n_files_proc = 55
-
-    pool_proc = 100
+    metars = MEP.get_metars(metars_file, verbose=verbose)
+    rwy_list = get_runway_list(airport)
 
     f_data = []
     pool = mp.Pool(processes=pool_proc)
@@ -80,9 +92,6 @@ def main(start_n, fidder, do_write):
         print("Processing batch starting with "
               + str(main_count + 1).zfill(5) + " of "
               + str(fli_len).zfill(5))
-        fidder.write("Processing batch starting with "
-                     + str(main_count + 1).zfill(5) + " of "
-                     + str(fli_len).zfill(5) + '\n')
 
         p_list = []
         # First we load several files at once
@@ -113,11 +122,12 @@ def main(start_n, fidder, do_write):
             if (flight.stop + timedelta(minutes=5) < end_time):
                 p_list.append(pool.apply_async(OSF.proc_fl,
                                                args=(flight,
-                                                     VABB.rwy_list,
+                                                     metars,
+                                                     rwy_list,
                                                      odirs,
                                                      colormap,
-                                                     True,
-                                                     False,)))
+                                                     do_plot,
+                                                     verbose,)))
             else:
                 f_data.append(flight)
 
@@ -196,12 +206,5 @@ def main(start_n, fidder, do_write):
         nogfid.close()
 
 
-# Use this to start processing from a given file number.
-# Can be helpful if processing fails at some point.
-init_num = 0
-
-fid = open('/home/proud/Desktop/log.log', 'w')
-
-main(init_num, fid, False)
-
-fid.close()
+if __name__ == '__main__':
+    main()

Metars_Get_Data.py

@@ -0,0 +1,30 @@
+import click
+import requests
+
+METAR_URL_TEMPLATE = 'https://mesonet.agron.iastate.edu/cgi-bin/request/asos.py?station={station}&data=metar&year1={year1}&month1={month1}&day1={day1}&year2={year2}&month2={month2}&day2={day2}&tz=Etc%2FUTC&format=onlycomma&latlon=no&missing=M&trace=T&direct=no&report_type=1&report_type=2'
+
+
+@click.command()
+@click.option('--station', default='VABB')
+@click.option('--start-dt', default='2019-08-10')
+@click.option('--end-dt', default='2019-08-21')
+@click.option('--outfile', default=None)
+def main(station, start_dt, end_dt, outfile):
+    if outfile is None:
+        outfile = f'{station}_METAR'
+    year1, month1, day1 = start_dt.split('-')
+    year2, month2, day2 = end_dt.split('-')
+    url = METAR_URL_TEMPLATE.format(
+        station=station,
+        year1=year1, month1=month1, day1=day1,
+        year2=year2, month2=month2, day2=day2,
+        )
+    resp = requests.get(url)
+    resp.raise_for_status()
+
+    with open(outfile, 'w') as fout:
+        fout.writelines(line + '\n' for line in resp.text.splitlines()[1:])
+
+
+if __name__ == '__main__':
+    main()

OS_Airports/RWY.py

@@ -1,3 +1,19 @@
+import importlib
+import os
+
+from numpy import (
+    arcsin,
+    arctan2,
+    cos,
+    degrees,
+    pi,
+    radians,
+    sin,
+    )
+import pandas
+import requests
+
+
 class rwy_data:
     ''' Defines a new runway for an airport. Takes the form:
     Name: Name of the runway, i.e: '01L'
@@ -15,12 +31,12 @@ class rwy_data:
     polynomial fit: y = f0 * x^6 + f1 * x^5 ... f6
     '''
     def __init__(self, name, mainhdg, heading, rwy, rwy2, gate,
-                 lons1, lonm, lonp1,
-                 lats1, latm, latp1,
-                 hdgs1, hdgm, hdgp1,
-                 gals1, galm, galp1,
-                 alts1, altm, altp1,
-                 rocs1, rocm, rocp1):
+                 lons1=None, lonm=None, lonp1=None,
+                 lats1=None, latm=None, latp1=None,
+                 hdgs1=None, hdgm=None, hdgp1=None,
+                 gals1=None, galm=None, galp1=None,
+                 alts1=None, altm=None, altp1=None,
+                 rocs1=None, rocm=None, rocp1=None):
 
         self.name = name
         self.mainhdg = mainhdg
@@ -47,3 +63,87 @@ def __init__(self, name, mainhdg, heading, rwy, rwy2, gate,
         self.rocm = rocm
         self.rocp1 = rocp1
 
+    def __eq__(self, other):
+        return vars(self) == vars(other)
+
+    def __repr__(self):
+        return (
+            'rwy_data(' +
+            ', '.join([
+                f'{key}={val!r}'
+                for key, val in vars(self).items() if val is not None
+                ]) +
+            ')'
+            )
+
+
+def gate(rwy, rwy2, dist=2.9e3, R=6371e3):
+    φ1, λ1 = radians(rwy[0]), radians(rwy[1])
+    φ2, λ2 = radians(rwy2[0]), radians(rwy2[1])
+    Δλ = λ2 - λ1
+    θ = arctan2(sin(Δλ) * cos(φ2), cos(φ1) * sin(φ2) -
+                sin(φ1) * cos(φ2) * cos(Δλ)) + pi
+    φ3 = arcsin(sin(φ1) * cos(dist/R) +
+                cos(φ1) * sin(dist/R) * cos(θ))
+    λ3 = λ1 + arctan2(sin(θ) * sin(dist/R) * cos(φ1),
+                      cos(dist/R) - sin(φ1) * sin(φ2))
+    return [degrees(φ3), degrees(λ3)]
+
+
+def retrieve_runway_list(name, runways_csv='runways.csv'):
+    result = []
+    if not os.path.exists(runways_csv):
+        text = requests.get('https://ourairports.com/data/runways.csv').text
+        with open(runways_csv, 'w') as f:
+            f.write(text)
+    runways = pandas.read_csv(runways_csv)
+    runways = runways.query(f'airport_ident == "{name}"')
+    for idx, row in runways.iterrows():
+        le_heading = row.le_heading_degT
+        he_heading = row.he_heading_degT
+        if le_heading > 180:
+            le_heading = -(360 - le_heading)
+        if he_heading > 180:
+            he_heading = -(360 - he_heading)
+        result.append(rwy_data(
+            name=row.le_ident,
+            mainhdg=row.le_heading_degT,
+            heading=[
+                le_heading-10,
+                le_heading,
+                le_heading,
+                le_heading+10,
+                ],
+            rwy=[row.le_latitude_deg, row.le_longitude_deg],
+            rwy2=[row.he_latitude_deg, row.he_longitude_deg],
+            gate=gate(
+                [row.le_latitude_deg, row.le_longitude_deg],
+                [row.he_latitude_deg, row.he_longitude_deg],
+                ),
+            ))
+        result.append(rwy_data(
+            name=row.he_ident,
+            mainhdg=row.he_heading_degT,
+            heading=[
+                he_heading-10,
+                he_heading,
+                he_heading,
+                he_heading+10,
+                ],
+            rwy=[row.he_latitude_deg, row.he_longitude_deg],
+            rwy2=[row.le_latitude_deg, row.le_longitude_deg],
+            gate=gate(
+                [row.he_latitude_deg, row.he_longitude_deg],
+                [row.le_latitude_deg, row.le_longitude_deg],
+                ),
+            ))
+    return result
+
+
+def get_runway_list(name):
+    try:
+        module = importlib.import_module(f'OS_Airports.{name}')
+    except ModuleNotFoundError:
+        return retrieve_runway_list(name)
+    else:
+        return getattr(module, 'rwy_list')

OS_Funcs.py

@@ -2,7 +2,6 @@
 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
@@ -11,10 +10,6 @@
 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):
     """Guess which runway a flight is attempting to land on.
 
@@ -255,11 +250,12 @@ def check_ga(fd, verbose, first_pos=-1):
     return ga_flag, bpt
 
 
-def proc_fl(flight, check_rwys, odirs, colormap, do_save, verbose):
+def proc_fl(flight, metars, 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 dict of METARS, each as a metobs class
         -   A list storing potential landing runways to check
         -   A 4-element list specifying various output directories:
             -   normal plot output
@@ -351,7 +347,7 @@ def proc_fl(flight, check_rwys, odirs, colormap, do_save, verbose):
     # Correct barometric altitudes
     t_alt = fd['alts']
     l_time = fd['strt'] + (fd['dura'] / 2)
-    l_time = pd.Timestamp(l_time, tz='UTC')
+    l_time = pd.Timestamp(l_time)
     bmet, tdiff = find_closest_metar(l_time, metars)
     if (bmet is not None):
         t_alt = correct_baro(t_alt, bmet.temp, bmet.pres)
@@ -363,16 +359,17 @@ def proc_fl(flight, check_rwys, odirs, colormap, do_save, verbose):
     # Now the actual go-around check
     ga_flag, gapt = check_ga(fd, True)
 
+    if (ga_flag):
+        odir_pl = odirs[1]
+        odir_np = odirs[3]
+    else:
+        odir_pl = odirs[0]
+        odir_np = odirs[2]
+
     # 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,
@@ -381,8 +378,7 @@ def proc_fl(flight, check_rwys, odirs, colormap, do_save, verbose):
                      bpos=None)
     if (ga_flag):
         ga_time = pd.Timestamp(fd['strt'] +
-                               pd.Timedelta(seconds=fd['time'][gapt]),
-                               tz='UTC')
+                               pd.Timedelta(seconds=fd['time'][gapt]))
     else:
         gapt = 0
         ga_time = fd['strt']
@@ -679,6 +675,7 @@ def do_labels(fd):
     try:
         labels = flph.fuzzylabels(fd['time'], fd['alts'],
                                   fd['spds'], fd['rocs'], twindow=15)
+        labels = np.array(labels)
     except Exception as e:
         print("Error creating spline", e, fd['call'])
         quit()