tcheck_stationarity.py - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
(HTM) git clone git://src.adamsgaard.dk/pism
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
tcheck_stationarity.py (6567B)
---
1 #!/usr/bin/env python3
2
3 # @package check_stationarity
4 # \author Andy Aschwanden, University of Alaska Fairbanks, USA
5 # \brief Script to evaluate stationarity of a variable.
6 # \details Given a time series of a variabale \f$ X \f$, it computes the
7 # time rate of change of the p-norm of \f$ X \f$.
8 # \f[\frac{d}{dt} || X ||_{p} =
9 # \frac{ \left [\sum_{i}^{m} \left( X_{i}^{n+1} - X_{i}^{n} \right)^{p} \right]
10 # ^{1/p} }{ t^{n+1} - t^{n+1} } \f] where \f$ X_{i}^{n} \f$ is the value at
11 # time \f$ n \f$ and coordinate \f$ i \f$.
12
13
14 try:
15 from netCDF4 import Dataset as CDF
16 except:
17 print("netCDF4 is not installed!")
18 sys.exit(1)
19
20 import numpy as np
21 import pylab as plt
22 from optparse import OptionParser
23
24 __author__ = "Andy Aschwanden"
25
26 # If no threshold is given, set it to 1e1
27 THRESHOLD = 1e1
28 # This is the default variable to calculate norm from
29 X = 'enthalpybase'
30 # Default norm is the euclidian norm, 2-norm
31 PNORM = float(2)
32
33 # Set up the option parser
34 parser = OptionParser()
35 parser.usage = "usage: %prog [options] FILE"
36 parser.description = '''Check stationarity of a variable in FILE by calculating
37 the rate of change of its p-norm. That is
38 d/dt || X ||_{p} =
39 (\sum_{i}^{m} (E_{i}^{n+1}-E_{i}^{n})^p)^(1/p)/(t^{n+1}-t^{n+1}),
40 where E_{i}^{n} is the value at time n and coordinate i.'''
41 parser.add_option("-p", "--pnorm", dest="pnorm", type='float',
42 help="use P norm (default p = 2)",
43 metavar="P", default=PNORM)
44 parser.add_option("-s", "--stride", dest="stride", type="int",
45 help="stride, plot only every stride value", default=1)
46 parser.add_option("-t", "--threshold", dest="threshold",
47 help="draws a line horizontal line at THRESHOLD",
48 metavar="THRESHOLD", default=THRESHOLD)
49 parser.add_option("-v", "--variable", dest="varname", type='string',
50 help="calculate from from variable X (default=enthalpybase)",
51 metavar="VAR", default=X)
52 # Run the option parser
53 (options, args) = parser.parse_args()
54
55 # Assign threshold
56 threshold = options.threshold
57 # Assign variable name
58 varname = options.varname
59 # Assign norm
60 p = options.pnorm
61 # Assign stride value
62 stride = options.stride
63
64 if len(args) == 1:
65 # Retrieve command line argument (name of input file)
66 infile = args[0]
67 else:
68 print('wrong number arguments, must be 1')
69 # If number of arguments is wrong, print help to give user some clues.
70 parser.print_help()
71 exit(0)
72
73
74 # Opens a netCDF file.
75 #
76 # Open netCDF file and check that time dimension exists.
77 # On success, a pointer to the file is returned, otherwise an error is issued.
78 def open_ncfile(infile):
79
80 try:
81 nc = CDF(infile, 'r')
82 try:
83 'time' in list(nc.variables.keys())
84 except:
85 print('Variable t not found in file %s' % infile)
86 except IOError:
87 pass
88
89 return nc
90
91
92 # Calculate rate of change.
93 #
94 # Calculate rate of change of p-norm of variable \f$ X \f$.
95 def getRateOfChange(t, X, p, varname):
96
97 # differenetiate time along time axis
98 dt = np.diff(t[::stride], axis=0)
99
100 Xdim = X.ndim
101
102 if Xdim == 1:
103 X = X[::stride]
104 dXp = np.diff(X)
105 dXpdt = dXp / dt
106 elif Xdim == 3:
107 if 'mask' in list(nc.variables.keys()):
108 mask = np.array(nc.variables['mask'][::stride, :, :]) # (time,y,x)
109 k = np.nonzero((mask == 1) ^ (mask == 2) ^ (mask == 3))
110 mask2 = np.ones_like(mask)
111 mask2[k] = 0
112 # use masked values (i.e. ignore ocean and ice-free land)
113 X = np.ma.array(data=X[::stride, :, :], mask=mask2)
114 else:
115 X = np.array(X[::stride, :, :])
116
117 dX = np.diff(X, axis=0)
118 nt, ny, nx = dX.shape
119 # convert (t,y,x) -> (t,y*x) so that np.nansum needs
120 # to be applied only once
121 dX = dX.reshape(nt, nx * ny)
122 dXp = (np.nansum(dX ** p, axis=1)) ** (1. / p)
123 dXpdt = dXp / dt
124 else:
125 print('error: dim n = %i of variable %s not supported, must be 1 or 3'
126 % (Xdim, varname))
127
128 return dXpdt
129
130
131 # Unit converter
132 def unit_converter(data, inunit, outunit):
133 '''
134 Unit converter. Takes an (numpy) array, valid udunits inunits and outunits
135 as strings, and returns the array in outunits.
136
137 Parameters
138 ----------
139 data : array_like
140 inunit : string
141 unit to convert from, must be UDUNITS-compatible string
142 outunit : string
143 unit to conver to, must be UDUNITS-compatible string
144
145 Returns
146 -------
147 out : array_like
148
149 Example
150 -------
151 >>> import numpy as np
152 >>> c = Converter("kg","Gt")
153 >>> out = c(np.array([1,2])*1e12)
154 >>> out = array([ 1., 2.])
155 '''
156
157 inunit = str(inunit)
158 outunit = str(outunit)
159
160 if not (inunit == outunit):
161 try:
162 from udunits import Converter
163 c = Converter(inunit, outunit)
164 outdata = c(data)
165 except:
166 print("No udunits module found, you're on your own.\n -> I am assuming input unit is m, and will convert to km.\n -> Installation of Constantine's awesome python wrapper for udunits is highly recommended.\n -> Download it from https://code.google.com/p/python-udunits/.")
167 c = 1. / 1e3
168 outdata = c * data
169 else:
170 outdata = data
171
172 return outdata
173
174
175 if __name__ == "__main__":
176
177 # Open netCDF file
178 nc = open_ncfile(infile)
179
180 # time variable t
181 t_units = nc.variables['time'].units
182 try:
183 date_prefix = 'since'
184 unit_array = t_units.split(date_prefix)
185 except:
186 date_prefix = 'before'
187 unit_array = t_units.split(date_prefix)
188 out_units = 'years ' + date_prefix + unit_array[1]
189
190 t = unit_converter(nc.variables['time'][:], t_units, out_units)
191 if varname in list(nc.variables.keys()):
192 var = nc.variables[varname]
193 else:
194 print("error: variable '%s' not found in %s" % (varname, infile))
195 exit(0)
196
197 # Calculate rate of change from time t, variable var,
198 # norm p and variable name varname
199 dVpdt = getRateOfChange(t, var, p, varname)
200
201 # Make plot with log-scale y-axis
202 plt.figure()
203 plt.hold(True)
204 plt.semilogy(t[1::stride], dVpdt, 'b', lw=2)
205 plt.plot([t[1:][0], t[1:][-1]], [threshold, threshold], 'r', lw=1)
206 plt.xlabel(('time [%s]' % out_units))
207 plt.ylabel(('d %s / dt [%s/a]' % (varname, var.units)))
208 plt.title(('rate of change of %s' % varname), fontsize=12)
209 plt.show()
210
211 # eventually, close the file
212 nc.close()