Hands-On Tutorial

[1]:

import numpy as np
import xarray as xr
import random
import matplotlib.pyplot as plt

np.random.seed(0)
random.seed(0)

Define some functions

[2]:

historical_time = xr.cftime_range("1971-01-01", "2000-12-31", freq="D", calendar="noleap")
future_time = xr.cftime_range("2001-01-01", "2030-12-31", freq="D", calendar="noleap")

get_hist_temp_for_lat = lambda val: 273.15 - (val * np.cos(2 * np.pi * historical_time.dayofyear / 365) + 2 * np.random.random_sample((historical_time.size,)) + 273.15 + .1 * (historical_time - historical_time[0]).days / 365)
get_rand = lambda: np.random.rand() if np.random.rand() > .5 else  -np.random.rand()

[3]:

latitudes = np.arange(23,27,1)
some_data = [get_hist_temp_for_lat(val) for val in latitudes]
data = np.array([some_data, np.array(some_data)+1])

Create dummy data

[4]:

attrs = {"units": "°C"}

[5]:

obsh = xr.DataArray(
    data,
    dims=("lon", "lat", "time"),
    coords={"time": historical_time, "lat": latitudes, "lon": [0,1]},
    attrs=attrs
).transpose("time","lat","lon").to_dataset(name="tas")#.to_netcdf("observations.nc", )

simh = xr.DataArray(
    data-2,
    dims=("lon", "lat", "time"),
    coords={"time": historical_time, "lat": latitudes, "lon": [0,1]},
    attrs=attrs
).transpose("time","lat","lon").to_dataset(name="tas")#.to_netcdf("control.nc", )

simp = xr.DataArray(
    data-1,
    dims=("lon", "lat", "time"),
    coords={"time": future_time, "lat": latitudes, "lon": [0,1]},
    attrs=attrs
).transpose("time","lat","lon").to_dataset(name="tas")#.to_netcdf("scenario.nc", )

obsp = xr.DataArray(
    data+1,
    dims=("lon", "lat", "time"),
    coords={"time": historical_time, "lat": latitudes, "lon": [0,1]},
    attrs=attrs
).transpose("time","lat","lon").to_dataset(name="tas")#.to_netcdf("observations_future.nc")

Plot created toy data

[6]:

plt.figure(figsize=(10,5),dpi=216)
simh["tas"].groupby("time.dayofyear").mean(...).plot(label="$T_{sim,h}$")
simp["tas"].groupby("time.dayofyear").mean(...).plot(label="$T_{sim,p}$")
obsh["tas"].groupby("time.dayofyear").mean(...).plot(label="$T_{obs,h}$")
plt.title("Historical modeled and obseved temperatures; and predicted temperatures")
plt.xlim(0,365)
plt.gca().grid(alpha=.3)
plt.legend();
_images/examples_9_0.png

  • Modeled historical temperatures are to warm in comparison to the observed temperatures.

Import module to adjust the data

[7]:

from cmethods import adjust

Apply QDM adjustment

[8]:

# to adjust a 3d dataset
qdm_result = adjust(
    method = "quantile_delta_mapping",
    obs = obsh["tas"],
    simh = simh["tas"],
    simp = simp["tas"],
    n_quantiles = 1000,
    kind = "+", # to calculate the relative rather than the absolute change, "*" can be used instead of "+" (this is prefered when adjusting precipitation)
)

[9]:

qdm_result

[9]:

<xarray.Dataset>
Dimensions:  (lat: 4, lon: 2, time: 10950)
Coordinates:
  * lat      (lat) int64 23 24 25 26
  * lon      (lon) int64 0 1
  * time     (time) object 2001-01-01 00:00:00 ... 2030-12-31 00:00:00
Data variables:
    tas      (time, lat, lon) float64 -23.09 -22.09 -24.46 ... -29.84 -28.84

Visualize QDM result

[10]:

plt.figure(figsize=(10,5),dpi=216)
simh["tas"].groupby("time.dayofyear").mean(...).plot(label="$T_{sim,h}$")
simp["tas"].groupby("time.dayofyear").mean(...).plot(label="$T_{sim,p}$")
obsh["tas"].groupby("time.dayofyear").mean(...).plot(label="$T_{obs,h}$")
qdm_result.groupby("time.dayofyear").mean(...).tas.plot(label="$T^{*QDM}_{sim,p}$")
plt.title("Historical and predicted modeled and adjusted temperatures")
plt.xlim(0,365)
plt.gca().grid(alpha=.3)
plt.legend();
_images/examples_17_0.png

  • After the adjustment, the predicted temperatures got warmer (\(\mu T^{*QDM}_{sim,p} > \mu T_{sim,p}\))

[11]:

ls_result = adjust(
    method="linear_scaling",
    obs=obsh["tas"],
    simh=simh["tas"],
    simp=simp["tas"],
    group="time.month",
    kind="+"
)

[12]:

vs_result = adjust(
    method="variance_scaling",
    obs=obsh["tas"],
    simh=simh["tas"],
    simp=simp["tas"],
    group="time.month",
    kind = "+"
)

[13]:

dm_result = adjust(
    method="delta_method",
    obs=obsh["tas"],
    simh=simh["tas"],
    simp=simp["tas"],
    group="time.month",
    kind="+"
)

[14]:

plt.figure(figsize=(10,5),dpi=216)
simh["tas"].groupby("time.dayofyear").mean(...).plot(label="$T_{sim,h}$")
simp["tas"].groupby("time.dayofyear").mean(...).plot(label="$T_{sim,p}$")
obsh["tas"].groupby("time.dayofyear").mean(...).plot(label="$T_{obs,h}$")
ls_result.groupby("time.dayofyear").mean(...).tas.plot(label="$T^{*LS}_{sim,p}$")
vs_result.groupby("time.dayofyear").mean(...).tas.plot(label="$T^{*VS}_{sim,p}$")
dm_result.groupby("time.dayofyear").mean(...).tas.plot(label="$T^{*DM}_{sim,p}$")
qdm_result.groupby("time.dayofyear").mean(...).tas.plot(label="$T^{*QDM}_{sim,p}$")
plt.title("Historical and predicted modeled and adjusted temperatures")
plt.xlim(0,365)
plt.gca().grid(alpha=.3)
plt.legend();
_images/examples_22_0.png

(… because of dummy data - all adjusted datasets seem to have the same result)

It is also possible to adjust individual time series:

[15]:

ls_result = adjust(
    method="linear_scaling",
    obs=obsh["tas"].sel(lat=23, lon=0, method="nearest"),
    simh=simh["tas"].sel(lat=23, lon=0, method="nearest"),
    simp=simp["tas"].sel(lat=23, lon=0, method="nearest"),
    kind="+",
    group="time.month"
)

[16]:

ls_result

[16]:

<xarray.Dataset>
Dimensions:  (time: 10950)
Coordinates:
  * time     (time) object 2001-01-01 00:00:00 ... 2030-12-31 00:00:00
    lat      int64 23
    lon      int64 0
Data variables:
    tas      (time) float64 -23.09 -23.42 -23.18 -23.04 ... -25.58 -26.33 -25.64

For an individual bias correction, ``None`` can be passed to the ``group`` parameter. This disables the month-dependand scaling and uses the whole teme-series as the basis. This enables to create custom timeframes that can be adjusted/corrected individually.