Run CoCiP over a flight

This tutorial walks through an example of running the Contrail Cirrus Predicition (CoCiP) model evaluation along a flight trajectory.

References

  • Schumann, U. “A Contrail Cirrus Prediction Model.” Geoscientific Model Development 5, no. 3 (May 3, 2012): 543–80. https://doi.org/10.5194/gmd-5-543-2012.

  • Schumann, U., B. Mayer, K. Graf, and H. Mannstein. “A Parametric Radiative Forcing Model for Contrail Cirrus.” Journal of Applied Meteorology and Climatology 51, no. 7 (July 2012): 1391–1406. https://doi.org/10.1175/JAMC-D-11-0242.1.

  • Teoh, Roger, Ulrich Schumann, Arnab Majumdar, and Marc E. J. Stettler. “Mitigating the Climate Forcing of Aircraft Contrails by Small-Scale Diversions and Technology Adoption.” Environmental Science & Technology 54, no. 5 (March 3, 2020): 2941–50. https://doi.org/10.1021/acs.est.9b05608.

  • Teoh, Roger, Ulrich Schumann, Edward Gryspeerdt, Marc Shapiro, Jarlath Molloy, George Koudis, Christiane Voigt, and Marc Stettler. “Aviation Contrail Climate Effects in the North Atlantic from 2016–2021.” Atmospheric Chemistry and Physics Discussions, March 30, 2022, 1–27. https://doi.org/10.5194/acp-2022-169.

[1]:

import pandas as pd
import xarray as xr

from pycontrails import Flight
from pycontrails.datalib.ecmwf import ERA5
from pycontrails.models.cocip import Cocip
from pycontrails.models.humidity_scaling import ConstantHumidityScaling
from pycontrails.physics import units

Load Flight

Load flight trajectory from dataset prepared by Roger Teoh in https://doi.org/10.5194/acp-2022-169

[2]:

# load flight waypoints
df_flight = pd.read_csv("data/flight-cocip.csv")
df_flight.head()

[2]:
Longitude (degrees) Latitude (degrees) Altitude (feet) UTC time True airspeed (m s-1) Mach Number Aircraft mass (kg) Fuel mass flow rate (kg s-1) Overall propulsion efficiency nvPM number emissions index (kg-1) ICAO Aircraft Type Wingspan (m)
0 -10.070 55.185 36000 1546651185 230.858 0.791 236479.000 1.654 0.4 1500000000000000 A359 64.75
1 -10.273 55.222 36000 1546651245 230.682 0.790 236379.755 1.657 0.4 1500000000000000 A359 64.75
2 -10.476 55.258 36000 1546651305 230.563 0.789 236280.355 1.659 0.4 1500000000000000 A359 64.75
3 -10.680 55.295 36000 1546651365 230.501 0.789 236180.791 1.661 0.4 1500000000000000 A359 64.75
4 -10.883 55.331 36000 1546651425 230.476 0.789 236081.128 1.662 0.4 1500000000000000 A359 64.75 
[3]:

# constant properties along the length of the flight
attrs = {
    "flight_id": "fid",
    "aircraft_type": df_flight["ICAO Aircraft Type"].values[0],
    "wingspan": df_flight["Wingspan (m)"].values[0],
}

Process the flight into a format expected by pycontrails. See pycontrails.Flight for interface details.

[4]:

# convert UTC timestamp to np.datetime64
df_flight["time"] = pd.to_datetime(df_flight["UTC time"], origin="unix", unit="s")

# set altitude in m
df_flight["altitude"] = units.ft_to_m(df_flight["Altitude (feet)"])

# rename a few columns for compatibility with `Flight` requirements
df_flight = df_flight.rename(
    columns={
        "Longitude (degrees)": "longitude",
        "Latitude (degrees)": "latitude",
        "True airspeed (m s-1)": "true_airspeed",
        "Mach Number": "mach_number",
        "Aircraft mass (kg)": "aircraft_mass",
        "Fuel mass flow rate (kg s-1)": "fuel_flow",
        "Overall propulsion efficiency": "engine_efficiency",
        "nvPM number emissions index (kg-1)": "nvpm_ei_n",
    }
)

