Reproducing Inventory Results

The Impact Inventory is the dataset backing our Impact Explorer. It contains contrail impact modeling for millions of flights. Each of the flights in the inventory was processed from raw ADS-B data through trajectory QA/QC processes, then run through CoCiP to model contrail formation and evolution. The process for ingesting and processing the ADS-B data, cleaning and sanitizing the trajectories, then modeling the flight contrail impacts with CoCiP are in our flights-pipeline. The results of this pipeline are the Impact Inventory and made available through the Impact Explorer.

In this notebook, we will examine how we configure and run CoCiP using two example flight trajectories as well as compare with results rom our research-oriented v0/trajectory/cocip endpoint, and results from the Impact Inventory itself (which are served from our v1/inventory endpoints).

# Start by importing dependencies
import os
from typing import List

import folium
import numpy as np
import pandas as pd
import requests
import xarray as xr
from dotenv import load_dotenv
from pycontrails import Flight, MetDataset
from pycontrails.models.cocip import Cocip
from pycontrails.models.humidity_scaling import (
    ExponentialBoostLatitudeCorrectionHumidityScaling,
)
from pycontrails.models.ps_model import PSFlight

# Load the .env file - it contains the api.contrails.org API key in the
# CONTRAILS_API_KEY environment variable
load_dotenv()

True

Let’s load two example flight trajectories. Here we’re using anonymized flight information, purely to illustrate our process for running CoCiP. These are the lat, lon, altitude trajectories for two flights that were sanitized from the raw ADS-B data by our flights-pipeline. These trajectories were then used as inputs to CoCiP to generate the flights’ entries in the Impact Inventory:

# read in files
example1 = pd.read_csv("example1_trajectory.csv")
example2 = pd.read_csv("example2_trajectory.csv")

example1.head()

	Unnamed: 0	longitude	latitude	time	altitude_ft	engine_uid	aircraft_type_icao	flight_id
0	0	-15.373030	27.958881	2024-01-20 18:20:00	1180.0	01P08CM105	A320	2XHCRwXWNF0nQem0AjILcQ==
1	1	-15.351441	28.011040	2024-01-20 18:21:00	2318.0	01P08CM105	A320	2XHCRwXWNF0nQem0AjILcQ==
2	2	-15.317189	28.072734	2024-01-20 18:22:00	4505.0	01P08CM105	A320	2XHCRwXWNF0nQem0AjILcQ==
3	3	-15.273730	28.135082	2024-01-20 18:23:00	7316.0	01P08CM105	A320	2XHCRwXWNF0nQem0AjILcQ==
4	4	-15.227630	28.201253	2024-01-20 18:24:00	9362.0	01P08CM105	A320	2XHCRwXWNF0nQem0AjILcQ==

print("Example 1 date range: {} - {}".format(example1["time"].min(), example1["time"].max()))
print("Example 2 date range: {} - {}".format(example2["time"].min(), example2["time"].max()))

Example 1 date range: 2024-01-20 18:20:00 - 2024-01-20 22:04:00
Example 2 date range: 2024-01-28 13:43:00 - 2024-01-28 17:51:00

Plot the flight trajectories

Let’s look at the flight trajectories we’re considering here.

# Some folium plotting functions for visualizing the flight paths


def add_points_to_map(map_obj, plot_df, color, correct_dateline_for_plot):
    points = []
    for i, (_, row) in enumerate(plot_df.iterrows()):
        lat = row["latitude"]
        lon = row["longitude"]
        alt = row["altitude_ft"]
        if correct_dateline_for_plot:
            if lon < 0:
                lon += 360
        else:
            speed = None
            point_color = color

        points.append((lat, lon))
        folium.CircleMarker(
            location=[lat, lon],
            tooltip=f'i={i}, timestamp={row["time"]}, point=({lat},{lon}), speed={speed}, altitude={alt}',
            radius=2,
            color=color,
        ).add_to(map_obj)

    folium.PolyLine(points, color=color, weight=1).add_to(map_obj)


def map_flight(plot_df1, correct_dateline_for_plot=True, bounds=None):
    map_obj = folium.Map()
    add_points_to_map(map_obj, plot_df1, "black", correct_dateline_for_plot)

    if bounds is None:
        bounds = map_obj.get_bounds()

    map_obj.fit_bounds(bounds)
    return map_obj

map_flight(example1, correct_dateline_for_plot=False)

Make this Notebook Trusted to load map: File -> Trust Notebook

