forked from LisaNeef/DARTpy
-
Notifications
You must be signed in to change notification settings - Fork 0
/
VORTEX.py
122 lines (96 loc) · 4.45 KB
/
VORTEX.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
119
120
121
122
## Python module for various diagnostics related to the stratospheric polar vortex
# load necessary modules
import DART as dart
import MJO as mjo
import numpy as np
import matplotlib.pyplot as plt
import brewer2mpl
import plot_tools
#--------------------------------------------------------------------
def plot_climate_indices(E,index_name,copies_to_plot,climatology_option = 'NODA',hostname='taurus',verbose=False):
"""
This subroutine computes a bunch of simple climate indices
for a given experiment, and plots them for all copies given by 'copies_to_plot'
INPUTS:
copies_to_plot: list containing keywords for what copies to plot. Here are the options:
+ any valid copystring in DART output data (e.g. "ensemble member 1")
+ 'ensemble' = plot the entire ensemble
+ 'ensemble mean' = plot the ensemble mean
+ 'operational' = plot the operational value of this index
"""
# create a list of copies to load
copy_list = []
if "copystring" in copies_to_plot:
copy_list.append(E['copystring'])
if ("ensemble" in copies_to_plot):
N = dart.get_ensemble_size_per_run(E['exp_name'])
for iens in np.arange(1,N+1):
if iens < 10:
spacing = ' '
else:
spacing = ' '
copy_list.append("ensemble member"+spacing+str(iens))
if ("ensemble mean" in copies_to_plot):
copy_list.append("ensemble mean")
# retrieve desired index for all the copies in the list
L = []
for copy in copy_list:
E['copystring'] = copy
x = compute_climate_indices(E,index_name,climatology_option,hostname,verbose)
L.append(x)
# plot it
for copy,climate_index in zip(copy_list,L):
# define a color for the ensemble mean - depending on experiment
lcolor = "#000000"
if E['exp_name'] == 'W0910_NODA':
bmap = brewer2mpl.get_map('Dark2', 'qualitative', 7)
lcolor = bmap.hex_colors[1]
if E['exp_name'] == 'W0910_GLOBAL':
# make this one black, since it's sort of the reference
lcolor = "#000000"
if E['exp_name'] == 'W0910_TROPICS':
bmap = brewer2mpl.get_map('Dark2', 'qualitative', 7)
lcolor = bmap.hex_colors[0]
# make the ensemble a lighter version of their original color
ensemble_color = plot_tools.colorscale(lcolor, 1.6)
# here is the plot
if (copy == 'ensemble mean'):
plt.plot(E['daterange'],climate_index,color=lcolor,linewidth=2)
if "ensemble member" in copy:
plt.plot(E['daterange'],climate_index,color=ensemble_color,linewidth=1)
def compute_climate_indices(E,index_name,climatology_option = 'NODA',hostname='taurus',verbose=False):
"""
This subroutine computes various simple climate indices for a dataset
defined by an experiment dictionary.
Currently supporting the following indices:
+ 'Aleutian Low' index of Garfinkel et al. (2010)
+ 'East European High' index of Garfinkel et al. (2010)
+ 'AO Proxy' -- Polar Cap GPH Anomaly at 500hPa -- it's a proxy for the AO suggested by Cohen et al. (2002)
* note however that we define the polar cap as everything north of 70N, I think Cohen et al do 60N
+ 'Vortex Strength' -- Polar Cap GPH Anomaly averaged 3-30hPa -- it's a measure of vortex strength suggested by Garfinkel et al. 2012
"""
# modify the experiment dictionary to retrieve the right index
EI = dart.climate_index_dictionaries(index_name)
E['levrange'] = EI['levrange']
E['latrange'] = EI['latrange']
E['lonrange'] = EI['lonrange']
E['variable'] = EI['variable']
# for all indices defined so far, compute the anomaly
# with respect to climatology
# this uses an anomaly subroutine from the MJO module
A,C,lat,lon,lev = mjo.ano(E,climatology_option = climatology_option,verbose=verbose)
# Aleutian Low and East European high indices are single points, so just return the anomaly
if (index_name == 'Aleutian Low') or (index_name == 'East European High'):
index_out = A
# for the Polar Cap GPH -based indices, average over latitude and longitude
# here can use a subroutine written for MJO stuff in the MJO module
if (index_name == 'AO Proxy') or (index_name == 'Vortex Strength'):
lat1,lon1,Aave = mjo.aave(E,A,lat,lon,season=None,variable_name=None,averaging_dimension='all')
# for the AO proxy, reverse the sign so that it's more intuitive -- a positive GPH anomaly is related to a negative AO
if (index_name == 'AO Proxy') or (index_name == 'Vortex Strength'):
index_out = -Aave
# for vortex strength, average between 3 and 30 hPa
if (index_name == 'Vortex Strength'):
index_out = np.nanmean(Aave,axis=0)
# return index over desired daterange
return index_out