Feature/load-factor-improvements

CHANGELOG.md

@@ -2,9 +2,14 @@
 
 ## v0.54.3 (unreleased)
 
+### Breaking changes
+
+- Update the default load factor from 70% to 83% to be consistent with historical data. This is used whenever an aircraft performance model is run without a specified load factor.
+
 ### Features
 
 - Create new function `ps_nominal_optimize_mach` which computes the optimal mach number given a set of operating conditions.
+- Add a new `jet.aircraft_load_factor` function to estimate aircraft (passenger/cargo) load factor based on historical monthly and regional load factors provided by IATA. This improves upon the default load factor assumption. Historical load factor databases will be continuously updated as new data is released.
 
 ### Fixes
 

pycontrails/core/aircraft_performance.py

@@ -26,7 +26,9 @@
 from pycontrails.utils.types import ArrayOrFloat
 
 #: Default load factor for aircraft performance models.
-DEFAULT_LOAD_FACTOR = 0.7
+#: See :func:`pycontrails.physics.jet.aircraft_load_factor`
+#: for a higher precision approach to estimating the load factor.
+DEFAULT_LOAD_FACTOR = 0.83
 
 
 # --------------------------------------
@@ -36,13 +38,13 @@
 
 @dataclasses.dataclass
 class CommonAircraftPerformanceParams:
-    """Params for :class:`AircraftPerformanceParams` and :class`AircraftPerformanceGridParams`."""
+    """Params for :class:`AircraftPerformanceParams` and :class:`AircraftPerformanceGridParams`."""
 
     #: Account for "in-service" engine deterioration between maintenance cycles.
     #: Default value is set to +2.5% increase in fuel consumption.
-    # Reference:
-    # Gurrola Arrieta, M.D.J., Botez, R.M. and Lasne, A., 2024. An Engine Deterioration Model for
-    # Predicting Fuel Consumption Impact in a Regional Aircraft. Aerospace, 11(6), p.426.
+    #: Reference:
+    #: Gurrola Arrieta, M.D.J., Botez, R.M. and Lasne, A., 2024. An Engine Deterioration Model for
+    #: Predicting Fuel Consumption Impact in a Regional Aircraft. Aerospace, 11(6), p.426.
     engine_deterioration_factor: float = 0.025
 
 

pycontrails/core/flight.py

@@ -121,6 +121,8 @@ class Flight(GeoVectorDataset):
           calculations with the ICAO Aircraft Emissions Databank (EDB).
         - ``max_mach_number``: Maximum Mach number at cruise altitude. Used by
           some aircraft performance models to clip true airspeed.
+        - ``load_factor``: The load factor used in determining the aircraft's
+          take-off weight. Used by some aircraft performance models.
 
         Numeric quantities that are constant over the entire flight trajectory
         should be included as attributes.

pycontrails/models/ps_model/ps_model.py

@@ -63,6 +63,10 @@ class PSFlight(AircraftPerformance):
     Poll & Schumann (2022). An estimation method for the fuel burn and other performance
     characteristics of civil transport aircraft. Part 3 Generalisation to cover climb,
     descent and holding. Aero. J., submitted.
+
+    See Also
+    --------
+    pycontrails.physics.jet.aircraft_load_factor
     """
 
     name = "PSFlight"

pycontrails/physics/jet.py

@@ -1,22 +1,28 @@
 """Jet aircraft trajectory and performance parameters.
 
 This module includes common functions to calculate jet aircraft trajectory
-and performance parameters, including fuel quantities, mass, thrust setting
-and propulsion efficiency.
+and performance parameters, including fuel quantities, mass, thrust setting,
+propulsion efficiency and load factors.
 """
 
 from __future__ import annotations
 
+import functools
 import logging
+import pathlib
 
 import numpy as np
 import numpy.typing as npt
+import pandas as pd
 
 from pycontrails.core import flight
 from pycontrails.physics import constants, units
 from pycontrails.utils.types import ArrayOrFloat, ArrayScalarLike
 
 logger = logging.getLogger(__name__)
+_path_to_static = pathlib.Path(__file__).parent / "static"
+PLF_PATH = _path_to_static / "iata-passenger-load-factors-20241115.csv"
+CLF_PATH = _path_to_static / "iata-cargo-load-factors-20241115.csv"
 
 
 # -------------------
@@ -43,7 +49,7 @@ def acceleration(
 
     See Also
     --------
-    :func:`flight.segment_duration`
+    pycontrails.Flight.segment_duration
     """
     dv_dt = np.empty_like(true_airspeed)
     dv_dt[:-1] = np.diff(true_airspeed) / segment_duration[:-1]