map_flight(example2, correct_dateline_for_plot=False)

Make this Notebook Trusted to load map: File -> Trust Notebook

Let’s run CoCiP as we do for our inventory

The inventory is generated with our flights-pipeline which runs pycontrails at version 0.60.3. Here is where that is run in the pipeline:

https://github.com/contrailcirrus/flights-pipeline/blob/3df5dcf2948350abf2b2db9a13ce411fbfaffddd/trajectory-worker/lib/handlers.py#L411

Note that in this process, we pull in our own ERA5 meteorological data from our own data stores in Google Cloud. You will have to supply your own source of meterological data conforming to one of the pycontrails meterological data source. We use ECMWF ERA5 for the Impact Inventory.

Important Note: The meterological data source is quite important. Having sufficiently fine vertical resolution is necessary so that the relatively thin contrail-forming regions aren’t missed. We interpolate our ERA5 data from model levels 66-93 to pressure levels [134, 141, 148, 155, 163, 171, 180, 188, 197, 207, 217, 227, 237, 248, 260, 272, 284, 297, 310, 323, 337, 352, 367, 383, 399, 416].

Load ERA5 data

We allow for up to 12 hours of contrail evolution. Since both flights end past noon on their respective days, we load the following day’s ERA5 data as well so we have appropriate meteorological data to simulate contrail evolution into the next day for long-lived contrails.

Modify the google cloud storage URLs to match your own datasets to run this example.

# era5 zarr stores for 2024-01-20--21
# includes following day to allow for 12-hour contrail evolution
zarr_pl_stores_ex1 = [
    "gs://contrails-301217-ecmwf-era5-zarr-v2/20240120_pl.zarr",
    "gs://contrails-301217-ecmwf-era5-zarr-v2/20240121_pl.zarr",
]
zarr_sl_stores_ex1 = [
    "gs://contrails-301217-ecmwf-era5-zarr-v2/20240120_sl.zarr",
    "gs://contrails-301217-ecmwf-era5-zarr-v2/20240121_sl.zarr",
]

# era5 zarr stores for 2024-01-28--29
# includes following day to allow for 12-hour contrail evolution
zarr_pl_stores_ex2 = [
    "gs://contrails-301217-ecmwf-era5-zarr-v2/20240128_pl.zarr",
    "gs://contrails-301217-ecmwf-era5-zarr-v2/20240129_pl.zarr",
]
zarr_sl_stores_ex2 = [
    "gs://contrails-301217-ecmwf-era5-zarr-v2/20240128_sl.zarr",
    "gs://contrails-301217-ecmwf-era5-zarr-v2/20240129_sl.zarr",
]


def load_era5_zarr_stores(era5_pl_src: List[str], era5_sl_src: List[str]):
    pl_datasets = []
    sl_datasets = []
    for pl_src in era5_pl_src:
        pl_datasets.append(xr.open_zarr(pl_src))
    for sl_src in era5_sl_src:
        sl_datasets.append(xr.open_zarr(sl_src))

    pl = xr.concat(pl_datasets, dim="time")
    sl = xr.concat(sl_datasets, dim="time")

    met = MetDataset(pl, provider="ECMWF", dataset="ERA5", product="reanalysis")
    variables = Cocip.ecmwf_met_variables()
    met = met.standardize_variables(variables)

    rad = MetDataset(sl, provider="ECMWF", dataset="ERA5", product="reanalysis")
    variables = Cocip.ecmwf_rad_variables()
    rad = rad.standardize_variables(variables)

    return met, rad


met_ex1, rad_ex1 = load_era5_zarr_stores(zarr_pl_stores_ex1, zarr_sl_stores_ex1)
met_ex2, rad_ex2 = load_era5_zarr_stores(zarr_pl_stores_ex2, zarr_sl_stores_ex2)

Run CoCiP

This is how we run CoCiP for our inventory.

First check that aircraft models are available in the Paul-Schumann PS model. In the flights-pipeline, we fall back to the BADA model if the PS model is unavailable for an aircraft, and eject it if BADA is also unavailable.

ps_model = PSFlight(
    fill_low_altitude_with_isa_temperature=True,
    fill_low_altitude_with_zero_wind=True,
)

print(
    "Example 1 aircraft type {} available: {}".format(
        example1["aircraft_type_icao"][0],
        ps_model.check_aircraft_type_availability(aircraft_type=example1["aircraft_type_icao"][0]),
    )
)

