Example Forecast Data

Inspect example contrail forecast data source.

Download data

# download example data file
file_id = "1HftWWr4Z9ho_mjm1kbaL35SnQS9pzyjG"
!wget "https://drive.usercontent.google.com/download?id={file_id}&export=download" -O forecast.nc

--2024-11-13 12:38:33--  https://drive.usercontent.google.com/download?id=1HftWWr4Z9ho_mjm1kbaL35SnQS9pzyjG&export=download
Resolving drive.usercontent.google.com (drive.usercontent.google.com)... 142.251.41.1
Connecting to drive.usercontent.google.com (drive.usercontent.google.com)|142.251.41.1|:443... connected.
HTTP request sent, awaiting response... 200 OK
Length: 66477996 (63M) [application/octet-stream]
Saving to: ‘forecast.nc’

forecast.nc         100%[===================>]  63.40M  8.97MB/s    in 6.6s

2024-11-13 12:38:43 (9.57 MB/s) - ‘forecast.nc’ saved [66477996/66477996]

Load in xarray

import xarray as xr

ds = xr.load_dataset("forecast.nc")
ds

<xarray.Dataset>
Dimensions:                  (longitude: 1440, latitude: 641, flight_level: 18,
                              time: 1)
Coordinates:
  * longitude                (longitude) float32 -180.0 -179.8 ... 179.5 179.8
  * latitude                 (latitude) float32 -80.0 -79.75 ... 79.75 80.0
  * flight_level             (flight_level) int32 270 280 290 ... 420 430 440
  * time                     (time) datetime64[ns] 2024-09-25T18:00:00
    forecast_reference_time  (time) datetime64[ns] 2024-09-25T12:00:00
Data variables:
    contrails                (longitude, latitude, flight_level, time) float32 ...
Attributes:
    aircraft_class:  default
    model:           v1

xarray.Dataset

• Dimensions:
  • longitude: 1440
  • latitude: 641
  • flight_level: 18
  • time: 1
• Coordinates: (5)
  • longitude
    (longitude)
    float32
    -180.0 -179.8 ... 179.5 179.8

    array([-180.  , -179.75, -179.5 , ...,  179.25,  179.5 ,  179.75],
          dtype=float32)

  • latitude
    (latitude)
    float32
    -80.0 -79.75 -79.5 ... 79.75 80.0

    array([-80.  , -79.75, -79.5 , ...,  79.5 ,  79.75,  80.  ], dtype=float32)

  • flight_level
    (flight_level)
    int32
    270 280 290 300 ... 410 420 430 440

    long_name :

    Flight level

    units :

    hectofeet

    array([270, 280, 290, 300, 310, 320, 330, 340, 350, 360, 370, 380, 390, 400,
           410, 420, 430, 440], dtype=int32)

  • time
    (time)
    datetime64[ns]
    2024-09-25T18:00:00

    array(['2024-09-25T18:00:00.000000000'], dtype='datetime64[ns]')

  • forecast_reference_time
    (time)
    datetime64[ns]
    2024-09-25T12:00:00

    array(['2024-09-25T12:00:00.000000000'], dtype='datetime64[ns]')

• Data variables: (1)
  • contrails
    (longitude, latitude, flight_level, time)
    float32
    0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0

    long_name :

    Contrail forcing index

    units :

    valid_min :

    0

    valid_max :

    4

    array([[[[0.        ],
             [0.        ],
             [0.        ],
             ...,
             [0.        ],
             [0.        ],
             [0.        ]],
    
            [[0.        ],
             [0.        ],
             [0.        ],
             ...,
             [0.        ],
             [0.        ],
             [0.        ]],
    
            [[0.03977207],
             [0.15373209],
             [0.31042677],
             ...,
    ...
             ...,
             [0.        ],
             [0.        ],
             [0.        ]],
    
            [[0.54221934],
             [0.8032166 ],
             [1.3094051 ],
             ...,
             [0.        ],
             [0.        ],
             [0.        ]],
    
            [[0.593475  ],
             [0.8938293 ],
             [1.4172134 ],
             ...,
             [0.        ],
             [0.        ],
             [0.        ]]]], dtype=float32)

• Indexes: (4)
  • longitude
    PandasIndex

    PandasIndex(Index([ -180.0, -179.75,  -179.5, -179.25,  -179.0, -178.75,  -178.5, -178.25,
            -178.0, -177.75,
           ...
             177.5,  177.75,   178.0,  178.25,   178.5,  178.75,   179.0,  179.25,
             179.5,  179.75],
          dtype='float32', name='longitude', length=1440))

  • latitude
    PandasIndex

    PandasIndex(Index([ -80.0, -79.75,  -79.5, -79.25,  -79.0, -78.75,  -78.5, -78.25,  -78.0,
           -77.75,
           ...
            77.75,   78.0,  78.25,   78.5,  78.75,   79.0,  79.25,   79.5,  79.75,
             80.0],
          dtype='float32', name='latitude', length=641))

  • flight_level
    PandasIndex

    PandasIndex(Index([270, 280, 290, 300, 310, 320, 330, 340, 350, 360, 370, 380, 390, 400,
           410, 420, 430, 440],
          dtype='int32', name='flight_level'))

  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2024-09-25 18:00:00'], dtype='datetime64[ns]', name='time', freq=None))

• Attributes: (2)

  aircraft_class :

  default

  model :

  v1

Plot continous data

# Plot continuous data at Flight Level 320 for the first timestep
ds["contrails"] \
    .isel(time=0) \
    .sel(flight_level=320) \
    .plot(x="longitude", y="latitude", cmap="Reds", figsize=(12, 8));

Interpret as polygons

Generate polygons for each contrail forcing index

from pycontrails.core import polygon  # pip install pycontrails

import numpy as np
from matplotlib import pyplot as plt

# Create 2D dataarray from one time and flight level
da = ds["contrails"].isel(time=0).sel(flight_level=320)

# Thresholds of interest
thresholds = [1, 2, 3, 4]

# Polygon options
longitude = da["longitude"].to_numpy()
latitude = da["latitude"].to_numpy()
min_area = 0
epsilon = 0

# Set up figure
fig, ax = plt.subplots(figsize=(12,8))
cmap = plt.get_cmap('Reds')
colors = cmap(np.linspace(0.3, 1, len(thresholds)))

for (c, th) in zip(colors, thresholds):

    # calculate polygons for threshold
    multipolygon = polygon.find_multipolygon(
        da.values.copy(),
        threshold=th,
        min_area=min_area,
        epsilon=epsilon,
        longitude=longitude,
        latitude=latitude,
    )

    # plot all polygons
    for (i, poly) in enumerate(multipolygon.geoms):
        if i == 0:
            ax.plot(*poly.exterior.xy, color=c, lw=1, label=th)
        else:
            ax.plot(*poly.exterior.xy, color=c, lw=1)

plt.legend();
plt.xlabel("Longitude");
plt.ylabel("Latitude");