Changeset 921 for trunk/UTIL/PYTHON/planetoplot_v2

Timestamp:

Mar 30, 2013, 11:21:01 AM (12 years ago)

Author:

aslmd

Message:

UTIL PYTHON planetoplot_v2. added meanarea possibility for averages ponderate with areas. see aire.py (and set_ppclass.txt for presets). corrected a bug with customplot. separated getfield from computations.

Location:

trunk/UTIL/PYTHON/planetoplot_v2

Files:

• examples/aire.py (added)
• pp.py (modified) (1 diff)
• ppclass.py (modified) (18 diffs)
• set_ppclass.txt (modified) (1 diff)

trunk/UTIL/PYTHON/planetoplot_v2/pp.py

@@ -36 +36 @@
 parser.add_option('-o','--output',action='store',dest='filename',type="string",default="myplot",help="name of output files")
 parser.add_option('-d','--directory',action='store',dest='folder',type="string",default="./",help="directory of output files")
+parser.add_option('-u','--compute',action='store',dest='compute',type="string",default="mean",help='')
 # plot --> upper case
 # -- generic

trunk/UTIL/PYTHON/planetoplot_v2/ppclass.py

@@ -33 +33 @@
 zefile = "set_ppclass.txt"
 glob_listx = [] ; glob_listy = [] ; glob_listz = [] ; glob_listt = []
+glob_listarea = []
 try: 
     f = open(whereset+zefile, 'r') ; lines = f.readlines()
@@ -39 +40 @@
     for stuff in lines[11].strip().split(';'): glob_listz.append(stuff)
     for stuff in lines[14].strip().split(';'): glob_listt.append(stuff)
+    for stuff in lines[17].strip().split(';'): glob_listarea.append(stuff)
 except IOError: 
     print "warning: "+zefile+" not in "+whereset+" ; no presets."
@@ -151 +153 @@
         self.nplot = 0
         self.p = None
+        self.customplot = False
         ## what could be defined by the user
         self.file = file
@@ -497 +500 @@
               obj.method_z = mz ; obj.method_t = mt           
               ### get index
-              obj.getindextime(st,t,self.stridet)
-              obj.getindexvert(sz,z,self.stridez)
-              obj.getindexhori(sx,sy,x,y,self.stridex,self.stridey)
+              obj.getindextime(dalist=st,ind=t,stride=self.stridet)
+              obj.getindexvert(dalist=sz,ind=z,stride=self.stridez)
+              obj.getindexhori(dalistx=sx,dalisty=sy,indx=x,indy=y,stridex=self.stridex,stridey=self.stridey)
         # change status
         self.status = "defined"
@@ -518 +521 @@
         ## first get fields
         ## ... only what is needed is extracted from the files
-        ## ... and averages are computed
+        ## ... and computations are performed
         for i in range(self.nfin):
          for j in range(self.nvin):
@@ -527 +530 @@
               obj = self.request[i][j][t][z][y][x]
               obj.getfield()
+              obj.computations()
         # change status
         self.status = "retrieved"
@@ -731 +735 @@
             ## adapted space for labels etc
             ## ... except for ortho because there is no label anyway
-            customplot = self.p[0].field.ndim == 2 \
+            self.customplot = self.p[0].field.ndim == 2 \
                         and self.p[0].mapmode == True \
                         and self.p[0].proj not in ["ortho"]
-            if customplot:
+            if self.customplot:
                 margin = 0.07
                 self.fig.subplots_adjust(left=margin,right=1-margin,bottom=margin,top=1-margin)
@@ -744 +748 @@
             self.n = self.plotin.n
             self.howmanyplots = self.plotin.howmanyplots
+            self.customplot = self.plotin.customplot
         # LOOP on all subplots
         # NB: cannot use 'for pl in self.p' if self.plotin not None
@@ -780 +785 @@
         print "**** Done step: makeplot"
         if (self.n == self.howmanyplots): 
-            ppplot.save(mode=self.out,filename=self.filename,folder=self.folder,custom=customplot)
+            ppplot.save(mode=self.out,filename=self.filename,folder=self.folder,custom=self.customplot)
             mpl.close()
         # SAVE A PICKLE FILE WITH THE self.p ARRAY OF OBJECTS
@@ -979 +984 @@
             print '!! ERROR !! File '+self.file+' does not contain variable: '+self.var
             print '..... try instead with ',self.f.variables.keys() ; exit()
+
+    # copy attributes from another existing object
+    # --------------------------------------------
+    def copy(self,source):
+        for k, v in vars(source).items():
+            setattr(self,k,v)
 
     # get x,y,z,t dimensions from NETCDF file
@@ -1077 +1088 @@
     ### TBD: il faudrait ne prendre que les indices qui correspondent a l interieur d un plot (dans all)
     # -------------------------------
-    def getindextime(self,dalist,ind,stride):
+    def getindextime(self,dalist=None,ind=None,stride=1):
         if self.method_t == "free": 
             self.index_t = np.arange(0,self.dim_t,stride)
@@ -1101 +1112 @@
     ### TBD: il faudrait ne prendre que les indices qui correspondent a l interieur d un plot (dans all)
     # -------------------------------
-    def getindexvert(self,dalist,ind,stride):
+    def getindexvert(self,dalist=None,ind=None,stride=1):
         if self.method_z == "free": 
             self.index_z = np.arange(0,self.dim_z,stride)
@@ -1129 +1140 @@
     ### TBD: il faudrait ne prendre que les indices qui correspondent a l interieur d un plot (dans all)
     # -------------------------------
