You can run this notebook in a live session or view it on Github.
Temporal and Spatial Smoothing#
This demo demonstrates climpred
’s capabilities to postprocess decadal prediction output before skill verification.
Here, we showcase a set of methods to smooth out noise in the spatial and temporal domain.
# linting
%load_ext nb_black
%load_ext lab_black
%matplotlib inline
from climpred import PerfectModelEnsemble, HindcastEnsemble
from climpred.tutorial import load_dataset
import matplotlib.pylab as plt
import xarray as xr
# Sea surface temperature
v = "tos"
initialized3d = load_dataset("MPI-PM-DP-3D")[v]
control3d = load_dataset("MPI-control-3D")[v]
climpred
requires that lead
dimension has an attribute called units
indicating what time units the lead
is assocated with. Options are: years
, seasons
, months
, weeks
, pentads
, days
, hours
, minutes
, seconds
. For the this data, the lead
units
are years
.
initialized3d["lead"].attrs = {"units": "years"}
Temporal smoothing#
In order to reduce temporal noise, decadal predictions are recommended to take multi-year averages [Goddard et al., 2013].
pm = PerfectModelEnsemble(initialized3d).add_control(control3d)
/Users/aaron.spring/Coding/climpred/climpred/utils.py:191: UserWarning: Assuming annual resolution starting Jan 1st due to numeric inits. Please change ``init`` to a datetime if it is another resolution. We recommend using xr.CFTimeIndex as ``init``, see https://climpred.readthedocs.io/en/stable/setting-up-data.html.
warnings.warn(
/Users/aaron.spring/Coding/climpred/climpred/utils.py:191: UserWarning: Assuming annual resolution starting Jan 1st due to numeric inits. Please change ``init`` to a datetime if it is another resolution. We recommend using xr.CFTimeIndex as ``init``, see https://climpred.readthedocs.io/en/stable/setting-up-data.html.
warnings.warn(
PerfectModelEnsemble.smooth()
({"time":x})
aggregates over x
timesteps in time dimensions lead
and time
. Here it does not matter whether you specify lead
or time
, temporal smoothing is applied to both time dimensions. Note that the time dimension labels are not changed by this temporal smoothing immediately.
pm_tsmoothed = pm.smooth({"lead": 3})
print("initialized", pm_tsmoothed.get_initialized().coords, "\n")
print("control", pm_tsmoothed.get_control().coords)
initialized Coordinates:
lon (y, x) float64 -47.25 -47.69 -48.12 -48.54 ... 131.3 132.5 133.8
lat (y, x) float64 76.36 76.3 76.24 76.17 ... -77.25 -77.39 -77.54
* lead (lead) int64 1 2 3
* init (init) object 3014-01-01 00:00:00 ... 3237-01-01 00:00:00
* member (member) int64 1 2 3 4
valid_time (lead, init) object 3015-01-01 00:00:00 ... 3240-01-01 00:00:00
control Coordinates:
lon (y, x) float64 -47.25 -47.69 -48.12 -48.54 ... 131.3 132.5 133.8
lat (y, x) float64 76.36 76.3 76.24 76.17 ... -77.25 -77.39 -77.54
* time (time) object 3000-01-01 00:00:00 ... 3047-01-01 00:00:00
But after calling PerfectModelEnsemble.verify()
, the correct time aggregation label is automatically set.
pm_tsmoothed.verify(metric="rmse", comparison="m2e", dim=["init", "member"])[v].plot(
col="lead", vmin=0, vmax=0.7, yincrease=False, x="x"
)
<xarray.plot.facetgrid.FacetGrid at 0x14b302df0>

Compare to the prediction skill without smoothing (with rmse
smaller is better):
pm.verify(metric="rmse", comparison="m2e", dim=["init", "member"])[v].plot(
col="lead", vmin=0, vmax=0.7, yincrease=False, x="x"
)
<xarray.plot.facetgrid.FacetGrid at 0x14b4f4970>

You can plot time-aggregated skill timeseries by specifying x=lead_center
while plotting. Lead centers are placed in the center of the lead-aggregation boundaries. Note how the correlation of forecasts is increasing for longer time-aggregations because noise is smoothed out.
pm = PerfectModelEnsemble(load_dataset("MPI-PM-DP-1D"))
pm = pm.add_control(load_dataset("MPI-control-1D"))
pm_NA = pm.sel(area="North_Atlantic_SPG", period="DJF")[['tos']]
/Users/aaron.spring/Coding/climpred/climpred/utils.py:191: UserWarning: Assuming annual resolution starting Jan 1st due to numeric inits. Please change ``init`` to a datetime if it is another resolution. We recommend using xr.CFTimeIndex as ``init``, see https://climpred.readthedocs.io/en/stable/setting-up-data.html.
warnings.warn(
/Users/aaron.spring/Coding/climpred/climpred/utils.py:191: UserWarning: Assuming annual resolution starting Jan 1st due to numeric inits. Please change ``init`` to a datetime if it is another resolution. We recommend using xr.CFTimeIndex as ``init``, see https://climpred.readthedocs.io/en/stable/setting-up-data.html.
warnings.warn(
for time_agg in [1, 2, 3, 4]:
pm_NA.smooth({"lead": time_agg}).verify(
metric="acc", comparison="m2e", dim=["init", "member"]
).tos.plot(x="lead_center", marker="o", label=f"{time_agg}-years smoothed")
plt.legend()
<matplotlib.legend.Legend at 0x14b28d940>

Now, we showcase these smoothing features for HindcastEnsemble.smooth()
.
v = "SST"
initialized = load_dataset("CESM-DP-SST-3D")[v]
reconstruction = load_dataset("FOSI-SST-3D")[v]
# Move reconstruction into proper anomaly space
reconstruction = reconstruction - reconstruction.sel(time=slice(1964, 2014)).mean(
"time"
)
hindcast = HindcastEnsemble(initialized)
hindcast = hindcast.add_observations(reconstruction)
/Users/aaron.spring/Coding/climpred/climpred/utils.py:191: UserWarning: Assuming annual resolution starting Jan 1st due to numeric inits. Please change ``init`` to a datetime if it is another resolution. We recommend using xr.CFTimeIndex as ``init``, see https://climpred.readthedocs.io/en/stable/setting-up-data.html.
warnings.warn(
/Users/aaron.spring/Coding/climpred/climpred/utils.py:191: UserWarning: Assuming annual resolution starting Jan 1st due to numeric inits. Please change ``init`` to a datetime if it is another resolution. We recommend using xr.CFTimeIndex as ``init``, see https://climpred.readthedocs.io/en/stable/setting-up-data.html.
warnings.warn(
hindcast.smooth({"lead": 5}).verify(
metric="rmse", comparison="e2r", dim="init", alignment="same_verif"
)[v].plot(col="lead", robust=True)
<xarray.plot.facetgrid.FacetGrid at 0x14b1e1eb0>

Spatial smoothing#
In order to reduce spatial noise, global decadal predictions are recommended to get regridded to a 5 degree longitude x 5 degree latitude grid as recommended [Goddard et al., 2013].
v = "tos"
pm_ssmoothed = pm.smooth({"lon": 5, "lat": 5})
pm_ssmoothed.get_initialized().coords
/Users/aaron.spring/anaconda3/envs/climpred-dev/lib/python3.9/site-packages/xarray/core/dataarray.py:784: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
return key in self.data
/Users/aaron.spring/anaconda3/envs/climpred-dev/lib/python3.9/site-packages/xesmf/backend.py:53: UserWarning: Latitude is outside of [-90, 90]
warnings.warn('Latitude is outside of [-90, 90]')
/Users/aaron.spring/anaconda3/envs/climpred-dev/lib/python3.9/site-packages/xarray/core/dataarray.py:784: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
return key in self.data
/Users/aaron.spring/anaconda3/envs/climpred-dev/lib/python3.9/site-packages/xesmf/backend.py:53: UserWarning: Latitude is outside of [-90, 90]
warnings.warn('Latitude is outside of [-90, 90]')
Coordinates:
* lead (lead) int64 1 2 3 4 5
* init (init) object 3014-01-01 00:00:00 ... 3237-01-01 00:00:00
* member (member) int64 1 2 3 4
valid_time (lead, init) object 3015-01-01 00:00:00 ... 3242-01-01 00:00:00
* lon (lon) float64 -180.0 -175.0 -170.0 -165.0 ... 170.0 175.0 180.0
* lat (lat) float64 -83.97 -78.97 -73.97 -68.97 ... 81.03 86.03 91.03
pm_ssmoothed.verify(metric="rmse", comparison="m2e", dim=["init", "member"])[v].T.plot(
col="lead", robust=True, yincrease=True,
)
<xarray.plot.facetgrid.FacetGrid at 0x149f67670>

Spatial smoothing guesses the names corresponding to lon
and lat
in the coordinates of the PredictionEnsemble
underlying datasets.
hindcast.get_initialized().coords
Coordinates:
TLAT (nlat, nlon) float64 -9.75 -9.75 -9.75 ... -0.1336 -0.1336
TLONG (nlat, nlon) float64 250.8 251.9 253.1 ... 276.7 277.8 278.9
* init (init) object 1954-01-01 00:00:00 ... 2017-01-01 00:00:00
* lead (lead) int32 1 2 3 4 5 6 7 8 9 10
TAREA (nlat, nlon) float64 3.661e+13 3.661e+13 ... 3.714e+13 3.714e+13
valid_time (lead, init) object 1955-01-01 00:00:00 ... 2027-01-01 00:00:00
hindcast.smooth({"lon": 1, "lat": 1}).get_initialized().coords
/Users/aaron.spring/anaconda3/envs/climpred-dev/lib/python3.9/site-packages/xarray/core/dataarray.py:784: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
return key in self.data
/Users/aaron.spring/anaconda3/envs/climpred-dev/lib/python3.9/site-packages/xarray/core/dataarray.py:784: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
return key in self.data
Coordinates:
* init (init) object 1954-01-01 00:00:00 ... 2017-01-01 00:00:00
* lead (lead) int32 1 2 3 4 5 6 7 8 9 10
valid_time (lead, init) object 1955-01-01 00:00:00 ... 2027-01-01 00:00:00
* lon (lon) float64 250.8 251.8 252.8 253.8 ... 277.8 278.8 279.8
* lat (lat) float64 -9.75 -8.75 -7.75 -6.75 ... -1.75 -0.7503 0.2497
HindcastEnsemble.smooth()
("goddard2013")
creates 4-year means and 5x5 degree regridding as suggested in [Goddard et al., 2013].
pm.smooth("goddard2013").verify(
metric="acc", comparison="m2e", dim=["init", "member"]
).coords
/Users/aaron.spring/anaconda3/envs/climpred-dev/lib/python3.9/site-packages/xarray/core/dataarray.py:784: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
return key in self.data
/Users/aaron.spring/anaconda3/envs/climpred-dev/lib/python3.9/site-packages/xesmf/backend.py:53: UserWarning: Latitude is outside of [-90, 90]
warnings.warn('Latitude is outside of [-90, 90]')
/Users/aaron.spring/anaconda3/envs/climpred-dev/lib/python3.9/site-packages/xarray/core/dataarray.py:784: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
return key in self.data
/Users/aaron.spring/anaconda3/envs/climpred-dev/lib/python3.9/site-packages/xesmf/backend.py:53: UserWarning: Latitude is outside of [-90, 90]
warnings.warn('Latitude is outside of [-90, 90]')
Coordinates:
* lead (lead) <U3 '1-4' '2-5'
* lon (lon) float64 -180.0 -175.0 -170.0 -165.0 ... 170.0 175.0 180.0
* lat (lat) float64 -83.97 -78.97 -73.97 -68.97 ... 81.03 86.03 91.03
lead_center (lead) float64 2.5 3.5
References#
L. Goddard, A. Kumar, A. Solomon, D. Smith, G. Boer, P. Gonzalez, V. Kharin, W. Merryfield, C. Deser, S. J. Mason, B. P. Kirtman, R. Msadek, R. Sutton, E. Hawkins, T. Fricker, G. Hegerl, C. a. T. Ferro, D. B. Stephenson, G. A. Meehl, T. Stockdale, R. Burgman, A. M. Greene, Y. Kushnir, M. Newman, J. Carton, I. Fukumori, and T. Delworth. A verification framework for interannual-to-decadal predictions experiments. Climate Dynamics, 40(1-2):245–272, January 2013. doi:10/f4jjvf.