ISSR

Model ice super-saturated regions (ISSR) of the atmosphere.

Met Data

• Requires account with Copernicus Data Portal and local ~/.cdsapirc file with credentials.

import matplotlib.pyplot as plt
import pandas as pd
from matplotlib.colors import ListedColormap

from pycontrails import Flight
from pycontrails.datalib.ecmwf import ERA5
from pycontrails.models.issr import ISSR

# ignore pycontrails warning about ECMWF humidity scaling
import warnings

warnings.filterwarnings(message=r"[\s\S]* humidity scaling [\s\S]*", action="ignore")

Get Data

time = ("2022-03-01 00:00:00", "2022-03-01 08:00:00")
pressure_levels = [300, 250, 200]
variables = ["t", "q"]

era5 = ERA5(time=time, variables=variables, pressure_levels=pressure_levels)
met = era5.open_metdataset()

Evaluate model

# run model for across full input domain
# outputs global ice super-saturated regions as 1.0, all other regions as 0.0
issr_mds = ISSR(met).eval()
issr = issr_mds["issr"]

# edge detection algorithm using differentiation to reduce the areas to lines
issr_edges = issr.find_edges()

da = issr.data.isel(time=0)
da.plot(x="longitude", y="latitude", row="level", cmap="Reds", figsize=(6, 12));

# plot issr edges for each pressure level
da = issr_edges.data.isel(time=0)
da.plot(x="longitude", y="latitude", row="level", cmap="Reds", figsize=(6, 12));

Interpolate

Run model along a flight path

# Load flight
df = pd.read_csv("data/flight.csv", parse_dates=["time"])
fl = Flight(data=df, flight_id="acdd1b", callsign="AAL1158")

fl

Flight [4 keys x 175 length, 3 attributes]

Attributes
time	[2022-03-01 00:50:00, 2022-03-01 03:47:00]
longitude	[-97.026, -77.036]
latitude	[32.931, 38.854]
altitude	[190.5, 11582.4]
flight_id	acdd1b
callsign	AAL1158
crs	EPSG:4326

	longitude	latitude	altitude	time
0	-77.035950	38.829315	236.22	2022-03-01 00:50:00
1	-77.038223	38.772675	708.66	2022-03-01 00:51:00
2	-77.114231	38.744568	9471.66	2022-03-01 00:52:00
3	-77.201965	38.739888	2019.30	2022-03-01 00:53:00
4	-77.286191	38.745117	3032.76	2022-03-01 00:54:00
...	...	...	...	...
170	-97.025925	32.931379	190.50	2022-03-01 03:43:00
171	-97.025922	32.930649	190.50	2022-03-01 03:44:00
172	-97.025922	32.930649	190.50	2022-03-01 03:45:00
173	-97.025922	32.930649	190.50	2022-03-01 03:46:00
174	-97.025922	32.930649	190.50	2022-03-01 03:47:00

175 rows × 4 columns

# run model for across full input domain
# outputs global ice super-saturated regions as 1, all other regions as 0
# np.nan is returned outside of the met domain
fl_out = ISSR(met=met).eval(source=fl)
fl_out["issr"]

array([nan, nan,  0., nan, nan, nan, nan, nan, nan, nan, nan, nan, nan,
        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.,  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.,  0.,  1., nan, nan, nan, nan,
       nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan,
       nan, nan, nan, nan, nan, nan])

# Get the length of the Flight in the ISSR region
fl_out.length_met("issr")

190812.60842238227

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

# Create colormap with red for ISSR and blue for non-ISSR
cmap = ListedColormap(["b", "r"])

ax.scatter(fl_out["longitude"], fl_out["latitude"], c=fl_out["issr"], cmap=cmap)


# Create legend
legend_elements = [
    plt.Line2D([0], [0], marker="o", color="w", label="ISSR", markerfacecolor="r", markersize=10),
    plt.Line2D(
        [0], [0], marker="o", color="w", label="non-ISSR", markerfacecolor="b", markersize=10
    ),
]
ax.legend(handles=legend_elements, loc="upper left")

ax.set(xlabel="longitude", ylabel="latitude");