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