Computing indices on weather forecasts

The PAVICS data catalog includes the latest weather forecast from the Global Ensemble Prediction System (GEPS) from Environment and Climate Change Canada. For the 20 members in the ensemble, two variables are available (precipitation and air temperature), every 3 hours for the first 8 days of the forecast, then every 6 hours for the following 8 days.

This notebook shows how to access the forecast data and compute climate indices using xclim. The first step is to open the catalog and get the URL to the data.

import warnings

import numba

warnings.simplefilter("ignore", category=numba.core.errors.NumbaDeprecationWarning)

import intake
import xarray as xr
import xclim

# fmt: off
forecast = ""  # TEST_USE_PROD_DATA
# fmt: on

cat = intake.open_esm_datastore(forecast)
license_type title institution member variable_name variable_long_name path
0 permissive Global Ensemble Prediction System (GEPS) - ECCC Canadian Meteorological Service - Montreal NaN ['pr', 'tas', 'member'] ['depth of water-equivalent precipitation', '2...
url = cat.df.path[0]
ds = xr.open_dataset(url)
<xarray.Dataset> Size: 4GB
Dimensions:  (lon: 720, lat: 361, time: 97, member: 20)
  * lon      (lon) float64 6kB 0.0 0.5 1.0 1.5 2.0 ... 358.0 358.5 359.0 359.5
  * lat      (lat) float64 3kB -90.0 -89.5 -89.0 -88.5 ... 88.5 89.0 89.5 90.0
    reftime  datetime64[ns] 8B ...
  * time     (time) datetime64[ns] 776B 2024-05-09 ... 2024-05-25
  * member   (member) float32 80B 1.0 2.0 3.0 4.0 5.0 ... 17.0 18.0 19.0 20.0
Data variables:
    pr       (member, time, lat, lon) float32 2GB ...
    tas      (member, time, lat, lon) float32 2GB ...
Attributes: (12/14)
    GRIB_edition:            2
    GRIB_centre:             cwao
    GRIB_centreDescription:  Canadian Meteorological Service - Montreal 
    GRIB_subCentre:          0
    Conventions:             CF-1.7
    institution:             Canadian Meteorological Service - Montreal 
    ...                      ...
    abstract:                Global ensemble forecasts are made twice a day u...
    dataset_id:              GEPS
    type:                    forecast
    license_type:            permissive

As the name suggests, GEPS is a global forecast, but to keep things simple, here we only pick one grid cell near Montréal. The following graph shows the time series of both precipitation and temperature. The member can be selected using the slider on the right of the graphics.

import hvplot.xarray

mtl = ds.sel(lon=45.5, lat=360 - 73.5, method="nearest")"time", width=300) + mtl.tas.hvplot(x="time", width=300)
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[3], line 1
----> 1 import hvplot.xarray
      3 mtl = ds.sel(lon=45.5, lat=360 - 73.5, method="nearest")
      4"time", width=300) + mtl.tas.hvplot(x="time", width=300)

ModuleNotFoundError: No module named 'hvplot'

Note that precipitation is cumulative starting from 0 at the beginning of the forecast. To get daily precipitation amounts, we need to first compute the 3- or 6-hourly precipitation, then sum daily.

from xclim.core.units import amount2rate

with xr.set_options(keep_attrs=True):
    pr_delta ="time")
    pr = pr_delta.resample(time="D").sum()
    # The units are mm, but for use with xclim, it's better to indicate this is a precipitation rate (mm/d).
    pr = amount2rate(pr, out_units="mm/d")
    # To avoid xclim warnings, we also fix some attributes to reflect the fact that precipitation is now a mean flux.
    pr.attrs["cell_methods"] = "time: mean"
    pr.attrs["standard_name"] = "precipitation_flux"

# It's then possible to compute xclim indicators, but of course, only on frequencies smaller than 16 days.
wetdays = xclim.atmos.wetdays(pr, freq="7D")

For temperature, here we’ll convert the time series to daily mean values. Once this is done, xclim indicators can be computed easily. Here, we compute heating and cooling degree days over periods of 7 days.

tas = mtl.tas.resample(time="D").mean(keep_attrs=True)
hdd = xclim.atmos.heating_degree_days(tas, freq="7D").dropna(dim="time")
cdd = xclim.atmos.cooling_degree_days(tas, freq="7D").dropna(dim="time")
hdd.hvplot.hist(groupby="time", legend=False, width=300) + cdd.hvplot.hist(
    groupby="time", legend=False, width=300