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")
# Define dates for forecast initialization and flight
issued_time = datetime(2024, 8, 26, 0, 0, 0)
flight_time = datetime(2024, 8, 27, 16, 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-27T10:15")
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': 'Bjorn Stevens',
'DropSondes': 'Helene Glöckner',
'HAMP': 'Clara Bayley',
'SMART/VELOX': 'Kevin Wolf',
'SpecMACS': 'Zekican Demiralay',
'WALES' : 'Georgios Dekoutsidis',
'Flight Documentation': 'Nicola Rochetin',
'Ground Support': 'Marius Rixen',
}
# 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} 14: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 = 13.0
lat_ec_south = 5.75
# latitude of circle centers
lat_c_south = 4.0
lat_c_north = 11.25
c_atr_nw = LatLon(17.433,-23.500, label = "c_atr")
c_atr_se = LatLon(14.500,-22.500, label = "c_atr")
meteor = LatLon(13.000,-22.500, label = "after-meteor")
xsonde = LatLon(meteor.lat, -23.500, label = "x-sonde")
atr_turn = LatLon(sal.lat-0.2, -22.500, label = "after-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_mid = c_south.towards(c_north).assign(label = "c_mid")
# extra way point on track
x_wp1 = c_mid.towards(c_north,fraction=2/5.).assign(label = "x_wp1")
# ec underpass
ec_under = c_north.towards(ec_north,fraction=1./3.).assign(label = "ec_under")
ec_under = LatLon(ec_under.lat, find_ec_lon(ec_under.lat, ec_lons, ec_lats), label = ec_under.label)
# Define flight track
outbound_legs = [
airport,
ec_north.assign(fl=410),
c_north.assign(fl=410),
c_mid.assign(fl=410),
]
ec_legs = [
ec_south.assign(fl=410),
IntoCircle(c_south.assign(fl=430), radius, -330,enter=220),
ec_south.assign(fl=430),
IntoCircle(c_mid.assign(fl=430), radius, 360),
x_wp1.assign(fl=450),
IntoCircle(c_north.assign(fl=450), radius, 360),
ec_under.assign(fl=450),
]
inbound_legs = [
ec_north.assign(fl=450),
xsonde.assign(fl=450),
meteor.assign(fl=450),
c_atr_se.assign(fl=350),
IntoCircle(c_atr_se.assign(fl=350), atr_radius, 360, enter=90),
atr_turn.assign(fl=350),
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)
# extra way points off track
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 south, 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',
'ec_under':f'EarthCARE at this point at {str(ec_time_at_lat(ec_track, ec_under.lat).values)[11:19]} UTC',
}
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')
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")
