Run CoCiP with FDR or QAR data#

References#

[1]:
import numpy as np
import xarray as xr
import pandas as pd
from matplotlib import pyplot as plt
import scipy

from pycontrails import Flight
from pycontrails.datalib.ecmwf import ERA5
from pycontrails.models.cocip import Cocip
from pycontrails.models.humidity_scaling import HistogramMatching
from pycontrails.models.ps_model import PSFlight
from pycontrails.models.emissions import Emissions

Download met data#

[2]:
time_bounds = ("2022-03-01 00:00:00", "2022-03-01 23:00:00")
pressure_levels = (300, 250, 200)  # 30,000 ft to 38,000 ft

era5pl = ERA5(
    time=time_bounds,
    variables=Cocip.met_variables + Cocip.optional_met_variables,
    pressure_levels=pressure_levels,
)
era5sl = ERA5(time=time_bounds, variables=Cocip.rad_variables)

# download data from ERA5 (or open from cache)
met = era5pl.open_metdataset()
rad = era5sl.open_metdataset()

Read in data and resample#

In this example, the sample data provides TAS, aircraft mass, and fuel flow per engine. Any values provided as input to the CoCiP model will not be overwritten when CoCiP is run, and so we should make sure we are saving the most accurate portions of the FDR data. Typically, thurst values provided by FDRs are noisy and innaccurate and so it is recommended to drop the thrust value if provided and recompute thrust through an aircraft performance model (below).

In this case, fuel flow is provided per engine, and so total fuel flow must be computed by summing the two values together.

Here we are also resampling the FDR data to a one minute sampling period, which is recommended when running CoCiP. Note that the resample_and_fill function will interpolate time, position, and altitude, but for the remaining columns, it will simply choose the nearsest value. There are much more accurate ways to resample fuel flow data, but this should generally be sufficient to estimate contrail impacts.

[3]:
attrs = {
    "flight_id": "test",
    "aircraft_type": "B77W",
    "engine_uid": "01P21GE217",  # 01P21GE217 -> GE90-115B
    # "n_engine":2 # This shouldn't be needed?
}

df = pd.read_csv("data/flight-fdr.csv")
df["fuel_flow"] = df["fuel_flow_1"] + df["fuel_flow_2"]  # Checked
fl = Flight(df, attrs=attrs)

fl = fl.resample_and_fill(freq="60s", drop=False)
fl
[3]:
Flight [10 keys x 136 length, 4 attributes]

Attributes
time[2022-03-01 00:15:00, 2022-03-01 02:30:00]
longitude[-39.926, -25.0]
latitude[34.0, 39.97]
altitude[10900.0, 10900.0]
flight_idtest
aircraft_typeB77W
engine_uid01P21GE217
crsEPSG:4326
longitude latitude altitude flight_id true_airspeed aircraft_mass fuel_flow_1 fuel_flow_2 fuel_flow time
0 -25.000000 34.000000 10900.0 test 215.628175 237650.498951 0.786975 0.786975 1.573950 2022-03-01 00:15:00
1 -25.111111 34.044444 10900.0 test 215.906118 237556.061693 0.786984 0.786984 1.573968 2022-03-01 00:16:00
2 -25.220370 34.088148 10900.0 test 216.176652 237463.196123 0.787013 0.787013 1.574027 2022-03-01 00:17:00
3 -25.331481 34.132593 10900.0 test 216.433880 237368.802459 0.786496 0.786496 1.572992 2022-03-01 00:18:00
4 -25.442593 34.177037 10900.0 test 216.684851 237274.417192 0.786600 0.786600 1.573200 2022-03-01 00:19:00
... ... ... ... ... ... ... ... ... ... ...
131 -39.483333 39.793333 10900.0 test 230.120409 225427.222884 0.775681 0.775681 1.551362 2022-03-01 02:26:00
132 -39.594444 39.837778 10900.0 test 230.055073 225334.320532 0.773240 0.773240 1.546479 2022-03-01 02:27:00
133 -39.703704 39.881481 10900.0 test 229.945472 225243.173069 0.771585 0.771585 1.543171 2022-03-01 02:28:00
134 -39.814815 39.925926 10900.0 test 229.772775 225150.752744 0.768882 0.768882 1.537763 2022-03-01 02:29:00
135 -39.925926 39.970370 10900.0 test 229.550583 225058.592246 0.767088 0.767088 1.534177 2022-03-01 02:30:00

