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...

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
)

Computing indices on weather forecasts¶

Previous topic

Next topic

This Page