# clean up a few columns before building Flight class
df_flight = df_flight.drop(
    columns=["ICAO Aircraft Type", "Wingspan (m)", "UTC time", "Altitude (feet)"]
)

fl = Flight(data=df_flight, attrs=attrs)
fl

[4]:
Flight [10 keys x 162 length, 4 attributes]

Attributes
time[2019-01-05 01:19:45, 2019-01-05 04:00:21]
longitude[-50.0, -10.07]
latitude[55.185, 61.089]
altitude[10972.8, 10972.8]
flight_idfid
aircraft_typeA359
wingspan64.75
crsEPSG:4326
longitude latitude true_airspeed mach_number aircraft_mass fuel_flow engine_efficiency nvpm_ei_n altitude time
0 -10.070 55.185 230.858 0.791 236479.000 1.654 0.4 1500000000000000 10972.8 2019-01-05 01:19:45
1 -10.273 55.222 230.682 0.790 236379.755 1.657 0.4 1500000000000000 10972.8 2019-01-05 01:20:45
2 -10.476 55.258 230.563 0.789 236280.355 1.659 0.4 1500000000000000 10972.8 2019-01-05 01:21:45
3 -10.680 55.295 230.501 0.789 236180.791 1.661 0.4 1500000000000000 10972.8 2019-01-05 01:22:45
4 -10.883 55.331 230.476 0.789 236081.128 1.662 0.4 1500000000000000 10972.8 2019-01-05 01:23:45
... ... ... ... ... ... ... ... ... ... ...
157 -49.002 61.027 254.806 0.860 220251.787 1.814 0.4 1500000000000000 10972.8 2019-01-05 03:56:41
158 -49.274 61.019 255.051 0.861 220142.920 1.830 0.4 1500000000000000 10972.8 2019-01-05 03:57:41
159 -49.546 61.010 255.317 0.862 220033.115 1.811 0.4 1500000000000000 10972.8 2019-01-05 03:58:41
160 -49.818 61.000 255.063 0.862 219924.426 0.309 0.4 1500000000000000 10972.8 2019-01-05 03:59:41
161 -50.000 61.000 235.015 0.794 219912.051 1.567 0.4 1500000000000000 10972.8 2019-01-05 04:00:21

162 rows × 10 columns

Load meteorology from ECMWF

[5]:

# get met domain from Flight
time = (
    pd.to_datetime(fl["time"][0]).floor("H"),
    pd.to_datetime(fl["time"][-1]).ceil("H") + pd.Timedelta("10H"),
)

# select pressure levels
pressure_levels = [
    400,
    350,
    300,
    250,
    225,
    200,
    175,
    150,
]

[6]:

# downloads met data from CDS
era5pl = ERA5(time=time, variables=Cocip.met_variables, pressure_levels=pressure_levels)
era5sl = ERA5(
    time=time,
    variables=Cocip.rad_variables,
)

[7]:

# create `MetDataset` from sources
met = era5pl.open_metdataset()
rad = era5sl.open_metdataset()

Set up model

[8]:

params = {
    "process_emissions": False,
    "verbose_outputs": True,
    "humidity_scaling": ConstantHumidityScaling(rhi_adj=0.98),
}
cocip = Cocip(met=met, rad=rad, params=params)

Run model

[9]:

fl_out = cocip.eval(source=fl)

