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> Size: 4GB Dimensions: (lon: 720, lat: 361, time: 97, member: 20) Coordinates: * 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 2025-02-12 ... 2025-02-28 * 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_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.groupby("time").mean().hvplot.hist(legend=False, width=300) + cdd.groupby(
"time"
).mean().hvplot.hist(legend=False, width=300)