-    def getindexhori(self,dalistx,dalisty,indx,indy,stridex,stridey):
+    def getindexhori(self,dalistx=None,dalisty=None,indx=None,indy=None,stridex=1,stridey=1):
         ## get what is the method over x and y axis
         test = self.method_x+self.method_y
@@ -1309 +1320 @@
              self.field = masked
              self.field.set_fill_value([np.NaN])
-        # now ready to compute [TBD?? we would like to have e.g. mean over x,y and min over t]
+
+    # perform computations
+    # -------------------------------
+    # available: mean, max, min, meanarea
+    # TB: integrals? for derivatives, define a function self.dx()
+    def computations(self):  
+        nt,nz,ny,nx = self.field.shape
+        # treat the case of mean on fields normalized with grid mesh area
+        # ... this is done in the .area() method. 
+        # after that self.field contains field*area/totarea
+        if "area" in self.compute: self.area()
+        # now ready to compute [TBD: we would like to have e.g. mean over x,y and min over t ??]
         if self.method_t == "comp":
             if self.verbose: print "**** OK. Computing over t axis."
-            if self.compute == "mean": self.field = ppcompute.mean(self.field,axis=0)
+            if "mean" in self.compute: self.field = ppcompute.mean(self.field,axis=0)
             elif self.compute == "min": self.field = ppcompute.min(self.field,axis=0)
             elif self.compute == "max": self.field = ppcompute.max(self.field,axis=0)
@@ -1319 +1341 @@
         if self.method_z == "comp": 
             if self.verbose: print "**** OK. Computing over z axis."
-            if self.compute == "mean": self.field = ppcompute.mean(self.field,axis=1)
+            if "mean" in self.compute: self.field = ppcompute.mean(self.field,axis=1)
             elif self.compute == "min": self.field = ppcompute.min(self.field,axis=1)
             elif self.compute == "max": self.field = ppcompute.max(self.field,axis=1)
+            else: print "!! ERROR !! operation not supported." ; exit()
             nz = 1 ; self.field = np.reshape(self.field,(nt,nz,ny,nx))
         if self.method_y == "comp": 
@@ -1328 +1351 @@
             elif self.compute == "min": self.field = ppcompute.min(self.field,axis=2)
             elif self.compute == "max": self.field = ppcompute.max(self.field,axis=2)
+            elif self.compute == "meanarea": self.field = ppcompute.sum(self.field,axis=2)
+            else: print "!! ERROR !! operation not supported." ; exit()
             ny = 1 ; self.field = np.reshape(self.field,(nt,nz,ny,nx))
             if self.field_x.ndim == 2: self.field_x = self.field_x[0,:] # TBD: this is OK for regular grid but not for irregular
@@ -1335 +1360 @@
             elif self.compute == "min": self.field = ppcompute.min(self.field,axis=3)
             elif self.compute == "max": self.field = ppcompute.max(self.field,axis=3)
+            elif self.compute == "meanarea": self.field = ppcompute.sum(self.field,axis=3)
+            else: print "!! ERROR !! operation not supported." ; exit()
             nx = 1 ; self.field = np.reshape(self.field,(nt,nz,ny,nx))
             if self.field_y.ndim == 2: self.field_y = self.field_y[:,0] # TBD: this is OK for regular grid but not for irregular
@@ -1343 +1370 @@
         if self.verbose: 
             print "**** OK. Final shape for "+self.var+" after averaging and squeezing",self.field.shape
+    
+    # get areas for computations and ponderate field by those areas
+    # -------------------------------------------------------------
+    def area(self):
+        if self.verbose: print "**** OK. Get area array for computations."
+        # create a request object for area
+        # ... and copy known attributes from self
+        aire = onerequest()
+        aire.copy(self)
+        # get area field name
+        aire.var = "nothing"
+        for c in glob_listarea:
+         if c in aire.f.variables.keys():
+            aire.var = c
+        # do not try to calculate areas automatically 
+        if aire.var == "nothing":
+            print "!! ERROR !! area variable not found... needs to be added in set_ppclass.txt?"
+            exit()
+        # define area request dimensions
+        aire.getdim()
+        # ensure this is a 2D horizontal request and define indexes
+        #    ... areas are not supposed to vary with time and height
+        aire.method_x = "free" ; aire.method_y = "free"
+        aire.getindexhori() ; aire.dimplot = 2
+        aire.method_z = "fixed" ; aire.field_z = np.array([0]) ; aire.index_z = np.array([0])
+        aire.method_t = "fixed" ; aire.field_t = np.array([0]) ; aire.index_t = np.array([0])
+        # read the 2D area array in netCDF file
+        aire.getfield()
+        # normalize by total area
+        aire.field = np.squeeze(aire.field)
+        totarea = ppcompute.sum(ppcompute.sum(aire.field,axis=1),axis=0)
+        if self.verbose: print "**** OK. Total area is: ",totarea
+        aire.field = aire.field / totarea
+        # reduce with self horizontal indexes
+        if "fixed" in self.method_x+self.method_y:
+            aire.field = aire.field[self.index_y,self.index_x]
+        # tile area array over self t and z axis so that area field could be multiplied with self.field
+        aire.field = np.tile(aire.field,(self.index_t.size,self.index_z.size,1,1))
+        if self.field.shape != aire.field.shape:
+            print "!! ERROR !! Problem in area(). Check array shapes." ; exit()
+        else:
+            self.field = self.field*aire.field
 
     # define coordinates for plot

trunk/UTIL/PYTHON/planetoplot_v2/set_ppclass.txt

@@ -14 +14 @@
   #
 Times;time;Time;time_counter
+  # area: possible names to search for. add in the line with a separating ;
+  #
+area;aire