-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy patheof_rotation_example.py
118 lines (86 loc) · 5.29 KB
/
eof_rotation_example.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
"""
Calculating OMI with rotation
Use this file to generate the rotated OMI as described in Weidman, S., Kleiner, N., & Kuang, Z. (2022). A rotation
procedure to improve seasonally varying empirical orthogonal function bases for MJO indices. Geophysical Research Letters,
49, e2022GL099998. https://doi.org/10.1029/2022GL099998.
** UPDATED JANUARY 2023 **
The rotation algorithm has now been fully integrated into the mjoindices package which
can be found on github at https://github.com/cghoffmann/mjoindices. The package description was originally published in
Hoffmann, C.G., Kiladis, G.N., Gehne, M. and von Savigny, C., 2021. A Python Package to Calculate the OLR-Based Index
of the Madden- Julian-Oscillation (OMI) in Climate Science and Weather Forecasting. Journal of Open Research Software, 9(1),
p.9. DOI: http://doi.org/10.5334/jors.331.
The example files in this folder now use the updated mjoindices package. Installation instructions are included in the
github link above. This file is simply a guide to using the rotation postprocessing option within that package. Example data
is also uploaded in that package, but you can use your own OLR dataset if desired. Please read the README and citation
information in the mjoindices documentation.
***************************
The original OMI algorithm is described in Kiladis, G.N., J. Dias, K.H. Straub, M.C. Wheeler, S.N. Tulich, K. Kikuchi,
K.M. Weickmann, and M.J. Ventrice, 2014: A Comparison of OLR and Circulation-Based Indices for Tracking the MJO.
Mon. Wea. Rev., 142, 1697–1715, https://doi.org/10.1175/MWR-D-13-00301.1
"""
from pathlib import Path
import os.path
import mjoindices.olr_handling as olr
import mjoindices.omi.omi_calculator as omi
import mjoindices.empirical_orthogonal_functions as eof
import mjoindices.principal_components as pc
import numpy as np
# create filename paths for OLR, EOFs, and PCs
OLR_path = Path(os.path.abspath('')) / 'observed_data' / 'olr.day.mean.nc'
eofs_dir = Path(os.path.abspath('')) / 'EOFs'
pc_dir = Path(os.path.abspath('')) / 'PCs'
# load in OLR dataset
raw_olr = olr.load_noaa_interpolated_olr(OLR_path)
# Restrict the OLR dataset to the dates used in the paper. This can be changed to use a longer time period,
# but EOF results will change slightly.
time_start = np.datetime64('1979-01-01')
time_end = np.datetime64('2012-12-31')
shorter_olr = olr.restrict_time_coverage(raw_olr, time_start, time_end)
# Interpolate the dataset to the same grid as used in the original OMI calcluation
# (not required, but testing has only been done using the interpolated grid).
interpolated_olr = olr.interpolate_spatial_grid_to_original(shorter_olr)
# Calculate the EOFs using the rotation algorithm. This may take an hour or more, especially on a personal computer.
# See mjoindices documentation for details on these parameters.
# set parameters for rotation postprocessing step
rotation_pp_params = {"sign_doy1reference": True}
# calculate EOFs using rotation postprocessing
rot_eofs = omi.calc_eofs_from_olr(interpolated_olr,
leap_year_treatment="original",
eofs_postprocessing_type="eof_rotation",
eofs_postprocessing_params=rotation_pp_params)
# Save EOFs in your defined EOF folder.
rot_eofs.save_all_eofs_to_npzfile(eofs_dir / 'EOFs_rotated.npz')
# Similarly, calculate the EOFs from the original OMI algorithm. This also may take up to an hour or more.
# Change interpolate_eofs to True to directly reproduce the method from Kiladis et al., 2014. This will interpolate
# the EOFs during the beginning of November.
kiladis_pp_params = {"sign_doy1reference": True,
"interpolate_eofs": False}
norot_eofs = omi.calc_eofs_from_olr(interpolated_olr,
leap_year_treatment="original",
eofs_postprocessing_type="kiladis2014",
eofs_postprocessing_params=kiladis_pp_params)
norot_eofs.save_all_eofs_to_npzfile(eofs_dir / 'EOFs_unrotated.npz')
# If you have already calculated the EOFs, uncomment and reload them from your EOF folder.
# eof_rot = eofs_dir / 'EOFs_rotated.npz'
# rot_eofs = eof.restore_all_eofs_from_npzfile(eof_rot)
# eof_norot = eofs_dir / 'EOFs_unrotated.npz'
# norot_eofs = eof.restore_all_eofs_from_npzfile(eof_norot)
# Now use the EOFs to calculate the principal components.
# The time variable is restricted to the same years as the EOF calculation here.
pcs_rot = omi.calculate_pcs_from_olr(raw_olr,
rot_eofs,
time_start,
time_end,
use_quick_temporal_filter=False)
pcs_norot = omi.calculate_pcs_from_olr(raw_olr,
norot_eofs,
time_start,
time_end,
use_quick_temporal_filter=False)
pcs_rot.save_pcs_to_txt_file(pc_dir / 'PCs_rotated.txt')
pcs_norot.save_pcs_to_txt_file(pc_dir / 'PCs_unrotated.txt')
# If you have already calculated PCs, uncomment and reload them here.
pc_rot_path = pc_dir / 'PCs_rotated.txt'
pcs_rot = pc.load_pcs_from_txt_file(pc_rot_path)
pc_norot_path = pc_dir / 'PCs_unrotated.txt'
pcs_norot = pc.load_pcs_from_txt_file(pc_norot_path)