Weather module [DRAFT]

src/weather/draw.py

@@ -0,0 +1,49 @@
+import matplotlib.pyplot as plt
+import cartopy.crs as ccrs
+import cartopy.feature as cfeature
+from pyproj import Geod
+import numpy as np
+
+
+def plot_flight_arc(mission):
+
+    # Set up map
+
+    fig = plt.figure(figsize=(10, 7))
+    ax = plt.axes(projection=ccrs.PlateCarree())
+    ax.set_extent([-130, -65, 22, 50], crs=ccrs.PlateCarree())
+
+    # Add map features
+    ax.add_feature(cfeature.STATES, linewidth=0.5)
+    ax.add_feature(cfeature.COASTLINE)
+    ax.add_feature(cfeature.BORDERS, linestyle=':')
+    ax.add_feature(cfeature.LAND, facecolor='lightgray')
+    ax.add_feature(cfeature.OCEAN, facecolor='lightblue')
+    ax.gridlines(draw_labels=True, linewidth=0.5, color='gray', alpha=0.5)
+
+    # Extract OD lat lon position
+    lon_dep, lat_dep, _ = mission["dep_location"]
+    lon_arr, lat_arr, _ = mission["arr_location"]
+
+    # Use WGS84 ellipse definition
+    geod = Geod(ellps="WGS84")
+
+    # Get 100 equally spaced points along the arc
+    points = geod.npts(lon_dep, lat_dep, lon_arr, lat_arr, 100)
+    lons = [lon_dep] + [pt[0] for pt in points] + [lon_arr]
+    lats = [lat_dep] + [pt[1] for pt in points] + [lat_arr]
+
+    ax.plot(lons, lats, 'k-', linewidth=0.8, alpha=0.7)
+
+    # Mark endpoints
+    ax.plot(lon_dep, lat_dep, 'go', markersize=6, label="Departure")
+    ax.plot(lon_arr, lat_arr, 'ro', markersize=6, label="Arrival")
+
+    # Add legend and title
+    ax.legend(loc='lower left')
+    dep_code = mission['dep_airport']
+    arr_code = mission['arr_airport']
+    ax.set_title(f"Flight Arc: {dep_code} → {arr_code}")
+
+    plt.tight_layout()
+    plt.savefig("sample.png", dpi=300)

src/weather/main.py

@@ -0,0 +1,92 @@
+import json
+import draw as dr
+import utils as util
+import trajectory as traj
+
+
+mission_path = "../missions/sample_missions_10.json"
+weather_data = "ERA5/sample.grib"
+
+# Load JSON data
+with open(mission_path, 'r') as file:
+    missions = json.load(file)
+
+
+# Select a dummy mission for now
+first_mission = missions[0]
+
+
+# Get mission points based on mission def
+trajectory = traj.create_dummy_traj(first_mission)
+
+# Unpack trajectory -->
+lons = trajectory["lons"]
+lats = trajectory["lats"]
+alts = trajectory["H"]
+tas = trajectory["TAS"]
+
+# Initialize lists -->
+gs_list, heading_list, u_list, v_list = [], [], [], []
+
+# Build wind interpolator ->
+u_interp, v_interp, meta = util.build_era5_interpolators(weather_data)
+
+# Loop over all mission points to compute heading and ground speed ->
+for i in range(len(lons)):
+    if i < len(lons) - 1:
+        lon_next, lat_next = lons[i + 1], lats[i + 1]
+    else:
+        # For the last point, repeat previous heading
+        lon_next, lat_next = lons[i - 1], lats[i - 1]
+
+    gs, heading, u, v = util.compute_ground_speed(
+        lon=lons[i],
+        lat=lats[i],
+        lon_next=lon_next,
+        lat_next=lat_next,
+        alt_ft=alts[i],
+        tas_kts=tas[i],
+        u_interp=u_interp,
+        v_interp=v_interp,
+    )
+
+    # Append to respective lists for sanity check ->
+    gs_list.append(gs)
+    heading_list.append(heading)
+    u_list.append(u)
+    v_list.append(v)
+
+# Append to trajectory ->
+trajectory["GS"] = gs_list
+trajectory["heading_rad"] = heading_list
+trajectory["u_wind"] = u_list
+trajectory["v_wind"] = v_list
+
+
+# Print the first few points
+print("\n")
+print("#--------------------------------SAMPLE TRAJECTORY (SNAPSHOT) -----------------------------#\n")
+for lon, lat, tas, gs, alt in zip(trajectory["lons"][:5], trajectory["lats"][:5], trajectory["TAS"][:5], trajectory["GS"][:5], trajectory["H"][:5]):
+    print(f"Lon: {lon:.2f}, Lat: {lat:.2f}, TAS: {tas:.2f} kt, GS: {gs:.2f} kt, Alt: {alt:.2f} ft")
+print("#------------------------------------------------------------------------------------------#\n")
+
+
+
+#track, heading, drift, tas, u, v, wind_mag = util.get_tas(trajectory, "era5.grib")
+
+#for i in range(5):
+#    print(f"Pt {i}: Track={track[i]:.1f}°, Heading={heading[i]:.1f}°, Drift={drift[i]:+.1f}°, "
+#          f"TAS={tas[i]:.1f} kt, U={u[i]:.1f}, V={v[i]:.1f}, Wind={wind_mag[i]:.1f} m/s")
+
+
+
+#dr.plot_flight_arc(first_mission)
+
+#util.get_flight_track(first_mission)
+
+
+
+
+
+
+

src/weather/trajectory.py

@@ -0,0 +1,132 @@
+from pyproj import Geod
+import numpy as np
+
+
+
+def get_mission_points(mission):
+    """
+    Generates a discretized set of latitude and longitude points along a great-circle path
+    between departure and arrival locations, assuming constant cruise altitude and ground speed.
+
+    Parameters
+    ----------
+    mission : dict
+        Dictionary containing origin and destination coordinates with the following keys:
+            - 'dep_location': tuple of (longitude, latitude, altitude) for the departure point [degrees, degrees, feet]
+            - 'arr_location': tuple of (longitude, latitude, altitude) for the arrival point [degrees, degrees, feet]
+
+    Returns
+    -------
+    dict
+        A dictionary containing:
+            - 'lons' : list of longitudes along the flight path [degrees]
+            - 'lats' : list of latitudes along the flight path [degrees]
+            - 'GS'   : list of assumed constant ground speed at each point [knots]
+            - 'H'    : list of assumed constant cruise altitude at each point [feet]
+
+    Notes
+    -----
+    - The path is discretized using 100 intermediate points (plus endpoints), resulting in 102 total waypoints.
+    - This function uses `pyproj.Geod` with the WGS84 ellipsoid to compute a geodesic path.
+    - All values for speed and altitude are placeholders; replace them with mission-specific data in actual use.
+    """
+
+    # Instantiate WGS84 ellipsoid
+    geod = Geod(ellps ="WGS84")
+
+    # Extract OD lat-lon 
+    lon_dep, lat_dep, _ = mission["dep_location"]
+    lon_arr, lat_arr, _ = mission["arr_location"]
+
+    # Assume 100 points for discretization (for demo only)
+    # Note: This will change when flying actual missions
+
+    points = geod.npts(lon_dep, lat_dep, lon_arr, lat_arr, 100)
+
+    lons = [lon_dep] + [pt[0] for pt in points] + [lon_arr]
+    lats = [lat_dep] + [pt[1] for pt in points] + [lat_arr]
+
+    # Assign a dummy ground speed
+    ground_speeds = [450] * len(lons)
+
+    # Assign a dummy cruise altitude
+    altitude_ft = [35000] * len(lons)
+
+    return {
+        "lons": lons,
+        "lats": lats,
+        "GS": ground_speeds,
+        "H": altitude_ft
+    }
+
+def create_dummy_traj(mission):
+    """
+    Generates a discretized trapezoidal trajectory profile (altitude and speed)
+    along a geodesic path between departure and arrival locations.
+
+    Parameters
+    ----------
+    mission : dict
+        Dictionary with:
+            - 'dep_location': (longitude, latitude, altitude) of departure [degrees, degrees, feet]
+            - 'arr_location': (longitude, latitude, altitude) of arrival [degrees, degrees, feet]
+
+    Returns
+    -------
+    dict
+        Dictionary containing:
+            - 'lons': list of longitudes [degrees]
+            - 'lats': list of latitudes [degrees]
+            - 'TAS'  : list of True Air Speed [knots]
+            - 'H'   : list of altitudes [feet]
+    """
+    # Instantiate GGS84 ellipsoid for geodesic calculation
+    geod = Geod(ellps="WGS84")
+
+    # Extract dept + Arrival lat-lon
+    lon_dep, lat_dep, _ = mission["dep_location"]
+    lon_arr, lat_arr, _ = mission["arr_location"]
+
+    # Assume 100 points for trajectory discretization (demo only)
+    n_total = 100
+    points = geod.npts(lon_dep, lat_dep, lon_arr, lat_arr, n_total)
+    lons = [lon_dep] + [pt[0] for pt in points] + [lon_arr]
+    lats = [lat_dep] + [pt[1] for pt in points] + [lat_arr]
+
+    # Define climb and descent length as 25 % of total trajectory
+    n_climb = n_descent = n_total // 4
+    n_cruise = len(lons) - n_climb - n_descent
+
+    # Define end points for speed and altitude 
+    alt_start, alt_cruise = 1000, 35000  # feet
+    spd_start, spd_cruise = 140, 450 # knots
+
+    # Define climb profile
+    climb_alt = np.linspace(alt_start, alt_cruise, n_climb)
+    climb_spd = np.linspace(spd_start, spd_cruise, n_climb)
+
+    # Define cruise profile
+    cruise_alt = np.full(n_cruise, alt_cruise)
+    cruise_spd = np.full(n_cruise, spd_cruise)
+
+    # Define descent profile
+    descent_alt = np.linspace(alt_cruise, alt_start, n_descent)
+    descent_spd = np.linspace(spd_cruise, spd_start, n_descent)
+
+
+
+    # Collect height and TAS
+    H = np.concatenate([climb_alt, cruise_alt, descent_alt])
+    TAS = np.concatenate([climb_spd, cruise_spd, descent_spd])
+
+    H = H[:len(lons)]
+    TAS = TAS[:len(lons)]
+
+    return {
+        "lons": lons,
+        "lats": lats,
+        "TAS": TAS,
+        "H": H
+    }
+
+