-
Notifications
You must be signed in to change notification settings - Fork 0
/
EquatorialAtlantic_Impact_Nino.py
111 lines (98 loc) · 5.49 KB
/
EquatorialAtlantic_Impact_Nino.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
import numpy as np
from spy4cast import Dataset, Region, Month
from spy4cast.spy4cast import Preprocess, MCA, Crossvalidation, Validation
## CONFIGURATION
# We will use sea surface temperature both for dataset and predictor, but
# with differnt regions
predictor = Dataset("HadISST_sst-1970_2020.nc", "./datasets").open("sst").slice(
Region(lat0=-30, latf=10,
lon0=-60, lonf=15,
month0=Month.JUN, monthf=Month.AUG,
year0=1970, yearf=2019),
skip=1
)
predictand = Dataset("HadISST_sst-1970_2020.nc", "./datasets").open("sst").slice(
Region(lat0=-30, latf=30,
lon0=-200, lonf=-60,
month0=Month.DEC, monthf=Month.FEB,
# year0 and yearf refer to monthf so the slice will span from DEC 1970 to FEB 2020
year0=1971, yearf=2020),
skip=1
)
# There is a lag of 6 months (from June to December)
## METHODOLOGY
# First step. Preprocess variables: anomaly and reshaping
predictor_preprocessed = Preprocess(predictor, period=8, 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")
# Second step. MCA: expansion coefficients and correlation and regression maps
nm = 3
alpha = 0.05
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)
# Third step. Crossvalidation: skill and hidcast evaluation and products
cross = Crossvalidation(predictor_preprocessed, predictand_preprocessed, nm, alpha)
cross.save("cross_", "./data-EquatorialAtalantic_Impact_Nino/")
# cross = Crossvalidation.load("cross_", "./data-EquatorialAtalantic_Impact_Nino/", dsy=predictor_preprocessed, dsz=predictand_preprocessed)
mca.plot(save_fig=True, name="mca.png", folder="./plots-EquatorialAtalantic_Impact_Nino/", ruy_ticks=[-1, -0.5, 0, 0.5, 1], ruz_ticks=[-1, -0.5, 0, 0.5, 1])
cross.plot(save_fig=True, name="cross.png", 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(-2, 2, 10), z_ticks=np.linspace(-2, 2, 5))
# VALIDATION
# It is important to say that this is **not** a good example of validation becuase it uses two periods
# where the relation is **non stationary**, so the results are not good.
# However, it shows how to apply *Validation* to any dataset with any configuration.
# To apply validation we first preprocess the training data
training_y = Preprocess(Dataset("HadISST_sst-1970_2020.nc", "./datasets").open("sst").slice(
Region(lat0=-30, latf=10,
lon0=-60, lonf=15,
month0=Month.JUN, monthf=Month.AUG,
year0=1970, yearf=2000),
skip=1
), period=8, order=4)
training_z = Preprocess(Dataset("HadISST_sst-1970_2020.nc", "./datasets").open("sst").slice(
Region(lat0=-30, latf=30,
lon0=-200, lonf=-60,
month0=Month.DEC, monthf=Month.FEB,
year0=1971, yearf=2001),
skip=1
))
# Optionally save so that we save computing time for next runs
training_y.save("training_y_", "./data-EquatorialAtalantic_Impact_Nino/")
# training_y = Preprocess.load("training_y_", "./data/-EquatorialAtalantic_Impact_Nino")
training_z.save("training_z_", "./data-EquatorialAtalantic_Impact_Nino/")
# training_z = Preprocess.load("training_z_", "./data/-EquatorialAtalantic_Impact_Nino")
# Train the prediction model
training_mca = MCA(training_y, training_z, nm=6, alpha=0.05)
training_mca.save("training_mca_", "./data-EquatorialAtalantic_Impact_Nino/")
# training_mca = MCA.load("training_mca_", "./data-EquatorialAtalantic_Impact_Nino/", dsy=training_y, dsz=training_z)
# We now validate agains the period from 2001 - 2020
validating_y = Preprocess(Dataset("HadISST_sst-1970_2020.nc", "./datasets").open("sst").slice(
Region(lat0=-30, latf=10,
lon0=-60, lonf=15,
month0=Month.JUN, monthf=Month.AUG,
year0=2010, yearf=2019),
skip=1
), period=8, order=4)
validating_z = Preprocess(Dataset("HadISST_sst-1970_2020.nc", "./datasets").open("sst").slice(
Region(lat0=-30, latf=30,
lon0=-200, lonf=-60,
month0=Month.DEC, monthf=Month.FEB,
year0=2011, yearf=2020),
skip=1
))
# Optionally save so that we save computing time for next runs
validating_y.save("validating_y_", "./data-EquatorialAtalantic_Impact_Nino/")
# validating_y = Preprocess.load("training_y_", "./data/-EquatorialAtalantic_Impact_Nino")
validating_z.save("validating_z_", "./data-EquatorialAtalantic_Impact_Nino/")
# validating_z = Preprocess.load("validating_z_", "./data/-EquatorialAtalantic_Impact_Nino")
validation = Validation(training_mca, validating_y, validating_z)
validation.save("validation_z_", "./data-EquatorialAtalantic_Impact_Nino/")
validation.plot(save_fig=True, folder="./plots-EquatorialAtalantic_Impact_Nino/",
name="validation.png", version='default')
validation.plot_zhat(2016, save_fig=True, folder="./plots-EquatorialAtalantic_Impact_Nino/",
name="zhat_2016_validation.png", z_levels=np.linspace(-4.1, 4.1, 10))