136 rows × 10 columns

(Optional) Run Aircraft Performance Model#

The CoCiP module will automatically run an aircraft performance model to compute any missing values needed. In this case, we still need to compute engine efficiency and estimated thrust force. For completeness, we show how this can be computed directly from the an aircraft performance model.

[4]:
perf = PSFlight(met=met)
fp = perf.eval(fl)
[5]:
fig, ax = plt.subplots()
ax2 = ax.twinx()
ax3 = ax2.twinx()
ax2.set_yticks([])
fp.dataframe.plot(ax=ax, x="time", y="engine_efficiency", style="r", legend=False)
fp.dataframe.plot(ax=ax2, x="time", y="thrust", style="b", legend=False)
fp.dataframe.plot(ax=ax3, x="time", y="fuel_flow", style="k", legend=False)
_ = ax3.legend(
    [ax.get_lines()[0], ax2.get_lines()[0], ax3.get_lines()[0]],
    ["Engine Efficiency", "Thrust", "Fuel Flow"],
    bbox_to_anchor=(1.15, 0.5),
)
../_images/notebooks_run-cocip-with-fdr_10_0.png

Compute Aircraft Emissions#

Given fuel flow data, met data, aircraft, and engine type, we have sufficient information to run the emissions module and compute nvPM estimates needed as input for CoCiP. We have not specified engine type, so in the case, the default engine type will be assumed, which for the B77W is the GE-90 115B. To specific a different engine type, add the key engine_uid to the attrs dict when creating the Flight object.

Note that the emissions module estimates the aircraft thrust setting by comparing the fuel flow to the maximum fuel flow. This is the perferred way to estimate emissions using the ICAO emissions inventory and is the default behavior of the emissions module even when aircraft thrust is provided in the input.

[6]:
emissions = Emissions(met=met, humidity_scaling=HistogramMatching())
fl = emissions.eval(fl)
fl
[6]:
Flight [31 keys x 136 length, 17 attributes]

