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, FlightPlan\
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 = 133e3
atr_radius = 72e3
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, 9, 5, 0, 0, 0)
flight_time = datetime(2024, 9, 6, 12, 0, 0)
flight_index = f"HALO-{flight_time.strftime('%Y%m%d')}a"
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': 'Julia Windmiller',
'DropSondes': 'Bjorn Stevens',
'HAMP': 'Yuting Wu',
'SMART/VELOX': 'Michael Schäfer',
'SpecMACS': 'Lea Volkmer',
'WALES' : 'Martin Wirth',
'Flight Documentation': 'Silke Gross',
'Ground Support': 'Brett McKim',
}
# 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 the 5th and 6th of September
flight_time_old = datetime(2024, 9, 5, 12, 0, 0)
ec_track = orcestra.sat.SattrackLoader("EARTHCARE", "2024-09-04", kind="PRE", roi="BARBADOS").get_track_for_day(f"{flight_time_old:%Y-%m-%d}")
ec_track = ec_track.sel(time=slice(f"{flight_time_old:%Y-%m-%d} 14:00", None))
ec_lons, ec_lats = ec_track.lon.values, ec_track.lat.values
load_ec = False
if load_ec:
ec_06 = orcestra.sat.earthcare_track_loader().get_track_for_day(f"{flight_time:%Y-%m-%d}").sel(time=slice(f"{flight_time:%Y-%m-%d} 14:00", f"{flight_time:%Y-%m-%d} 17:00"))
ec_06 = ec_06.where((ec_06.lat >-5)&(ec_06.lat <25)&(ec_06.lon >-60)&(ec_06.lon <-20), drop = True)
# Latitudes where we enter, underfly, and leave the ec track (visually estimated)
lat_ec_north = bco.lat
lat_ec_under = 10.0
lat_ec_south = 8.0
# 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")
ec_south2 = LatLon(lat_ec_south+1, find_ec_lon(lat_ec_south+1, ec_lons, ec_lats), label = "ec_south_2")
if load_ec:
ec_06_crossing = LatLon(lat_ec_south, find_ec_lon(lat_ec_south, ec_06.lon, ec_06.lat), label = "ec_crossing")
# 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")
# create circles
lat_circ = lat_ec_south
c_track_on = LatLon(lat_circ, -32.0, label = "c_track_on")
c_east = LatLon(lat_circ, -35.0, label = "c_east")
c_center = LatLon(lat_circ, -40.0, label = "c_center")
lat_c_west = 12.0
c_west = LatLon(lat_c_west, find_ec_lon(lat_c_west, ec_lons, ec_lats), label = "c_west")
# extra waypoints
lat_circ_alt = lat_ec_south
meteor = LatLon(lat_circ_alt, -38.0, label = "meteor")
c_center_s = LatLon(lat_circ_alt, -40.0, label = "c_center_s")
# Define flight track
outbound_legs = [
sal.assign(time='2024-09-06T10:00Z'),
c_track_on.assign(fl=410),
]
circle_legs = [
meteor.assign(fl=430),
# ec_06_crossing.assign(fl=430),
]
ec_legs = [
ec_south.assign(fl=450),
ec_south2.assign(fl=450),
ec_under.assign(fl=450),
IntoCircle(c_west.assign(fl=450), radius, -360, enter = 90),
#ec_north.assign(fl=450),
]
inbound_legs = [
bco,
]
waypoints = outbound_legs + circle_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")
extra_waypoints = [meteor]
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, ax, levels = [40.0, 50.0, 52.0, 54.0, 56.0, 58.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")
if load_ec:
plt.plot(ec_06.lon, ec_06.lat, "k--", linewidth = 1)
plt.scatter(ec_06_crossing.lon, ec_06_crossing.lat, c = 'k')
