Changeset 960 for trunk/UTIL


Ignore:
Timestamp:
May 15, 2013, 8:29:27 PM (12 years ago)
Author:
aslmd
Message:

UTIL PYTHON planetoplot_v2. MAJOR: created an easy way to load a variable from a netcdf file, see easy_get_field.py [note: now to avoid confusion with .f attributes containing field, .f() method is named .func()] ; MINOR: generic units keyword, quiet mode, dummy xy axis when not given in ppplot, bug fixes for one-value 0D requests.

Location:
trunk/UTIL/PYTHON/planetoplot_v2
Files:
1 added
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/UTIL/PYTHON/planetoplot_v2/examples/easy_get_field.py

    r931 r960  
    11#! /usr/bin/env python
    22from ppclass import pp
    3 import numpy as np
    4 import matplotlib.pyplot as mpl
    53
     4########
     5## ppclass allows for getting fields very easily from a netcdf file
     6## ... this is contained in the "f" attributes of the pp object
     7## ... and the "l" attributes of the pp object contains useful information
     8## ... two methods getf() and getfl() are helpful shortcuts
     9###########################################
     10
     11## 1 --> an unique and straightforward request
     12##   --> very easy ! see also minimal_field.py
     13##############################################
     14icetot = pp(file="/home/aymeric/Big_Data/DATAPLOT/diagfired.nc",var="icetot",t=0.5).getf()
     15print "icetot", icetot[10,10]
     16
     17## 2 --> simple multiple request, no labelling
     18##############################################
     19test = pp(file="/home/aymeric/Big_Data/DATAPLOT/diagfired.nc",t=0.5)
     20test.var = ["mtot","icetot"]
     21allf = test.getf() # or allf = test.get().f
     22##
     23mtot = allf[0]
     24icetot = allf[1]
     25print "mtot", mtot[10,10]
     26print "icetot", icetot[10,10]
     27
     28## 3 --> complex multiple requests and labelling
     29################################################
     30test = pp(file="/home/aymeric/Big_Data/DATAPLOT/diagfired.nc")
     31test.var = ["mtot","icetot"]
     32test.t = [0.4,0.5]
     33allf,lab = test.getfl()
     34##
     35icetot04 = allf[lab.index("_v=icetot_t=0.4_")]
     36mtot04 = allf[lab.index("_v=mtot_t=0.4_")]
     37icetot05 = allf[lab.index("_v=icetot_t=0.5_")]
     38mtot05 = allf[lab.index("_v=mtot_t=0.5_")]
     39print "mtot04", mtot04[10,10]
     40print "icetot04", icetot04[10,10]
     41print "mtot05", mtot05[10,10]
     42print "icetot05", icetot05[10,10]
     43
     44## 4 --> an example of complete labelling
     45## .... a rather unlikely example ....
     46## .... but shows label ordering ....
     47#########################################
    648test = pp()
    7 test.file = [ \
    8   "/planeto/jllmd/Earth/Nmix/F1525/output/diagfi4.nc", \
    9   "/planeto/jllmd/Earth/Nmix/F1550/output/diagfi4.nc", \
    10   "/planeto/jllmd/Earth/Nmix/F1500/output/diagfi4.nc" ]
    11 test.var = ["tsurf","ASR"]
    12 test.x = "-180.,175."
    13 test.y = "-90.,90."
    14 test.t = "0,1000"
    15 #test.verbose = True
    16 test.compute = "meanarea"
    17 test.get()
    18 
    19 # -----------------------------
    20 # self.allfield contains data !
    21 # -----------------------------
    22 # index ordering:
    23 # - file
    24 # - var
    25 # - request t
    26 # - request z
    27 # - request y
    28 # - request x
    29 # - field dimensions
    30 # -----------------------------
    31 
    32 
    33 for ii in range(test.nfin):
    34     print "tsurf",test.allfield[ii][0][0][0][0][0]
    35     print "asr",test.allfield[ii][1][0][0][0][0]
    36 
    37 #mpl.plot(tsurf,asr)
    38 #mpl.show()
     49test.file = ["/home/aymeric/Big_Data/DATAPLOT/diagfired.nc","/home/aymeric/Big_Data/DATAPLOT/diagfired.nc"]
     50test.var = ["u","v"]
     51test.x = [10.,20.]
     52test.y = [10.,20.]
     53test.z = [10.,20.]
     54test.t = [0.4,0.5]
     55print "... please wait. this one is a bit stupid..."
     56allf,lab = test.getfl()
     57## note label ordering: file -- var -- x -- y -- z -- t
     58l1 = "_f=#2_v=u_x=10.0_y=10.0_z=20.0_t=0.4_"
     59l2 = "_f=#2_v=v_x=10.0_y=10.0_z=20.0_t=0.4_"
     60u_example = allf[lab.index(l1)]
     61v_example = allf[lab.index(l2)]
     62print l1, u_example
     63print l2, v_example

  • trunk/UTIL/PYTHON/planetoplot_v2/examples/function.py

    r933 r960  
    1919var2 = var2 / 3.72
    2020
    21 S = var2.f(var1)
     21S = var2.func(var1)
    2222
    2323S.p[0].marker = 'o'

  • trunk/UTIL/PYTHON/planetoplot_v2/examples/hodograph.py

    r923 r960  
    1616
    1717# u as a function of v
    18 hodo = u.f(v)
     18hodo = u.func(v)
    1919hodo.filename = "hodograph"
    2020hodo.makeplot()
    2121
    2222# v as a function of u
    23 hodo2 = v.f(u)
     23hodo2 = v.func(u)
    2424hodo2.makeplot()

  • trunk/UTIL/PYTHON/planetoplot_v2/examples/meso_profile.py

    r923 r960  
    3838
    3939# define potential temperature as a function of height
    40 S = tpot.f(z)
     40S = tpot.func(z)
    4141
    4242# change a few plot settings

  • trunk/UTIL/PYTHON/planetoplot_v2/examples/scatter.py

    r923 r960  
    1515ps.getdefineplot()
    1616
    17 S = ps.f(tsurf)
     17S = ps.func(tsurf)
    1818S.p[0].lstyle=""
    1919S.p[0].marker="h"
    … …  
    2626icetot.getdefineplot()
    2727
    28 S2 = icetot.f(tsurf)
     28S2 = icetot.func(tsurf)
    2929S2.p[0].lstyle=""
    3030S2.p[0].marker="D"
    … …  
    4545wind = u**2 + v**2
    4646wind = wind**0.5
    47 S3 = wind.f(ps)
     47S3 = wind.func(ps)
    4848S3.p[0].lstyle=""
    4949S3.p[0].marker="o"
    … …  
    7070ps.getdefineplot()
    7171
    72 S = ps.f(tsurf)
     72S = ps.func(tsurf)
    7373S.p[0].lstyle=""
    7474S.p[0].marker="h"

  • trunk/UTIL/PYTHON/planetoplot_v2/pp_reload.py

    r942 r960  
    2626
    2727    yeah = pp()
     28    yeah.quiet = True
    2829    yeah.defineplot(loadfile=files)
    2930    yeah.out = opt.out

  • trunk/UTIL/PYTHON/planetoplot_v2/ppclass.py

    r942 r960  
    2525###################################
    2626# where settings files are located...
     27# ... this can be hardcoded here
    2728whereset = None
    2829whereset = ppcompute.findset(whereset)
    … …  
    3940    for stuff in lines[17].strip().split(';'): glob_listarea.append(stuff)
    4041except IOError:
    41     print "warning: "+zefile+" not in "+whereset+" ; no presets."
     42    print "PPCLASS warning: "+zefile+" not in "+whereset+" ; no presets."
    4243
    4344##################################
    … …  
    135136                      stridez=1,stridet=1,\
    136137                      compute="mean",\
    137                       verbose=False,noproj=False,\
     138                      verbose=False,\
     139                      quiet=False,\
     140                      noproj=False,\
    138141                      superpose=False,\
    139142                      plotin=None,\
    … …  
    154157                      label=None,\
    155158                      changetime=None,\
     159                      units=None,\
    156160                      title=None):
    157161        self.request = None
     162        self.nrequest = 0
    158163        self.nfin = 0 ; self.nvin = 0
    159164        self.nplotx = None ; self.nploty = None
    … …  
    165170        self.p = None
    166171        self.customplot = False
    167         self.allfield = None
     172        self.f = None
     173        self.l = None
    168174        ## what could be defined by the user
    169175        self.file = file
    … …  
    177183        self.compute = compute
    178184        self.verbose = verbose
     185        self.quiet = quiet
    179186        self.noproj = noproj
    180187        self.plotin = plotin
    … …  
    198205        self.color = color
    199206        self.label = label
     207        self.units = units
    200208        self.title = title
    201209
    202210    # print status
    203211    def printstatus(self):
     212      if not self.quiet:
    204213        if self.filename == "THIS_IS_A_CLONE":
    205214            pass
    … …  
    243252            self.color = other.color
    244253            self.label = other.label
     254            self.units = other.units
    245255            self.title = other.title
    246256            self.includedate = other.includedate
    … …  
    319329              if not isnum:   
    320330                  ope = other.request[i][j][t][z][y][x].field
    321                   if obj_ref.field.shape != ope.shape:
    322                     print "!! ERROR !! The two fields for operation do not have the same shape.",self.field.shape,other.field.shape
     331                  if ope.ndim == 0:
     332                    ope = float(ope) # if no dimension then this means that ope is a single value (not to be kept as an array)
     333                  elif obj_ref.field.shape != ope.shape:
     334                    print "!! ERROR !! The two fields for operation do not have the same shape.",obj_ref.field.shape,ope.shape
    323335                    exit()
    324336              else:           
    … …  
    347359              if not isnum:
    348360                  ope = other.request[i][j][t][z][y][x].field
    349                   if obj_ref.field.shape != ope.shape:
    350                     print "!! ERROR !! The two fields for operation do not have the same shape.",self.field.shape,other.field.shape
     361                  if ope.ndim == 0:
     362                    ope = float(ope) # if no dimension then this means that ope is a single value (not to be kept as an array)
     363                  elif obj_ref.field.shape != ope.shape:
     364                    print "!! ERROR !! The two fields for operation do not have the same shape.",obj_ref.field.shape,ope.shape
    351365                    exit()
    352366              else:
    … …  
    375389              if not isnum:
    376390                  ope = other.request[i][j][t][z][y][x].field
    377                   if obj_ref.field.shape != ope.shape:
    378                     print "!! ERROR !! The two fields for operation do not have the same shape.",self.field.shape,other.field.shape
     391                  if ope.ndim == 0:
     392                    ope = float(ope) # if no dimension then this means that ope is a single value (not to be kept as an array)
     393                  elif obj_ref.field.shape != ope.shape:
     394                    print "!! ERROR !! The two fields for operation do not have the same shape.",obj_ref.field.shape,ope.shape
    379395                    exit()
    380396              else:
    … …  
    403419              if not isnum:
    404420                  ope = other.request[i][j][t][z][y][x].field
    405                   if obj_ref.field.shape != ope.shape:
    406                     print "!! ERROR !! The two fields for operation do not have the same shape.",self.field.shape,other.field.shape
     421                  if ope.ndim == 0:
     422                    ope = float(ope) # if no dimension then this means that ope is a single value (not to be kept as an array)
     423                  elif obj_ref.field.shape != ope.shape:
     424                    print "!! ERROR !! The two fields for operation do not have the same shape.",obj_ref.field.shape,ope.shape
    407425                    exit()
    408426              else:
    … …  
    526544                       for z in range(self.nplotz)] for t in range(self.nplott)] \
    527545                       for j in range(self.nvin)]   for i in range(self.nfin)]
     546        # store how many onerequest() objects are in self.request
     547        self.nrequest = self.nfin*self.nvin*self.nplotx*self.nploty*self.nplotz*self.nplott
    528548        # loop on onerequest() objects
    529549        for i in range(self.nfin):
    … …  
    571591        # check if things were done OK before
    572592        if self.status != "defined": print "!! ERROR !! Please use .define() to define your pp object." ; exit()
    573         # create the list of allfield() objects
    574         # --> so that the user can easily access values
    575         self.allfield = [[[[[[ \
    576                         [] \
    577                         for x in range(self.nplotx)] for y in range(self.nploty)] \
    578                         for z in range(self.nplotz)] for t in range(self.nplott)] \
    579                         for j in range(self.nvin)]   for i in range(self.nfin)]
     593        ## create the list of f() and l() objects
     594        ## --> so that the user can easily access values (and labels for easy exploration)
     595        ## --> see example easy_get_field
     596        self.f = [ [] for iii in range(self.nrequest) ]
     597        self.l = [ [] for iii in range(self.nrequest) ]
     598        count = 0
    580599        ## first get fields
    581600        ## ... only what is needed is extracted from the files
    … …  
    590609              obj.getfield()
    591610              obj.computations()
    592               self.allfield[i][j][t][z][y][x] = obj.field
     611              # save fields in self.f for the user
     612              self.f[count] = obj.field
     613              # save a label in self.l for the user
     614              self.l[count] = "_"
     615              if self.nfin > 1:   self.l[count] = self.l[count] + "f=#"+str(int(i+1))+'_'
     616              if self.nvin > 1:   self.l[count] = self.l[count] + "v="+obj.var+'_'
     617              if self.nplotx > 1: self.l[count] = self.l[count] + "x="+str(self.x[x])+'_'
     618              if self.nploty > 1: self.l[count] = self.l[count] + "y="+str(self.y[y])+'_'
     619              if self.nplotz > 1: self.l[count] = self.l[count] + "z="+str(self.z[z])+'_'
     620              if self.nplott > 1: self.l[count] = self.l[count] + "t="+str(self.t[t])+'_'
     621              count = count + 1
     622        ## make it simple: self.f is simply the data array if self.nrequest=1
     623        if self.nrequest == 1: self.f = self.f[0]
    593624        # change status
    594625        self.status = "retrieved"
    … …  
    602633        self.retrieve()
    603634        return self 
     635
     636    ###########################################################
     637    # getf: a shortcut method for the define + retrieve chain #
     638    #       ... in which the output is self.f                 #
     639    #       ... and the ppclass is kept quiet                 #
     640    ###########################################################
     641    def getf(self):
     642        self.quiet = True
     643        self.get()
     644        return self.f
     645
     646    ############################################################
     647    # getfl: a shortcut method for the define + retrieve chain #
     648    #       ... in which the output is self.f, self.l          #
     649    #       ... and the ppclass is kept quiet                 #
     650    ############################################################
     651    def getfl(self):
     652        self.quiet = True
     653        self.get()
     654        return self.f,self.l
    604655
    605656    ########################################
    … …  
    730781                    if self.ycoeff is not None: plobj.ycoeff = self.ycoeff
    731782                    if self.title is not None: plobj.title = self.title
     783                    if self.units is not None: plobj.units = self.units
    732784                    # -- 1D specific
    733785                    if dp == 1:
    … …  
    799851            if self.verbose: print "**** OK. expect %i plots" % (self.howmanyplots)
    800852        else:
    801             exit() # because this means that we only had 0D values !
     853            pass # because this means that we only had 0D values !
    802854        # final status
    803855        self.status = "definedplot"
    … …  
    813865    ##############################################################################################
    814866    def makeplot(self):
     867      if self.howmanyplots > 0:
    815868        self.printstatus()
    816869        # a few initial operations
    … …  
    888941            pickle.dump(self.p, filehandler)
    889942        except IOError:
    890             print "!! WARNING !! Saved object file not written. Probably do not have permission to write here."
     943            if self.verbose: print "!! WARNING !! Saved object file not written. Probably do not have permission to write here."
    891944        return self
    892945
    … …  
    915968        return self
    916969
    917     ##############################################################
    918     # f: operation on two pp objects being on status 'definedplot'
     970    #################################################################
     971    # func: operation on two pp objects being on status 'definedplot'
    919972    # this allows for one field being function of another one
    920     # e.g. u.f(v) means u will be displayed as a function of v
    921     # ... no need to do defineplot after u.f(v), makeplot directly
    922     ##############################################################
    923     def f(self,other):
     973    # e.g. u.func(v) means u will be displayed as a function of v
     974    # ... no need to do defineplot after u.func(v), makeplot directly
     975    #################################################################
     976    def func(self,other):
    924977        # preamble: for this operation to work, defineplot() must have been done
    925978        if self.status != "definedplot":
    … …  
    11491202          if "grid points" not in self.name_x:
    11501203            if self.field_x.all() == self.field_x[0,0]:
    1151                print "!! WARNING !! xy axis look undefined. creating non-dummy ones."
     1204               if self.verbose: print "!! WARNING !! xy axis look undefined. creating non-dummy ones."
    11521205               self.field_x = np.array(range(self.dim_x)) ; self.name_x = "x grid points"
    11531206               self.field_y = np.array(range(self.dim_y)) ; self.name_y = "y grid points"
    … …  
    12421295                            self.field_t = self.field_t % 24 # so that the user is not mistaken!
    12431296            else:
    1244                 print "!! WARNING !! This time change is not implemented. Nothing is done."
     1297                if self.verbose: print "!! WARNING !! This time change is not implemented. Nothing is done."
    12451298            if self.verbose: print "**** OK. new t axis values [%5.1f,%5.1f]" % (self.field_t.min(),self.field_t.max())
    12461299
    … …  
    13961449    def getfield(self):
    13971450        ## first tell what is to be done
    1398         if self.dimplot > 2:                       print "**** !! ERROR !! "+str(self.dimplot)+"D plots not supported!" ; exit()
    1399         elif self.dimplot == 0 and self.verbose:   print "**** OK. 0D value requested."
    1400         elif self.dimplot == 1 and self.verbose:   print "**** OK. 1D plot requested."
    1401         elif self.verbose:                         print "**** OK. 2D section requested."
     1451        if self.verbose:
     1452          if self.dimplot > 2:                       print "**** !! WARNING !! "+str(self.dimplot)+"D plots will not be supported!"
     1453          elif self.dimplot == 0 and self.verbose:   print "**** OK. 0D value requested."
     1454          elif self.dimplot == 1 and self.verbose:   print "**** OK. 1D plot requested."
     1455          elif self.verbose:                         print "**** OK. 2D section requested."
    14021456        # well, now get field from netcdf file
    14031457        # part below is necessary otherwise there is an index error below

  • trunk/UTIL/PYTHON/planetoplot_v2/ppplot.py

    r947 r960  
    534534            print "!! WARNING !! Transposing axes"
    535535            self.field = np.transpose(self.field)
     536        # set dummy xy axis if not defined
     537        if self.absc is None:
     538            self.absc = np.array(range(self.field.shape[0]))
     539            self.mapmode = False
     540            print "!! WARNING !! dummy coordinates on x axis"
     541        if self.ordi is None:
     542            self.ordi = np.array(range(self.field.shape[1]))
     543            self.mapmode = False
     544            print "!! WARNING !! dummy coordinates on y axis"
    536545        # bound field
    537546        zevmin, zevmax = calculate_bounds(self.field,vmin=self.vmin,vmax=self.vmax)
Note: See TracChangeset for help on using the changeset viewer.