-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathfigura_comparacion_PC1-PDO.py
113 lines (92 loc) · 3.67 KB
/
figura_comparacion_PC1-PDO.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
"""
En esta rutina ploteo las series temporales de PC1 de SSTa,
y PDO en el SWA y SEOP
Dani Risaro
Enero 2020
"""
import warnings
import os
warnings.filterwarnings('ignore')
import numpy as np
import pandas as pd
import xarray as xr
import matplotlib.pyplot as plt
from eofs.standard import Eof
warnings.filterwarnings('ignore')
os.chdir('/home/daniu/Documentos/rutinas/')
from edof_mpc_py import meanun
from indices_climaticos_mensuales import carga_indices
time_tot = pd.date_range('1982-01-01','2017-12-31', freq='MS')
# load index data
ind = ['AAO','SAM','ENSO 3.4','SOI','PDO','IPO','NAO']
ind_df = pd.DataFrame(index=time_tot, columns=ind)
AAO, SAM, NINO34, SOI, PDO, IPO, NAO = carga_indices()
ind_df['AAO'] = AAO.loc['1982':'2017']
ind_df['SAM'] = SAM.loc['1982':'2017']
ind_df['ENSO 3.4'] = NINO34.loc['1982':'2017']
ind_df['SOI'] = SOI.loc['1982':'2017']
ind_df['PDO'] = PDO.loc['1982':'2017']
ind_df['IPO'] = IPO.loc['1982':'2017']
ind_df['NAO'] = NAO.loc['1982':'2017']
ind_df_filt = ind_df.rolling(window=36, center=True).mean().dropna()
ONI = ind_df['ENSO 3.4'].rolling(window=3, center=True).mean().dropna()
# load SST date already filtered
archivo = '/home/daniu/Documentos/datos_reynolds/output/filt36_anom_sst_monthly_reynolds_1982-2017_swa_sep_detrended.nc'
data = xr.open_dataset(archivo)
lon_sst = data.lon.values
lat_sst = data.lat.values
time_sst = data.time.values
# EOF calculator
coslat = np.cos(np.deg2rad(data.coords['lat'].values))
wgts = np.sqrt(coslat)[..., np.newaxis]
solver_sst = Eof(data.sst.values, weights=wgts)
cant_modos = 3
scaling_pc = 1 # * *0* : Un-scaled EOFs (default).
scaling_eof = 2 # * *1* : EOFs are divided by the square-root of their eigenvalues.
# * *2* : EOFs are multiplied by the square-root of their eigenvalues.
eof_sst = solver_sst.eofs(neofs=cant_modos, eofscaling=scaling_eof)
pc_sst = solver_sst.pcs(npcs=cant_modos, pcscaling=scaling_pc)
varfrac_sst = solver_sst.varianceFraction()
lambdas_sst = solver_sst.eigenvalues()
time = data.time.values
PCs = ['PC 1','PC 2','PC 3']
PC_df = pd.DataFrame(pc_sst, index=ind_df_filt.index[0:-1], columns=PCs)
PC_df.loc['2016-07-01'] = PC_df.loc['2016-06-01']
oni = ONI.loc['1983-07-01':'2016-07-01']
oni = (oni - oni.mean())/oni.std()
pdo = (ind_df_filt['PDO'] - ind_df_filt['PDO'].mean())/ind_df_filt['PDO'].std()
enso = (ind_df_filt['ENSO 3.4'] - ind_df_filt['ENSO 3.4'].mean())/ind_df_filt['ENSO 3.4'].std()
pc1 = (PC_df['PC 1'] - PC_df['PC 1'].mean())/PC_df['PC 1'].std()
time = enso.index.values
idx, = np.where(np.abs(oni)>1)
idx = list(idx)
from operator import itemgetter
from itertools import *
groups = []
for k, g in groupby(enumerate(idx), lambda x: x[0]-x[1]):
groups.append(list(map(itemgetter(1), g)))
figname = 'pdo-pc1-enso'
figprops = dict(figsize=(5, 2.7), dpi=72)
fig = plt.figure(**figprops)
ax = plt.axes([0.05, 0.05, 0.9, 0.9])
ax.plot(pdo, label='PDO', lw=.5, color='k', linestyle='--')
ax.plot(enso, label='ENSO 3.4', lw=.5, color='saddlebrown', linestyle='--')
ax.plot(pc1, label='PC1', lw=.5, color='b')
ax.axhline(y=0, lw=.5, color='k')
for g in groups:
if len(g)>3:
xi = g[0]
xf = g[-1]
ti = time[xi]
tf = time[xf]
print(ti, tf)
if oni[g[0]]>0:
ax.axvspan(ti, tf, facecolor='lightsalmon', alpha=0.15)
else:
ax.axvspan(ti, tf, facecolor='lightblue', alpha=0.15)
ax.legend(fontsize=6)
ax.tick_params('both', labelsize=6)
ax.set_xlabel('Time [Years]', fontsize=6)
ax.set_ylim([-2.5, 2.5])
fig.savefig('/home/daniu/Documentos/figuras/' + figname + '.pdf', bbox_inches='tight')
fig.savefig('/home/daniu/Documentos/figuras/' + figname, dpi=300, bbox_inches='tight')