@@ -71,8 +77,8 @@ def climb_descent_angle(
 
     See Also
     --------
-    :func:`flight.segment_rocd`
-    :func:`flight.segment_true_airspeed`
+    pycontrails.Flight.segment_rocd
+    pycontrails.Flight.segment_true_airspeed
     """
     rocd_ms = units.ft_to_m(rocd) / 60.0
     sin_theta = rocd_ms / true_airspeed
@@ -319,8 +325,8 @@ def reserve_fuel_requirements(
 
     See Also
     --------
-    :func:`flight.segment_phase`
-    :func:`fuel_burn`
+    pycontrails.Flight.segment_phase
+    fuel_burn
     """
     segment_phase = flight.segment_phase(rocd, altitude_ft)
 
@@ -345,6 +351,128 @@ def reserve_fuel_requirements(
 # -------------
 
 
+@functools.cache
+def _historical_regional_load_factor(path: pathlib.Path) -> pd.DataFrame:
+    """Load the historical regional load factor database.
+
+    Daily load factors are estimated from linearly interpolating the monthly statistics.
+
+    Returns
+    -------
+    pd.DataFrame
+        Historical regional load factor for each day.
+
+    Notes
+    -----
+    The monthly **passenger load factor** for each region is compiled from IATA's monthly
+    publication of the Air Passenger Market Analysis, where the static file will be continuously
+    updated. The report estimates the regional passenger load factor by dividing the revenue
+    passenger-km (RPK) by the available seat-km (ASK).
+
+    The monthly **cargo load factor** for each region is compiled from IATA's monthly publication
+    of the Air Cargo Market Analysis, where the static file will be continuously updated.
+    The report estimates the regional cargo load factor by dividing the freight tonne-km (FTK)
+    by the available freight tonne-km (AFTK).
+    """
+    df = pd.read_csv(path, index_col="Date", parse_dates=True, date_format="%d/%m/%Y")
+    return df.resample("D").interpolate()
+
+
+AIRPORT_TO_REGION = {
+    "A": "Asia Pacific",
+    "B": "Europe",
+    "C": "North America",
+    "D": "Africa",
+    "E": "Europe",
+    "F": "Africa",
+    "G": "Africa",
+    "H": "Africa",
+    "K": "North America",
+    "L": "Europe",
+    "M": "Latin America",
+    "N": "Asia Pacific",
+    "O": "Middle East",
+    "P": "Asia Pacific",
+    "R": "Asia Pacific",
+    "S": "Latin America",
+    "T": "Latin America",
+    "U": "Asia Pacific",
+    "V": "Asia Pacific",
+    "W": "Asia Pacific",
+    "Y": "Asia Pacific",
+    "Z": "Asia Pacific",
+}
+
+
+def aircraft_load_factor(
+    origin_airport_icao: str | None = None,
+    first_waypoint_time: pd.Timestamp | None = None,
+    *,
+    freighter: bool = False,
+) -> float:
+    """
+    Estimate passenger/cargo load factor based on historical data.
+
+    Accounts for regional and seasonal differences.
+
+    Parameters
+    ----------
+    origin_airport_icao : str | None
+        ICAO code of origin airport. If None is provided, then globally averaged values will be
+        assumed at `first_waypoint_time`.
+    first_waypoint_time : pd.Timestamp | None
+        First waypoint UTC time. If None is provided, then regionally or globally averaged values
+        from the trailing twelve months will be used.
+    freighter: bool
+        Historical cargo load factor will be used if true, otherwise use passenger load factor.
+
+    Returns
+    -------
+    float
+        Passenger/cargo load factor [0 - 1], unitless
+    """
+    # If origin airport is provided, use regional load factor
+    if origin_airport_icao is not None:
+        first_letter = origin_airport_icao[0]
+        region = AIRPORT_TO_REGION.get(first_letter, "Global")
+    else:
+        region = "Global"
+
+    # Use passenger or cargo database
+    if freighter:
+        lf_database = _historical_regional_load_factor(CLF_PATH)
+    else:
+        lf_database = _historical_regional_load_factor(PLF_PATH)
+
+    # If `first_waypoint_time` is None, global/regional averages for the trailing twelve months
+    # will be assumed.
+    if first_waypoint_time is None:
+        t1 = lf_database.index[-1]
+        t0 = t1 - pd.DateOffset(months=12) + pd.DateOffset(days=1)
+        return lf_database.loc[t0:t1, region].mean().item()
+
+    date = first_waypoint_time.floor("D")
+
+    # If `date` is more recent than the historical data, then use most recent load factors
+    # from trailing twelve months as seasonal values are stable except in COVID years (2020-22).
+    if date > lf_database.index[-1]:
+        if date.month == 2 and date.day == 29:  # remove any leap day
+            date = date.replace(day=28)
+
+        filt = (lf_database.index.month == date.month) & (lf_database.index.day == date.day)
+        date = lf_database.index[filt][-1]
+
+    # (2) If `date` is before the historical data, then use 2019 load factors.
+    elif date < lf_database.index[0]:
+        if date.month == 2 and date.day == 29:  # remove any leap day
+            date = date.replace(day=28)
+
+        filt = (lf_database.index.month == date.month) & (lf_database.index.day == date.day)
+        date = lf_database.index[filt][0]
+
+    return lf_database.at[date, region].item()
+
+
 def aircraft_weight(aircraft_mass: ArrayOrFloat) -> ArrayOrFloat:
     """Calculate the aircraft weight at each waypoint.
 
@@ -413,7 +541,8 @@ def initial_aircraft_mass(
 
     See Also
     --------
-    :func:`reserve_fuel_requirements`
+    reserve_fuel_requirements
+    aircraft_load_factor
     """
     tom = operating_empty_weight + load_factor * max_payload + total_fuel_burn + total_reserve_fuel
     return min(tom, max_takeoff_weight)
@@ -463,9 +592,10 @@ def update_aircraft_mass(
 
     See Also
     --------
-    :func:`fuel_burn`
-    :func:`reserve_fuel_requirements`
-    :func:`initial_aircraft_mass`
+    fuel_burn
+    reserve_fuel_requirements
+    initial_aircraft_mass
+    aircraft_load_factor
     """
     if takeoff_mass is None:
         takeoff_mass = initial_aircraft_mass(

pycontrails/physics/static/iata-cargo-load-factors-20241115.csv

@@ -0,0 +1,71 @@
+Date,Global,Africa,Asia Pacific,Europe,Latin America,Middle East,North America
+15/12/2018,0.488,0.381,0.54,0.567,0.291,0.488,0.414
+15/1/2019,0.451,0.354,0.501,0.501,0.299,0.421,0.4
+15/2/2019,0.447,0.363,0.473,0.53,0.297,0.466,0.379
+15/3/2019,0.495,0.384,0.556,0.56,0.323,0.488,0.416
+15/4/2019,0.463,0.374,0.518,0.496,0.325,0.458,0.405
+15/5/2019,0.468,0.386,0.52,0.513,0.353,0.469,0.398
+15/6/2019,0.454,0.324,0.522,0.498,0.337,0.44,0.382
+15/7/2019,0.45,0.323,0.519,0.485,0.354,0.453,0.373
+15/8/2019,0.446,0.302,0.516,0.477,0.372,0.435,0.377
+15/9/2019,0.464,0.329,0.539,0.501,0.379,0.459,0.381
+15/10/2019,0.477,0.361,0.539,0.533,0.364,0.477,0.394
+15/11/2019,0.496,0.404,0.538,0.569,0.403,0.497,0.413
+15/12/2019,0.467,0.368,0.519,0.53,0.3,0.47,0.395
+15/1/2020,0.45,0.356,0.474,0.501,0.311,0.426,0.424
+15/2/2020,0.464,0.368,0.543,0.531,0.342,0.461,0.372
+15/3/2020,0.545,0.425,0.656,0.63,0.411,0.532,0.429
+15/4/2020,0.58,0.486,0.691,0.648,0.554,0.525,0.487
+15/5/2020,0.576,0.612,0.643,0.625,0.561,0.483,0.526
+15/6/2020,0.573,0.547,0.645,0.62,0.512,0.494,0.521
+15/7/2020,0.564,0.489,0.639,0.594,0.464,0.53,0.506
+15/8/2020,0.548,0.502,0.616,0.568,0.478,0.535,0.489
+15/9/2020,0.569,0.507,0.642,0.62,0.456,0.579,0.484
+15/10/2020,0.576,0.502,0.617,0.651,0.443,0.606,0.496
+15/11/2020,0.582,0.496,0.631,0.655,0.436,0.6,0.5
+15/12/2020,0.573,0.51,0.639,0.653,0.367,0.597,0.482
+15/1/2021,0.589,0.48,0.665,0.627,0.39,0.569,0.532
+15/2/2021,0.575,0.476,0.692,0.641,0.429,0.598,0.453
+15/3/2021,0.588,0.499,0.661,0.685,0.453,0.613,0.472
+15/4/2021,0.578,0.504,0.633,0.681,0.457,0.598,0.473
+15/5/2021,0.572,0.502,0.646,0.656,0.423,0.589,0.469
+15/6/2021,0.565,0.48,0.676,0.626,0.381,0.581,0.458
+15/7/2021,0.544,0.455,0.654,0.598,0.387,0.536,0.443
+15/8/2021,0.542,0.43,0.698,0.575,0.404,0.529,0.437
+15/9/2021,0.553,0.428,0.68,0.604,0.37,0.558,0.447
+15/10/2021,0.561,0.45,0.661,0.626,0.421,0.572,0.449
+15/11/2021,0.559,0.434,0.654,0.631,0.446,0.572,0.444
+15/12/2021,0.542,0.502,0.634,0.623,0.413,0.556,0.43
+15/1/2022,0.541,0.492,0.609,0.584,0.417,0.513,0.474
+15/2/2022,0.532,0.502,0.592,0.636,0.476,0.529,0.429
+15/3/2022,0.549,0.494,0.638,0.671,0.448,0.526,0.442
+15/4/2022,0.516,0.49,0.631,0.578,0.419,0.504,0.419
+15/5/2022,0.505,0.495,0.627,0.548,0.387,0.487,0.411
+15/6/2022,0.492,0.447,0.608,0.507,0.383,0.488,0.404
+15/7/2022,0.472,0.452,0.563,0.493,0.374,0.469,0.398
+15/8/2022,0.467,0.418,0.547,0.502,0.374,0.466,0.393
+15/9/2022,0.481,0.451,0.572,0.528,0.381,0.478,0.396
+15/10/2022,0.487,0.437,0.561,0.558,0.384,0.48,0.401
+15/11/2022,0.491,0.458,0.545,0.569,0.382,0.475,0.419
+15/12/2022,0.472,0.432,0.528,0.559,0.322,0.454,0.406
+15/1/2023,0.448,0.439,0.452,0.541,0.325,0.411,0.423
+15/2/2023,0.456,0.468,0.464,0.574,0.361,0.445,0.4
+15/3/2023,0.462,0.489,0.485,0.57,0.366,0.456,0.393
+15/4/2023,0.427,0.482,0.442,0.497,0.364,0.431,0.373
+15/5/2023,0.415,0.448,0.422,0.489,0.333,0.41,0.373
+15/6/2023,0.432,0.446,0.468,0.476,0.337,0.446,0.374
+15/7/2023,0.421,0.417,0.457,0.472,0.322,0.411,0.37
+15/8/2023,0.42,0.388,0.443,0.484,0.326,0.407,0.377
+15/9/2023,0.438,0.436,0.466,0.5,0.319,0.424,0.392
+15/10/2023,0.452,0.416,0.472,0.53,0.354,0.46,0.392
+15/11/2023,0.467,0.421,0.479,0.57,0.363,0.469,0.408
+15/12/2023,0.459,0.41,0.479,0.562,0.316,0.455,0.403
+15/1/2024,0.457,0.431,0.446,0.555,0.344,0.439,0.435
+15/2/2024,0.451,0.451,0.432,0.584,0.376,0.463,0.396
+15/3/2024,0.473,0.473,0.475,0.581,0.402,0.496,0.404
+15/4/2024,0.439,0.429,0.445,0.515,0.387,0.447,0.387
+15/5/2024,0.446,0.438,0.453,0.518,0.362,0.461,0.397
+15/6/2024,0.458,0.385,0.496,0.507,0.336,0.473,0.388
+15/7/2024,0.444,0.4,0.48,0.496,0.338,0.458,0.382
+15/8/2024,0.44,0.378,0.466,0.501,0.359,0.445,0.387
+15/9/2024,0.456,0.392,0.485,0.525,0.368,0.474,0.389