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()