ValueError when running CoCiP on a Flight object

[murakrishn]

Description

I am trying to run CoCiP as mentioned in this tutorial using pycontrails (version 0.54.2).

import numpy as np
import pandas as pd
from matplotlib import pyplot as plt

from pycontrails import Flight, MetDataset
from pycontrails.datalib.ecmwf import ERA5
from pycontrails.models.cocip import Cocip
from pycontrails.models.humidity_scaling import ConstantHumidityScaling

plt.rcParams["figure.figsize"] = (10, 6)

time_bounds = ("2019-03-01 00:00:00", "2019-03-01 23:00:00")
pressure_levels = (300, 250, 200)

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()

# demo synthetic flight
flight_attrs = {
    "flight_id": "test",
    # set constants along flight path
    "true_airspeed": 226.099920796651,  # true airspeed, m/s
    "thrust": 0.22,  # thrust_setting
    "nvpm_ei_n": 1.897462e15,  # non-volatile emissions index
    "aircraft_type": "E190",
    "wingspan": 48,  # m
    "n_engine": 2,
}

# Example flight
df = pd.DataFrame()
df["longitude"] = np.linspace(-25, -40, 100)
df["latitude"] = np.linspace(34, 40, 100)
df["altitude"] = np.linspace(10900, 10900, 100)
df["engine_efficiency"] = np.linspace(0.34, 0.35, 100)
df["fuel_flow"] = np.linspace(2.1, 2.4, 100)  # kg/s
df["aircraft_mass"] = np.linspace(154445, 154345, 100)  # kg
df["time"] = pd.date_range("2019-03-01T00:15:00", "2019-03-01T02:30:00", periods=100)

flight = Flight(df, attrs=flight_attrs)
flight

params = {
    "dt_integration": np.timedelta64(10, "m"),
    # The humidity_scaling parameter is only used for ECMWF ERA5 data
    # based on Teoh 2020 and Teoh 2022 - https://acp.copernicus.org/preprints/acp-2022-169/acp-2022-169.pdf
    # Here we use an example of constantly scaling the humidity value by 0.99
    "humidity_scaling": ConstantHumidityScaling(rhi_adj=0.99),
}
cocip = Cocip(met=met, rad=rad, params=params)

However, I am not able to instantiate a CoCiP as it errors out with the following trace:

Traceback (most recent call last):
  File "C:\Users\...\run_cocip_on_flight.py", line 61, in <module>
    cocip = Cocip(met=met, rad=rad, params=params)
  File "C:\Users\...\lib\site-packages\pycontrails\models\cocip\cocip.py", li
ne 317, in __init__
    self.met, self.rad = process_met_datasets(met, rad, compute_tau_cirrus)
  File "C:\Users\...\lib\site-packages\pycontrails\models\cocip\cocip.py", li
ne 1468, in process_met_datasets
    met = add_tau_cirrus(met)
  File "C:\Users\...\lib\site-packages\pycontrails\models\cocip\cocip.py", li
ne 1492, in add_tau_cirrus
    met.data["tau_cirrus"] = tau_cirrus.tau_cirrus(met)
  File "C:\Users\...\lib\site-packages\pycontrails\models\tau_cirrus.py", lin
e 83, in tau_cirrus
    dz = -dask.array.gradient(geopotential_height, axis=geopotential_height.get_axis_num("level"))
  File "C:\Users\...\lib\site-packages\dask\array\routines.py", line 701, in
gradient
    raise ValueError(
ValueError: Chunk size must be larger than edge_order + 1. Minimum chunk for axis 2 is 1. Rechunk to proceed.

Can anyone provide me with some guidance on what I could do to resolve this error?

Any help is much appreciated! TIA!

Details

• Version: 0.54.2
• Modules being used in the example: pycontrails.core.Flight, pycontrails.datalib.ecmwf, pycontrails.models.cocip, pycontrails.models.humidity_scaling`
• OS: Windows 11

[zebengberg]

Thanks for alerting us! This may have caused by a change to the netcdf chunks disseminated by ECMWF, but I'm not completely sure.

For now, you could use a workaround like this:

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()
met.data = met.data.chunk(level=-1)  # <---------------- RECHUNK HERE
rad = era5sl.open_metdataset()

I'll try to patch pycontrails to avoid this.

[murakrishn]

That fixes the error for me! Thank you so much!