from dataclasses import asdict
from datetime import datetime
import cartopy.crs as ccrs
import easygems.healpix as egh
import intake
import matplotlib.pyplot as plt
import numpy as np
import orcestra
import orcestra.flightplan as fp
import orcestra.sat
from orcestra.flightplan import LatLon, IntoCircle, bco, sal, mindelo, find_ec_lon, vertical_preview, to_kml\
def ec_time_at_lat(ec_track, lat):
e = np.datetime64("2024-08-01")
s = np.timedelta64(1, "ns")
return (((ec_track.swap_dims({"time":"lat"}).time - e) / s).interp(lat=lat) * s + e)
# Global coordinates and definitions that should not change from flight to flight
lon_min, lon_max, lat_min, lat_max = -65, -5, -5, 25
radius = 130e3
atr_radius = 72e3
band = "east"
airport = sal if band == "east" else bco
natal = LatLon(-5 - 47/60. - 42.00/3600.,-35 - 12/60. - 33.98/3600., label = "natal")
# Basic information
lon_min, lon_max, lat_min, lat_max = -65, -5, -5, 25
# Define dates for forecast initialization and flight
issued_time = datetime(2024, 8, 27, 0, 0, 0)
flight_time = datetime(2024, 8, 29, 12, 0, 0)
flight_index = f"HALO-{flight_time.strftime('%Y%m%d')}a"
# adjust takeoff time to match EC overpass
takeoff_time = np.datetime64("2025-08-29T12:20:00")
print(
f"Initalization date of IFS forecast: {issued_time}\n"
f"Flight date: {flight_time:%Y-%m-%d}\n"
f"Flight index: {flight_index}"
)
crew = {'Mission PI': 'Silke Gross',
'DropSondes': 'Nina Robbins',
'HAMP': 'Christian Heske',
'SMART/VELOX': 'Michael Schäfer',
'SpecMACS': 'Veronika Pörtge',
'WALES' : 'Georgios Dekoutsidis',
'Flight Documentation': 'Basile Poujol',
'Ground Support': 'Julia Windmiller',
}
# Load forecast data
cat = intake.open_catalog("https://tcodata.mpimet.mpg.de/internal.yaml")
ds = cat.HIFS(datetime = issued_time).to_dask().pipe(egh.attach_coords)
# Load ec satellite track for
ec_track = orcestra.sat.SattrackLoader("EARTHCARE", "2024-08-26", kind="PRE").get_track_for_day(f"{flight_time:%Y-%m-%d}")
ec_track = ec_track.sel(time=slice(f"{flight_time:%Y-%m-%d} 06:00", None))
ec_lons, ec_lats = ec_track.lon.values, ec_track.lat.values
# Latitudes where we enter and leave the ec track (visually estimated)
lat_ec_north = 15.0
lat_ec_south = 2.5
# latitude of circle centers
lat_c_south_s = 3.5
lat_c_south = 4.5
lat_c_south_n = 5.5
lat_c_north = 13.0
lat_c_mid_s = 7.5
lat_c_mid = 8.5
lat_c_mid_n = 9.5
lat_ec_under = 5.0
c_atr_nw = LatLon(18.58125000,-24.27616667, label = "c_atr")
c_atr_se = LatLon(15.79318333,-24.82891944, label = "c_atr")
# create ec track
ec_north = LatLon(lat_ec_north, find_ec_lon(lat_ec_north, ec_lons, ec_lats), label = "ec_north")
ec_south = LatLon(lat_ec_south, find_ec_lon(lat_ec_south, ec_lons, ec_lats), label = "ec_south")
# create circles
c_north = LatLon(lat_c_north, find_ec_lon(lat_c_north, ec_lons, ec_lats), label = "c_north")
c_south = LatLon(lat_c_south, find_ec_lon(lat_c_south, ec_lons, ec_lats), label = "c_south")
c_south_s = LatLon(lat_c_south_s, find_ec_lon(lat_c_south_s, ec_lons, ec_lats), label = "c_south_s")
c_south_n = LatLon(lat_c_south_n, find_ec_lon(lat_c_south_n, ec_lons, ec_lats), label = "c_south_n")
c_mid = LatLon(lat_c_mid, find_ec_lon(lat_c_mid, ec_lons, ec_lats), label = "c_mid")
c_mid_s = LatLon(lat_c_mid_s, find_ec_lon(lat_c_mid_s, ec_lons, ec_lats), label = "c_mid_s")
c_mid_n = LatLon(lat_c_mid_n, find_ec_lon(lat_c_mid_n, ec_lons, ec_lats), label = "c_mid_n")
# ec underpass
ec_under = LatLon(lat_ec_under, find_ec_lon(lat_ec_under, ec_lons, ec_lats), label = "ec_under")
ec_under = ec_under.assign(time=str(ec_time_at_lat(ec_track, ec_under.lat).values)+"Z")
# Define flight track
outbound_legs = [
airport,
mindelo,
ec_north.assign(fl=410),
]
ec_legs = [
IntoCircle(c_south.assign(fl=430), radius, 360),
ec_south.assign(fl=410),
ec_under.assign(fl=450),
IntoCircle(c_mid.assign(fl=430), radius, 360),
IntoCircle(c_north.assign(fl=450), radius, 360),
]
inbound_legs = [
ec_north.assign(fl=450),
IntoCircle(c_atr_nw.assign(fl=350), atr_radius, 360),
IntoCircle(c_atr_se.assign(fl=350), atr_radius, 360),
airport,
]
waypoints = outbound_legs + ec_legs + inbound_legs
waypoint_centers = []
for point in waypoints:
if isinstance(point, IntoCircle):
point = point.center
waypoint_centers.append(point)
path = fp.expand_path(waypoints, dx=10e3)
plan = path.isel(distance = path.waypoint_indices).to_dataframe().set_index("waypoint_labels")
xwp_2 = LatLon(lat_c_south-1, find_ec_lon(lat_c_south-1, ec_lons, ec_lats), label = "xwp2")
xwp_3 = LatLon(c_atr_nw.lat,c_atr_nw.lon, label = "xwp3")
extra_waypoints = [xwp_2,xwp_3]
notes = {'c_south_in':f' {radius/1852:2.0f} nm circle centered at {c_south.format_pilot()}, enter from north, CCW',
'c_mid_in':f' {radius/1852:2.0f} nm circle centered at {c_mid.format_pilot()}, enter from north, CCW',
'c_north_in':f' {radius/1852:2.0f} nm circle centered at {c_north.format_pilot()}, enter from south, CCW',
'c_atr_in':f' {atr_radius/1852:2.0f} nm circle centered at {c_atr_se.format_pilot()}, enter from west, CW',
'xwp2':'Alternative center for c_south',
'xwp3':'Alternative center for c_atr',
}
plt.figure(figsize = (14, 8))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_extent([lon_min, lon_max, lat_min, lat_max], crs=ccrs.PlateCarree())
ax.coastlines(alpha=1.0)
ax.gridlines(draw_labels=True, dms=True, x_inline=False, y_inline=False, alpha = 0.25)
cwv_flight_time = ds["tcwv"].sel(time=flight_time, method = "nearest")
fp.plot_cwv(cwv_flight_time, levels = [50.0, 60.0])
plt.title(f"{flight_time}\n(CWV forecast issued on {issued_time})")
plt.plot(ec_lons, ec_lats, c='k', ls='dotted')
if (False):
plt.plot([natal.lon,sal.lon], [natal.lat,sal.lat], c='purple', ls='dashed')
for wp in waypoint_centers:
plt.scatter(wp.lon,wp.lat,s=10.,color='k')
for wp in extra_waypoints:
plt.scatter(wp.lon,wp.lat,s=10.,color='r',marker='o')
fp.plot_path(path, ax, color="C1")
