-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathcreate_ly_series.py
118 lines (96 loc) · 5.58 KB
/
create_ly_series.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
111
112
113
114
115
116
117
118
#script to plot all wanted plots
from leadyear import ly_series
import config as cfg
from leadyear import calculate_leadyear
from preprocessing import get_variable
import preprocessing as pp
from decorrelation_time import decorrelation_time
from plots import correlation_plot, plot_variable_mask
from plots import bias_plot
from residuals import residual
import matplotlib
import h5py as h5
import numpy as np
import matplotlib.pyplot as plt
matplotlib.use('Agg')
cfg.set_args()
#plot correlation for a specific lead year
#lys = calculate_leadyear(cfg.start_year, cfg.end_year, lead_year='2 9')
#lys.plot()
#plot leadyear correlation for all lead years
#ly = ly_series(cfg.start_year, cfg.end_year)
#ly.ly_series()
#plot decorrelation time for HadISSTs
#define start and end years, threshold for decorrelation mask
start_year = 1960
end_year = 2015
#normal monthly hadissts
HadIsst = get_variable(cfg.observation_path, name='HadIsst', start_year=start_year, end_year=end_year)
HadIsst = HadIsst.__getitem__()
#save continent mask
continent_mask = np.where(HadIsst[0]==0.1, np.nan, 1)
f = h5.File('tmp/continent_mask.hdf5', 'w')
f.create_dataset('contintent_mask', continent_mask.shape, dtype = 'float32',data = continent_mask)
f.close()
#historical model run, annual mean
hist = get_variable(cfg.historical_path, name=cfg.hist_name, ensemble_members=cfg.ensemble_member, start_year=start_year,
end_year=end_year, start_month='01', variable='tos', ensemble=True, mean='annual', mode='hist')
hist = hist.__getitem__()
#annual hadisst timeseries
HadIsst_annual = get_variable(path = cfg.observation_path, name = 'HadI_annual', start_year = start_year, end_year = end_year, mean='annual')
HadIsst_annual = HadIsst_annual.__getitem__()
#create residual observation timeseries
residual_dataset = residual(lead_year=1, start_year=start_year)
HadIsst_res = residual_dataset.obs_res(HadIsst_annual, hist) * continent_mask
#plot residual and normal annual decorrelation times
decor = decorrelation_time(HadIsst_annual, del_t=1, threshold=2, name='HadIsst_annual')
dc, mask_decor_1 = decor.__getitem__()
decor.plot()
decor = decorrelation_time(HadIsst_annual, del_t=4, threshold=4, name='HadIsst_annual')
dc, mask_decor_4 = decor.__getitem__()
decor.plot()
#decor = decorrelation_time(HadIsst_annual, del_t=8, threshold=threshold, name='HadIsst_annual')
#dc, mask_decor_8 = decor.__getitem__()
#decor.plot()
decor = decorrelation_time(HadIsst_res, del_t=1, threshold=2, name='HadIsst_annual_res')
dc, mask_decor_1_res = decor.__getitem__()
decor.plot()
decor = decorrelation_time(HadIsst_res, del_t=4, threshold=4, name='HadIsst_annual_res')
dc, mask_decor_4_res = decor.__getitem__()
decor.plot()
#decor = decorrelation_time(HadIsst_res, del_t=8, threshold=threshold, name='HadIsst_annual_res')
#dc, mask_decor_8_res = decor.__getitem__()
#decor.plot()
#plot ssh bias and correlation
Aviso_ssh = get_variable(path=cfg.tmp_path + 'Aviso_Ssh_full/Aviso_Ssh_full_2000_2018.nc', name='Aviso_ssh', start_year=2000, end_year=2017, variable='var', mean='monthly')
Aviso_ssh = Aviso_ssh.__getitem__()
Assi_ssh = get_variable(path=cfg.assi_path, name='_dcppA-assim_r', ensemble_members=cfg.ensemble_member, mod_year=cfg.ssh_mod, start_year=2000, end_year=2017, start_month='01', start_year_file=1950, end_year_file=2017, variable='zos', ensemble=True, time_edit=True, mean='monthly')
Assi_ssh = Assi_ssh.__getitem__()
mask_ssh_8 = correlation_plot(Aviso_ssh, Assi_ssh, del_t=8, name_1='Aviso_ssh', name_2='Assimilation_ssh')
mask_ssh_4 = correlation_plot(Aviso_ssh, Assi_ssh, del_t=4, name_1='Aviso_ssh', name_2='Assimilation_ssh')
mask_ssh_1 = correlation_plot(Aviso_ssh, Assi_ssh, del_t=1, name_1='Aviso_ssh', name_2='Assimilation_ssh')
bias_plot(Aviso_ssh, Assi_ssh, name_1='Aviso_ssh', name_2='Assimilation_ssh')
#plot correlation between ocean heat content and ssts
IAP_Ohc = get_variable(path = cfg.ohc_path, name='IAP_Ohc', start_year=196100, end_year=201600, variable='heatcontent', time_edit=False)
IAP_Ohc = IAP_Ohc.__getitem__()
mask_ohc_8 = correlation_plot(HadIsst, IAP_Ohc, del_t=8, name_1='Residual_HadIsst', name_2='IAP_Ohc')
mask_ohc_4 = correlation_plot(HadIsst, IAP_Ohc, del_t=4, name_1='Residual_HadIsst', name_2='IAP_Ohc')
mask_ohc_1 = correlation_plot(HadIsst, IAP_Ohc, del_t=1, name_1='Residual_HadIsst', name_2='IAP_Ohc')
#final_mask_8 = mask_decor_8 * mask_ohc_8 * mask_ssh_8
final_mask_4 = mask_decor_4 * mask_ohc_4 * mask_ssh_4
final_mask_1 = mask_decor_1 * mask_ohc_1 * mask_ssh_1
#final_mask_8_res = mask_decor_8_res * mask_ohc_8 * mask_ssh_8
final_mask_4_res = mask_decor_4_res * mask_ohc_4 * mask_ssh_4
final_mask_1_res = mask_decor_1_res * mask_ohc_1 * mask_ssh_1
f4 = h5.File(cfg.tmp_path + 'correlation/correlation' + str(cfg.start_year) + '_' + str(cfg.end_year) + '_' + str(12) + '.hdf5', 'r')
f8 = h5.File(cfg.tmp_path + 'correlation/correlation' + str(cfg.start_year) + '_' + str(cfg.end_year) + '_' + str(20) + '.hdf5', 'r')
f1 = h5.File(cfg.tmp_path + 'correlation/correlation' + str(cfg.start_year) + '_' + str(cfg.end_year) + '_' + str(1) + '.hdf5', 'r')
res_hind_4 = f4.get('res_hind_corr')
res_hind_8 = f8.get('res_hind_corr')
res_hind_1 = f1.get('res_hind_corr')
plot_variable_mask(res_hind_1, final_mask_1_res, 'Residual_res_hindcast_leadyear_1')
plot_variable_mask(res_hind_4, final_mask_4_res, 'Residual_res_hindcast_leadyear_2-5')
#plot_variable_mask(res_hind_8, final_mask_8_res, 'Residual_res_hindcast_leadyear_2-9')
plot_variable_mask(res_hind_1, final_mask_1, 'Residual_hindcast_leadyear_1')
plot_variable_mask(res_hind_4, final_mask_4, 'Residual_hindcast_leadyear_2-5')
#plot_variable_mask(res_hind_8, final_mask_8, 'Residual_hindcast_leadyear_2-9')