Attributes
time[2022-03-01 00:15:00, 2022-03-01 02:30:00]
longitude[-39.926, -25.0]
latitude[34.0, 39.97]
altitude[10900.0, 10900.0]
flight_idtest
aircraft_typeB77W
engine_uid01P21GE217
crsEPSG:4326
n_engine2
gaseous_data_sourceFFM2
nvpm_data_sourceICAO EDB
total_co239390.76441156896
total_h2o15337.3346711712
total_so214.963253337727998
total_sulphates0.30537251709648977
total_oc0.2493875556288
total_nox192.979281298076
total_co2.915210297723201
total_hc0.6230983842461807
total_nvpm_mass0.18164446517898938
total_nvpm_number3.54753797881968e+18
longitude latitude altitude flight_id true_airspeed aircraft_mass fuel_flow_1 fuel_flow_2 fuel_flow time ... co2 h2o so2 sulphates oc nox co hc nvpm_mass nvpm_number
0 -25.000000 34.000000 10900.0 test 215.628175 237650.498951 0.786975 0.786975 1.573950 2022-03-01 00:15:00 ... 298.326415 116.157483 0.113324 0.002313 0.001889 1.450262 0.021463 0.004587 0.001376 2.686732e+16
1 -25.111111 34.044444 10900.0 test 215.906118 237556.061693 0.786984 0.786984 1.573968 2022-03-01 00:16:00 ... 298.329927 116.158851 0.113326 0.002313 0.001889 1.450494 0.021460 0.004587 0.001376 2.686764e+16
2 -25.220370 34.088148 10900.0 test 216.176652 237463.196123 0.787013 0.787013 1.574027 2022-03-01 00:17:00 ... 298.340990 116.163158 0.113330 0.002313 0.001889 1.450803 0.021459 0.004587 0.001376 2.686863e+16
3 -25.331481 34.132593 10900.0 test 216.433880 237368.802459 0.786496 0.786496 1.572992 2022-03-01 00:18:00 ... 298.144830 116.086781 0.113255 0.002311 0.001888 1.449472 0.021447 0.004584 0.001375 2.685097e+16
4 -25.442593 34.177037 10900.0 test 216.684851 237274.417192 0.786600 0.786600 1.573200 2022-03-01 00:19:00 ... 298.184382 116.102181 0.113270 0.002312 0.001888 1.450112 0.021455 0.004586 0.001375 2.685453e+16
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
131 -39.483333 39.793333 10900.0 test 230.120409 225427.222884 0.775681 0.775681 1.551362 2022-03-01 02:26:00 ... 294.045243 114.490550 0.111698 0.002280 0.001862 1.444279 0.022017 0.004706 0.001356 2.648176e+16
132 -39.594444 39.837778 10900.0 test 230.055073 225334.320532 0.773240 0.773240 1.546479 2022-03-01 02:27:00 ... 293.119693 114.130175 0.111347 0.002272 0.001856 1.436100 0.021918 0.004685 0.001352 2.639840e+16
133 -39.703704 39.881481 10900.0 test 229.945472 225243.173069 0.771585 0.771585 1.543171 2022-03-01 02:28:00 ... 292.492624 113.886017 0.111108 0.002268 0.001852 1.430317 0.021840 0.004668 0.001349 2.634193e+16
134 -39.814815 39.925926 10900.0 test 229.772775 225150.752744 0.768882 0.768882 1.537763 2022-03-01 02:29:00 ... 291.467679 113.486941 0.110719 0.002260 0.001845 1.421186 0.021731 0.004645 0.001344 2.624962e+16
135 -39.925926 39.970370 10900.0 test 229.550583 225058.592246 0.767088 0.767088 1.534177 2022-03-01 02:30:00 ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN

136 rows × 31 columns

[7]:
fig, ax = plt.subplots()
ax2 = ax.twinx()
fl.dataframe.plot(ax=ax, x="time", y="nvpm_mass", style="r", legend=False)
fl.dataframe.plot(ax=ax2, x="time", y="fuel_flow", style="b", legend=False)
ax2.legend(
    [
        ax.get_lines()[0],
        ax2.get_lines()[0],
    ],
    ["nvPM Mass", "Fuel Flow"],
    bbox_to_anchor=(1.15, 0.5),
)
[7]:
<matplotlib.legend.Legend at 0x176d15010>
../_images/notebooks_run-cocip-with-fdr_13_1.png

Run CoCiP over the flight#

In order to predict contrail impact, we still need an estimate of engine efficiency, which is needed to determine if the Schmidt-Appleman Criteria is satisfied. If we provide the CoCiP module with an aircraft performance model, in this case the Poll-Schumann model, then this value will be estimated for us.

In the code below, the CoCiP module will run the Poll-Schmann model over the flight using the provided aircraft mass in order to estimate the thrust force required for the plane to fly the specified trajectory. This value will then be used along with the provided fuel flow data to estimate engine efficiency. Fuel flow and emissions numbers will not be recomputed here, and the warning from aircraft_performance.py can be ignored — three iterations of the performance model are only needed to estimate aircraft mass, which in this case is taken from the FDR data.

[8]:
cocip = Cocip(
    met=met, rad=rad, aircraft_performance=PSFlight(), humidity_scaling=HistogramMatching()
)
fl = cocip.eval(fl)
fl
[8]:
Flight [67 keys x 136 length, 25 attributes]

