| |
| """ |
| Compute ground truth for Table 2 test cases. |
| Dataset: step 630 = 16:30 UTC, replay steps 600-719. |
| """ |
| import pandas as pd |
| import json |
| import math |
| from pathlib import Path |
|
|
| ROOT = Path(__file__).resolve().parent.parent |
| PARQUET = ROOT / "data/opensky_snapshot/states.parquet" |
| ROUTES = ROOT / "data/opensky_snapshot/routes.json" |
| AIRPORTS = ROOT / "airports.csv" |
| STEP = 630 |
|
|
|
|
| def get_route_for_step(routes, icao24, step): |
| r = routes.get(icao24, {}) |
| segs = r.get('segments', []) |
| for s in segs: |
| if s['startStep'] <= step <= s['endStep']: |
| org = s.get('origin', {}).get('iata', '?') if s.get('origin') else '?' |
| dst = s.get('destination', {}).get('iata', '?') if s.get('destination') else '?' |
| org_icao = s.get('origin', {}).get('icao', '?') if s.get('origin') else '?' |
| dst_icao = s.get('destination', {}).get('icao', '?') if s.get('destination') else '?' |
| return org, dst, org_icao, dst_icao, s.get('routeType', '?') |
| return '?', '?', '?', '?', '?' |
|
|
|
|
| def haversine_km(lat1, lon1, lat2, lon2): |
| R = 6371 |
| phi1, phi2 = math.radians(lat1), math.radians(lat2) |
| dphi = math.radians(lat2 - lat1) |
| dlam = math.radians(lon2 - lon1) |
| a = math.sin(dphi/2)**2 + math.cos(phi1)*math.cos(phi2)*math.sin(dlam/2)**2 |
| c = 2*math.atan2(math.sqrt(a), math.sqrt(1-a)) |
| return R * c |
|
|
|
|
| def main(): |
| df = pd.read_parquet(PARQUET) |
| with open(ROUTES) as f: |
| routes_data = json.load(f) |
| routes = routes_data.get('routes', routes_data) |
| airports = pd.read_csv(AIRPORTS) |
|
|
| s630 = df[df['step'] == STEP] |
| route_icaos = set() |
| for icao24, r in routes.items(): |
| for s in r.get('segments', []): |
| if s['startStep'] <= STEP <= s['endStep'] and (s.get('origin') or s.get('destination')): |
| route_icaos.add(icao24) |
| break |
| route_aircraft = s630[s630['icao24'].isin(route_icaos)] |
|
|
| |
| icao24_tc1 = '0100dc' |
| row = s630[s630['icao24'] == icao24_tc1].iloc[0] |
| org, dst, org_icao, dst_icao, _ = get_route_for_step(routes, icao24_tc1, STEP) |
| baro_m = row['baroaltitude'] |
| baro_ft = baro_m * 3.28084 |
| ap = airports[airports['iata_code'] == dst].iloc[0] |
| dist_km = haversine_km(row['lat'], row['lon'], ap['latitude_deg'], ap['longitude_deg']) |
| eta_min = dist_km / (row['velocity'] * 3.6) * 60 |
|
|
| print("=" * 60) |
| print("TABLE 2 GROUND TRUTH (step=630, 16:30 UTC)") |
| print("=" * 60) |
| print() |
| print("TC1: Callsign lookup → read route") |
| print(" callsign: MSR788 | icao24: 0100dc") |
| print(" origin: MUC (EDDM) | destination: CAI (HECA)") |
| print(" Verification: get_route_for_step(routes, '0100dc', 630)") |
| print() |
| print("TC2: Callsign lookup → read parameters") |
| print(" lat: {:.6f} | lon: {:.6f}".format(row['lat'], row['lon'])) |
| print(" baro_alt: {:.2f} m ({:.0f} ft)".format(baro_m, baro_ft)) |
| print(" speed: {:.2f} m/s | heading: {:.1f} deg".format(row['velocity'], row['heading'])) |
| print(" Verification: states.parquet step=630, icao24=0100dc") |
| print() |
| print("TC3: altitude_ft > 35000") |
| print(" Result: PASS" if baro_ft > 35000 else " Result: FAIL") |
| print(" Verification: baro_ft = {:.0f} > 35000".format(baro_ft)) |
| print() |
| print("TC4: ETA") |
| print(" dist_km: {:.1f} | speed_mps: {:.2f}".format(dist_km, row['velocity'])) |
| print(" ETA_min: {:.1f}".format(eta_min)) |
| print(" Verification: ETA = dist_km / (speed_mps * 3.6) * 60") |
| print() |
|
|
| |
| ap_drs = airports[airports['iata_code'] == 'DRS'].iloc[0] |
| ap_lat, ap_lon = ap_drs['latitude_deg'], ap_drs['longitude_deg'] |
| dists = route_aircraft.apply(lambda r: haversine_km(ap_lat, ap_lon, r['lat'], r['lon']), axis=1) |
| within80 = route_aircraft[dists <= 80] |
| within30 = route_aircraft[dists <= 30] |
|
|
| inbound = [] |
| outbound = [] |
| for _, r in within80.iterrows(): |
| o, d, _, _ = get_route_for_step(routes, r['icao24'], STEP)[:4] |
| if d == 'DRS': |
| inbound.append(r['callsign']) |
| if o == 'DRS': |
| outbound.append(r['callsign']) |
|
|
| print("TC5: Airport search → count route aircraft (DRS, 80km)") |
| print(" count: {} | callsigns: {}".format(len(within80), sorted(within80['callsign'].tolist()))) |
| print(" Verification: route aircraft with dist <= 80km from DRS") |
| print() |
| print("TC6: Airport search → MIN altitude (DRS, 80km)") |
| min_row = within80.loc[within80['baroaltitude'].idxmin()] |
| print(" callsign: {} | baro_alt: {:.2f} m ({:.0f} ft)".format(min_row['callsign'], min_row['baroaltitude'], min_row['baroaltitude']*3.28084)) |
| print(" Verification: min(baroaltitude) among route aircraft within 80km") |
| print() |
| print("TC7: Airport search → first inbound/outbound (DRS, 80km)") |
| print(" inbound (dest=DRS): {} | outbound (origin=DRS): {}".format( |
| sorted(inbound)[0] if inbound else 'NONE', |
| sorted(outbound)[0] if outbound else 'NONE')) |
| print(" Verification: lexicographic min of callsigns with dest=DRS, origin=DRS") |
| print() |
| print("TC8: Airport search → MAX groundspeed (DRS, 80km)") |
| max_row = within80.loc[within80['velocity'].idxmax()] |
| print(" callsign: {} | speed: {:.2f} m/s".format(max_row['callsign'], max_row['velocity'])) |
| print(" Verification: max(velocity) among route aircraft within 80km") |
| print() |
|
|
| |
| dest_to_callsigns = {} |
| for _, r in route_aircraft.iterrows(): |
| d = get_route_for_step(routes, r['icao24'], STEP)[1] |
| if d and d != '?': |
| dest_to_callsigns.setdefault(d, []).append(r['callsign']) |
| tc9_dst = next(d for d, c in dest_to_callsigns.items() if len(c) >= 2) |
| tc9_calls = dest_to_callsigns[tc9_dst][:2] |
| print("TC9: Two callsigns → same destination") |
| print(" callsign1: {} | callsign2: {} | destination: {} | same_destination: 1".format(tc9_calls[0], tc9_calls[1], tc9_dst)) |
| print(" Verification: both have destination={} in routes at step 630".format(tc9_dst)) |
| print() |
|
|
| |
| ap_ycd = airports[airports['iata_code'] == 'YCD'].iloc[0] |
| ap_ycd_lat, ap_ycd_lon = ap_ycd['latitude_deg'], ap_ycd['longitude_deg'] |
| inbound_ycd = [] |
| for _, r in route_aircraft.iterrows(): |
| d = get_route_for_step(routes, r['icao24'], STEP)[1] |
| if d == 'YCD': |
| dist = haversine_km(r['lat'], r['lon'], ap_ycd_lat, ap_ycd_lon) |
| if 20 < dist < 200: |
| eta = (dist - 10) / (r['velocity'] * 3.6) * 60 |
| inbound_ycd.append({'callsign': r['callsign'], 'dist': dist, 'speed': r['velocity'], 'eta': eta}) |
| inbound_ycd.sort(key=lambda x: x['eta']) |
| print("TC10: Geofence prediction (YCD, 10km circle)") |
| print(" airport: YCD (CYCD)") |
| print(" flight1: {} | dist={:.1f}km | ETA ~{:.1f} min".format(inbound_ycd[0]['callsign'], inbound_ycd[0]['dist'], inbound_ycd[0]['eta'])) |
| print(" flight2: {} | dist={:.1f}km | ETA ~{:.1f} min".format(inbound_ycd[1]['callsign'], inbound_ycd[1]['dist'], inbound_ycd[1]['eta'])) |
| print(" enters_first: {}".format(inbound_ycd[0]['callsign'])) |
| print(" Verification: (dist_km - 10) / (speed_mps * 3.6) * 60") |
| print() |
|
|
| |
| print("TC11: Zoom 30km (DRS)") |
| print(" airport: DRS (EDDC) | radius: 30km") |
| print(" count: {} | callsigns: {}".format(len(within30), sorted(within30['callsign'].tolist()))) |
| print(" Verification: route aircraft within 30km of DRS") |
| print() |
|
|
| |
| ap_prn = airports[airports['iata_code'] == 'PRN'].iloc[0] |
| ap_tia = airports[airports['iata_code'] == 'TIA'].iloc[0] |
| dists_tia = route_aircraft.apply(lambda r: haversine_km(ap_tia['latitude_deg'], ap_tia['longitude_deg'], r['lat'], r['lon']), axis=1) |
| within20_tia = route_aircraft[dists_tia <= 20] |
| max_alt_row = within20_tia.loc[within20_tia['baroaltitude'].idxmax()] |
| print("TC12: Pan-required (PRN → TIA)") |
| print(" start: PRN (BKPR) | pan_to: TIA (LATI)") |
| print(" within 20km of TIA: callsign={} | altitude={:.0f} m".format(max_alt_row['callsign'], max_alt_row['baroaltitude'])) |
| print(" Verification: max(baroaltitude) among route aircraft within 20km of TIA") |
| print() |
|
|
| |
| FRANCE = {'lat_min': 41.3, 'lat_max': 51.1, 'lon_min': -5.1, 'lon_max': 9.6} |
| df_before = df[(df['step'] < STEP) & (df['icao24'].isin(route_icaos))] |
| in_france = df_before[ |
| (df_before['lat'] >= FRANCE['lat_min']) & (df_before['lat'] <= FRANCE['lat_max']) & |
| (df_before['lon'] >= FRANCE['lon_min']) & (df_before['lon'] <= FRANCE['lon_max']) |
| ] |
| flew_over = in_france['icao24'].unique() |
| at_630 = set(s630['icao24'].unique()) |
| tc13_icao = next(x for x in flew_over if x in at_630) |
| tc13_callsign = s630[s630['icao24'] == tc13_icao].iloc[0]['callsign'] |
| print("TC13: Common sense - flew over country") |
| print(" flight: {} | country: France".format(tc13_callsign)) |
| print(" bbox: lat {}-{}, lon {}-{}".format(FRANCE['lat_min'], FRANCE['lat_max'], FRANCE['lon_min'], FRANCE['lon_max'])) |
| print(" Verification: trajectory has points in bbox at steps < 630") |
| print() |
|
|
| |
| print("TC14: Multi-path - inbound to airport") |
| print(" callsign: MSR788 | destination: CAI") |
| print(" Verification: route destination=CAI at step 630") |
| print() |
|
|
| |
| print("TC15: Holding pattern") |
| print(" callsign: TB2 T303 | airport: BJV (LTFE)") |
| print(" Verification: trajectory has many points in small area (looping) near BJV at step 630") |
| print() |
|
|
|
|
| if __name__ == "__main__": |
| main() |
|
|