/Users/marcshapiro/computing/contrailcirrus/pycontrails/pycontrails/models/cocip/cocip.py:2189: UserWarning: At time 2019-01-05T15:30:00.000000, the contrail has no intersection with the met data. This is likely due to the contrail being advected outside the met domain.
  warnings.warn(

Review output

The output flight has the original flight data with many new variables added from the evaluation.

[10]:

fl_out

[10]:
Flight [66 keys x 162 length, 9 attributes]

Attributes
time[2019-01-05 01:19:45, 2019-01-05 04:00:21]
longitude[-50.0, -10.07]
latitude[55.185, 61.089]
altitude[10972.8, 10972.8]
flight_idfid
aircraft_typeA359
wingspan64.75
crsEPSG:4326
rhi_adj0.98
humidity_scaling_nameconstant_scale
humidity_scaling_formularhi -> rhi / rhi_adj
humidity_scaling_rhi_adj0.98
pycontrails_version0.49.3.dev50
waypoint longitude latitude true_airspeed mach_number aircraft_mass fuel_flow engine_efficiency nvpm_ei_n altitude ... rf_sw_mean rf_sw_min rf_sw_max rf_lw_mean rf_lw_min rf_lw_max rf_net_mean rf_net_min rf_net_max cocip
0 0 -10.070 55.185 230.858 0.791 236479.000 1.654 0.4 1500000000000000 10972.8 ... NaN NaN NaN NaN NaN NaN NaN NaN NaN 0.0
1 1 -10.273 55.222 230.682 0.790 236379.755 1.657 0.4 1500000000000000 10972.8 ... NaN NaN NaN NaN NaN NaN NaN NaN NaN 0.0
2 2 -10.476 55.258 230.563 0.789 236280.355 1.659 0.4 1500000000000000 10972.8 ... NaN NaN NaN NaN NaN NaN NaN NaN NaN 0.0
3 3 -10.680 55.295 230.501 0.789 236180.791 1.661 0.4 1500000000000000 10972.8 ... NaN NaN NaN NaN NaN NaN NaN NaN NaN 0.0
4 4 -10.883 55.331 230.476 0.789 236081.128 1.662 0.4 1500000000000000 10972.8 ... NaN NaN NaN NaN NaN NaN NaN NaN NaN 0.0
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
157 157 -49.002 61.027 254.806 0.860 220251.787 1.814 0.4 1500000000000000 10972.8 ... NaN NaN NaN NaN NaN NaN NaN NaN NaN 0.0
158 158 -49.274 61.019 255.051 0.861 220142.920 1.830 0.4 1500000000000000 10972.8 ... NaN NaN NaN NaN NaN NaN NaN NaN NaN 0.0
159 159 -49.546 61.010 255.317 0.862 220033.115 1.811 0.4 1500000000000000 10972.8 ... NaN NaN NaN NaN NaN NaN NaN NaN NaN 0.0
160 160 -49.818 61.000 255.063 0.862 219924.426 0.309 0.4 1500000000000000 10972.8 ... NaN NaN NaN NaN NaN NaN NaN NaN NaN 0.0
161 161 -50.000 61.000 235.015 0.794 219912.051 1.567 0.4 1500000000000000 10972.8 ... NaN NaN NaN NaN NaN NaN NaN NaN NaN 0.0

162 rows × 66 columns

[11]:

fl_out.dataframe.columns

[11]:

Index(['waypoint', 'longitude', 'latitude', 'true_airspeed', 'mach_number',
       'aircraft_mass', 'fuel_flow', 'engine_efficiency', 'nvpm_ei_n',
       'altitude', 'time', 'flight_id', 'level', 'air_pressure',
       'air_temperature', 'specific_humidity', 'u_wind', 'v_wind', 'rhi',
       'tau_cirrus', 'specific_cloud_ice_water_content', 'rho_air', 'sdr',
       'segment_length', 'sin_a', 'cos_a', 'G', 'T_sat_liquid', 'rh',
       'rh_critical_sac', 'sac', 'T_critical_sac', 'width', 'depth', 'rhi_1',
       'air_temperature_1', 'specific_humidity_1', 'altitude_1',
       'persistent_1', 'dT_dz', 'ds_dz', 'dz_max', 'rho_air_1', 'iwc_1',
       'n_ice_per_m_1', 'ef', 'contrail_age', 'sdr_mean', 'sdr_min', 'sdr_max',
       'rsr_mean', 'rsr_min', 'rsr_max', 'olr_mean', 'olr_min', 'olr_max',
       'rf_sw_mean', 'rf_sw_min', 'rf_sw_max', 'rf_lw_mean', 'rf_lw_min',
       'rf_lw_max', 'rf_net_mean', 'rf_net_min', 'rf_net_max', 'cocip'],
      dtype='object')

The cocip variable describes where persistent contrails form. It can take on values:

  • 1: Persistent contrails form

  • 0: No persistent contrails form

[12]:

fl_out["cocip"]

[12]:

array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0.])

The model class contains information about the contrail created:

  • cocip.source the original input flight

  • cocip.contrail will be defined as a pandas DataFrame if a contrail is created.

  • cocip.contrail_dataset is the same data but formatted as an xarray Dataset.

