ARCO ERA5

This notebook demonstrates how to load ARCO ERA5 data from Google Cloud Storage through the pycontrails ERA5ARCO interface.

import dask.config
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import xarray as xr

from pycontrails import Flight, MetDataset
from pycontrails.datalib.ecmwf import ERA5ARCO
from pycontrails.models.cocip import Cocip
from pycontrails.models.humidity_scaling import ConstantHumidityScaling
from pycontrails.models.issr import ISSR
from pycontrails.models.ps_model import PSFlight
from pycontrails.physics import units

Load 13 hours of pressure level and single level data from the ARCO ERA5 dataset. Select all the variables needed for the Cocip model.

Remove the cachestore=None keyword parameter to cache a copy of the data locally. This will make re-running the notebook much faster.

time = ("2019-01-01T00", "2019-01-01T12")

era5_pl = ERA5ARCO(
    time=time,
    variables=(*Cocip.met_variables, *Cocip.optional_met_variables),
    cachestore=None,
)
met = era5_pl.open_metdataset()


era5_sl = ERA5ARCO(
    time=time,
    variables=Cocip.rad_variables,
    pressure_levels=-1,
    cachestore=None,
)
rad = era5_sl.open_metdataset()

# Materialize the data into memory.
# We run the computation in single-threaded mode to use a bit less memory
# at the expense of a slower computation time.
# Under the hood, this computes the model level to pressure level conversion
# on one zarr chunk at a time rather than computing several across threads

with dask.config.set(scheduler="single-threaded"):
    met.data.load()
    rad.data.load()

ISSRs

Visualize ISSRs over common cruise altitudes over the continental US.

issr = ISSR(met=met, humidity_scaling=ConstantHumidityScaling(rhi_adj=0.95))
source = MetDataset(
    xr.Dataset(
        coords={
            "time": ["2019-01-01T00"],
            "longitude": np.arange(-125, -60, 0.25),
            "latitude": np.arange(25, 50, 0.25),
            "level": units.ft_to_pl(np.arange(27000, 45000, 1000)),
        }
    )
)
da = issr.eval(source).data["issr"]

# Use altitude_ft as the vertical coordinate for plotting
da["altitude_ft"] = units.pl_to_ft(da["level"]).round().astype(int)
da = da.swap_dims(level="altitude_ft")
da = da.sel(altitude_ft=da["altitude_ft"].values[::-1])
da = da.squeeze()

da.plot(x="longitude", y="latitude", col="altitude_ft", col_wrap=4, add_colorbar=False);

CoCiP

Run CoCiP on a collection of synthetic flights at different cruise altitudes over the continental US.

The flights constructed below are in no way realistic. They are simply used to showcase running CoCiP with model level met data.

jfk = 40.64, -73.78
lax = 33.94, -118.40

fl_list = [
    Flight(
        longitude=[lax[1], jfk[1]],
        latitude=[lax[0], jfk[0]],
        altitude_ft=[altitude_ft, altitude_ft],
        time=[np.datetime64("2019-01-01T00"), np.datetime64("2019-01-01T04:30")],
        aircraft_type="B738",
        flight_id=f"cruise {altitude_ft} ft",
        origin="LAX",
        destination="JFK",
    ).resample_and_fill()
    for altitude_ft in range(27000, 45000, 1000)
]

cocip = Cocip(
    met=met,
    rad=rad,
    dt_integration="1min",
    max_age="8h",
    aircraft_performance=PSFlight(),
    humidity_scaling=ConstantHumidityScaling(rhi_adj=0.95),
)
fl_list = cocip.eval(fl_list)

Visualize CoCiP results

Visualize the results of CoCiP on the synthetic flights.

summary = {fl.attrs["flight_id"]: fl["ef"].sum() / 1e12 for fl in fl_list}
pd.Series(summary).to_frame("Energy Forcing [TJ]").astype(int)

	Energy Forcing [TJ]
cruise 27000 ft	0
cruise 28000 ft	0
cruise 29000 ft	44
cruise 30000 ft	18
cruise 31000 ft	0
cruise 32000 ft	0
cruise 33000 ft	0
cruise 34000 ft	0
cruise 35000 ft	51
cruise 36000 ft	92
cruise 37000 ft	216
cruise 38000 ft	193
cruise 39000 ft	156
cruise 40000 ft	186
cruise 41000 ft	42
cruise 42000 ft	8
cruise 43000 ft	0
cruise 44000 ft	0

fig, ax = plt.subplots(figsize=(10, 6))

for fl in fl_list:
    x = fl["time"]
    y = fl["ef"]
    ax.plot(x, y, label=fl.attrs["flight_id"])
    ax.fill_between(x, y, alpha=0.1)

ax.set_ylabel("Energy forcing [J]")
ax.set_xlabel("Waypoint time")

for tick in ax.get_xticklabels():
    tick.set_rotation(45)

ax.set_yscale("symlog", linthresh=1e10)
ax.set_ylim(-1e9, 3e13)

ax.legend();