Attributes
time[2022-03-01 00:15:00, 2022-03-01 02:30:00]
longitude[-39.926, -25.0]
latitude[34.0, 39.97]
altitude[10900.0, 10900.0]
flight_idtest
aircraft_typeB77W
engine_uid01P21GE217
crsEPSG:4326
n_engine2
gaseous_data_sourceFFM2
nvpm_data_sourceICAO EDB
total_co239390.76441156896
total_h2o15337.3346711712
total_so214.963253337727998
total_sulphates0.30537251709648977
total_oc0.2493875556288
total_nox192.979281298076
total_co2.915210297723201
total_hc0.6230983842461807
total_nvpm_mass0.18164446517898938
total_nvpm_number3.54753797881968e+18
aircraft_performance_modelPSFlight
wingspan64.8
max_mach0.89
max_altitude13136.880000000001
total_fuel_burn12469.37778144
humidity_scaling_namehistogram_matching
humidity_scaling_formulaera5_quantiles -> iagos_quantiles
pycontrails_version0.49.3.dev50
waypoint longitude latitude altitude flight_id true_airspeed aircraft_mass fuel_flow_1 fuel_flow_2 fuel_flow ... n_ice_per_m_1 ef contrail_age sdr_mean rsr_mean olr_mean rf_sw_mean rf_lw_mean rf_net_mean cocip
0 0 -25.000000 34.000000 10900.0 test 215.628175 237650.498951 0.786975 0.786975 1.573950 ... NaN 0.0 0 days NaN NaN NaN NaN NaN NaN 0.0
1 1 -25.111111 34.044444 10900.0 test 215.906118 237556.061693 0.786984 0.786984 1.573968 ... NaN 0.0 0 days NaN NaN NaN NaN NaN NaN 0.0
2 2 -25.220370 34.088148 10900.0 test 216.176652 237463.196123 0.787013 0.787013 1.574027 ... NaN 0.0 0 days NaN NaN NaN NaN NaN NaN 0.0
3 3 -25.331481 34.132593 10900.0 test 216.433880 237368.802459 0.786496 0.786496 1.572992 ... NaN 0.0 0 days NaN NaN NaN NaN NaN NaN 0.0
4 4 -25.442593 34.177037 10900.0 test 216.684851 237274.417192 0.786600 0.786600 1.573200 ... 2.453534e+10 0.0 0 days NaN NaN NaN NaN NaN NaN 0.0
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
131 131 -39.483333 39.793333 10900.0 test 230.120409 225427.222884 0.775681 0.775681 1.551362 ... NaN 0.0 0 days NaN NaN NaN NaN NaN NaN 0.0
132 132 -39.594444 39.837778 10900.0 test 230.055073 225334.320532 0.773240 0.773240 1.546479 ... NaN 0.0 0 days NaN NaN NaN NaN NaN NaN 0.0
133 133 -39.703704 39.881481 10900.0 test 229.945472 225243.173069 0.771585 0.771585 1.543171 ... NaN 0.0 0 days NaN NaN NaN NaN NaN NaN 0.0
134 134 -39.814815 39.925926 10900.0 test 229.772775 225150.752744 0.768882 0.768882 1.537763 ... NaN 0.0 0 days NaN NaN NaN NaN NaN NaN 0.0
135 135 -39.925926 39.970370 10900.0 test 229.550583 225058.592246 0.767088 0.767088 1.534177 ... NaN 0.0 0 days NaN NaN NaN NaN NaN NaN 0.0

136 rows × 67 columns

Visualize Contrail Impact#

First, we plot cumulative EF by flight waypoint (one minute sample period)

[9]:
fl.dataframe.plot.scatter(
    x="longitude",
    y="latitude",
    c="ef",
    cmap="coolwarm",
    vmin=-1e13,
    vmax=1e13,
    title="EF generated by flight waypoint",
);
../_images/notebooks_run-cocip-with-fdr_17_0.png

Next, we plot the evolution of the contrail as it advects along each flight segment

[10]:
ax = plt.axes()

cocip.source.dataframe.plot(
    "longitude",
    "latitude",
    color="k",
    ax=ax,
    label="Flight path",
)
cocip.contrail.plot.scatter(
    "longitude",
    "latitude",
    c="rf_lw",
    cmap="Reds",
    ax=ax,
    label="Contrail LW RF",  # Contrail age?
)
ax.legend();
../_images/notebooks_run-cocip-with-fdr_19_0.png