[13]:

cocip.contrail

[13]:
waypoint flight_id formation_time time age longitude latitude altitude level continuous ... dn_dt_agg dn_dt_turb rf_sw rf_lw rf_net persistent ef timestep age_hours dt_integration
index
0 29 fid 2019-01-05 01:48:44 2019-01-05 02:00:00 0 days 00:11:16 -16.010456 56.334371 10929.200527 228.855145 True ... 1.567664e-20 0.000025 0.0 0.092565 0.092565 True 2.118821e+09 0 0.187778 0 days 00:11:16
1 30 fid 2019-01-05 01:49:44 2019-01-05 02:00:00 0 days 00:10:16 -16.243723 56.380602 10926.894766 228.938191 True ... 3.985451e-19 0.000027 0.0 1.043019 1.043019 True 5.918572e+09 0 0.171111 0 days 00:10:16
2 31 fid 2019-01-05 01:50:44 2019-01-05 02:00:00 0 days 00:09:16 -16.478084 56.425542 10924.667135 229.018447 True ... 6.609394e-19 0.000029 0.0 1.955334 1.955334 True 8.176894e+09 0 0.154444 0 days 00:09:16
3 32 fid 2019-01-05 01:51:44 2019-01-05 02:00:00 0 days 00:08:16 -16.714148 56.469048 10922.745458 229.087698 True ... 8.881788e-19 0.000030 0.0 2.946739 2.946739 True 9.090643e+09 0 0.137778 0 days 00:08:16
4 33 fid 2019-01-05 01:52:44 2019-01-05 02:00:00 0 days 00:07:16 -16.951516 56.510025 10921.141514 229.145512 True ... 9.416210e-19 0.000031 0.0 3.009561 3.009561 True 6.849343e+09 0 0.121111 0 days 00:07:16
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
0 55 fid 2019-01-05 02:14:48 2019-01-05 13:00:00 0 days 00:00:00 -9.957466 74.162396 10938.604760 228.516685 False ... 1.933380e-18 0.000002 0.0 3.184593 3.184593 True 0.000000e+00 22 0.000000 0 days 00:30:00
0 55 fid 2019-01-05 02:14:48 2019-01-05 13:30:00 0 days 00:00:00 -7.452062 74.581339 10914.404661 229.388471 False ... 1.674621e-18 0.000001 0.0 3.146509 3.146509 True 0.000000e+00 23 0.000000 0 days 00:30:00
0 55 fid 2019-01-05 02:14:48 2019-01-05 14:00:00 0 days 00:00:00 -4.724895 74.903333 10874.310228 230.838762 False ... 1.218293e-18 0.000001 0.0 2.435764 2.435764 True 0.000000e+00 24 0.000000 0 days 00:30:00
0 55 fid 2019-01-05 02:14:48 2019-01-05 14:30:00 0 days 00:00:00 -1.807077 75.115034 10854.571581 231.555468 False ... 1.016172e-18 0.000001 0.0 2.323749 2.323749 True 0.000000e+00 25 0.000000 0 days 00:30:00
0 55 fid 2019-01-05 02:14:48 2019-01-05 15:00:00 0 days 00:00:00 1.202617 75.218719 10856.300880 231.492606 False ... 1.791175e-18 0.000001 0.0 NaN NaN True 0.000000e+00 26 0.000000 0 days 00:30:00

230 rows × 57 columns

We can visualize the contrail on top of the original flight trajectory using pandas plotting capabilities

[14]:

ax = cocip.source.dataframe.plot(
    "longitude", "latitude", color="k", label=fl.attrs["flight_id"], figsize=(12, 8)
)
cocip.contrail.plot.scatter("longitude", "latitude", c="rf_lw", cmap="Reds", ax=ax);
../_images/notebooks_run-cocip-on-flight_24_0.png 
[15]:

ax = cocip.source.dataframe.plot(
    "longitude", "latitude", color="k", label=fl.attrs["flight_id"], figsize=(12, 8)
)
cocip.contrail.plot.scatter(
    "longitude", "latitude", c="ef", cmap="coolwarm", vmin=-1e12, vmax=1e12, ax=ax
);
../_images/notebooks_run-cocip-on-flight_25_0.png