Jupyter Notebook Version for calculating Equatorial Impact of the Niño¶
⚠️ The data used in this notebook is from HadISST and can be downloaded using this download link
After download, decompress the file and save it in ./datasets under the name “HadISST_sst.nc”
[1]:
from spy4cast import Dataset, Region, Month
from spy4cast.spy4cast import Preprocess, MCA, Crossvalidation
import numpy as np
Configuration¶
[2]:
# We will use sea surface temperature both for dataset and predictor, but
# with differnt regions
predictor = Dataset("HadISST_sst.nc", "./datasets").open("sst").slice(
Region(lat0=-30, latf=10,
lon0=-60, lonf=15,
month0=Month.JUN, monthf=Month.AUG,
year0=1890, yearf=2019),
)
predictand = Dataset("HadISST_sst.nc", "./datasets").open("sst").slice(
Region(lat0=-30, latf=30,
lon0=-200, lonf=-60,
month0=Month.DEC, monthf=Month.FEB,
# year0, yearf refer to monthf -> will span from DEC 1970 to FEB 2020
year0=1891, yearf=2020),
)
# There is a lag of 6 months (from June to December)
Methodology¶
Preprocessing¶
[3]:
# First step. Preprocess variables: anomaly and reshaping
predictor_preprocessed = Preprocess(predictor, period=7, order=4)
predictor_preprocessed.save("y_", "./data-EquatorialAtalantic_Impact_Nino/")
# Save matrices as .npy for fast loading. To load use:
# predictor_preprocessed = Preprocess.load("y_",
# "./data-EquatorialAtalantic_Impact_Nino/")
predictand_preprocessed = Preprocess(predictand)
predictand_preprocessed.save("z_", "./data-EquatorialAtalantic_Impact_Nino/")
# predictand_preprocessed = Preprocess.load("z_",
# "./data/-EquatorialAtalantic_Impact_Nino")
[INFO] Preprocessing data for variable sst took: 0.431 seconds
[INFO] Saving Preprocess data in `./data-EquatorialAtalantic_Impact_Nino/y_*.npy`
[INFO] Preprocessing data for variable sst took: 1.643 seconds
[INFO] Saving Preprocess data in `./data-EquatorialAtalantic_Impact_Nino/z_*.npy`
MCA¶
[6]:
# Second step. MCA: expansion coefficients and correlation and regression maps
nm = 3
alpha = 0.1
mca = MCA(predictor_preprocessed, predictand_preprocessed, nm, alpha)
mca.save("mca_", "./data-EquatorialAtalantic_Impact_Nino/")
# mca = MCA.load("mca_", "./data-EquatorialAtalantic_Impact_Nino/",
# dsy=predictor_preprocessed, dsz=predictand_preprocessed)
[INFO] Applying MCA
Shapes: Z(8400, 130)
Y(3000, 130)
Regions: Z DJF (30.00ºS, 30.00ºN - 160.00ºE, 60.00ºW)
Y JJA (30.00ºS, 10.00ºN - 60.00ºW, 15.00ºE)
Took: 9.537 seconds
[INFO] Saving MCA data in `./data-EquatorialAtalantic_Impact_Nino/mca_*.npy`
Plot MCA¶
[7]:
mca.plot(save_fig=True, name="mca.png",
variable="s",
folder="./plots-EquatorialAtalantic_Impact_Nino/",
y_levels=np.arange(-.51, .51, .1),
z_levels=np.arange(-.51, .51, .1),
width_ratios=[1, 1, 1])
[INFO] Saving plot with path ./plots-EquatorialAtalantic_Impact_Nino/mca.png
[7]:
((<Figure size 1200x800 with 11 Axes>,),
(<AxesSubplot:title={'center':'Us Vs mode 1'}>,
<AxesSubplot:title={'center':'Us Vs mode 2'}>,
<AxesSubplot:title={'center':'Us Vs mode 3'}>,
<GeoAxesSubplot:title={'center':'SUY mode 1. SCF=90.3%'}>,
<GeoAxesSubplot:title={'center':'SUY mode 2. SCF=4.6%'}>,
<GeoAxesSubplot:title={'center':'SUY mode 3. SCF=2.3%'}>,
<GeoAxesSubplot:title={'center':'SUZ mode 1. SCF=90.3%'}>,
<GeoAxesSubplot:title={'center':'SUZ mode 2. SCF=4.6%'}>,
<GeoAxesSubplot:title={'center':'SUZ mode 3. SCF=2.3%'}>))
[8]:
# Plot on a larger map using global regression
map_y = Preprocess(Dataset("HadISST_sst.nc", "./datasets").open("sst").slice(
Region(lat0=-60, latf=60, lon0=-180, lonf=50, month0=Month.JUN, monthf=Month.AUG,
year0=1890, yearf=2019)))
map_z = Preprocess(Dataset("HadISST_sst.nc", "./datasets").open("sst").slice(
Region(lat0=-60, latf=60, lon0=-300, lonf=-40, month0=Month.DEC, monthf=Month.FEB,
year0=1891
, yearf=2020)))
[INFO] Preprocessing data for variable sst took: 0.302 seconds
[INFO] Preprocessing data for variable sst took: 4.325 seconds
[9]:
mca.plot(save_fig=True,
name="mca_global.png",
folder="./plots-EquatorialAtalantic_Impact_Nino/",
map_y=map_y, map_z=map_z, figsize=(11, 6),
variable="s",
width_ratios=[1, 1, 1],
y_levels=np.arange(-.5, .51, .1),
y_ticks=np.arange(-.5, .51, .25),
z_levels=np.arange(-.51, .51, .1),
z_ticks=np.arange(-.5, .51, .25), )
[INFO] Saving plot with path ./plots-EquatorialAtalantic_Impact_Nino/mca_global.png
[9]:
((<Figure size 1100x600 with 11 Axes>,),
(<AxesSubplot:title={'center':'Us Vs mode 1'}>,
<AxesSubplot:title={'center':'Us Vs mode 2'}>,
<AxesSubplot:title={'center':'Us Vs mode 3'}>,
<GeoAxesSubplot:title={'center':'SUY mode 1. SCF=90.3%'}>,
<GeoAxesSubplot:title={'center':'SUY mode 2. SCF=4.6%'}>,
<GeoAxesSubplot:title={'center':'SUY mode 3. SCF=2.3%'}>,
<GeoAxesSubplot:title={'center':'SUZ mode 1. SCF=90.3%'}>,
<GeoAxesSubplot:title={'center':'SUZ mode 2. SCF=4.6%'}>,
<GeoAxesSubplot:title={'center':'SUZ mode 3. SCF=2.3%'}>))
Only first mode will be used as it offers 50% of the predictability
Crossvalidation¶
[22]:
# Third step. Crossvalidation: skill and hidcast evaluation and products
# num_svdvals=1 to speed up calculation (we don't need precise scf)
# cross = Crossvalidation(predictor_preprocessed, predictand_preprocessed, nm=1, alpha=alpha, num_svdvals=1)
# cross.save("cross_", "./data-EquatorialAtalantic_Impact_Nino/")
cross = Crossvalidation.load("cross_", "./data-EquatorialAtalantic_Impact_Nino/",
dsy=predictor_preprocessed, dsz=predictand_preprocessed)
[INFO] Loading Crossvalidation data from `./data-EquatorialAtalantic_Impact_Nino/cross_*` took 0.028 seconds
Plotting crossvalidation¶
[33]:
cross.plot(
save_fig=True, name="cross.png", plot_type="pcolor",
figsize=(12, 8),
nm=1,
cmap="plasma",
map_levels=np.arange(0, 0.41, .1),
map_ticks=np.arange(0, 0.41, .1),
folder="./plots-EquatorialAtalantic_Impact_Nino/",
)
# cross.plot_zhat(1998, figsize=(12, 10), save_fig=True, name="zhat_1998.png",
# folder="./plots-EquatorialAtalantic_Impact_Nino/",
# z_levels=np.linspace(0, 2, 10), z_ticks=np.linspace(-2, 2, 5))
[INFO] Saving plot with path ./plots-EquatorialAtalantic_Impact_Nino/cross.png
[33]:
((<Figure size 1200x800 with 6 Axes>,),
(<GeoAxesSubplot:title={'center':'ACC map'}>,
<AxesSubplot:title={'center':'ACC time series'}>,
<GeoAxesSubplot:title={'center':'RMSE map'}>,
<AxesSubplot:title={'center':'RMSE time series'}>))
Analyse predictability¶
[6]:
from scipy import stats
[7]:
win = 20
nt = mca.Us.shape[1]
r_uv = np.empty(nt - win + 1)
p_uv = np.empty(nt - win + 1)
for i in range(nt - win + 1):
r_uv[i], p_uv[i] = stats.pearsonr(mca.Us[0, i:i+win], mca.Vs[0, i:i+win])
[8]:
import matplotlib.pyplot as plt
[12]:
fig = plt.figure(figsize=(7, 4))
ax = fig.add_subplot()
ytime = mca.dsy.time
ax.plot(ytime[:nt-win+1], r_uv, ".", color="black", label=r"Correlation $(U_s,\,V_s)$")
ax.scatter(mca.dsy.time[:nt-win+1][p_uv < alpha], r_uv[p_uv < alpha], color="red", alpha=0.5, label="$p_{val} < 0.1$")
ax.grid()
ax.xaxis.set_ticklabels([rf"${int(year)} \rightarrow {int(year+win-1)}$" for year in ax.get_xticks()], rotation=15)
ax.set_xlabel("$Y$ time")
ax.set_ylabel("Correlation")
fig.tight_layout()
ax.legend()
fig.savefig("./plots-EquatorialAtalantic_Impact_Nino/correlation_u_v.png")
## Validation
Let’s use the period where the correlation is high and significant (around 1970 to 1995) to train the model and validate in the years after
[34]:
from spy4cast.spy4cast import Validation
[35]:
# To apply validation we first preprocess the training data
training_y = Preprocess(Dataset("HadISST_sst.nc", "./datasets").open("sst").slice(
Region(lat0=-30, latf=10, lon0=-60, lonf=15, month0=Month.JUN, monthf=Month.AUG,
year0=1970, yearf=1995)), period=7, order=4)
training_z = Preprocess(Dataset("HadISST_sst.nc", "./datasets").open("sst").slice(
Region(lat0=-30, latf=30, lon0=-200, lonf=-60, month0=Month.DEC, monthf=Month.FEB,
year0=1971, yearf=1996)))
[INFO] Preprocessing data for variable sst took: 0.185 seconds
[INFO] Preprocessing data for variable sst took: 0.251 seconds
[36]:
training_mca = MCA(training_y, training_z, nm=1, alpha=alpha)
training_mca.plot()
[INFO] Applying MCA
Shapes: Z(8400, 26)
Y(3000, 26)
Regions: Z DJF (30.00ºS, 30.00ºN - 160.00ºE, 60.00ºW)
Y JJA (30.00ºS, 10.00ºN - 60.00ºW, 15.00ºE)
Took: 7.445 seconds
[36]:
((<Figure size 1700x566.667 with 5 Axes>,),
(<AxesSubplot:title={'center':'Us Vs mode 1'}>,
<GeoAxesSubplot:title={'center':'RUY mode 1. SCF=47.6%'}>,
<GeoAxesSubplot:title={'center':'RUZ mode 1. SCF=47.6%'}>))
[ ]:
training_cross = Crossvalidation(training_y, training_z, nm=1, alpha=alpha, num_svdvals=1)
training_cross.save("cross_training_", "./data-EquatorialAtalantic_Impact_Nino/")
# training_cross = Crossvalidation.load("cross_training_", "./data-EquatorialAtalantic_Impact_Nino/", dsy=training_y, dsz=training_z)
[58]:
training_cross.plot(
save_fig=True,
name="cross_training.png",
plot_type="pcolor",
cmap="plasma",
map_levels=np.arange(0, 1.1, .1),
map_ticks=np.arange(0, 1.1, .2),
folder="./plots-EquatorialAtalantic_Impact_Nino/")
[INFO] Saving plot with path ./plots-EquatorialAtalantic_Impact_Nino/cross_high.png
[58]:
((<Figure size 1600x800 with 6 Axes>,),
(<GeoAxesSubplot:title={'center':'ACC map'}>,
<AxesSubplot:title={'center':'ACC time series'}>,
<GeoAxesSubplot:title={'center':'RMSE map'}>,
<AxesSubplot:title={'center':'RMSE time series'}>))
[9]:
from spy4cast.meteo import Anom
[42]:
validating_y1 = Preprocess(Dataset("HadISST_sst.nc", "./datasets").open("sst").slice(
Region(lat0=-30, latf=10, lon0=-60, lonf=15, month0=Month.JUN, monthf=Month.AUG,
year0=1996, yearf=2019)))
validating_z1 = Preprocess(Dataset("HadISST_sst.nc", "./datasets").open("sst").slice(
Region(lat0=-30, latf=30, lon0=-200, lonf=-60, month0=Month.DEC, monthf=Month.FEB,
year0=1997, yearf=2020)))
validation1 = Validation(training_mca, validating_y1, validating_z1)
validation1.plot(
save_fig=True,
figsize=(12, 8),
folder="./plots-EquatorialAtalantic_Impact_Nino/",
name="validation19962019.png",
cmap="magma",
version='default',
map_levels=np.arange(0, 1.01, .1),
map_ticks=np.arange(0, 1.01, .2),
)
[INFO] Preprocessing data for variable sst took: 0.013 seconds
[INFO] Preprocessing data for variable sst took: 0.245 seconds
[INFO] Applying Validation
Training
----------
Shapes: Z(8400, 26)
Y(3000, 26)
Regions: Z DJF (30.00ºS, 30.00ºN - 160.00ºE, 60.00ºW)
Y JJA (30.00ºS, 10.00ºN - 60.00ºW, 15.00ºE)
Validation
----------
Shapes: Z(8400, 24)
Y(3000, 24)
Regions: Z DJF (30.00ºS, 30.00ºN - 160.00ºE, 60.00ºW)
Y JJA (30.00ºS, 10.00ºN - 60.00ºW, 15.00ºE)
Took: 0.590 seconds
[INFO] Saving plot with path ./plots-EquatorialAtalantic_Impact_Nino/validation19962019.png
[42]:
((<Figure size 1200x800 with 6 Axes>,),
(<GeoAxesSubplot:title={'center':'Correlation in space between z and zhat'}>,
<AxesSubplot:title={'center':'Correlation in space between z and zhat'}>,
<GeoAxesSubplot:title={'center':'RMSE'}>,
<AxesSubplot:title={'center':'RMSE time series'}>))
[46]:
validation1.plot_zhat(
[1997, 1998, 2009],
save_fig=True,
figsize=(11, 7),
name="validated_zhat.png",
folder="./plots-EquatorialAtalantic_Impact_Nino/",
y_levels=np.arange(-1.1, 1.11, .1),
z_levels=np.arange(-1.1, 1.11, .1))
[INFO] Saving plot with path ./plots-EquatorialAtalantic_Impact_Nino/validated_zhat.png
[46]:
(<Figure size 1100x700 with 15 Axes>,
(<GeoAxesSubplot:title={'center':'$Y(2008)$'}>,
<GeoAxesSubplot:title={'center':'$\\hat Z(2009)$'}>,
<GeoAxesSubplot:title={'center':'Z(2009)'}>))
[ ]:
[61]:
atlantic_nino_anomaly = Anom(
Dataset("HadISST_sst.nc", "./datasets").open("sst").slice(Region(-3, 3, -20, 0, Month.JUN, Month.AUG, 1970, 2019)), "ts")
fig = plt.figure(figsize=(10, 5))
ax = fig.add_subplot()
ax.plot(atlantic_nino_anomaly.time.values, atlantic_nino_anomaly.data.values, label="Atlantic Niño index", color="black", alpha=.5)
ax.legend()
ax.grid()
[ ]:
[ ]: