How to show 'Potential Contrail region' in plot?

[ashk12]

I'm plotting contrails for a flight in jupyter notebook in my local machine. As you can see in the left of this image, I also want to show 'Potential Contrail region'. Can you please help me how can I do that? Screenshot is from map.contrails.org

[zebengberg]

You're interested in plotting the vertical "profile" view, similar to your screenshot?

It's a little tricky to do this natively in pycontrails. We have an example notebook in the Contrails API with some better explanations.

Below is a quick and crude way to do this in pycontrails.

import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
from pycontrails import Flight, GeoVectorDataset
from pycontrails.models.pcr import PCR
from pycontrails.datalib.ecmwf import ERA5
from pycontrails.models.humidity_scaling import ConstantHumidityScaling

# Create some synthetic toy flight
fl = Flight(longitude=[0, 20], latitude=[0, 20], level=[280, 250], time=["2023-10-24T06", "2023-10-24T10"])
fl = fl.resample_and_fill("1T", nominal_rocd=1)

# Download some weather data around the flight
era5 = ERA5(
    time=("2023-10-24T06", "2023-10-24T10"),
    grid=1,
    variables=["t", "q"],
    pressure_levels=[200, 225, 250, 300],
)
met = era5.open_metdataset()

# Define a vector with the same footprint as fl by taking the product with a few flight levels of interest
lon = fl["longitude"]
lat = fl["latitude"]
time = fl["time"]
altitude_ft = np.arange(31000, 39000, 1000, dtype=float)

lon2d, lat2d, time2d, altitude_ft2d = np.broadcast_arrays(
    lon[np.newaxis, :], lat[np.newaxis, :], time[np.newaxis, :], altitude_ft[:, np.newaxis]
)

vector = GeoVectorDataset(
    longitude=lon2d.ravel(),
    latitude=lat2d.ravel(),
    time=time2d.ravel(),
    altitude_ft=altitude_ft2d.ravel(),
)

# Run the PCR model, reshape back to a 2D array
model = PCR(met=met, humidity_scaling=ConstantHumidityScaling(rhi_adj=0.6))
pcr = model.eval(vector)["pcr"]
pcr = pcr.reshape(lon2d.shape)
da = xr.DataArray(pcr, dims=("altitude_ft", "time"), coords={"altitude_ft": altitude_ft, "time": time})

# Plot both the 2D array and the flight together
fig, ax = plt.subplots()
da.plot(x="time", y="altitude_ft", ax=ax, cmap="coolwarm", alpha=0.5)
ax.plot(fl["time"], fl.altitude_ft, color="black")