print(
    "Example 2 aircraft type {} available: {}".format(
        example2["aircraft_type_icao"][0],
        ps_model.check_aircraft_type_availability(aircraft_type=example2["aircraft_type_icao"][0]),
    )
)

Example 1 aircraft type A320 available: True
Example 2 aircraft type B38M available: True

Set up params and helper functions to run CoCiP:

STATIC_PARAMS = dict(
    humidity_scaling=ExponentialBoostLatitudeCorrectionHumidityScaling(),
    dt_integration="5min",
    max_altitude_m=None,
    min_altitude_m=None,
    interpolation_use_indices=True,
    interpolation_bounds_error=False,
    filter_sac=True,
    copy_source=True,
    met_longitude_buffer=(
        10.0,
        10.0,
    ),  # default; potential perf gains fomr reducing
    met_latitude_buffer=(10.0, 10.0),  # default; potential perf gains from reducing
    met_level_buffer=(20, 20),  # reduced to same buffer used in api preprocessor
    max_age=np.timedelta64(12, "h"),
)


def create_flight(dataset) -> Flight:
    """
    Create Flight from trajectory.
    """
    engine_uid = dataset["engine_uid"][0]
    flight_id = dataset["flight_id"][0]
    aircraft_type = dataset["aircraft_type_icao"][0]

    return Flight(
        longitude=list(dataset["longitude"]),
        latitude=list(dataset["latitude"]),
        altitude_ft=list(dataset["altitude_ft"]),
        time=list(pd.to_datetime(dataset["time"])),
        attrs=dict(
            flight_id=flight_id,
            aircraft_type=aircraft_type,
            engine_uid=engine_uid,
        ),
    )


def run(flight: Flight, met_dataset: MetDataset, rad_dataset: MetDataset) -> Flight:
    """
    Run the cocip trajectory model.
    """
    perf_model = PSFlight(
        fill_low_altitude_with_isa_temperature=True,
        fill_low_altitude_with_zero_wind=True,
    )

    model = Cocip(
        met=met_dataset, rad=rad_dataset, aircraft_performance=perf_model, **STATIC_PARAMS
    )

    result = model.eval(flight)
    return result, model

flight_ex1 = create_flight(example1)
result_ex1, _ = run(flight_ex1, met_ex1, rad_ex1)
flight_ex2 = create_flight(example2)
result_ex2, _ = run(flight_ex2, met_ex2, rad_ex2)

print(f"CoCiP energy forcing for example 1: {result_ex1['ef'].sum():e} J")
print(f"CoCiP energy forcing for example 2: {result_ex2['ef'].sum():e} J")

CoCiP energy forcing for example 1: 1.628444e+15 J
CoCiP energy forcing for example 2: 2.908131e+13 J

Compare with results from the impact inventory

Let’s pull in some (lightly edited/anonymized) impact inventory results. These files are the result of pulling impact inventory results from the v1/inventory/segments API endpoint, and select out only results matching the example flights above. The files were lightly sanitized to remove some identifying information, but not in a way that affects the trajectory or results.

example1_inventory = pd.read_csv("example1_inventory.csv")
example2_inventory = pd.read_csv("example2_inventory.csv")
print(f"Result from example 1 inventory: {example1_inventory['sum_ef_mj'].sum()*10**6:e} J")
print(f"Result from example 2 inventory: {example2_inventory['sum_ef_mj'].sum()*10**6:e} J")

Result from example 1 inventory: 1.628445e+15 J
Result from example 2 inventory: 2.908129e+13 J

These are within rounding errors of the same results. The impact inventory energy forcing results are stored as integer megajoule quantities for each 1-minute segment, so summing these over a whole flight produces small errors compared with the values directly from CoCiP which are not truncated to round numbers of MJ.

Compare to using inventory segments trajectory to run CoCiP

Results from the v1/inventory/segments endpoint include full trajectories. Each 1-minute segment is an output from CoCiP and has a lat_start, lat_end, lon_start, lon_end, time_start, time_end, and mean_altitude_ft. CoCiP breaks segments up in a forward-looking way. That is, the 0th element represents the segment from the 0th to 1st elements of the input. The final element is dropped, so CoCiP outputs are one element shorter than the trajectory inputs.

The key difference here is that we have stored the mean altitude in feet for each segment, rather than starting and ending altitudes. This means the altitude does not exactly match what is in the input trajectory. This ultimately causes small differences in CoCiP outputs using the trajectory included in the inventory to run CoCiP again. We will see how much of a difference that makes in these two example flights.

