Flight plan - HALO-20240928a

Crew

Job	Name
PI	Julia Windmiller
WALES	Konstantin Krüger
HAMP	Christian Heske
Dropsondes	Allison Wing
Smart/VELOX	Patrizia Schoch
SpecMACS	Tobias Zinner
Flight Documentation	Daniel Klocke
Ground contact	Luca Schmidt, Sebastian Ortega

Flight plan

Show code cell source Hide code cell source

# Define HALO flight and crew
from datetime import datetime

aircraft = "HALO"
flight_time = datetime(2024, 9, 28, 12, 30)
flight_id = f"{aircraft}-{flight_time.strftime('%Y%m%d')}a"
crew = {
    "Mission PI": "Julia Windmiller",
    "DropSondes": "Allison Wing",
    "HAMP": "Christian Heske",
    "SMART/VELOX": "Patrizia Schoch",
    "SpecMACS": "Tobias Zinner",
    "WALES": "Konstantin Krüger",
    "Flight Documentation": "Daniel Klocke",
    "Ground Support": "Luca Schmidt, Sebastian Ortega",
}

Show code cell source Hide code cell source

# Load forecasts
import easygems.healpix as egh
import intake
import numpy as np
import orcestra.sat
from orcestra.flightplan import tbpb, LatLon

radius = 72e3 * 1.852
halo_speed = 238.5  # at FL450 [m/s]

sat_fcst_date = "2024-09-24"  # date to get satelite forecast(s) from
ifs_fcst_time = "2024-09-24"  # date to get IFS forecast(s) from
ifs_fcst_time = np.datetime64(ifs_fcst_time + "T12:00:00")


# Load satellite track
print(
    f"SATELITE TRACK FORECAST FROM: {sat_fcst_date} FOR FLIGHT DAY: {flight_time:%Y-%m-%d}"
)
ec_track_full = orcestra.sat.SattrackLoader(
    "EARTHCARE", sat_fcst_date, kind="PRE", roi="BARBADOS"
).get_track_for_day(f"{flight_time:%Y-%m-%d}")

ec_track = ec_track_full.where(
    (ec_track_full.lat>-10)&
    (ec_track_full.lat<25)&
    (ec_track_full.lon>-60)&
    (ec_track_full.lon<-35), 
    drop = True).sel(time = slice(
    f"{flight_time:%Y-%m-%d} 10:00", f"{flight_time:%Y-%m-%d} 22:00"))

ec_lons, ec_lats = ec_track.lon.values, ec_track.lat.values

pace_track_full = orcestra.sat.pace_track_loader() \
    .get_track_for_day(f"{flight_time:%Y-%m-%d}")

pace_track = pace_track_full.where(
    (pace_track_full.lat>-10)&
    (pace_track_full.lat<25)&
    (pace_track_full.lon>-60)&
    (pace_track_full.lon<-35), 
    drop = True).sel(time = slice(
    f"{flight_time:%Y-%m-%d} 10:00", f"{flight_time:%Y-%m-%d} 22:00"))

# Load IFS forecast
cat = "https://tcodata.mpimet.mpg.de/internal.yaml"
ifs_ds = (
    intake.open_catalog(cat)
    .HIFS(datetime=ifs_fcst_time)
    .to_dask()
    .pipe(egh.attach_coords)
)

SATELITE TRACK FORECAST FROM: 2024-09-24 FOR FLIGHT DAY: 2024-09-28

Show code cell source Hide code cell source

