| 1 | ############################################### |
|---|
| 2 | ## PLANETOPLOT ## |
|---|
| 3 | ## --> PPCLASS ## |
|---|
| 4 | ## A generic and versatile Python module ## |
|---|
| 5 | ## ... to read netCDF files and plot ## |
|---|
| 6 | ############################################### |
|---|
| 7 | ## Author: Aymeric Spiga. 02-03/2013 ## |
|---|
| 8 | ############################################### |
|---|
| 9 | # python built-in librairies |
|---|
| 10 | import os |
|---|
| 11 | import time as timelib |
|---|
| 12 | import pickle |
|---|
| 13 | # added librairies |
|---|
| 14 | import numpy as np |
|---|
| 15 | import netCDF4 |
|---|
| 16 | import matplotlib.pyplot as mpl |
|---|
| 17 | # personal librairies |
|---|
| 18 | import ppplot |
|---|
| 19 | import ppcompute |
|---|
| 20 | ############################################### |
|---|
| 21 | |
|---|
| 22 | ################################### |
|---|
| 23 | #### HEADER ## |
|---|
| 24 | #### ... executed when imported ## |
|---|
| 25 | ################################### |
|---|
| 26 | # where settings files are located... |
|---|
| 27 | whereset = None |
|---|
| 28 | whereset = ppcompute.findset(whereset) |
|---|
| 29 | # ... we load user-defined automatic settings from set_ppclass.txt |
|---|
| 30 | zefile = "set_ppclass.txt" |
|---|
| 31 | glob_listx = [] ; glob_listy = [] ; glob_listz = [] ; glob_listt = [] |
|---|
| 32 | glob_listarea = [] |
|---|
| 33 | try: |
|---|
| 34 | f = open(whereset+zefile, 'r') ; lines = f.readlines() |
|---|
| 35 | for stuff in lines[5].strip().split(';'): glob_listx.append(stuff) |
|---|
| 36 | for stuff in lines[8].strip().split(';'): glob_listy.append(stuff) |
|---|
| 37 | for stuff in lines[11].strip().split(';'): glob_listz.append(stuff) |
|---|
| 38 | for stuff in lines[14].strip().split(';'): glob_listt.append(stuff) |
|---|
| 39 | for stuff in lines[17].strip().split(';'): glob_listarea.append(stuff) |
|---|
| 40 | except IOError: |
|---|
| 41 | print "warning: "+zefile+" not in "+whereset+" ; no presets." |
|---|
| 42 | |
|---|
| 43 | ################################## |
|---|
| 44 | #### USEFUL GENERIC FUNCTIONS #### |
|---|
| 45 | ################################## |
|---|
| 46 | |
|---|
| 47 | # inspect variables and dimensions in a netCDF file |
|---|
| 48 | def inspect(filename): |
|---|
| 49 | print "**** INSPECT FILE",filename |
|---|
| 50 | test = netCDF4.Dataset(filename) |
|---|
| 51 | print "**** VARIABLES: ",test.variables.keys() |
|---|
| 52 | for dim in test.dimensions.keys(): |
|---|
| 53 | output = "**** DIMENSION: "+str(dim)+" "+str(len(test.dimensions[dim])) |
|---|
| 54 | try: output = output + " ----> "+str(test.variables[dim][0])+" "+str(test.variables[dim][-1]) |
|---|
| 55 | except: pass |
|---|
| 56 | print output ; output = "" |
|---|
| 57 | |
|---|
| 58 | # check a tab and exit if wrong. if just one string make it a list. |
|---|
| 59 | # (if allownumber, convert this into a string). |
|---|
| 60 | def checktab(tab,mess="",allownone=False,allownumber=False): |
|---|
| 61 | if tab is None: |
|---|
| 62 | if not allownone: print "pp.define: no "+mess ; exit() |
|---|
| 63 | else: pass |
|---|
| 64 | else: |
|---|
| 65 | if not isinstance(tab, list): |
|---|
| 66 | if isinstance(tab, str): |
|---|
| 67 | tab = [tab] |
|---|
| 68 | elif (isinstance(tab, int) or isinstance(tab, float)) and allownumber: |
|---|
| 69 | tab = [str(tab)] |
|---|
| 70 | else: |
|---|
| 71 | print "pp.define: "+mess+" should be either a string or a list of strings!" ; exit() |
|---|
| 72 | elif isinstance(tab, list): |
|---|
| 73 | if isinstance(tab[0],str): |
|---|
| 74 | pass |
|---|
| 75 | elif (isinstance(tab[0], int) or isinstance(tab[0], float)) and allownumber: |
|---|
| 76 | for iii in range(len(tab)): tab[iii] = str(tab[iii]) |
|---|
| 77 | else: |
|---|
| 78 | print "pp.define: "+mess+" should be either a string or a list of strings!" ; exit() |
|---|
| 79 | return tab |
|---|
| 80 | |
|---|
| 81 | # determine which method is to be applied to a given dimension |
|---|
| 82 | def findmethod(tab): |
|---|
| 83 | if tab is None: output = "free" |
|---|
| 84 | elif tab[0,0] != tab[0,1]: output = "comp" |
|---|
| 85 | else: output = "fixed" |
|---|
| 86 | return output |
|---|
| 87 | |
|---|
| 88 | # read what is given by the user (version of T. Navarro) |
|---|
| 89 | def readslices(saxis): |
|---|
| 90 | if saxis == None: |
|---|
| 91 | zesaxis = None |
|---|
| 92 | else: |
|---|
| 93 | zesaxis = np.empty((len(saxis),2)) |
|---|
| 94 | for i in range(len(saxis)): |
|---|
| 95 | a = separatenames(saxis[i]) |
|---|
| 96 | if len(a) == 1: |
|---|
| 97 | zesaxis[i,:] = float(a[0]) |
|---|
| 98 | else: |
|---|
| 99 | zesaxis[i,0] = float(a[0]) |
|---|
| 100 | zesaxis[i,1] = float(a[1]) |
|---|
| 101 | return zesaxis |
|---|
| 102 | |
|---|
| 103 | # look for comas in the input name to separate different names (files, variables,etc ..) |
|---|
| 104 | # (needed by readslices) |
|---|
| 105 | def separatenames (name): |
|---|
| 106 | if name is None: names = None |
|---|
| 107 | else: |
|---|
| 108 | names = [] ; stop = 0 ; currentname = name |
|---|
| 109 | while stop == 0: |
|---|
| 110 | indexvir = currentname.find(',') |
|---|
| 111 | if indexvir == -1: stop = 1 ; name1 = currentname |
|---|
| 112 | else: name1 = currentname[0:indexvir] |
|---|
| 113 | names = np.concatenate((names,[name1])) |
|---|
| 114 | currentname = currentname[indexvir+1:len(currentname)] |
|---|
| 115 | return names |
|---|
| 116 | |
|---|
| 117 | ####################### |
|---|
| 118 | ### THE MAIN OBJECT ### |
|---|
| 119 | ####################### |
|---|
| 120 | class pp(): |
|---|
| 121 | |
|---|
| 122 | # print out a help string when help is invoked on the object |
|---|
| 123 | def __repr__(self): |
|---|
| 124 | whatprint = 'pp object. \"help(pp)\" for more information\n' |
|---|
| 125 | return whatprint |
|---|
| 126 | |
|---|
| 127 | # default settings |
|---|
| 128 | # -- user can define settings by two methods. |
|---|
| 129 | # -- 1. yeah = pp(file="file.nc") |
|---|
| 130 | # -- 2. yeah = pp() ; yeah.file = "file.nc" |
|---|
| 131 | def __init__(self,file=None,var="notset",\ |
|---|
| 132 | filegoal=None,vargoal=None,\ |
|---|
| 133 | x=None,y=None,z=None,t=None,\ |
|---|
| 134 | stridex=1,stridey=1,\ |
|---|
| 135 | stridez=1,stridet=1,\ |
|---|
| 136 | compute="mean",\ |
|---|
| 137 | verbose=False,noproj=False,\ |
|---|
| 138 | superpose=False,\ |
|---|
| 139 | plotin=None,\ |
|---|
| 140 | forcedimplot=-1,\ |
|---|
| 141 | out="gui",\ |
|---|
| 142 | filename="myplot",\ |
|---|
| 143 | folder="./",\ |
|---|
| 144 | includedate=True,\ |
|---|
| 145 | xlabel=None,ylabel=None,\ |
|---|
| 146 | xcoeff=None,ycoeff=None,\ |
|---|
| 147 | proj=None,\ |
|---|
| 148 | vmin=None,vmax=None,\ |
|---|
| 149 | div=None,\ |
|---|
| 150 | colorb=None,\ |
|---|
| 151 | lstyle=None,\ |
|---|
| 152 | marker=None,\ |
|---|
| 153 | color=None,\ |
|---|
| 154 | label=None,\ |
|---|
| 155 | changetime=None,\ |
|---|
| 156 | title=None): |
|---|
| 157 | self.request = None |
|---|
| 158 | self.nfin = 0 ; self.nvin = 0 |
|---|
| 159 | self.nplotx = None ; self.nploty = None |
|---|
| 160 | self.nplotz = None ; self.nplott = None |
|---|
| 161 | self.status = "init" |
|---|
| 162 | self.fig = None ; self.subv = None ; self.subh = None |
|---|
| 163 | self.n = 0 ; self.howmanyplots = 0 |
|---|
| 164 | self.nplot = 0 |
|---|
| 165 | self.p = None |
|---|
| 166 | self.customplot = False |
|---|
| 167 | self.allfield = None |
|---|
| 168 | ## what could be defined by the user |
|---|
| 169 | self.file = file |
|---|
| 170 | self.var = var |
|---|
| 171 | self.filegoal = filegoal |
|---|
| 172 | self.vargoal = vargoal |
|---|
| 173 | self.x = x ; self.y = y ## if None, free dimension |
|---|
| 174 | self.z = z ; self.t = t ## if None, free dimension |
|---|
| 175 | self.stridex = stridex ; self.stridey = stridey |
|---|
| 176 | self.stridez = stridez ; self.stridet = stridet |
|---|
| 177 | self.compute = compute |
|---|
| 178 | self.verbose = verbose |
|---|
| 179 | self.noproj = noproj |
|---|
| 180 | self.plotin = plotin |
|---|
| 181 | self.superpose = superpose |
|---|
| 182 | self.forcedimplot = forcedimplot |
|---|
| 183 | self.out = out |
|---|
| 184 | self.filename = filename |
|---|
| 185 | self.folder = folder |
|---|
| 186 | self.includedate = includedate |
|---|
| 187 | self.changetime = changetime |
|---|
| 188 | ## here are user-defined plot settings |
|---|
| 189 | ## -- if not None, valid on all plots in the pp() objects |
|---|
| 190 | self.xlabel = xlabel ; self.xcoeff = xcoeff |
|---|
| 191 | self.ylabel = ylabel ; self.ycoeff = ycoeff |
|---|
| 192 | self.proj = proj |
|---|
| 193 | self.vmin = vmin ; self.vmax = vmax |
|---|
| 194 | self.div = div |
|---|
| 195 | self.colorb = colorb |
|---|
| 196 | self.lstyle = lstyle |
|---|
| 197 | self.marker = marker |
|---|
| 198 | self.color = color |
|---|
| 199 | self.label = label |
|---|
| 200 | self.title = title |
|---|
| 201 | |
|---|
| 202 | # print status |
|---|
| 203 | def printstatus(self): |
|---|
| 204 | if self.filename == "THIS_IS_A_CLONE": |
|---|
| 205 | pass |
|---|
| 206 | else: |
|---|
| 207 | print "**** PPCLASS. Done step: " + self.status |
|---|
| 208 | |
|---|
| 209 | ##################################################### |
|---|
| 210 | # EMULATE OPERATORS + - * / ** << FOR PP() OBJECTS # |
|---|
| 211 | ##################################################### |
|---|
| 212 | |
|---|
| 213 | # define the operation << |
|---|
| 214 | # ... e.g. obj2 << obj1 |
|---|
| 215 | # ... means: get init for pp object obj2 from another pp object obj1 |
|---|
| 216 | # ... (this should solve the affectation trap obj2 = obj1) |
|---|
| 217 | def __lshift__(self,other): |
|---|
| 218 | if other.__class__.__name__ == "pp": |
|---|
| 219 | self.file = other.file |
|---|
| 220 | self.var = other.var |
|---|
| 221 | self.filegoal = other.filegoal |
|---|
| 222 | self.vargoal = other.vargoal |
|---|
| 223 | self.x = other.x ; self.y = other.y ## if None, free dimension |
|---|
| 224 | self.z = other.z ; self.t = other.t ## if None, free dimension |
|---|
| 225 | self.stridex = other.stridex ; self.stridey = other.stridey |
|---|
| 226 | self.stridez = other.stridez ; self.stridet = other.stridet |
|---|
| 227 | self.verbose = other.verbose |
|---|
| 228 | self.noproj = other.noproj |
|---|
| 229 | self.plotin = other.plotin |
|---|
| 230 | self.superpose = other.superpose |
|---|
| 231 | self.forcedimplot = other.forcedimplot |
|---|
| 232 | self.out = other.out |
|---|
| 233 | self.filename = other.filename |
|---|
| 234 | self.folder = other.folder |
|---|
| 235 | self.xlabel = other.xlabel ; self.xcoeff = other.xcoeff |
|---|
| 236 | self.ylabel = other.ylabel ; self.ycoeff = other.ycoeff |
|---|
| 237 | self.proj = other.proj |
|---|
| 238 | self.vmin = other.vmin ; self.vmax = other.vmax |
|---|
| 239 | self.div = other.div |
|---|
| 240 | self.colorb = other.colorb |
|---|
| 241 | self.lstyle = other.lstyle |
|---|
| 242 | self.marker = other.marker |
|---|
| 243 | self.color = other.color |
|---|
| 244 | self.label = other.label |
|---|
| 245 | self.title = other.title |
|---|
| 246 | self.includedate = other.includedate |
|---|
| 247 | self.changetime = other.changetime |
|---|
| 248 | else: |
|---|
| 249 | print "!! ERROR !! argument must be a pp object." ; exit() |
|---|
| 250 | |
|---|
| 251 | # check the compatibility of two objects for operations |
|---|
| 252 | # --> if other is a pp class, test sizes and return isnum = False |
|---|
| 253 | # --> if other is an int or a float, return isnum = True |
|---|
| 254 | # --> otherwise, just print an error and exit |
|---|
| 255 | def checktwo(self,other): |
|---|
| 256 | if other.__class__.__name__ == "pp": |
|---|
| 257 | isnum = False |
|---|
| 258 | if self.status in ["init","defined"] or other.status in ["init","define"]: |
|---|
| 259 | print "!! ERROR !! Please use .retrieve to get fields for plots with one of your pp operands." ; exit() |
|---|
| 260 | if self.nfin != other.nfin or \ |
|---|
| 261 | self.nvin != other.nvin or \ |
|---|
| 262 | self.nplott != other.nplott or \ |
|---|
| 263 | self.nplotz != other.nploty or \ |
|---|
| 264 | self.nploty != other.nploty or \ |
|---|
| 265 | self.nplotx != other.nplotx : |
|---|
| 266 | print "!! ERROR !! The two operands do not have the same number of files, variables, t z y x requests." |
|---|
| 267 | exit() |
|---|
| 268 | elif isinstance(other,int) or isinstance(other,float): |
|---|
| 269 | isnum = True |
|---|
| 270 | else: |
|---|
| 271 | print "!! ERROR !! The operand is neither a pp class nor an integer or a float." ; exit() |
|---|
| 272 | return isnum |
|---|
| 273 | |
|---|
| 274 | # define a selective copy of a pp() object for operations |
|---|
| 275 | # ... copy.copy() is not conservative (still acts like a pointer) |
|---|
| 276 | # ... copy.deepcopy() does not work with netCDF objects |
|---|
| 277 | # so what is done here is a copy of everything except |
|---|
| 278 | # (to avoid sharing with self and therefore modifying self through operations) |
|---|
| 279 | # - request attribute of pp() object |
|---|
| 280 | # - field attribute of the onerequest() objects |
|---|
| 281 | def selective_copy(self): |
|---|
| 282 | if self.status in ["init","defined"]: |
|---|
| 283 | print "!! ERROR !! Please use .retrieve to get fields for the object you want to copy from." ; exit() |
|---|
| 284 | the_clone = pp() |
|---|
| 285 | for k, v in vars(self).items(): |
|---|
| 286 | if k != "request": |
|---|
| 287 | setattr(the_clone,k,v) |
|---|
| 288 | the_clone.verbose = False |
|---|
| 289 | the_clone.filename = "THIS_IS_A_CLONE" # trick to avoid additional outputs |
|---|
| 290 | the_clone.define() |
|---|
| 291 | for i in range(self.nfin): |
|---|
| 292 | for j in range(self.nvin): |
|---|
| 293 | for t in range(self.nplott): |
|---|
| 294 | for z in range(self.nplotz): |
|---|
| 295 | for y in range(self.nploty): |
|---|
| 296 | for x in range(self.nplotx): |
|---|
| 297 | obj_ref = self.request[i][j][t][z][y][x] |
|---|
| 298 | obj = the_clone.request[i][j][t][z][y][x] |
|---|
| 299 | for k, v in vars(obj_ref).items(): |
|---|
| 300 | if k != "field": |
|---|
| 301 | setattr(obj,k,v) |
|---|
| 302 | the_clone.status = "retrieved" |
|---|
| 303 | the_clone.filename = self.filename |
|---|
| 304 | return the_clone |
|---|
| 305 | |
|---|
| 306 | # define the operation + on two objects. or with an int/float. |
|---|
| 307 | # ... with selective_copy the self object is not modified. |
|---|
| 308 | def __add__(self,other): |
|---|
| 309 | isnum = self.checktwo(other) |
|---|
| 310 | the_clone = self.selective_copy() |
|---|
| 311 | for i in range(self.nfin): |
|---|
| 312 | for j in range(self.nvin): |
|---|
| 313 | for t in range(self.nplott): |
|---|
| 314 | for z in range(self.nplotz): |
|---|
| 315 | for y in range(self.nploty): |
|---|
| 316 | for x in range(self.nplotx): |
|---|
| 317 | obj = the_clone.request[i][j][t][z][y][x] |
|---|
| 318 | obj_ref = self.request[i][j][t][z][y][x] |
|---|
| 319 | if not isnum: |
|---|
| 320 | ope = other.request[i][j][t][z][y][x].field |
|---|
| 321 | if obj_ref.field.shape != ope.shape: |
|---|
| 322 | print "!! ERROR !! The two fields for operation do not have the same shape.",self.field.shape,other.field.shape |
|---|
| 323 | exit() |
|---|
| 324 | else: |
|---|
| 325 | ope = other |
|---|
| 326 | goal = self.vargoal[j] + self.filegoal[i] |
|---|
| 327 | if ("vector" in goal) or ("contour" in goal): |
|---|
| 328 | if self.verbose: print "!! WARNING !! No operation was made on contours and vectors. This can be debatted actually." |
|---|
| 329 | obj.field = obj_ref.field |
|---|
| 330 | else: |
|---|
| 331 | obj.field = obj_ref.field + ope |
|---|
| 332 | return the_clone |
|---|
| 333 | |
|---|
| 334 | # define the operation - on two objects. or with an int/float. |
|---|
| 335 | # ... with selective_copy the self object is not modified. |
|---|
| 336 | def __sub__(self,other): |
|---|
| 337 | isnum = self.checktwo(other) |
|---|
| 338 | the_clone = self.selective_copy() |
|---|
| 339 | for i in range(self.nfin): |
|---|
| 340 | for j in range(self.nvin): |
|---|
| 341 | for t in range(self.nplott): |
|---|
| 342 | for z in range(self.nplotz): |
|---|
| 343 | for y in range(self.nploty): |
|---|
| 344 | for x in range(self.nplotx): |
|---|
| 345 | obj = the_clone.request[i][j][t][z][y][x] |
|---|
| 346 | obj_ref = self.request[i][j][t][z][y][x] |
|---|
| 347 | if not isnum: |
|---|
| 348 | ope = other.request[i][j][t][z][y][x].field |
|---|
| 349 | if obj_ref.field.shape != ope.shape: |
|---|
| 350 | print "!! ERROR !! The two fields for operation do not have the same shape.",self.field.shape,other.field.shape |
|---|
| 351 | exit() |
|---|
| 352 | else: |
|---|
| 353 | ope = other |
|---|
| 354 | goal = self.vargoal[j] + self.filegoal[i] |
|---|
| 355 | if ("vector" in goal) or ("contour" in goal): |
|---|
| 356 | if self.verbose: print "!! WARNING !! No operation was made on contours and vectors. This can be debatted actually." |
|---|
| 357 | obj.field = obj_ref.field |
|---|
| 358 | else: |
|---|
| 359 | obj.field = obj_ref.field - ope |
|---|
| 360 | return the_clone |
|---|
| 361 | |
|---|
| 362 | # define the operation * on two objects. or with an int/float. |
|---|
| 363 | # ... with selective_copy the self object is not modified. |
|---|
| 364 | def __mul__(self,other): |
|---|
| 365 | isnum = self.checktwo(other) |
|---|
| 366 | the_clone = self.selective_copy() |
|---|
| 367 | for i in range(self.nfin): |
|---|
| 368 | for j in range(self.nvin): |
|---|
| 369 | for t in range(self.nplott): |
|---|
| 370 | for z in range(self.nplotz): |
|---|
| 371 | for y in range(self.nploty): |
|---|
| 372 | for x in range(self.nplotx): |
|---|
| 373 | obj = the_clone.request[i][j][t][z][y][x] |
|---|
| 374 | obj_ref = self.request[i][j][t][z][y][x] |
|---|
| 375 | if not isnum: |
|---|
| 376 | ope = other.request[i][j][t][z][y][x].field |
|---|
| 377 | if obj_ref.field.shape != ope.shape: |
|---|
| 378 | print "!! ERROR !! The two fields for operation do not have the same shape.",self.field.shape,other.field.shape |
|---|
| 379 | exit() |
|---|
| 380 | else: |
|---|
| 381 | ope = other |
|---|
| 382 | goal = self.vargoal[j] + self.filegoal[i] |
|---|
| 383 | if ("vector" in goal) or ("contour" in goal): |
|---|
| 384 | if self.verbose: print "!! WARNING !! No operation was made on contours and vectors. This can be debatted actually." |
|---|
| 385 | obj.field = obj_ref.field |
|---|
| 386 | else: |
|---|
| 387 | obj.field = obj_ref.field * ope |
|---|
| 388 | return the_clone |
|---|
| 389 | |
|---|
| 390 | # define the operation / on two objects. or with an int/float. |
|---|
| 391 | # ... with selective_copy the self object is not modified. |
|---|
| 392 | def __div__(self,other): |
|---|
| 393 | isnum = self.checktwo(other) |
|---|
| 394 | the_clone = self.selective_copy() |
|---|
| 395 | for i in range(self.nfin): |
|---|
| 396 | for j in range(self.nvin): |
|---|
| 397 | for t in range(self.nplott): |
|---|
| 398 | for z in range(self.nplotz): |
|---|
| 399 | for y in range(self.nploty): |
|---|
| 400 | for x in range(self.nplotx): |
|---|
| 401 | obj = the_clone.request[i][j][t][z][y][x] |
|---|
| 402 | obj_ref = self.request[i][j][t][z][y][x] |
|---|
| 403 | if not isnum: |
|---|
| 404 | ope = other.request[i][j][t][z][y][x].field |
|---|
| 405 | if obj_ref.field.shape != ope.shape: |
|---|
| 406 | print "!! ERROR !! The two fields for operation do not have the same shape.",self.field.shape,other.field.shape |
|---|
| 407 | exit() |
|---|
| 408 | else: |
|---|
| 409 | ope = other |
|---|
| 410 | goal = self.vargoal[j] + self.filegoal[i] |
|---|
| 411 | if ("vector" in goal) or ("contour" in goal): |
|---|
| 412 | if self.verbose: print "!! WARNING !! No operation was made on contours and vectors. This can be debatted actually." |
|---|
| 413 | obj.field = obj_ref.field |
|---|
| 414 | else: |
|---|
| 415 | obj.field = obj_ref.field / ope |
|---|
| 416 | return the_clone |
|---|
| 417 | |
|---|
| 418 | # define the reverse operation float/int + object |
|---|
| 419 | def __radd__(self,other): |
|---|
| 420 | isnum = self.checktwo(other) |
|---|
| 421 | if not isnum: print "!! ERROR !! Operand should be a number" ; exit() |
|---|
| 422 | return self.__add__(other) |
|---|
| 423 | |
|---|
| 424 | # define the reverse operation float/int - object |
|---|
| 425 | def __rsub__(self,other): |
|---|
| 426 | isnum = self.checktwo(other) |
|---|
| 427 | if not isnum: print "!! ERROR !! Operand should be a number" ; exit() |
|---|
| 428 | return self.__sub__(other) |
|---|
| 429 | |
|---|
| 430 | # define the reverse operation float/int * object |
|---|
| 431 | def __rmul__(self,other): |
|---|
| 432 | isnum = self.checktwo(other) |
|---|
| 433 | if not isnum: print "!! ERROR !! Operand should be a number" ; exit() |
|---|
| 434 | return self.__mul__(other) |
|---|
| 435 | |
|---|
| 436 | # define the reverse operation float/int / object |
|---|
| 437 | def __rdiv__(self,other): |
|---|
| 438 | isnum = self.checktwo(other) |
|---|
| 439 | if not isnum: print "!! ERROR !! Operand should be a number" ; exit() |
|---|
| 440 | return self.__div__(other) |
|---|
| 441 | |
|---|
| 442 | # define the operation ** on one object. |
|---|
| 443 | # ... with selective_copy the self object is not modified. |
|---|
| 444 | def __pow__(self,num): |
|---|
| 445 | the_clone = self.selective_copy() |
|---|
| 446 | if isinstance(num,int) or isinstance(num,float): |
|---|
| 447 | for i in range(self.nfin): |
|---|
| 448 | for j in range(self.nvin): |
|---|
| 449 | for t in range(self.nplott): |
|---|
| 450 | for z in range(self.nplotz): |
|---|
| 451 | for y in range(self.nploty): |
|---|
| 452 | for x in range(self.nplotx): |
|---|
| 453 | obj = the_clone.request[i][j][t][z][y][x] |
|---|
| 454 | obj_ref = self.request[i][j][t][z][y][x] |
|---|
| 455 | goal = self.vargoal[j] + self.filegoal[i] |
|---|
| 456 | if ("vector" in goal) or ("contour" in goal): |
|---|
| 457 | if self.verbose: print "!! WARNING !! No operation was made on contours and vectors. This can be debatted actually." |
|---|
| 458 | obj.field = obj_ref.field |
|---|
| 459 | else: |
|---|
| 460 | obj.field = obj_ref.field ** num |
|---|
| 461 | else: |
|---|
| 462 | print "!! ERROR !! To define a power, either an int or a float is needed." ; exit() |
|---|
| 463 | return the_clone |
|---|
| 464 | |
|---|
| 465 | ### TBD: reverse power? for exponentials? |
|---|
| 466 | |
|---|
| 467 | ############################################################################################## |
|---|
| 468 | # define method |
|---|
| 469 | # --------- |
|---|
| 470 | # ... (file and var are either one string or a vector of strings) |
|---|
| 471 | # ... the goal of define is to define a 2D array of onerequest() objects (see class below) |
|---|
| 472 | # given the number of file, var, x, y, z, t asked by the user |
|---|
| 473 | # ... objectives for file or var are given through filegoal and vargoal |
|---|
| 474 | # --> possible values: main contour vector |
|---|
| 475 | # --------- |
|---|
| 476 | # ... then onerequest() objects are being defined more precisely |
|---|
| 477 | # by getting index_x index_y index_z index_t |
|---|
| 478 | # and setting method_x method_y method_z method_t to either |
|---|
| 479 | # - "free" for free dimensions (plot dimensions) |
|---|
| 480 | # - "comp" for averages, max, min |
|---|
| 481 | # - "fixed" for fixed dimensions (possibly several i.e. multislice) |
|---|
| 482 | ############################################################################################## |
|---|
| 483 | def define(self): |
|---|
| 484 | self.printstatus() |
|---|
| 485 | # initial check and get dimensions |
|---|
| 486 | self.file = checktab(self.file,mess="file") |
|---|
| 487 | self.nfin = len(self.file) |
|---|
| 488 | if self.verbose: |
|---|
| 489 | for i in range(self.nfin): inspect(self.file[i]) |
|---|
| 490 | self.var = checktab(self.var,mess="var") |
|---|
| 491 | self.nvin = len(self.var) |
|---|
| 492 | # check goal tabs for files and variables |
|---|
| 493 | # ... default is to plot everything |
|---|
| 494 | if self.filegoal is None: self.filegoal = ["main"]*self.nfin |
|---|
| 495 | if self.vargoal is None: self.vargoal = ["main"]*self.nvin |
|---|
| 496 | self.filegoal = checktab(self.filegoal, mess="filegoal") |
|---|
| 497 | self.vargoal = checktab(self.vargoal, mess="vargoal") |
|---|
| 498 | if len(self.filegoal) != self.nfin: print "!! ERROR !! filegoal must be the same size as file." ; exit() |
|---|
| 499 | if len(self.vargoal) != self.nvin: print "!! ERROR !! vargoal must be the same size as var." ; exit() |
|---|
| 500 | # variables: initial check |
|---|
| 501 | self.x = checktab(self.x,mess="x",allownone=True,allownumber=True) |
|---|
| 502 | self.y = checktab(self.y,mess="y",allownone=True,allownumber=True) |
|---|
| 503 | self.z = checktab(self.z,mess="z",allownone=True,allownumber=True) |
|---|
| 504 | self.t = checktab(self.t,mess="t",allownone=True,allownumber=True) |
|---|
| 505 | # for the moment not var- nor file- dependent. |
|---|
| 506 | # but this could be the case. |
|---|
| 507 | sx = readslices(self.x) ; sy = readslices(self.y) |
|---|
| 508 | sz = readslices(self.z) ; st = readslices(self.t) |
|---|
| 509 | # get methods |
|---|
| 510 | mx = findmethod(sx) ; my = findmethod(sy) |
|---|
| 511 | mz = findmethod(sz) ; mt = findmethod(st) |
|---|
| 512 | # get number of plots to be done |
|---|
| 513 | if mx in ["fixed","comp"]: self.nplotx = sx.size/2 |
|---|
| 514 | else: self.nplotx = 1 |
|---|
| 515 | if my in ["fixed","comp"]: self.nploty = sy.size/2 |
|---|
| 516 | else: self.nploty = 1 |
|---|
| 517 | if mz in ["fixed","comp"]: self.nplotz = sz.size/2 |
|---|
| 518 | else: self.nplotz = 1 |
|---|
| 519 | if mt in ["fixed","comp"]: self.nplott = st.size/2 |
|---|
| 520 | else: self.nplott = 1 |
|---|
| 521 | if self.verbose: print "**** OK. Plots over x,y,z,t -->",self.nplotx,self.nploty,self.nplotz,self.nplott |
|---|
| 522 | # create the list of onerequest() objects |
|---|
| 523 | self.request = [[[[[[ \ |
|---|
| 524 | onerequest() \ |
|---|
| 525 | for x in range(self.nplotx)] for y in range(self.nploty)] \ |
|---|
| 526 | for z in range(self.nplotz)] for t in range(self.nplott)] \ |
|---|
| 527 | for j in range(self.nvin)] for i in range(self.nfin)] |
|---|
| 528 | # loop on onerequest() objects |
|---|
| 529 | for i in range(self.nfin): |
|---|
| 530 | for j in range(self.nvin): |
|---|
| 531 | for t in range(self.nplott): |
|---|
| 532 | for z in range(self.nplotz): |
|---|
| 533 | for y in range(self.nploty): |
|---|
| 534 | for x in range(self.nplotx): |
|---|
| 535 | obj = self.request[i][j][t][z][y][x] |
|---|
| 536 | # fill in names for files and variables |
|---|
| 537 | obj.verbose = self.verbose |
|---|
| 538 | obj.file = self.file[i] |
|---|
| 539 | obj.var = self.var[j] |
|---|
| 540 | # indicate the computation method |
|---|
| 541 | obj.compute = self.compute |
|---|
| 542 | # open the files (the same file might be opened several times but this is cheap) |
|---|
| 543 | obj.openfile() |
|---|
| 544 | ### get x,y,z,t dimensions from file |
|---|
| 545 | obj.getdim() |
|---|
| 546 | ### possible time axis change |
|---|
| 547 | obj.changetime = self.changetime |
|---|
| 548 | obj.performtimechange() |
|---|
| 549 | ### get methods |
|---|
| 550 | obj.method_x = mx ; obj.method_y = my |
|---|
| 551 | obj.method_z = mz ; obj.method_t = mt |
|---|
| 552 | ### get index |
|---|
| 553 | obj.getindextime(dalist=st,ind=t,stride=self.stridet) |
|---|
| 554 | obj.getindexvert(dalist=sz,ind=z,stride=self.stridez) |
|---|
| 555 | obj.getindexhori(dalistx=sx,dalisty=sy,indx=x,indy=y,stridex=self.stridex,stridey=self.stridey) |
|---|
| 556 | # change status |
|---|
| 557 | self.status = "defined" |
|---|
| 558 | return self |
|---|
| 559 | |
|---|
| 560 | ############################################################################################## |
|---|
| 561 | # retrieve method |
|---|
| 562 | # --> for each element onerequest() in the array, get field .var from .f file |
|---|
| 563 | # --> see below the onerequest() class: |
|---|
| 564 | # - only get what is needed for computing and plotting |
|---|
| 565 | # - averages etc... are computed here |
|---|
| 566 | # --> RESULT: each onerequest() object has now its attribute .field filled |
|---|
| 567 | # --> if one wants to perform operations on fields, this should be done after retrieve() |
|---|
| 568 | ############################################################################################## |
|---|
| 569 | def retrieve(self): |
|---|
| 570 | self.printstatus() |
|---|
| 571 | # check if things were done OK before |
|---|
| 572 | if self.status != "defined": print "!! ERROR !! Please use .define() to define your pp object." ; exit() |
|---|
| 573 | # create the list of allfield() objects |
|---|
| 574 | # --> so that the user can easily access values |
|---|
| 575 | self.allfield = [[[[[[ \ |
|---|
| 576 | [] \ |
|---|
| 577 | for x in range(self.nplotx)] for y in range(self.nploty)] \ |
|---|
| 578 | for z in range(self.nplotz)] for t in range(self.nplott)] \ |
|---|
| 579 | for j in range(self.nvin)] for i in range(self.nfin)] |
|---|
| 580 | ## first get fields |
|---|
| 581 | ## ... only what is needed is extracted from the files |
|---|
| 582 | ## ... and computations are performed |
|---|
| 583 | for i in range(self.nfin): |
|---|
| 584 | for j in range(self.nvin): |
|---|
| 585 | for t in range(self.nplott): |
|---|
| 586 | for z in range(self.nplotz): |
|---|
| 587 | for y in range(self.nploty): |
|---|
| 588 | for x in range(self.nplotx): |
|---|
| 589 | obj = self.request[i][j][t][z][y][x] |
|---|
| 590 | obj.getfield() |
|---|
| 591 | obj.computations() |
|---|
| 592 | self.allfield[i][j][t][z][y][x] = obj.field |
|---|
| 593 | # change status |
|---|
| 594 | self.status = "retrieved" |
|---|
| 595 | return self |
|---|
| 596 | |
|---|
| 597 | ########################################################## |
|---|
| 598 | # get: a shortcut method for the define + retrieve chain # |
|---|
| 599 | ########################################################## |
|---|
| 600 | def get(self): |
|---|
| 601 | self.define() |
|---|
| 602 | self.retrieve() |
|---|
| 603 | return self |
|---|
| 604 | |
|---|
| 605 | ######################################## |
|---|
| 606 | # smooth: smooth the field in 1D or 2D # |
|---|
| 607 | ######################################## |
|---|
| 608 | ## TBD: smooth not OK with masked array in the end of retrieve() |
|---|
| 609 | def smooth(self,window): |
|---|
| 610 | if self.verbose: |
|---|
| 611 | print "!! WARNING !! Performing a smoothing with a window size",window |
|---|
| 612 | print "!! WARNING !! To come back to unsmoothed file, use .get() again" |
|---|
| 613 | for i in range(self.nfin): |
|---|
| 614 | for j in range(self.nvin): |
|---|
| 615 | for t in range(self.nplott): |
|---|
| 616 | for z in range(self.nplotz): |
|---|
| 617 | for y in range(self.nploty): |
|---|
| 618 | for x in range(self.nplotx): |
|---|
| 619 | obj = self.request[i][j][t][z][y][x] |
|---|
| 620 | if obj.field.ndim == 1: |
|---|
| 621 | print "!! ERROR !! 1D smoothing not supported yet because reduces array sizes." |
|---|
| 622 | exit() |
|---|
| 623 | # TBD TBD TBD |
|---|
| 624 | #obj.field = ppcompute.smooth1d(obj.field,window=window) |
|---|
| 625 | elif obj.field.ndim == 2: |
|---|
| 626 | obj.field = ppcompute.smooth2d(obj.field,window=window) |
|---|
| 627 | |
|---|
| 628 | ############################################################################################## |
|---|
| 629 | # defineplot method |
|---|
| 630 | # --> defineplot first defines arrays of plot objects and set each of them |
|---|
| 631 | # ... simple looping except cases where goal is not main (e.g. contour or vector) |
|---|
| 632 | # --> principle: each onerequest() object gives birth to a subplot |
|---|
| 633 | # --> defineplot vs. makeplot: defining plot and actually plotting it are clearly separated |
|---|
| 634 | # --> THE KEY OUPUT OF defineplot IS AN ARRAY self.p OF PLOT OBJECTS |
|---|
| 635 | # optional arguments |
|---|
| 636 | # --> extraplot: to indicate a number of plots to be added afterwards (use self.plotin) |
|---|
| 637 | # --> loadfile: to use self.p from a previously saved file |
|---|
| 638 | ############################################################################################## |
|---|
| 639 | def defineplot(self,extraplot=0,loadfile=None): |
|---|
| 640 | # ----------------------------------------------------- |
|---|
| 641 | # LOAD MODE: load a self.p object. count plots from it. |
|---|
| 642 | # ----------------------------------------------------- |
|---|
| 643 | if loadfile is not None: |
|---|
| 644 | try: filehandler = open(loadfile, 'r') ; self.p = pickle.load(filehandler) |
|---|
| 645 | except IOError: print "!! ERROR !! Cannot find object file to load." ; exit() |
|---|
| 646 | self.status = "definedplot" ; self.plotin = None |
|---|
| 647 | self.nplot = len(self.p) ; self.howmanyplots = self.nplot |
|---|
| 648 | return |
|---|
| 649 | # ----------------------------------------------------- |
|---|
| 650 | # REGULAR MODE |
|---|
| 651 | # ----------------------------------------------------- |
|---|
| 652 | self.printstatus() |
|---|
| 653 | # check if things were done OK before |
|---|
| 654 | if self.status in ["init","defined"]: |
|---|
| 655 | print "!! ERROR !! Please use .retrieve() to get fields for plots with your pp object." ; exit() |
|---|
| 656 | # check self.plotin (an existing fig on which to add plots afterwards) |
|---|
| 657 | if self.plotin.__class__.__name__ == "pp": |
|---|
| 658 | if self.plotin.fig is None: |
|---|
| 659 | self.plotin = None # this is an additional security in case |
|---|
| 660 | # a pp object is given without figure opened yet. |
|---|
| 661 | elif self.plotin is not None: |
|---|
| 662 | print "!! ERROR !! plotin argument must be a pp object." ; exit() |
|---|
| 663 | # initialize the array of subplot objects |
|---|
| 664 | # either something new or attributes coming from plotin object |
|---|
| 665 | if self.plotin is None: self.p = [ ] |
|---|
| 666 | else: self.p = self.plotin.p |
|---|
| 667 | # create an array of subplot objects |
|---|
| 668 | # ... in theory the order of looping can be changed without any harm |
|---|
| 669 | # ... the only important thing is to keep i,j,t,z,y,x resp. for file,var,t,z,y,x |
|---|
| 670 | count = 0 |
|---|
| 671 | for i in range(self.nfin): |
|---|
| 672 | if self.filegoal[i] == "main": |
|---|
| 673 | for j in range(self.nvin): |
|---|
| 674 | if self.vargoal[j] == "main": |
|---|
| 675 | for t in range(self.nplott): |
|---|
| 676 | for z in range(self.nplotz): |
|---|
| 677 | for y in range(self.nploty): |
|---|
| 678 | for x in range(self.nplotx): |
|---|
| 679 | # look at dimension and append the right plot object |
|---|
| 680 | obj = self.request[i][j][t][z][y][x] |
|---|
| 681 | dp = obj.dimplot |
|---|
| 682 | if dp == 1 or self.forcedimplot == 1: plobj = ppplot.plot1d() |
|---|
| 683 | elif dp == 2 or self.forcedimplot == 2: plobj = ppplot.plot2d() |
|---|
| 684 | elif dp == 0: print "**** OK. VALUES VALUES VALUES",obj.field |
|---|
| 685 | else: print "!! ERROR !! 3D or 4D plots not supported" ; exit() |
|---|
| 686 | # load abscissa and ordinate in obj |
|---|
| 687 | obj.definecoord() |
|---|
| 688 | # we start to define things here before appending |
|---|
| 689 | # (convenient: could be overridden by the user before makeplot) |
|---|
| 690 | # ... the if loop is necessary so that we can loop above on the dp=0 case |
|---|
| 691 | if dp in [1,2]: |
|---|
| 692 | # and define what to do in plobj |
|---|
| 693 | plobj.invert = obj.invert_axes |
|---|
| 694 | plobj.swap = obj.swap_axes |
|---|
| 695 | # axis labels |
|---|
| 696 | plobj.xlabel = obj.absclab ; plobj.ylabel = obj.ordilab |
|---|
| 697 | # superpose or not (this is mostly for saving purpose) |
|---|
| 698 | plobj.superpose = self.superpose |
|---|
| 699 | # get title, colormaps, labels, etc.. from var |
|---|
| 700 | plobj.var = obj.var |
|---|
| 701 | plobj.define_from_var() |
|---|
| 702 | # generic 1D/2D: load field and coord in plot object |
|---|
| 703 | plobj.field = obj.field # field to be plotted |
|---|
| 704 | plobj.absc = obj.absc # abscissa (or longitude) |
|---|
| 705 | plobj.ordi = obj.ordi # ordinate (or latitude) |
|---|
| 706 | # -- useless in 1D but not used anyway |
|---|
| 707 | # specific 1D plot stuff |
|---|
| 708 | if dp == 1: |
|---|
| 709 | # -- a default label |
|---|
| 710 | plobj.label = "" |
|---|
| 711 | if self.nfin > 1: plobj.label = plobj.label + " file #"+str(i+1) |
|---|
| 712 | if self.nvin > 1: plobj.label = plobj.label + " var "+plobj.var |
|---|
| 713 | if self.nplott > 1: plobj.label = plobj.label + " t="+str(self.t[t]) |
|---|
| 714 | if self.nplotz > 1: plobj.label = plobj.label + " z="+str(self.z[z]) |
|---|
| 715 | if self.nploty > 1: plobj.label = plobj.label + " y="+str(self.y[y]) |
|---|
| 716 | if self.nplotx > 1: plobj.label = plobj.label + " x="+str(self.x[x]) |
|---|
| 717 | # specific 2d plot stuff |
|---|
| 718 | if dp == 2: |
|---|
| 719 | # -- light grey background for missing values |
|---|
| 720 | if type(plobj.field).__name__ in 'MaskedArray': plobj.axisbg = '0.75' |
|---|
| 721 | # -- define if it is a map or a plot |
|---|
| 722 | plobj.mapmode = ( obj.method_x+obj.method_y == "freefree" \ |
|---|
| 723 | and "grid points" not in obj.name_x \ |
|---|
| 724 | and not self.noproj ) |
|---|
| 725 | # possible user-defined plot settings shared by all plots |
|---|
| 726 | if self.div is not None: plobj.div = self.div |
|---|
| 727 | if self.xlabel is not None: plobj.xlabel = self.xlabel |
|---|
| 728 | if self.xcoeff is not None: plobj.xcoeff = self.xcoeff |
|---|
| 729 | if self.ylabel is not None: plobj.ylabel = self.ylabel |
|---|
| 730 | if self.ycoeff is not None: plobj.ycoeff = self.ycoeff |
|---|
| 731 | if self.title is not None: plobj.title = self.title |
|---|
| 732 | # -- 1D specific |
|---|
| 733 | if dp == 1: |
|---|
| 734 | if self.lstyle is not None: plobj.lstyle = self.lstyle |
|---|
| 735 | if self.marker is not None: plobj.marker = self.marker |
|---|
| 736 | if self.color is not None: plobj.color = self.color |
|---|
| 737 | if self.label is not None: plobj.label = self.label |
|---|
| 738 | # -- 2D specific |
|---|
| 739 | elif dp == 2: |
|---|
| 740 | if self.proj is not None and not self.noproj: plobj.proj = self.proj |
|---|
| 741 | if self.vmin is not None: plobj.vmin = self.vmin |
|---|
| 742 | if self.vmax is not None: plobj.vmax = self.vmax |
|---|
| 743 | if self.colorb is not None: plobj.colorb = self.colorb |
|---|
| 744 | # finally append plot object |
|---|
| 745 | self.p.append(plobj) |
|---|
| 746 | count = count + 1 |
|---|
| 747 | # self.nplot is number of plot to be defined in this call to defineplot() |
|---|
| 748 | # (because of self.plotin this might less than length of self.p) |
|---|
| 749 | self.nplot = count |
|---|
| 750 | # --- superimposed contours and vectors --- |
|---|
| 751 | # we have to start another loop because we need forward information |
|---|
| 752 | # TBD: there is probably a more flexible way to do that |
|---|
| 753 | count = 0 |
|---|
| 754 | for i in range(self.nfin): |
|---|
| 755 | for j in range(self.nvin): |
|---|
| 756 | for t in range(self.nplott): |
|---|
| 757 | for z in range(self.nplotz): |
|---|
| 758 | for y in range(self.nploty): |
|---|
| 759 | for x in range(self.nplotx): |
|---|
| 760 | goal = self.vargoal[j] + self.filegoal[i] |
|---|
| 761 | obj = self.request[i][j][t][z][y][x] |
|---|
| 762 | if "mainmain" in goal and obj.dimplot == 2: |
|---|
| 763 | # the plot object we consider in the loop |
|---|
| 764 | pl = self.p[count] |
|---|
| 765 | # -- see if there is a contour requested... |
|---|
| 766 | # (we use try because we might be at the end of the list) |
|---|
| 767 | found = 0 |
|---|
| 768 | try: condvar = self.vargoal[j+1] |
|---|
| 769 | except: condvar = "itisok" |
|---|
| 770 | try: condfile = self.filegoal[i+1] |
|---|
| 771 | except: condfile = "itisok" |
|---|
| 772 | # ... get contour |
|---|
| 773 | ########################################## |
|---|
| 774 | # NB: contour is expected to be right after main otherwise it is not displayed |
|---|
| 775 | ########################################## |
|---|
| 776 | if condvar == "contour": |
|---|
| 777 | plobj.addcontour = self.request[i][j+1][t][z][y][x].field ; found += 1 |
|---|
| 778 | if condfile == "contour": |
|---|
| 779 | plobj.addcontour = self.request[i+1][j][t][z][y][x].field ; found += 1 |
|---|
| 780 | # see if there is a vector requested... |
|---|
| 781 | # (we use try because we might be at the end of the list) |
|---|
| 782 | try: condvar = self.vargoal[j+found+1]+self.vargoal[j+found+2] |
|---|
| 783 | except: condvar = "itisok" |
|---|
| 784 | try: condfile = self.filegoal[i+found+1]+self.filegoal[i+found+2] |
|---|
| 785 | except: condfile = "itisok" |
|---|
| 786 | # ... get vector and go directly to the next iteration |
|---|
| 787 | # (in some cases we would do this twice but this is cheap) |
|---|
| 788 | if "vector" in condvar: |
|---|
| 789 | plobj.addvecx = self.request[i][j+found+1][t][z][y][x].field |
|---|
| 790 | plobj.addvecy = self.request[i][j+found+2][t][z][y][x].field |
|---|
| 791 | if "vector" in condfile: |
|---|
| 792 | plobj.addvecx = self.request[i+found+1][j][t][z][y][x].field |
|---|
| 793 | plobj.addvecy = self.request[i+found+2][j][t][z][y][x].field |
|---|
| 794 | count = count + 1 |
|---|
| 795 | # COUNT PLOTS. if 0 just exit. |
|---|
| 796 | # self.howmanyplots is self.nplot + possible extraplots |
|---|
| 797 | self.howmanyplots = self.nplot + extraplot |
|---|
| 798 | if self.howmanyplots > 0: |
|---|
| 799 | if self.verbose: print "**** OK. expect %i plots" % (self.howmanyplots) |
|---|
| 800 | else: |
|---|
| 801 | exit() # because this means that we only had 0D values ! |
|---|
| 802 | # final status |
|---|
| 803 | self.status = "definedplot" |
|---|
| 804 | return self |
|---|
| 805 | |
|---|
| 806 | ############################################################################################## |
|---|
| 807 | # makeplot method |
|---|
| 808 | # --> after defineplot and before makeplot, user-defined plot settings can be easily given |
|---|
| 809 | # simply by modifying the attributes of each elements of self.p |
|---|
| 810 | # --> to change only one plot setting, no need to call defineplot again before makeplot |
|---|
| 811 | # --> in the end, the array self.p of plot objects is saved for easy and convenient replotting |
|---|
| 812 | # --> see practical examples in the folder 'examples' |
|---|
| 813 | ############################################################################################## |
|---|
| 814 | def makeplot(self): |
|---|
| 815 | self.printstatus() |
|---|
| 816 | # a few initial operations |
|---|
| 817 | # ------------------------ |
|---|
| 818 | if "definedplot" not in self.status: |
|---|
| 819 | print "!! ERROR !! Please use .defineplot() before .makeplot() can be used with your pp object." ; exit() |
|---|
| 820 | # open a figure and define subplots |
|---|
| 821 | # --------------------------------- |
|---|
| 822 | if self.plotin is None: |
|---|
| 823 | # start from scratch |
|---|
| 824 | self.fig = mpl.figure(figsize=(16,8)) |
|---|
| 825 | self.subv,self.subh = ppplot.definesubplot(self.howmanyplots,self.fig) |
|---|
| 826 | self.n = 0 |
|---|
| 827 | ## adapted space for labels etc |
|---|
| 828 | ## ... except for ortho because there is no label anyway |
|---|
| 829 | self.customplot = self.p[0].field.ndim == 2 \ |
|---|
| 830 | and self.p[0].mapmode == True \ |
|---|
| 831 | and self.p[0].proj not in ["ortho"] |
|---|
| 832 | if self.customplot: |
|---|
| 833 | margin = 0.07 |
|---|
| 834 | self.fig.subplots_adjust(left=margin,right=1-margin,bottom=margin,top=1-margin) |
|---|
| 835 | else: |
|---|
| 836 | # start from an existing figure. |
|---|
| 837 | # extraplot must have been set in the call to the previous figure. |
|---|
| 838 | self.fig = self.plotin.fig |
|---|
| 839 | self.subv,self.subh = self.plotin.subv,self.plotin.subh |
|---|
| 840 | self.n = self.plotin.n |
|---|
| 841 | self.howmanyplots = self.plotin.howmanyplots |
|---|
| 842 | self.customplot = self.plotin.customplot |
|---|
| 843 | # LOOP on all subplots |
|---|
| 844 | # NB: cannot use 'for pl in self.p' if self.plotin not None |
|---|
| 845 | # -------------------- |
|---|
| 846 | for count in range(self.nplot): |
|---|
| 847 | # the plot object we consider in the loop |
|---|
| 848 | pl = self.p[self.n] |
|---|
| 849 | # before making the plot, create a subplot. the first one is numbered 1 not 0. |
|---|
| 850 | # ... if pl.superpose, we use only one and only figure |
|---|
| 851 | # ... (and we have to be careful with not doing things several times) |
|---|
| 852 | if pl.superpose: |
|---|
| 853 | if self.n == 0: |
|---|
| 854 | self.fig.add_subplot(1,1,1,axisbg=pl.axisbg) # define one subplot (still needed for user-defined font sizes) |
|---|
| 855 | sav = pl.xlabel,pl.ylabel,pl.xcoeff,pl.ycoeff,pl.title,pl.swaplab # save titles and labels |
|---|
| 856 | else: |
|---|
| 857 | pl.invert = False ; pl.lstyle = None # don't invert again axis |
|---|
| 858 | # set saved titles and labels |
|---|
| 859 | if self.plotin is None: |
|---|
| 860 | pl.xlabel,pl.ylabel,pl.xcoeff,pl.ycoeff,pl.title,pl.swaplab = sav |
|---|
| 861 | else: |
|---|
| 862 | prev_plot = self.plotin.p[self.n-1] |
|---|
| 863 | pl.xlabel = prev_plot.xlabel |
|---|
| 864 | pl.ylabel = prev_plot.ylabel |
|---|
| 865 | pl.xcoeff = prev_plot.xcoeff |
|---|
| 866 | pl.ycoeff = prev_plot.ycoeff |
|---|
| 867 | pl.title = prev_plot.title |
|---|
| 868 | pl.swaplab = prev_plot.swaplab |
|---|
| 869 | else: |
|---|
| 870 | self.fig.add_subplot(self.subv,self.subh,self.n+1,axisbg=pl.axisbg) |
|---|
| 871 | if self.verbose: print "**** Done subplot %i / %i " %( self.n+1,self.howmanyplots ) |
|---|
| 872 | # finally make the plot |
|---|
| 873 | pl.make() |
|---|
| 874 | # increment plot count (and propagate this in plotin) |
|---|
| 875 | self.n = self.n+1 |
|---|
| 876 | if self.plotin is not None: self.plotin.n = self.n |
|---|
| 877 | # once completed show the plot (cannot show intermediate plotin) |
|---|
| 878 | # ... added a fix (customplot=True) for the label problem in basemap |
|---|
| 879 | print "**** PPCLASS. Done step: makeplot" |
|---|
| 880 | if (self.n == self.howmanyplots): |
|---|
| 881 | ppplot.save(mode=self.out,filename=self.filename,folder=self.folder,custom=self.customplot,includedate=self.includedate) |
|---|
| 882 | mpl.close() |
|---|
| 883 | # SAVE A PICKLE FILE WITH THE self.p ARRAY OF OBJECTS |
|---|
| 884 | if self.verbose: print "**** Saving session in "+self.filename + ".ppobj" |
|---|
| 885 | savfile = self.folder + "/" + self.filename + ".ppobj" |
|---|
| 886 | try: |
|---|
| 887 | filehandler = open(savfile, 'w') |
|---|
| 888 | pickle.dump(self.p, filehandler) |
|---|
| 889 | except IOError: |
|---|
| 890 | print "!! WARNING !! Saved object file not written. Probably do not have permission to write here." |
|---|
| 891 | return self |
|---|
| 892 | |
|---|
| 893 | ########################################################### |
|---|
| 894 | # plot: a shortcut method for the defineplot + plot chain # |
|---|
| 895 | ########################################################### |
|---|
| 896 | def plot(self,extraplot=0): |
|---|
| 897 | self.defineplot(extraplot=extraplot) |
|---|
| 898 | self.makeplot() |
|---|
| 899 | return self |
|---|
| 900 | |
|---|
| 901 | ####################################################### |
|---|
| 902 | # getplot: a shortcut method for the get + plot chain # |
|---|
| 903 | ####################################################### |
|---|
| 904 | def getplot(self,extraplot=0): |
|---|
| 905 | self.get() |
|---|
| 906 | self.plot(extraplot=extraplot) |
|---|
| 907 | return self |
|---|
| 908 | |
|---|
| 909 | ################################################################### |
|---|
| 910 | # getdefineplot: a shortcut method for the get + defineplot chain # |
|---|
| 911 | ################################################################### |
|---|
| 912 | def getdefineplot(self,extraplot=0): |
|---|
| 913 | self.get() |
|---|
| 914 | self.defineplot(extraplot=extraplot) |
|---|
| 915 | return self |
|---|
| 916 | |
|---|
| 917 | ############################################################## |
|---|
| 918 | # f: operation on two pp objects being on status 'definedplot' |
|---|
| 919 | # this allows for one field being function of another one |
|---|
| 920 | # e.g. u.f(v) means u will be displayed as a function of v |
|---|
| 921 | # ... no need to do defineplot after u.f(v), makeplot directly |
|---|
| 922 | ############################################################## |
|---|
| 923 | def f(self,other): |
|---|
| 924 | # preamble: for this operation to work, defineplot() must have been done |
|---|
| 925 | if self.status != "definedplot": |
|---|
| 926 | if self.verbose: print "!! WARNING !! performing defineplot on operand" |
|---|
| 927 | self.defineplot() |
|---|
| 928 | if other.status != "definedplot": |
|---|
| 929 | if self.verbose: print "!! WARNING !! performing defineplot on operand" |
|---|
| 930 | other.defineplot() |
|---|
| 931 | # check total number of plots |
|---|
| 932 | if self.howmanyplots != other.howmanyplots: |
|---|
| 933 | print "!! ERROR !! The two operands do not have the same number of subplots." |
|---|
| 934 | exit() |
|---|
| 935 | # and now operation. |
|---|
| 936 | count = 0 |
|---|
| 937 | while count < self.howmanyplots: |
|---|
| 938 | sobj = self.p[count] ; oobj = other.p[count] |
|---|
| 939 | if sobj.field.ndim !=1 or oobj.field.ndim !=1: |
|---|
| 940 | if self.verbose: print "!! WARNING !! Flattening arrays because more than one-dimensional." |
|---|
| 941 | sobj.field = np.ravel(sobj.field) |
|---|
| 942 | oobj.field = np.ravel(oobj.field) |
|---|
| 943 | sobj.absc = oobj.field |
|---|
| 944 | sobj.xlabel = oobj.ylabel |
|---|
| 945 | if sobj.absc.size > sobj.field.size: |
|---|
| 946 | if self.verbose: |
|---|
| 947 | print "!! WARNING !! Trying to define y=f(x) with x and y not at the same size.",sobj.absc.size,sobj.field.size |
|---|
| 948 | print "!! WARNING !! Modifying x to fit y size but please check." |
|---|
| 949 | sobj.absc = sobj.absc[0:sobj.field.size] |
|---|
| 950 | count = count + 1 |
|---|
| 951 | return self |
|---|
| 952 | |
|---|
| 953 | ########################################################### |
|---|
| 954 | # copyopt: get options from e.g. a parser |
|---|
| 955 | # ... allow for simple scripting and user-defined settings |
|---|
| 956 | # ... must be called between defineplot and makeplot |
|---|
| 957 | # REQUIRED: attributes of opt must be the same as in the pp object |
|---|
| 958 | ########################################################### |
|---|
| 959 | def getopt(self,opt): |
|---|
| 960 | # -- if only one, or less than the number of plots --> we take the first one |
|---|
| 961 | # -- if as many as number of plots --> OK, each plot has its own setting |
|---|
| 962 | # (except a few cases such as trans) |
|---|
| 963 | for iii in range(self.howmanyplots): |
|---|
| 964 | ### |
|---|
| 965 | try: self.p[iii].trans = opt.trans |
|---|
| 966 | except: pass |
|---|
| 967 | ### |
|---|
| 968 | try: self.p[iii].div = opt.div |
|---|
| 969 | except: pass |
|---|
| 970 | ### |
|---|
| 971 | try: self.p[iii].logy = opt.logy |
|---|
| 972 | except: pass |
|---|
| 973 | ### |
|---|
| 974 | try: self.p[iii].colorb = opt.colorb[iii] |
|---|
| 975 | except: |
|---|
| 976 | try: self.p[iii].colorb = opt.colorb[0] |
|---|
| 977 | except: pass |
|---|
| 978 | ### |
|---|
| 979 | try: self.p[iii].title = opt.title[iii] |
|---|
| 980 | except: |
|---|
| 981 | try: self.p[iii].title = opt.title[0] |
|---|
| 982 | except: pass |
|---|
| 983 | ### |
|---|
| 984 | try: self.p[iii].xlabel = opt.xlabel[iii] |
|---|
| 985 | except: |
|---|
| 986 | try: self.p[iii].xlabel = opt.xlabel[0] |
|---|
| 987 | except: pass |
|---|
| 988 | ### |
|---|
| 989 | try: self.p[iii].ylabel = opt.ylabel[iii] |
|---|
| 990 | except: |
|---|
| 991 | try: self.p[iii].ylabel = opt.ylabel[0] |
|---|
| 992 | except: pass |
|---|
| 993 | ### |
|---|
| 994 | try: self.p[iii].lstyle = opt.lstyle[iii] |
|---|
| 995 | except: |
|---|
| 996 | try: self.p[iii].lstyle = opt.lstyle[0] |
|---|
| 997 | except: pass |
|---|
| 998 | ### |
|---|
| 999 | try: self.p[iii].color = opt.color[iii] |
|---|
| 1000 | except: |
|---|
| 1001 | try: self.p[iii].color = opt.color[0] |
|---|
| 1002 | except: pass |
|---|
| 1003 | ### |
|---|
| 1004 | try: self.p[iii].marker = opt.marker[iii] |
|---|
| 1005 | except: |
|---|
| 1006 | try: self.p[iii].marker = opt.marker[0] |
|---|
| 1007 | except: pass |
|---|
| 1008 | ### |
|---|
| 1009 | try: self.p[iii].label = opt.label[iii] |
|---|
| 1010 | except: |
|---|
| 1011 | try: self.p[iii].label = opt.label[0] |
|---|
| 1012 | except: pass |
|---|
| 1013 | ### |
|---|
| 1014 | try: self.p[iii].proj = opt.proj[iii] |
|---|
| 1015 | except: |
|---|
| 1016 | try: self.p[iii].proj = opt.proj[0] |
|---|
| 1017 | except: pass |
|---|
| 1018 | ### |
|---|
| 1019 | try: self.p[iii].back = opt.back[iii] |
|---|
| 1020 | except: |
|---|
| 1021 | try: self.p[iii].back = opt.back[0] |
|---|
| 1022 | except: pass |
|---|
| 1023 | ### |
|---|
| 1024 | try: self.p[iii].area = opt.area[iii] |
|---|
| 1025 | except: |
|---|
| 1026 | try: self.p[iii].area = opt.area[0] |
|---|
| 1027 | except: pass |
|---|
| 1028 | ### |
|---|
| 1029 | try: self.p[iii].blon = opt.blon[iii] |
|---|
| 1030 | except: |
|---|
| 1031 | try: self.p[iii].blon = opt.blon[0] |
|---|
| 1032 | except: pass |
|---|
| 1033 | ### |
|---|
| 1034 | try: self.p[iii].blat = opt.blat[iii] |
|---|
| 1035 | except: |
|---|
| 1036 | try: self.p[iii].blat = opt.blat[0] |
|---|
| 1037 | except: pass |
|---|
| 1038 | ### |
|---|
| 1039 | try: self.p[iii].vmin = opt.vmin[iii] |
|---|
| 1040 | except: |
|---|
| 1041 | try: self.p[iii].vmin = opt.vmin[0] |
|---|
| 1042 | except: pass |
|---|
| 1043 | ### |
|---|
| 1044 | try: self.p[iii].vmax = opt.vmax[iii] |
|---|
| 1045 | except: |
|---|
| 1046 | try: self.p[iii].vmax = opt.vmax[0] |
|---|
| 1047 | except: pass |
|---|
| 1048 | |
|---|
| 1049 | ########################################################## |
|---|
| 1050 | ### THE ONEREQUEST SUBOBJECT TO PP (ON WHICH IT LOOPS) ### |
|---|
| 1051 | ########################################################## |
|---|
| 1052 | class onerequest(): |
|---|
| 1053 | |
|---|
| 1054 | # default settings. mostly initialized to diagnose problem, except dimplot, nplot, verbose, swap_axes, invert_axes |
|---|
| 1055 | # ------------------------------- |
|---|
| 1056 | def __init__(self): |
|---|
| 1057 | self.file = '!! file: I am not set, damned !!' |
|---|
| 1058 | self.f = None |
|---|
| 1059 | self.dim = None |
|---|
| 1060 | self.var = '!! var: I am not set, damned !!' |
|---|
| 1061 | self.index_x = [] ; self.index_y = [] ; self.index_z = [] ; self.index_t = [] |
|---|
| 1062 | self.index_x2d = [] ; self.index_y2d = [] |
|---|
| 1063 | self.method_x = '!! method_x: I am not set, damned !!' |
|---|
| 1064 | self.method_y = '!! method_y: I am not set, damned !!' |
|---|
| 1065 | self.method_z = '!! method_z: I am not set, damned !!' |
|---|
| 1066 | self.method_t = '!! method_t: I am not set, damned !!' |
|---|
| 1067 | self.field = None |
|---|
| 1068 | self.name_x = None ; self.name_y = None ; self.name_z = None ; self.name_t = None |
|---|
| 1069 | self.dim_x = None ; self.dim_y = None ; self.dim_z = None ; self.dim_t = None |
|---|
| 1070 | self.field_x = None ; self.field_y = None ; self.field_z = None ; self.field_t = None |
|---|
| 1071 | self.dimplot = 0 |
|---|
| 1072 | self.nplot = 1 |
|---|
| 1073 | self.absc = None ; self.ordi = None ; self.absclab = None ; self.ordilab = None |
|---|
| 1074 | self.verbose = True |
|---|
| 1075 | self.swap_axes = False ; self.invert_axes = False |
|---|
| 1076 | self.compute = None |
|---|
| 1077 | self.changetime = None |
|---|
| 1078 | |
|---|
| 1079 | # open a file. for now it is netcdf. TBD for other formats. |
|---|
| 1080 | # check that self.var is inside. |
|---|
| 1081 | # ------------------------------- |
|---|
| 1082 | def openfile(self): |
|---|
| 1083 | if not os.path.exists(self.file): print '!! ERROR !! I could not find the following file: '+self.file ; exit() |
|---|
| 1084 | if not os.path.isfile(self.file): print '!! ERROR !! This does not appear to be a file: '+self.file ; exit() |
|---|
| 1085 | self.f = netCDF4.Dataset(self.file) |
|---|
| 1086 | if self.verbose: print "**** OK. Opened file "+self.file |
|---|
| 1087 | if self.var not in self.f.variables.keys(): |
|---|
| 1088 | print '!! ERROR !! File '+self.file+' does not contain variable: '+self.var |
|---|
| 1089 | print '..... try instead with ',self.f.variables.keys() ; exit() |
|---|
| 1090 | |
|---|
| 1091 | # copy attributes from another existing object |
|---|
| 1092 | # -------------------------------------------- |
|---|
| 1093 | def copy(self,source): |
|---|
| 1094 | for k, v in vars(source).items(): |
|---|
| 1095 | setattr(self,k,v) |
|---|
| 1096 | |
|---|
| 1097 | # get x,y,z,t dimensions from NETCDF file |
|---|
| 1098 | # TBD: user could request for a specific altitude dimension |
|---|
| 1099 | # TBD: staggered variables could request specific dimensions |
|---|
| 1100 | # ------------------------------- |
|---|
| 1101 | def getdim(self): |
|---|
| 1102 | # GET SIZES OF EACH DIMENSION |
|---|
| 1103 | if self.verbose: print "**** OK. Found variable "+self.var |
|---|
| 1104 | shape = self.f.variables[self.var].shape |
|---|
| 1105 | self.dim = len(shape) |
|---|
| 1106 | if self.dim == 1: |
|---|
| 1107 | if self.verbose: print "**** OK. 1D field. I assume this varies with time." |
|---|
| 1108 | self.dim_x = 1 ; self.dim_y = 1 ; self.dim_z = 1 ; self.dim_t = shape[0] |
|---|
| 1109 | elif self.dim == 2: |
|---|
| 1110 | if self.verbose: print "**** OK. 2D field. I assume this is not-time-varying lat-lon map." |
|---|
| 1111 | self.dim_x = shape[1] ; self.dim_y = shape[0] ; self.dim_z = 1 ; self.dim_t = 1 |
|---|
| 1112 | elif self.dim == 3: |
|---|
| 1113 | if self.verbose: print "**** OK. 3D field. I assume this is time-varying lat-lon map." |
|---|
| 1114 | self.dim_x = shape[2] ; self.dim_y = shape[1] ; self.dim_z = 1 ; self.dim_t = shape[0] |
|---|
| 1115 | elif self.dim == 4: |
|---|
| 1116 | if self.verbose: print "**** OK. 4D field." |
|---|
| 1117 | self.dim_x = shape[3] ; self.dim_y = shape[2] ; self.dim_z = shape[1] ; self.dim_t = shape[0] |
|---|
| 1118 | # LONGITUDE. Try preset fields. If not present set grid points axis. |
|---|
| 1119 | self.name_x = "nothing" |
|---|
| 1120 | for c in glob_listx: |
|---|
| 1121 | if c in self.f.variables.keys(): |
|---|
| 1122 | self.name_x = c |
|---|
| 1123 | if self.name_x == "nothing": |
|---|
| 1124 | self.field_x = np.array(range(self.dim_x)) |
|---|
| 1125 | self.name_x = "x grid points" |
|---|
| 1126 | else: |
|---|
| 1127 | self.field_x = self.f.variables[self.name_x] |
|---|
| 1128 | # LATITUDE. Try preset fields. If not present set grid points axis. |
|---|
| 1129 | self.name_y = "nothing" |
|---|
| 1130 | for c in glob_listy: |
|---|
| 1131 | if c in self.f.variables.keys(): |
|---|
| 1132 | self.name_y = c |
|---|
| 1133 | if self.name_y == "nothing": |
|---|
| 1134 | self.field_y = np.array(range(self.dim_y)) |
|---|
| 1135 | self.name_y = "y grid points" |
|---|
| 1136 | else: |
|---|
| 1137 | self.field_y = self.f.variables[self.name_y] |
|---|
| 1138 | # ensure that lon and lat are 2D fields |
|---|
| 1139 | # 1. simple 1D case (not time-varying) |
|---|
| 1140 | if len(self.field_x.shape)*len(self.field_y.shape) == 1: |
|---|
| 1141 | if self.verbose: print "**** OK. recasting lon and lat as 2D fields." |
|---|
| 1142 | [self.field_x,self.field_y] = np.meshgrid(self.field_x,self.field_y) |
|---|
| 1143 | # 2. complex 3D case (time-varying, actually just copied over time axis) |
|---|
| 1144 | elif len(self.field_x.shape)*len(self.field_y.shape) == 9: |
|---|
| 1145 | if self.verbose: print "**** OK. reducing lon and lat as 2D fields. get rid of time." |
|---|
| 1146 | self.field_x = self.field_x[0,:,:] |
|---|
| 1147 | self.field_y = self.field_y[0,:,:] |
|---|
| 1148 | # if xy axis are apparently undefined, set 2D grid points axis. |
|---|
| 1149 | if "grid points" not in self.name_x: |
|---|
| 1150 | if self.field_x.all() == self.field_x[0,0]: |
|---|
| 1151 | print "!! WARNING !! xy axis look undefined. creating non-dummy ones." |
|---|
| 1152 | self.field_x = np.array(range(self.dim_x)) ; self.name_x = "x grid points" |
|---|
| 1153 | self.field_y = np.array(range(self.dim_y)) ; self.name_y = "y grid points" |
|---|
| 1154 | [self.field_x,self.field_y] = np.meshgrid(self.field_x,self.field_y) |
|---|
| 1155 | if self.dim_x > 1: |
|---|
| 1156 | if self.verbose: print "**** OK. x axis %4.0f values [%5.1f,%5.1f]" % (self.dim_x,self.field_x.min(),self.field_x.max()) |
|---|
| 1157 | if self.dim_y > 1: |
|---|
| 1158 | if self.verbose: print "**** OK. y axis %4.0f values [%5.1f,%5.1f]" % (self.dim_y,self.field_y.min(),self.field_y.max()) |
|---|
| 1159 | # ALTITUDE. Try preset fields. If not present set grid points axis. |
|---|
| 1160 | # WARNING: how do we do if several are available? |
|---|
| 1161 | self.name_z = "nothing" |
|---|
| 1162 | for c in glob_listz: |
|---|
| 1163 | if c in self.f.variables.keys(): |
|---|
| 1164 | self.name_z = c |
|---|
| 1165 | if self.name_z == "nothing": |
|---|
| 1166 | self.field_z = np.array(range(self.dim_z)) |
|---|
| 1167 | self.name_z = "z grid points" |
|---|
| 1168 | else: |
|---|
| 1169 | self.field_z = self.f.variables[self.name_z][:] # specify dimension |
|---|
| 1170 | # TBD: have to check that this is not a 3D field |
|---|
| 1171 | if self.dim_z > 1: |
|---|
| 1172 | if self.verbose: print "**** OK. z axis %4.0f values [%5.1f,%5.1f]" % (self.dim_z,self.field_z.min(),self.field_z.max()) |
|---|
| 1173 | # TIME. Try preset fields. |
|---|
| 1174 | self.name_t = "nothing" |
|---|
| 1175 | for c in glob_listt: |
|---|
| 1176 | if c in self.f.variables.keys(): |
|---|
| 1177 | self.name_t = c |
|---|
| 1178 | if self.verbose: print "**** OK. Found time variable: ",c |
|---|
| 1179 | try: |
|---|
| 1180 | # speed up: only get first value, last one. |
|---|
| 1181 | tabtime = self.f.variables[self.name_t] |
|---|
| 1182 | # (consider the case where tabtime is not dim 1) |
|---|
| 1183 | # (time is most often the first dimension) |
|---|
| 1184 | if tabtime.ndim == 2: tabtime = tabtime[:,0] |
|---|
| 1185 | elif tabtime.ndim == 3: tabtime = tabtime[:,0,0] |
|---|
| 1186 | elif tabtime.ndim == 4: tabtime = tabtime[:,0,0,0] |
|---|
| 1187 | # (now proceed) |
|---|
| 1188 | dafirst = tabtime[0] |
|---|
| 1189 | if self.dim_t == 1: |
|---|
| 1190 | self.field_t = np.array([dafirst]) |
|---|
| 1191 | else: |
|---|
| 1192 | daint = tabtime[1] - dafirst |
|---|
| 1193 | dalast = dafirst + (self.dim_t-1)*daint |
|---|
| 1194 | if dalast != tabtime[self.dim_t-1] and self.verbose: |
|---|
| 1195 | print "!! WARNING !! Time axis has been recast to be monotonic",dalast,tabtime[self.dim_t-1] |
|---|
| 1196 | self.field_t = np.linspace(dafirst,dalast,num=self.dim_t) |
|---|
| 1197 | except: |
|---|
| 1198 | # ... or if a problem encountered, define a simple time axis |
|---|
| 1199 | self.field_t = np.array(range(self.dim_t)) |
|---|
| 1200 | self.name_t = "t grid points" |
|---|
| 1201 | if self.dim_t > 1: |
|---|
| 1202 | if self.verbose: print "**** OK. t axis %4.0f values [%5.1f,%5.1f]" % (self.dim_t,self.field_t.min(),self.field_t.max()) |
|---|
| 1203 | |
|---|
| 1204 | # change time axis |
|---|
| 1205 | # ... add your options here! |
|---|
| 1206 | # -------------------------- |
|---|
| 1207 | def performtimechange(self): |
|---|
| 1208 | if self.changetime is not None \ |
|---|
| 1209 | and self.name_t != "t grid points": |
|---|
| 1210 | if self.verbose: print "**** OK. Converting time axis:",self.changetime |
|---|
| 1211 | if self.changetime == "mars_sol2ls": |
|---|
| 1212 | self.field_t = ppcompute.mars_sol2ls(self.field_t) |
|---|
| 1213 | else: |
|---|
| 1214 | print "!! WARNING !! This time change is not implemented. Nothing is done." |
|---|
| 1215 | |
|---|
| 1216 | # get list of index to be retrieved for time axis |
|---|
| 1217 | ### TBD: il faudrait ne prendre que les indices qui correspondent a l interieur d un plot (dans all) |
|---|
| 1218 | # ------------------------------- |
|---|
| 1219 | def getindextime(self,dalist=None,ind=None,stride=1): |
|---|
| 1220 | if self.method_t == "free": |
|---|
| 1221 | self.index_t = np.arange(0,self.dim_t,stride) |
|---|
| 1222 | if self.dim_t > 1: |
|---|
| 1223 | self.dimplot = self.dimplot + 1 |
|---|
| 1224 | if self.verbose: print "**** OK. t values. all." |
|---|
| 1225 | else: |
|---|
| 1226 | self.method_t = "fixed" |
|---|
| 1227 | if self.verbose: print "**** OK. no t dimension." |
|---|
| 1228 | elif self.method_t == "comp": |
|---|
| 1229 | start = np.argmin( np.abs( self.field_t - dalist[ind][0] ) ) |
|---|
| 1230 | stop = np.argmin( np.abs( self.field_t - dalist[ind][1] ) ) |
|---|
| 1231 | self.index_t = np.arange(start,stop,stride) |
|---|
| 1232 | if self.verbose: print "**** OK. t values. comp over interval ",self.field_t[start],self.field_t[stop]," nvalues=",self.index_t.size |
|---|
| 1233 | elif self.method_t == "fixed": |
|---|
| 1234 | self.index_t.append( np.argmin( np.abs( self.field_t - dalist[ind][0] ) )) |
|---|
| 1235 | if self.verbose: print "**** OK. t values",self.field_t[self.index_t] |
|---|
| 1236 | else: |
|---|
| 1237 | print "!! ERROR !! method "+self.method_t+" not supported" |
|---|
| 1238 | self.index_t = np.array(self.index_t) |
|---|
| 1239 | |
|---|
| 1240 | # get list of index to be retrieved for vertical axis |
|---|
| 1241 | ### TBD: il faudrait ne prendre que les indices qui correspondent a l interieur d un plot (dans all) |
|---|
| 1242 | # ------------------------------- |
|---|
| 1243 | def getindexvert(self,dalist=None,ind=None,stride=1): |
|---|
| 1244 | if self.method_z == "free": |
|---|
| 1245 | self.index_z = np.arange(0,self.dim_z,stride) |
|---|
| 1246 | if self.dim_z > 1: |
|---|
| 1247 | self.dimplot = self.dimplot + 1 |
|---|
| 1248 | if self.verbose: print "**** OK. z values. all." |
|---|
| 1249 | else: |
|---|
| 1250 | self.method_z = "fixed" |
|---|
| 1251 | if self.verbose: print "**** OK. no z dimension." |
|---|
| 1252 | elif self.method_z == "comp": |
|---|
| 1253 | start = np.argmin( np.abs( self.field_z - dalist[ind][0] ) ) |
|---|
| 1254 | stop = np.argmin( np.abs( self.field_z - dalist[ind][1] ) ) |
|---|
| 1255 | self.index_z = np.arange(start,stop,stride) |
|---|
| 1256 | if self.verbose: print "**** OK. z values. comp over interval",self.field_z[start],self.field_z[stop]," nvalues=",self.index_z.size |
|---|
| 1257 | elif self.method_z == "fixed": |
|---|
| 1258 | self.index_z.append( np.argmin( np.abs( self.field_z - dalist[ind][0] ) )) |
|---|
| 1259 | if self.verbose: print "**** OK. z values",self.field_z[self.index_z] |
|---|
| 1260 | else: |
|---|
| 1261 | if self.verbose: print "!! ERROR !! method "+self.method_z+" not supported" |
|---|
| 1262 | self.index_z = np.array(self.index_z) |
|---|
| 1263 | |
|---|
| 1264 | # get list of index to be retrieved for horizontal grid |
|---|
| 1265 | # --> index_x and index_y are slices to be retrieved from NETCDF files |
|---|
| 1266 | # --> index_x2d and index_y2d are the actual (x,y) coordinates corresponding to each relevant point |
|---|
| 1267 | # [this is slightly more complicated because 2D arrays for lat-lon projection possibly irregular] |
|---|
| 1268 | # NB: to append index we use lists (the most convenient) then we convert into a numpy.array |
|---|
| 1269 | ### TBD: il faudrait ne prendre que les indices qui correspondent a l interieur d un plot (dans all) |
|---|
| 1270 | # ------------------------------- |
|---|
| 1271 | def getindexhori(self,dalistx=None,dalisty=None,indx=None,indy=None,stridex=1,stridey=1): |
|---|
| 1272 | ## get what is the method over x and y axis |
|---|
| 1273 | test = self.method_x+self.method_y |
|---|
| 1274 | ## CASE 0, EASY CASES: |
|---|
| 1275 | ## - LAT IS FREE (we do here what must be done whatever LON case is) |
|---|
| 1276 | ## - LON IS FREE (we do here what must be done whatever LAT case is) |
|---|
| 1277 | ## - LAT IS COMP AND LON IS FREE |
|---|
| 1278 | ## - LON IS COMP AND LAT IS FREE |
|---|
| 1279 | if self.method_x == "free" or test in ["compfree","compcomp"]: |
|---|
| 1280 | self.index_x = range(0,self.dim_x,stridex) |
|---|
| 1281 | if self.dim_x > 1: |
|---|
| 1282 | if self.method_x == "free": self.dimplot = self.dimplot + 1 |
|---|
| 1283 | if self.verbose: print "**** OK. x values. all." |
|---|
| 1284 | else: |
|---|
| 1285 | self.method_x = "fixed" |
|---|
| 1286 | if self.verbose: print "**** OK. no x dimension." |
|---|
| 1287 | if self.method_y == "free" or test in ["freecomp","compcomp"]: |
|---|
| 1288 | self.index_y = range(0,self.dim_y,stridey) |
|---|
| 1289 | if self.dim_y > 1: |
|---|
| 1290 | if self.method_y == "free": self.dimplot = self.dimplot + 1 |
|---|
| 1291 | if self.verbose: print "**** OK. y values. all." |
|---|
| 1292 | else: |
|---|
| 1293 | self.method_y = "fixed" |
|---|
| 1294 | if self.verbose: print "**** OK. no y dimension." |
|---|
| 1295 | ## CASE 0 above, this is just for continuity for free. |
|---|
| 1296 | ## ... for comp we have to select bounds. |
|---|
| 1297 | ## ... TBD: take the bool array strategy for what follows! |
|---|
| 1298 | if self.method_x in ["free","comp"] and self.method_y in ["free","comp"]: |
|---|
| 1299 | ### ref1_dirty_hack |
|---|
| 1300 | ### ... for the moment this is a hack. but actually this is more powerful. |
|---|
| 1301 | if self.method_x == "comp": |
|---|
| 1302 | yeah = (self.field_x >= dalistx[indx][0])*(self.field_x <= dalistx[indx][1]) |
|---|
| 1303 | self.index_x = yeah[0,:] |
|---|
| 1304 | if self.method_y == "comp": |
|---|
| 1305 | yeah = (self.field_y >= dalisty[indy][0]) * (self.field_y <= dalisty[indy][1]) |
|---|
| 1306 | self.index_y = yeah[:,0] |
|---|
| 1307 | self.index_x2d = self.index_x |
|---|
| 1308 | self.index_y2d = self.index_y |
|---|
| 1309 | ## AND NOW THE LITTLE BIT MORE COMPLICATED CASES |
|---|
| 1310 | ## CASE 1 LAT AND LON ARE FIXED |
|---|
| 1311 | elif test == "fixedfixed": |
|---|
| 1312 | idy,idx = np.unravel_index( np.argmin( ( self.field_x - dalistx[indx][0])**2 + (self.field_y - dalisty[indy][0])**2 ), self.field_x.shape ) |
|---|
| 1313 | #TBD: pb with staggered coord |
|---|
| 1314 | if idx not in self.index_x: self.index_x.append(idx) |
|---|
| 1315 | if idy not in self.index_y: self.index_y.append(idy) |
|---|
| 1316 | self.index_x2d.append(idx) |
|---|
| 1317 | self.index_y2d.append(idy) |
|---|
| 1318 | ## CASE 2 LON IS FIXED BUT NOT LAT |
|---|
| 1319 | elif test in ["fixedfree","fixedcomp"]: |
|---|
| 1320 | # find where are requested x values for each y on the free dimension |
|---|
| 1321 | # NB: this does not work for non-bijective cases e.g. polar stereographic |
|---|
| 1322 | for iy in range(self.dim_y): |
|---|
| 1323 | idx = np.argmin( np.abs( self.field_x[iy,:] - dalistx[indx][0] ) ) |
|---|
| 1324 | # if comp is requested we select only indexes which yield values between requested min and max |
|---|
| 1325 | storeval = (self.method_y == "comp") and (self.field_y[iy,idx] > dalisty[indy][0]) and (self.field_y[iy,idx] < dalisty[indy][1]) |
|---|
| 1326 | storeval = storeval or (self.method_y == "free") |
|---|
| 1327 | if storeval: |
|---|
| 1328 | if idx not in self.index_x: self.index_x.append(idx) |
|---|
| 1329 | if iy not in self.index_y and self.method_y == "comp": self.index_y.append(iy) |
|---|
| 1330 | if idx not in self.index_x2d or iy not in self.index_y2d: |
|---|
| 1331 | self.index_x2d.append(idx) |
|---|
| 1332 | self.index_y2d.append(iy) |
|---|
| 1333 | ## CASE 3 LAT IS FIXED BUT NOT LON |
|---|
| 1334 | elif test in ["freefixed","compfixed"]: |
|---|
| 1335 | # find where are requested y values for each x on the free dimension |
|---|
| 1336 | # NB: this does not work for non-bijective cases e.g. polar stereographic |
|---|
| 1337 | for ix in range(self.dim_x): |
|---|
| 1338 | idy = np.argmin( np.abs( self.field_y[:,ix] - dalisty[indy][0] ) ) |
|---|
| 1339 | # if comp is requested we select only indexes which yield values between requested min and max |
|---|
| 1340 | storeval = (self.method_x == "comp") and (self.field_x[idy,ix] > dalistx[indx][0]) and (self.field_x[idy,ix] < dalistx[indx][1]) |
|---|
| 1341 | storeval = storeval or (self.method_x == "free") |
|---|
| 1342 | if storeval: |
|---|
| 1343 | if idy not in self.index_y: self.index_y.append(idy) |
|---|
| 1344 | if ix not in self.index_x and self.method_x == "comp": self.index_x.append(ix) |
|---|
| 1345 | if ix not in self.index_x2d or idy not in self.index_y2d: |
|---|
| 1346 | self.index_x2d.append(ix) |
|---|
| 1347 | self.index_y2d.append(idy) |
|---|
| 1348 | ## check index tab |
|---|
| 1349 | if len(self.index_x) == 0 or len(self.index_y) == 0: |
|---|
| 1350 | print "!! ERROR !! no indices found. check prescribed latitudes or longitudes" ; exit() |
|---|
| 1351 | ## ensure the array is a numpy array for getfield to work |
|---|
| 1352 | self.index_x = np.array(self.index_x) |
|---|
| 1353 | self.index_y = np.array(self.index_y) |
|---|
| 1354 | self.index_x2d = np.array(self.index_x2d) |
|---|
| 1355 | self.index_y2d = np.array(self.index_y2d) |
|---|
| 1356 | ### print extrema |
|---|
| 1357 | printx = self.field_x[np.ix_(self.index_y2d, self.index_x2d)] |
|---|
| 1358 | printy = self.field_y[np.ix_(self.index_y2d, self.index_x2d)] |
|---|
| 1359 | if self.verbose: |
|---|
| 1360 | print "**** OK. x values (min,max).", printx.min(),printx.max() |
|---|
| 1361 | print "**** OK. y values (min,max).", printy.min(),printy.max() |
|---|
| 1362 | |
|---|
| 1363 | # get the field from the NETCDF file and perform averages |
|---|
| 1364 | # ------------------------------- |
|---|
| 1365 | def getfield(self): |
|---|
| 1366 | ## first tell what is to be done |
|---|
| 1367 | if self.dimplot > 2: print "**** !! ERROR !! "+str(self.dimplot)+"D plots not supported!" ; exit() |
|---|
| 1368 | elif self.dimplot == 0 and self.verbose: print "**** OK. 0D value requested." |
|---|
| 1369 | elif self.dimplot == 1 and self.verbose: print "**** OK. 1D plot requested." |
|---|
| 1370 | elif self.verbose: print "**** OK. 2D section requested." |
|---|
| 1371 | # well, now get field from netcdf file |
|---|
| 1372 | # part below is necessary otherwise there is an index error below |
|---|
| 1373 | if self.index_x.size == 1: self.index_x = self.index_x[0] |
|---|
| 1374 | if self.index_y.size == 1: self.index_y = self.index_y[0] |
|---|
| 1375 | if self.index_z.size == 1: self.index_z = self.index_z[0] |
|---|
| 1376 | if self.index_t.size == 1: self.index_t = self.index_t[0] |
|---|
| 1377 | # then retrieve what is requested by user |
|---|
| 1378 | # each self.dim case corresponds to tests in the beginning of getdim. |
|---|
| 1379 | time0 = timelib.time() |
|---|
| 1380 | if self.verbose: print "**** OK. I am getting values from files. Please wait." |
|---|
| 1381 | if self.dim == 1: |
|---|
| 1382 | nt = self.index_t.size ; nz = 1 ; ny = 1 ; nx = 1 |
|---|
| 1383 | self.field = self.f.variables[self.var][self.index_t] |
|---|
| 1384 | elif self.dim == 2: |
|---|
| 1385 | nt = 1 ; nz = 1 ; ny = self.index_y.size ; nx = self.index_x.size |
|---|
| 1386 | self.field = self.f.variables[self.var][self.index_y,self.index_x] |
|---|
| 1387 | elif self.dim == 3: |
|---|
| 1388 | nt = self.index_t.size ; nz = 1 ; ny = self.index_y.size ; nx = self.index_x.size |
|---|
| 1389 | self.field = self.f.variables[self.var][self.index_t,self.index_y,self.index_x] |
|---|
| 1390 | # this is far faster than retrieving each term with a loop |
|---|
| 1391 | elif self.dim == 4: |
|---|
| 1392 | nt = self.index_t.size ; nz = self.index_z.size ; ny = self.index_y.size ; nx = self.index_x.size |
|---|
| 1393 | self.field = self.f.variables[self.var][self.index_t,self.index_z,self.index_y,self.index_x] |
|---|
| 1394 | else: |
|---|
| 1395 | print "!! ERROR !! field would have more than four dimensions ?" ; exit() |
|---|
| 1396 | # dirty hack (AS) ref1_dirty_hack |
|---|
| 1397 | # waiting for more fundamental modifications. case when self.index_y is a bool array. |
|---|
| 1398 | # ... be careful if no point... |
|---|
| 1399 | try: |
|---|
| 1400 | if type(self.index_x[0]) == np.bool_: nx = np.sum(self.index_x) ## gives the size of the True part! |
|---|
| 1401 | if type(self.index_y[0]) == np.bool_: ny = np.sum(self.index_y) ## gives the size of the True part! |
|---|
| 1402 | except: |
|---|
| 1403 | pass |
|---|
| 1404 | # NB: ... always 4D array but possibly with "size 1" dimensions |
|---|
| 1405 | # ... if one dimension is missing because 1D 2D or 3D requests, make it appear again |
|---|
| 1406 | self.field = np.reshape(self.field,(nt,nz,ny,nx)) |
|---|
| 1407 | if self.verbose: print "**** OK. I got %7.1e values. This took me %6.4f seconds" % (nx*ny*nz*nt,timelib.time() - time0) |
|---|
| 1408 | if self.verbose: print "**** OK. I got var "+self.var+" with shape",self.field.shape |
|---|
| 1409 | # reduce coordinates to useful points |
|---|
| 1410 | # ... TBD: this should be ordered in the case of non-regular projections |
|---|
| 1411 | if self.method_x in ["free","comp"] and self.method_y in ["free","comp"]: |
|---|
| 1412 | # we need 2D coordinates (free) or we get broadcast problem (comp) so we use np.ix |
|---|
| 1413 | self.field_x = self.field_x[np.ix_(self.index_y2d, self.index_x2d)] |
|---|
| 1414 | self.field_y = self.field_y[np.ix_(self.index_y2d, self.index_x2d)] |
|---|
| 1415 | else: |
|---|
| 1416 | # we are OK with 1D coordinates |
|---|
| 1417 | self.field_x = self.field_x[self.index_y2d, self.index_x2d] |
|---|
| 1418 | self.field_y = self.field_y[self.index_y2d, self.index_x2d] |
|---|
| 1419 | self.field_z = self.field_z[self.index_z] |
|---|
| 1420 | self.field_t = self.field_t[self.index_t] |
|---|
| 1421 | # extract relevant horizontal points |
|---|
| 1422 | # TBD: is compfree OK with computing on irregular grid? |
|---|
| 1423 | test = self.method_x + self.method_y |
|---|
| 1424 | if test in ["fixedfixed","freefree"]: |
|---|
| 1425 | pass |
|---|
| 1426 | elif test in ["fixedfree","fixedcomp"] or test in ["freefixed","compfixed"]: |
|---|
| 1427 | time0 = timelib.time() |
|---|
| 1428 | # now have to obtain the new indexes which correspond to the extracted self.field |
|---|
| 1429 | # for it to work with unique index, ensure that any index_* is a numpy array |
|---|
| 1430 | if not isinstance(self.index_x, np.ndarray): self.index_x = np.array([self.index_x]) |
|---|
| 1431 | if not isinstance(self.index_y, np.ndarray): self.index_y = np.array([self.index_y]) |
|---|
| 1432 | if not isinstance(self.index_z, np.ndarray): self.index_z = np.array([self.index_z]) |
|---|
| 1433 | if not isinstance(self.index_t, np.ndarray): self.index_t = np.array([self.index_t]) |
|---|
| 1434 | for val in self.index_x: self.index_x2d[np.where(self.index_x2d == val)] = np.where(self.index_x == val)[0] |
|---|
| 1435 | for val in self.index_y: self.index_y2d[np.where(self.index_y2d == val)] = np.where(self.index_y == val)[0] |
|---|
| 1436 | ##### VERY EXPENSIVE |
|---|
| 1437 | ## recast self.field with 2D horizontal arrays because we might have extracted |
|---|
| 1438 | ## more than what is to be actually plot or computed, in particular for comps on 2D lat,lon coordinates |
|---|
| 1439 | #self.field = self.field[np.ix_(self.index_t,self.index_z,self.index_y2d,self.index_x2d)] |
|---|
| 1440 | #(nt,nz,ny,nx) = self.field.shape |
|---|
| 1441 | # prepare the loop on all relevant horizontal points |
|---|
| 1442 | if self.method_x in ["comp","free"]: |
|---|
| 1443 | nnn = self.index_x2d.shape[0] ; what_I_am_supposed_to_do = "keepx" |
|---|
| 1444 | elif self.method_y in ["comp","free"]: |
|---|
| 1445 | nnn = self.index_y2d.shape[0] ; what_I_am_supposed_to_do = "keepy" |
|---|
| 1446 | # LOOP to extract only useful values over horizontal dimensions |
|---|
| 1447 | # only take diagonal terms, do not loop on all self.index_x2d*self.index_y2d |
|---|
| 1448 | # ... this method is fast enough, perhaps there is a faster way though |
|---|
| 1449 | # ... (for sure the method with np.diag is much slower) |
|---|
| 1450 | for iii in range(nnn): |
|---|
| 1451 | ix = self.index_x2d[iii] ; iy = self.index_y2d[iii] |
|---|
| 1452 | for iz in range(self.index_z.size): |
|---|
| 1453 | for it in range(self.index_t.size): |
|---|
| 1454 | if what_I_am_supposed_to_do == "keepx": self.field[it,iz,0,ix] = self.field[it,iz,iy,ix] |
|---|
| 1455 | elif what_I_am_supposed_to_do == "keepy": self.field[it,iz,iy,0] = self.field[it,iz,iy,ix] |
|---|
| 1456 | if self.verbose: print "**** OK. I got to pick the right values for your request. This took me %6.4f seconds" % (timelib.time() - time0) |
|---|
| 1457 | # we only keep the one value that was modified on the dimension which is not free |
|---|
| 1458 | if what_I_am_supposed_to_do == "keepx": self.field = self.field[:,:,0,:] ; ny = 1 ; self.field = np.reshape(self.field,(nt,nz,ny,nx)) |
|---|
| 1459 | elif what_I_am_supposed_to_do == "keepy": self.field = self.field[:,:,:,0] ; nx = 1 ; self.field = np.reshape(self.field,(nt,nz,ny,nx)) |
|---|
| 1460 | # make a mask in case there are non-NaN missing values. (what about NaN missing values?) |
|---|
| 1461 | # ... this is important for computations below (see ppcompute) |
|---|
| 1462 | masked = np.ma.masked_where(np.abs(self.field) > 1e25,self.field) |
|---|
| 1463 | if masked.mask.any() == True: |
|---|
| 1464 | if self.verbose: print "!! WARNING !! Values over +-1e25 are considered missing values." |
|---|
| 1465 | self.field = masked |
|---|
| 1466 | self.field.set_fill_value([np.NaN]) |
|---|
| 1467 | |
|---|
| 1468 | # perform computations |
|---|
| 1469 | # ------------------------------- |
|---|
| 1470 | # available: mean, max, min, meanarea |
|---|
| 1471 | # TB: integrals? for derivatives, define a function self.dx() |
|---|
| 1472 | def computations(self): |
|---|
| 1473 | nt,nz,ny,nx = self.field.shape |
|---|
| 1474 | # treat the case of mean on fields normalized with grid mesh area |
|---|
| 1475 | # ... this is done in the .area() method. |
|---|
| 1476 | # after that self.field contains field*area/totarea |
|---|
| 1477 | if "area" in self.compute: |
|---|
| 1478 | if "comp" in self.method_x+self.method_y: |
|---|
| 1479 | self.area() |
|---|
| 1480 | else: |
|---|
| 1481 | if self.verbose: print "!! WARNING !! No area accounted for (computing on t and/or z axis)." |
|---|
| 1482 | # now ready to compute [TBD: we would like to have e.g. mean over x,y and min over t ??] |
|---|
| 1483 | if self.method_t == "comp": |
|---|
| 1484 | if self.verbose: print "**** OK. Computing over t axis." |
|---|
| 1485 | if "mean" in self.compute: self.field = ppcompute.mean(self.field,axis=0) |
|---|
| 1486 | elif self.compute == "min": self.field = ppcompute.min(self.field,axis=0) |
|---|
| 1487 | elif self.compute == "max": self.field = ppcompute.max(self.field,axis=0) |
|---|
| 1488 | else: print "!! ERROR !! operation not supported." ; exit() |
|---|
| 1489 | nt = 1 ; self.field = np.reshape(self.field,(nt,nz,ny,nx)) |
|---|
| 1490 | if self.method_z == "comp": |
|---|
| 1491 | if self.verbose: print "**** OK. Computing over z axis." |
|---|
| 1492 | if "mean" in self.compute: self.field = ppcompute.mean(self.field,axis=1) |
|---|
| 1493 | elif self.compute == "min": self.field = ppcompute.min(self.field,axis=1) |
|---|
| 1494 | elif self.compute == "max": self.field = ppcompute.max(self.field,axis=1) |
|---|
| 1495 | else: print "!! ERROR !! operation not supported." ; exit() |
|---|
| 1496 | nz = 1 ; self.field = np.reshape(self.field,(nt,nz,ny,nx)) |
|---|
| 1497 | if self.method_y == "comp": |
|---|
| 1498 | if self.verbose: print "**** OK. Computing over y axis." |
|---|
| 1499 | if self.compute == "mean": self.field = ppcompute.mean(self.field,axis=2) |
|---|
| 1500 | elif self.compute == "min": self.field = ppcompute.min(self.field,axis=2) |
|---|
| 1501 | elif self.compute == "max": self.field = ppcompute.max(self.field,axis=2) |
|---|
| 1502 | elif self.compute == "meanarea": self.field = ppcompute.sum(self.field,axis=2) |
|---|
| 1503 | else: print "!! ERROR !! operation not supported." ; exit() |
|---|
| 1504 | ny = 1 ; self.field = np.reshape(self.field,(nt,nz,ny,nx)) |
|---|
| 1505 | if self.field_x.ndim == 2: self.field_x = self.field_x[0,:] # TBD: this is OK for regular grid but not for irregular |
|---|
| 1506 | if self.method_x == "comp": |
|---|
| 1507 | if self.verbose: print "**** OK. Computing over x axis." |
|---|
| 1508 | if self.compute == "mean": self.field = ppcompute.mean(self.field,axis=3) |
|---|
| 1509 | elif self.compute == "min": self.field = ppcompute.min(self.field,axis=3) |
|---|
| 1510 | elif self.compute == "max": self.field = ppcompute.max(self.field,axis=3) |
|---|
| 1511 | elif self.compute == "meanarea": self.field = ppcompute.sum(self.field,axis=3) |
|---|
| 1512 | else: print "!! ERROR !! operation not supported." ; exit() |
|---|
| 1513 | nx = 1 ; self.field = np.reshape(self.field,(nt,nz,ny,nx)) |
|---|
| 1514 | if self.field_y.ndim == 2: self.field_y = self.field_y[:,0] # TBD: this is OK for regular grid but not for irregular |
|---|
| 1515 | # remove all dimensions with size 1 to prepare plot (and check the resulting dimension with dimplot) |
|---|
| 1516 | self.field = np.squeeze(self.field) |
|---|
| 1517 | if self.field.ndim != self.dimplot: |
|---|
| 1518 | print "!! ERROR !! Problem: self.field is different than plot dimensions", self.field.ndim, self.dimplot ; exit() |
|---|
| 1519 | if self.verbose: |
|---|
| 1520 | print "**** OK. Final shape for "+self.var+" after averaging and squeezing",self.field.shape |
|---|
| 1521 | |
|---|
| 1522 | # get areas for computations and ponderate self.field by area/totarea |
|---|
| 1523 | # ------------------------------------------------------------------- |
|---|
| 1524 | def area(self): |
|---|
| 1525 | if self.verbose: print "**** OK. Get area array for computations." |
|---|
| 1526 | # create a request object for area |
|---|
| 1527 | # ... and copy known attributes from self |
|---|
| 1528 | aire = onerequest() |
|---|
| 1529 | aire.copy(self) |
|---|
| 1530 | # get area field name |
|---|
| 1531 | aire.var = "nothing" |
|---|
| 1532 | for c in glob_listarea: |
|---|
| 1533 | if c in aire.f.variables.keys(): |
|---|
| 1534 | aire.var = c |
|---|
| 1535 | # do not try to calculate areas automatically |
|---|
| 1536 | if aire.var == "nothing": |
|---|
| 1537 | print "!! ERROR !! area variable not found... needs to be added in set_ppclass.txt?" |
|---|
| 1538 | exit() |
|---|
| 1539 | # define area request dimensions |
|---|
| 1540 | aire.getdim() |
|---|
| 1541 | # ensure this is a 2D horizontal request and define indexes |
|---|
| 1542 | # ... areas are not supposed to vary with time and height |
|---|
| 1543 | aire.method_x = "free" ; aire.method_y = "free" |
|---|
| 1544 | aire.getindexhori() ; aire.dimplot = 2 |
|---|
| 1545 | aire.method_z = "fixed" ; aire.field_z = np.array([0]) ; aire.index_z = np.array([0]) |
|---|
| 1546 | aire.method_t = "fixed" ; aire.field_t = np.array([0]) ; aire.index_t = np.array([0]) |
|---|
| 1547 | # read the 2D area array in netCDF file |
|---|
| 1548 | aire.getfield() |
|---|
| 1549 | aire.field = np.squeeze(aire.field) |
|---|
| 1550 | # reduce with self horizontal indexes |
|---|
| 1551 | if "fixed" in self.method_x+self.method_y: |
|---|
| 1552 | aire.field = aire.field[self.index_y,self.index_x] |
|---|
| 1553 | # calculate total area |
|---|
| 1554 | # ... 2D comp is easy. 1D comp is a bit less easy but simple array manipulation. |
|---|
| 1555 | if "free" in self.method_x+self.method_y: |
|---|
| 1556 | if self.method_x == "free": |
|---|
| 1557 | totarea = ppcompute.sum(aire.field,axis=0) |
|---|
| 1558 | totarea = np.reshape(totarea,(1,totarea.size)) |
|---|
| 1559 | totarea = np.tile(totarea,(1,self.index_x)) |
|---|
| 1560 | elif self.method_y == "free": |
|---|
| 1561 | totarea = ppcompute.sum(aire.field,axis=1) |
|---|
| 1562 | totarea = np.reshape(totarea,(totarea.size,1)) |
|---|
| 1563 | totarea = np.tile(totarea,(1,self.index_x.size)) |
|---|
| 1564 | elif self.method_x == "comp" and self.method_y == "comp": |
|---|
| 1565 | aire.field = aire.field[np.ix_(self.index_y, self.index_x)] # reduce to requested indexes only |
|---|
| 1566 | totarea = ppcompute.sum(ppcompute.sum(aire.field,axis=1),axis=0) |
|---|
| 1567 | else: |
|---|
| 1568 | if self.verbose: print "!! WARNING !! Not account for areas. Only averaging over z and/or t axis." |
|---|
| 1569 | # normalize by total area |
|---|
| 1570 | print "**** OK. I can now normalize field by areas." |
|---|
| 1571 | aire.field = aire.field / totarea |
|---|
| 1572 | # tile area array over self t and z axis so that area field could be multiplied with self.field |
|---|
| 1573 | aire.field = np.tile(aire.field,(self.index_t.size,self.index_z.size,1,1)) |
|---|
| 1574 | if self.field.shape != aire.field.shape: |
|---|
| 1575 | print "!! ERROR !! Problem in area(). Check array shapes." |
|---|
| 1576 | print "Field vs. aire:",self.field.shape,aire.field.shape ; exit() |
|---|
| 1577 | else: |
|---|
| 1578 | self.field = self.field*aire.field |
|---|
| 1579 | |
|---|
| 1580 | # define coordinates for plot |
|---|
| 1581 | # ------------------------------- |
|---|
| 1582 | def definecoord(self): |
|---|
| 1583 | I_got_abs = False ; I_got_ord = False |
|---|
| 1584 | # here is the thing. time is usually taken as an abscissa so we start with time. |
|---|
| 1585 | if self.method_t == "free": |
|---|
| 1586 | self.absc = self.field_t ; self.absclab = self.name_t |
|---|
| 1587 | I_got_abs = True |
|---|
| 1588 | # then we usually have x as an abscissa. |
|---|
| 1589 | if self.method_x == "free": |
|---|
| 1590 | if I_got_abs: |
|---|
| 1591 | self.ordi = self.field_x ; self.ordilab = self.name_x |
|---|
| 1592 | I_got_ord = True |
|---|
| 1593 | else: |
|---|
| 1594 | self.absc = self.field_x ; self.absclab = self.name_x |
|---|
| 1595 | I_got_abs = True |
|---|
| 1596 | # ... or we have y |
|---|
| 1597 | if self.method_y == "free": |
|---|
| 1598 | if I_got_abs: |
|---|
| 1599 | self.ordi = self.field_y ; self.ordilab = self.name_y |
|---|
| 1600 | I_got_ord = True |
|---|
| 1601 | else: |
|---|
| 1602 | self.absc = self.field_y ; self.absclab = self.name_y |
|---|
| 1603 | I_got_abs = True |
|---|
| 1604 | # ... and we end with z because it is usually not an abscissa (profiles). |
|---|
| 1605 | if self.method_z == "free": |
|---|
| 1606 | if self.field_z[0] > self.field_z[1]: |
|---|
| 1607 | self.invert_axes = True # the axis will be turned upside-down |
|---|
| 1608 | if I_got_abs: |
|---|
| 1609 | self.ordi = self.field_z ; self.ordilab = self.name_z |
|---|
| 1610 | I_got_ord = True |
|---|
| 1611 | else: |
|---|
| 1612 | self.absc = self.field_z ; self.absclab = self.name_z |
|---|
| 1613 | I_got_abs = True |
|---|
| 1614 | self.swap_axes = True # says that altitude is not supposed to remain as an abscissa |
|---|
| 1615 | if I_got_abs and self.verbose: print "**** OK. abscissa:",self.absclab, self.absc.shape |
|---|
| 1616 | if I_got_ord and self.verbose: print "**** OK. ordinate:",self.ordilab, self.ordi.shape |
|---|
