import intake
import easygems.healpix as egh
from datetime import datetime, timedelta
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import orcestra.sat
import orcestra.flightplan as fp
from orcestra.flightplan import LatLon, IntoCircle, bco, sal, mindelo, find_ec_lon, vertical_preview, to_kml
# Define dates for forecast initialization and flight
issued_time = datetime(2024, 8, 24, 0, 0, 0)
flight_time = datetime(2024, 8, 25, 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}"
)
# 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-24", 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
# Domain definition
lon_min, lon_max, lat_min, lat_max = -65, -5, -10, 25
## Setting lat/lon coordinates
# Mass flux circle radius (m)
radius = 130e3
atr_radius = 70e3
# Setting region (Cabo Verde vs. Barbados)
airport = sal #bco
# Latitudes where we enter and leave the ec track (visually estimated)
lat_ec_north_out = 16.0
lat_ec_north_in = 12.0
lat_ec_south = 2.5
# ITCZ edges visually estimated from iwv contours
lat_c_south = 4.0
lat_c_north = 10.0
# Extra Points along the flight track
lat_c_south_alt_s = 3.25
lat_c_south_alt_n = 4.75
lat_c_north_alt_s = 9.25
lat_c_north_alt_n = 10.75
# Extra Points off the flight track
lat_ec_south_alt = 2.0
# Setting lat/lon coordinates
# Points where we get on ec track
north_ec_in = LatLon(lat_ec_north_in, find_ec_lon(lat_ec_north_in, ec_lons, ec_lats), label = "north_ec_in")
north_ec_out = LatLon(lat_ec_north_out, find_ec_lon(lat_ec_north_out, ec_lons, ec_lats), label = "north_ec_out")
south_ec = LatLon(lat_ec_south, find_ec_lon(lat_ec_south, ec_lons, ec_lats), label = "south_ec")
# Intersection of ITCZ edges with ec track
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")
# Center of middle circle
c_mid = c_south.towards(c_north).assign(label = "c_mid")
# EarthCARE underpass
ec_under = c_north.towards(north_ec_out)
ec_under = LatLon(ec_under.lat, find_ec_lon(ec_under.lat, ec_lons, ec_lats), label = "ec_under")
# ATR circle
atr_circ_sar = LatLon(lat=15.5, lon=-22.1, label="atr_sar")
# Extra Points along the way
c_south_alt_s = LatLon(lat_c_south_alt_s, find_ec_lon(lat_c_south_alt_s, ec_lons, ec_lats), label = "c_south_alt_s")
c_south_alt_n = LatLon(lat_c_south_alt_n, find_ec_lon(lat_c_south_alt_n, ec_lons, ec_lats), label = "c_south_alt_n")
c_north_alt_s = LatLon(lat_c_north_alt_s, find_ec_lon(lat_c_north_alt_s, ec_lons, ec_lats), label = "c_north_alt_s")
c_north_alt_n = LatLon(lat_c_north_alt_n, find_ec_lon(lat_c_north_alt_n, ec_lons, ec_lats), label = "c_north_alt_n")
# Extra Points
ec_south_alt = LatLon(lat_ec_south_alt, find_ec_lon(lat_ec_south_alt, ec_lons, ec_lats), label = "south_ec_alt")
# Define flight track, can be split into different legs
take_off_time = '2024-08-25T09:15:00Z'
waypoints = [
airport.assign(time=take_off_time),
north_ec_in.assign(fl=400),
c_north.assign(fl=400),
c_mid.assign(fl=400),
c_south.assign(fl=400),
south_ec.assign(fl=430),
c_south_alt_s.assign(fl=430),
IntoCircle(c_south.assign(fl=430), radius, 360),
c_south_alt_n.assign(fl=430),
IntoCircle(c_mid.assign(fl=430), radius, 360),
c_north_alt_s.assign(fl=430),
IntoCircle(c_north.assign(fl=450), radius, 360),
c_north_alt_n.assign(fl=450),
ec_under.assign(fl=450),
north_ec_out.assign(fl=450),
mindelo.assign(fl=450),
IntoCircle(atr_circ_sar.assign(fl=350), atr_radius, 360),
airport
]
path = fp.expand_path(waypoints, dx=10e3)
# Plot flight path
plt.figure(figsize = (12, 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 = [45.0, 48.0, 50.0, 52.0, 54.0])
plt.title(f"{flight_time}\n(CWV forecast issued on {issued_time})")
plt.plot(ec_lons, ec_lats)
fp.plot_path(path, ax, color="C1")