def find_intersect(ab: tuple[LatLon, LatLon], cd: tuple[LatLon, LatLon], res=1000):
  import pyproj
  import numpy as np

  geod = pyproj.Geod(ellps="WGS84")
  
  def points_on_line(a, b, num_points):
    lons = np.linspace(a.lon, b.lon, num_points)
    lats = np.linspace(a.lat, b.lat, num_points)
    return lons, lats

  a, b = ab
  c, d = cd
      
  ab_dist = geod.inv(a.lon, a.lat, b.lon, b.lat)[2]
  cd_dist = geod.inv(c.lon, c.lat, d.lon, d.lat)[2]
  n = int(ab_dist // res)
  m = int(cd_dist // res)

  lons1, lats1  = points_on_line(ab[0], ab[1], n)
  ab_points = np.asarray(list(zip(lons1, lats1)))
  lons2, lats2 = points_on_line(cd[0], cd[1], m)
  cd_points = np.asarray(list(zip(lons2, lats2)))

  dist = []
  for lon, lat in zip(lons1, lats1):
    lon = np.full(len(lons2), lon)
    lat = np.full(len(lats2), lat)
    dist.append(geod.inv(lon, lat, lons2, lats2)[2])
  dist = np.asarray(dist)

  idx = np.argwhere(dist == np.nanmin(dist))[0]
  assert len(idx) == 2 and "one and only one minimum distance must be found"
  lon = np.mean([lons1[idx[0]], lons2[idx[1]]])
  lat = np.mean([lats1[idx[0]], lats2[idx[1]]])

  return LatLon(lon = lon, lat = lat)

# Create Points of flight plan
from orcestra.flightplan import (
    LatLon,
    FlightPlan,
    point_on_track,
    IntoCircle,
)

radius = 72e3 * 1.852 # Mass flux circle radius (m)
radius_small = 39e3 * 1.852 # Mass flux circle radius ATR sized (m)
speed_halo = 240  # m/s

# EC Track

ref_east = LatLon(lat = 10, lon = -44.1, label = "reference value")
ec_south = point_on_track(ec_track, lat=11.0)
ec_north = point_on_track(ec_track, lat=13.0)

ec_start = find_intersect((tbpb, ref_east), (ec_south, ec_north)).assign(label = "ec_start")
ec_under = point_on_track(ec_track, lat=float(ec_start.lat - 1.5), with_time = True).assign(label = "ec_under", note = "meet EarthCARE")#, time = "2024-09-28T17:42Z")
ec_end = point_on_track(ec_track, lat=ec_start.lat - 3.5).assign(label = "ec_end")

# Circles

circle_west = ec_start.towards(tbpb, distance=radius).assign(label="west", note=f"PACE overpass at {pace_track.time.mean().values}")

fac_rad = 2.80
circle_mid_west = circle_west.towards(ref_east, distance=radius*fac_rad).assign(label="mid_west")
circle_mid_east = circle_west.towards(ref_east, distance=radius*fac_rad*2.0).assign(label="mid_east")
circle_east = circle_west.towards(ref_east, distance=radius*fac_rad*3.0).assign(label="east")
circ_last = ec_end.towards(tbpb, fraction= 0.5).assign(label = "circ_last")

# Create plan 

fl1 = 430
fl2 = 450

waypoints = [
    tbpb,
    IntoCircle(circle_east.assign(fl = fl1), radius = radius, angle = 360),
    IntoCircle(circle_mid_east.assign(fl = fl1), radius = radius, angle = 360, enter = -90),
    IntoCircle(circle_mid_west.assign(fl = fl2), radius = radius, angle = 360, enter = 90),
    IntoCircle(circle_west.assign(fl = fl2), radius = radius, angle = 360, enter = -90),
    ec_start.assign(fl = fl2),
    ec_under.assign(fl = fl2),
    ec_end.assign(fl = fl2),
    IntoCircle(circ_last.assign(fl = fl2), radius = radius_small, angle = 360),
    tbpb,
]

extra_waypoints = []

# FlightPlan and print short statement
plan = FlightPlan(
    path=waypoints,
    flight_id=flight_id,
    extra_waypoints=extra_waypoints,
    crew=crew,
    aircraft=aircraft,
)

msg = (
    f"Flight ID: {plan.flight_id}\n"
    + f"Take-off: {plan.takeoff_time:%H:%M %Z}\n"
    + f"Landing:  {plan.landing_time:%H:%M %Z}\n"
    + f"Duration: {plan.duration}"
)
print(msg)

Flight ID: HALO-20240928a
Take-off: 10:28 UTC
Landing:  19:28 UTC
Duration: 8:59:55.203430

Show code cell source Hide code cell source

# Plot CWV forecast
from HALO_20240919a_plan import plot_flight_plan_satellite_forecast
from orcestra.flightplan import plot_cwv

figsize = (14, 8)
lon_min, lon_max, lat_min, lat_max = -65, -5, -5, 25
domain_lonlat = [lon_min, lon_max, lat_min, lat_max]

def plot_ifs_cwv_forecast(fig, ax, ifs_ds, flight_time, levels=None):
    cwv_flight_time = ifs_ds["tcwv"].sel(time=flight_time, method="nearest")
    plot_cwv(cwv_flight_time, ax=ax, levels=levels)

plot_cwv_kwargs = {"flight_time": flight_time, "levels": [48, 50, 52, 54, 56]}

plot_flight_plan_satellite_forecast(
    figsize,
    flight_time,
    plan,
    domain_lonlat,
    is_ec_track=True,
    ec_track=ec_track,
    is_pace_track=True,
    pace_track=pace_track,
    forecast_overlay=True,
    ifs_ds=ifs_ds,
    ifs_fcst_time=ifs_fcst_time,
    forecast_title_label="CWV",
    plot_forecast_func=plot_ifs_cwv_forecast,
    plot_forecast_kwargs=plot_cwv_kwargs,
    atc_zones=True,
    worldfirs_json="worldfirs.json",
    is_meteor=False,
)

# Plot 10m winds forecast
import HALO_20240919a_sfc_winds as sfc_winds

figsize = (16, 9)

def plot_ifs_sfc_winds_forecast(fig, ax, ifs_ds, domain_lonlat, flight_time):
    u10m = ifs_ds["10u"].sel(time=flight_time, method="nearest")
    v10m = ifs_ds["10v"].sel(time=flight_time, method="nearest")
    windspeed_10m = np.sqrt(u10m**2 + v10m**2)
    sfc_winds._windspeed_plot(windspeed_10m, fig, ax)
    sfc_winds._wind_direction_plot(u10m, v10m, ax, domain_lonlat)
    sfc_winds._windspeed_contour(windspeed_10m, ax)
    sfc_winds._draw_confluence_contour(v10m, ax)

plot_ifs_sfc_winds_kwargs = {
    "domain_lonlat": domain_lonlat,
    "flight_time": flight_time,
}

fig, ax = plot_flight_plan_satellite_forecast(
    figsize,
    flight_time,
    plan,
    domain_lonlat,
    is_ec_track=True,
    ec_track=ec_track,
    is_pace_track=True,
    pace_track=pace_track,
    forecast_overlay=True,
    ifs_ds=ifs_ds,
    ifs_fcst_time=ifs_fcst_time,
    forecast_title_label="10m Winds",
    plot_forecast_func=plot_ifs_sfc_winds_forecast,
    plot_forecast_kwargs=plot_ifs_sfc_winds_kwargs,
    atc_zones=True,
    worldfirs_json="worldfirs.json",
    is_meteor=False,
)

/home/runner/miniconda3/envs/orcestra_book/lib/python3.12/site-packages/cartopy/io/__init__.py:241: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/50m_physical/ne_50m_coastline.zip
  warnings.warn(f'Downloading: {url}', DownloadWarning)

Show code cell source Hide code cell source

from orcestra.flightplan import vertical_preview

vertical_preview(waypoints)

Show code cell source Hide code cell source

plan.show_details()
plan.export()

Detailed Overview:
              TBPB         N13 04.48, W059 29.55, FL000, 10:28:53 UTC, 
circle around east         N10 09.86, W044 51.15, FL430, 12:34:55 UTC - 13:33:56 UTC, radius: 72 nm, 360° CW, enter from west
circle around mid_east     N10 52.20, W048 11.29, FL430, 13:41:27 UTC - 14:40:28 UTC, radius: 72 nm, 360° CW, enter from east
circle around mid_west     N11 32.36, W051 32.34, FL450, 15:25:23 UTC - 16:23:56 UTC, radius: 72 nm, 360° CW, enter from west, fly through circle before entering
circle around west         N12 10.18, W054 54.32, FL450, 16:31:23 UTC - 17:29:56 UTC, radius: 72 nm, 360° CW, enter from east, PACE overpass at 2024-09-28T16:35:35.000000000
to            ec_start     N11 55.17, W053 42.44, FL450, 17:30:10 UTC, 
to            ec_under     N10 25.17, W053 59.63, FL450, 17:41:58 UTC, meet EarthCARE
to            ec_end       N08 25.17, W054 22.39, FL450, 17:57:42 UTC, 
circle around circ_last    N10 45.47, W056 54.79, FL450, 18:19:14 UTC - 18:50:57 UTC, radius: 39 nm, 360° CW, enter from south east
to            TBPB         N13 04.48, W059 29.55, FL000, 19:28:48 UTC,