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 = "https://pavics.ouranos.ca/catalog/forecast.json"  # TEST_USE_PROD_DATA
# fmt: on

cat = intake.open_esm_datastore(forecast)
cat.df.head()
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... https://pavics.ouranos.ca/twitcher/ows/proxy/t...
# NBVAL_IGNORE_OUTPUT
url = cat.df.path[0]
ds = xr.open_dataset(url)
ds
<xarray.Dataset>
Dimensions:  (lon: 720, lat: 361, time: 97, member: 20)
Coordinates:
  * lon      (lon) float64 0.0 0.5 1.0 1.5 2.0 ... 357.5 358.0 358.5 359.0 359.5
  * lat      (lat) float64 -90.0 -89.5 -89.0 -88.5 -88.0 ... 88.5 89.0 89.5 90.0
    reftime  datetime64[ns] ...
  * time     (time) datetime64[ns] 2024-01-19 2024-01-19T03:00:00 ... 2024-02-04
  * member   (member) float32 1.0 2.0 3.0 4.0 5.0 ... 16.0 17.0 18.0 19.0 20.0
Data variables:
    pr       (member, time, lat, lon) float32 ...
    tas      (member, time, lat, lon) float32 ...
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_description:     https://weather.gc.ca/grib/grib2_ens_geps_e.html
    dataset_id:              GEPS
    type:                    forecast
    license_type:            permissive
    license:                 https://open.canada.ca/en/open-government-licenc...

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")
mtl.pr.hvplot(x="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 mtl.pr.hvplot(x="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 = mtl.pr.diff(dim="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
)