# Set up Flights objects from impact inventory dataframes
def create_flight_from_result(dataset) -> Flight:
    """
    Create a Flight from our impact dataset.
    """
    engine_uid = dataset["engine_uid"][0]
    flight_id = dataset["flight_id"][0]
    aircraft_type = dataset["aircraft_type_icao"][0]

    return Flight(
        longitude=list(dataset["lon_start"]),
        latitude=list(dataset["lat_start"]),
        altitude_ft=list(dataset["median_altitude_ft"]),
        time=list(pd.to_datetime(dataset["time_start"])),
        attrs=dict(
            flight_id=flight_id,
            aircraft_type=aircraft_type,
            engine_uid=engine_uid,
        ),
    )


# Create Flight objects
fl_ex1_inv = create_flight_from_result(example1_inventory)
fl_ex2_inv = create_flight_from_result(example2_inventory)

# Run CoCiP - same met data, slightly different flight trajectory data
result_ex1_inv, _ = run(fl_ex1_inv, met_ex1, rad_ex1)
result_ex2_inv, _ = run(fl_ex2_inv, met_ex2, rad_ex2)

print(f"Result from example 1 inventory trajectory CoCiP run: {result_ex1_inv['ef'].sum():e} J")
print(f"Result from example 2 inventory trajectory CoCiP run: {result_ex2_inv['ef'].sum():e} J")

Result from example 1 inventory trajectory CoCiP run: 1.624297e+15 J
Result from example 2 inventory trajectory CoCiP run: 2.935056e+13 J

We see that while the results are close, the differences in altitudes while climbing and descending during the route have accumulated to a 0.25% increase in EF for example 1, and 0.92% decrease in EF for example 2.

Compare with v0/trajectory/cocip API result

Let’s compare the above methods with the output of our research-focused v0/trajectory/cocip API method. This method will run CoCiP against a full flight trajectory, but it makes default assumptions about the aircraft type and engine.

# Define credentials
URL = "https://api.contrails.org"
API_KEY = os.environ["CONTRAILS_API_KEY"]  # put in your API key here or in an .env file
HEADERS = {"x-api-key": API_KEY}

Run the trajectory through v0/trajectory/cocip

We set up the JSON objects to post to this endpoint below.

# Set up flight json objects for use with the API
fl_ex1 = example1.copy()

ex1_payload = example1.assign(time=example1["time"].astype(str)).to_dict("list")
ex1_payload["aircraft_type"] = example1["aircraft_type_icao"][0]

fl_ex2 = example2.copy()

ex2_payload = example2.assign(time=example2["time"].astype(str)).to_dict("list")
ex2_payload["aircraft_type"] = example2["aircraft_type_icao"][0]

# Make the request for example 1 flight
r = requests.post(f"{URL}/v0/trajectory/cocip", json=ex1_payload, headers=HEADERS)
print(f"HTTP Response Code: {r.status_code} {r.reason}")
r_json = r.json()
ef_traj = r_json["energy_forcing"]

ex1_ef = np.array(ef_traj).sum()
print(f"v0/trajectory/cocip energy forcing for Example 1: {ex1_ef:e} J")

HTTP Response Code: 200 OK
v0/trajectory/cocip energy forcing for Example 1: 1.253250e+15 J

# Make the request for example 2 flight
r = requests.post(f"{URL}/v0/trajectory/cocip", json=ex2_payload, headers=HEADERS)
print(f"HTTP Response Code: {r.status_code} {r.reason}")
r_json = r.json()
ef_traj = r_json["energy_forcing"]

ex2_ef = np.array(ef_traj).sum()
print(f"v0/trajectory/cocip energy forcing for Example 2: {ex2_ef:e} J")

HTTP Response Code: 200 OK
v0/trajectory/cocip energy forcing for Example 2: 2.593000e+13 J

print(
    f"v0/trajectory/copcip EF percentage of inventory example 1 EF: {100*ex1_ef/(example1_inventory['sum_ef_mj'].sum()*10**6):.2f}%"
)
print(
    f"v0/trajectory/copcip EF percentage of inventory example 2 EF: {100*ex2_ef/(example2_inventory['sum_ef_mj'].sum()*10**6):.2f}%"
)

v0/trajectory/copcip EF percentage of inventory example 1 EF: 76.96%
v0/trajectory/copcip EF percentage of inventory example 2 EF: 89.16%