Changeset 409 for trunk/UTIL

Timestamp:

Nov 22, 2011, 3:14:32 PM (14 years ago)

Author:

acolaitis

Message:

PYTHON. 1/ mcs.py now works on a list of 3D and 4D variables given by -v or --var (it will also put aps and bps in the diagfi). 2/ added option --intas which is a template for the -i 2 gcm interpolation. -i 2 --intas mcs will interpolate the field on the same pressure levels as mcs. Works also with --intas tes. 3/ Corrected a small bug in make_gcm_netcdf

Location:

trunk/UTIL/PYTHON

Files:

• make_netcdf.py (modified) (1 diff)
• mcs.py (modified) (9 diffs)
• myplot.py (modified) (1 diff)
• myscript.py (modified) (1 diff)
• pp.py (modified) (2 diffs)
• zrecast_wrapper.py (modified) (2 diffs)

trunk/UTIL/PYTHON/make_netcdf.py

@@ -28 +28 @@
       zznames=[[]]*nvar
       for zvar in zvariables:
-          zznames[i]="var"+np.string(i+1)
+          zznames[i]="var"+str(i+1)
           i=i+1
    

trunk/UTIL/PYTHON/mcs.py

@@ -20 +20 @@
    import numpy as np
    from mymath import find_nearest 
-   from myplot import getfield
+   from myplot import getfield,separatenames
    from make_netcdf import make_gcm_netcdf
+   from zrecast_wrapper import call_zrecast
    parser = OptionParser() 
 
@@ -28 +29 @@
    parser.add_option('-f', '--file',   action='store',dest='file',     type="string",  default=None,  help='[NEEDED] filename.')
    parser.add_option('-m', '--mfile',  action='store',dest='mcsfile',     type="string",  default=None,  help='[NEEDED] filename for MCS comparison.')
+   parser.add_option('-v', '--var',    action='append',dest='var',      type="string",  default=None,  help='[NEEDED] Variables to process. (coma-separated list. aps and bps are always included.)')
 
    #############################
@@ -75 +77 @@
    ntimemax.set_fill_value([np.NaN])
 
-   # Variables
-
-   temp=getfield(nc,"temp")
-   dimensions = np.array(temp).shape
-   print "temp"
-   print dimensions
-   #aps=nc.variables["aps"][:]
-   #bps=nc.variables["bps"][:]
-
+   # Variables to treat
+
+   varz=[]
+   varznames=separatenames(opt.var[0])
+   n=0
+   for zn in varznames:
+       varz.append(getfield(nc,zn))
+       print "Found: "+zn+" with dimensions: "
+       print np.array(varz[n]).shape
+       n=n+1
+
+   nzvar=len(varz)
+   dimensions={}
+   vv=0
+   for var in varz:
+       a=len(np.array(var).shape)
+       if a == 3: dimensions[vv]=a
+       elif a == 4: dimensions[vv]=a
+       else:
+          print "Warning, only 3d and 4d variables are supported for time-recasting"
+          exit()
+       vv=vv+1
+
+   # Variables to save but not treated (only along z, or phisinit-like files)
+
+   aps=nc.variables["aps"][:]
+   bps=nc.variables["bps"][:]
+   fullnames=["aps","bps"]
+   for name in varznames:
+       fullnames.append("d"+name)
+       fullnames.append("n"+name)
+   print "Will output: "
+   print fullnames
    #############################
    ### Building
@@ -92 +118 @@
    varday=np.zeros([len(timemcs),nz,ny,nx])
    varnight=np.zeros([len(timemcs),nz,ny,nx])
-   vardayout=np.ma.masked_invalid(varday)
-   varnightout=np.ma.masked_invalid(varnight)
+   vardayout=np.zeros([nzvar,len(timemcs),nz,ny,nx])
+   varnightout=np.zeros([nzvar,len(timemcs),nz,ny,nx])
+   vardayout=np.ma.masked_invalid(vardayout)
+   varnightout=np.ma.masked_invalid(varnightout)
    i=0
    for ls in timemcs:
@@ -100 +128 @@
        istart=find_nearest(lstime,lsstart,strict=True)
        istop=find_nearest(lstime,lsstop,strict=True)
+       varchk=0
        if ((istart is np.NaN) or (istop is np.NaN)):
-          vardayout[i,:,:,:]=np.NaN
-          varnightout[i,:,:,:]=np.NaN
+          vardayout[:,i,:,:,:]=np.NaN
+          varnightout[:,i,:,:,:]=np.NaN
           print "Time interval skipped. Ls MCS: (",lsstart,';',lsstop,')',"// Ls diagfi: (",lstime.min(),';',lstime.max(),')'
           i=i+1
           continue
        print "--->> Processing Data. Ls MCS: (",lsstart,';',lsstop,')',"// Ls diagfi: (",lstime.min(),';',lstime.max(),')'
-
-       varchk=temp[istart:istop,:,:,:]
-       varchktime=time[istart:istop]
+       # warning, python's convention is that the second index of array[a:b] is the array index of element after the one being picked last.
+       # for that reason, array[0:0] is nan and array[0:1] is only one value. Hence, len(array[a:b+1]) is b-a+1 and not b-a+2
+       print "     .initialisation."
+       
+       
+       varchk=np.zeros([nzvar,istop-istart+1,nz,ny,nx])
+       vv=0
+       for variable in varz:
+           if dimensions[vv] is 3:
+              varchk[vv,:,0,:,:]=variable[istart:istop+1,:,:]
+           else:
+              varchk[vv,:,:,:,:]=variable[istart:istop+1,:,:,:]
+           vv=vv+1
+       varchktime=time[istart:istop+1]
        ndays=np.floor(varchktime[len(varchktime)-1])-np.floor(varchktime[0])
        dtmichk=dtimemin[i,:,:]
@@ -129 +169 @@
    ### (yea it's complicated)
 
-       vartmpnight=np.zeros([ndays,nz,ny,nx])
-       vartmpday=np.zeros([ndays,nz,ny,nx])
+       vartmpnight=np.zeros([nzvar,ndays,nz,ny,nx])
+       vartmpday=np.zeros([nzvar,ndays,nz,ny,nx])
        nd=0
-       
-
+       print "     .time indices MCS grid -> diagfi grid."
        while nd < ndays:
 
@@ -152 +191 @@
                 nitntmichk=find_nearest(varchktime_lon,ntmichk[niy,nix])
                 nitntmachk=find_nearest(varchktime_lon,ntmachk[niy,nix])
-                if ((nitdtmichk is np.NaN) or (nitdtmachk is np.NaN)):
-                    vartmpday[nd,:,iy,ix]=np.NaN
-                elif nitdtmichk > nitdtmachk:
-                    vartmpday[nd,:,iy,ix]=(np.mean(varchk[daystart+nitdtmichk:daystop+1,:,iy,ix],axis=0)+np.mean(varchk[daystart:daystart+nitdtmachk+1,:,iy,ix],axis=0))/2.
-                else:
-                    vartmpday[nd,:,iy,ix]=np.mean(varchk[daystart+nitdtmichk:daystart+nitdtmachk+1,:,iy,ix],axis=0)
-                if ((nitntmichk is np.NaN) or (nitntmachk is np.NaN)):
-                    vartmpnight[nd,:,iy,ix]=np.NaN
-                elif nitntmichk > nitntmachk:
-                    vartmpnight[nd,:,iy,ix]=(np.mean(varchk[daystart+nitntmichk:daystop+1,:,iy,ix],axis=0)+np.mean(varchk[daystart:daystart+nitntmachk+1,:,iy,ix],axis=0))/2.
-                else:                                                            
-                    vartmpnight[nd,:,iy,ix]=np.mean(varchk[daystart+nitntmichk:daystart+nitntmachk+1,:,iy,ix],axis=0)
-
+                for vv in np.arange(nzvar):
+                   if ((nitdtmichk is np.NaN) or (nitdtmachk is np.NaN)):
+                       vartmpday[vv,nd,:,iy,ix]=np.NaN
+                   elif nitdtmichk > nitdtmachk:
+                       vartmpday[vv,nd,:,iy,ix]=(np.mean(varchk[vv,daystart+nitdtmichk:daystop+1,:,iy,ix],axis=0)+np.mean(varchk[vv,daystart:daystart+nitdtmachk+1,:,iy,ix],axis=0))/2.
+                   else:
+                       vartmpday[vv,nd,:,iy,ix]=np.mean(varchk[vv,daystart+nitdtmichk:daystart+nitdtmachk+1,:,iy,ix],axis=0)
+                   if ((nitntmichk is np.NaN) or (nitntmachk is np.NaN)):
+                       vartmpnight[vv,nd,:,iy,ix]=np.NaN
+                   elif nitntmichk > nitntmachk:
+                       vartmpnight[vv,nd,:,iy,ix]=(np.mean(varchk[vv,daystart+nitntmichk:daystop+1,:,iy,ix],axis=0)+np.mean(varchk[vv,daystart:daystart+nitntmachk+1,:,iy,ix],axis=0))/2.
+                   else:                                                            
+                       vartmpnight[vv,nd,:,iy,ix]=np.mean(varchk[vv,daystart+nitntmichk:daystart+nitntmachk+1,:,iy,ix],axis=0)
                 iy=iy+1
              ix=ix+1
           nd=nd+1
 
+       print "     .creating bins."
+
        vartmpdaymasked=np.ma.masked_invalid(vartmpday)
        vartmpnightmasked=np.ma.masked_invalid(vartmpnight)
-
-       vardayout[i,:,:,:]=np.ma.mean(vartmpdaymasked[:,:,:,:],axis=0)
-       varnightout[i,:,:,:]=np.ma.mean(vartmpnightmasked[:,:,:,:],axis=0)
+       for vv in np.arange(nzvar):
+          vardayout[vv,i,:,:,:]=np.ma.mean(vartmpdaymasked[vv,:,:,:,:],axis=0)
+          varnightout[vv,i,:,:,:]=np.ma.mean(vartmpnightmasked[vv,:,:,:,:],axis=0)
+          print "          ."+varznames[vv]+".done"
        i=i+1
+
+   print "--->> Preparing Data for ncdf. Missing value is NaN."
 
    vardayout=np.ma.masked_invalid(vardayout)
@@ -181 +225 @@
    varnightout.set_fill_value([np.NaN])
 
-   vardayncdf=vardayout.filled()
-   varnightncdf=varnightout.filled()
+   all=[aps,bps]
+   for vv in np.arange(nzvar):
+       if dimensions[vv] == 3:
+          all.append(vardayout[vv,:,0,:,:].filled())
+          all.append(varnightout[vv,:,0,:,:].filled())
+       elif dimensions[vv] == 4:
+          all.append(vardayout[vv,:,:,:,:].filled())
+          all.append(varnightout[vv,:,:,:,:].filled())
 
    make_gcm_netcdf (zfilename="diagfi_MCS.nc", \
@@ -190 +240 @@
                         zalt=alt, \
                         ztime=timemcs, \
-                        zvariables=[vardayncdf,varnightncdf], \
-                        znames=["dtemp","ntemp"])
+                        zvariables=all, \
+                        znames=fullnames)

trunk/UTIL/PYTHON/myplot.py

@@ -664 +664 @@
     fmtvar    =     { \
              "TK":           "%.0f",\
+# Variables from TES ncdf format
              "T_NADIR_DAY":  "%.0f",\
              "T_NADIR_NIT":  "%.0f",\
+# Variables from tes.py ncdf format
              "TEMP_DAY":     "%.0f",\
              "TEMP_NIGHT":   "%.0f",\
+# Variables from MCS and mcs.py ncdf format
+             "DTEMP":     "%.0f",\
+             "NTEMP":   "%.0f",\
              "TPOT":         "%.0f",\
              "TSURF":        "%.0f",\

trunk/UTIL/PYTHON/myscript.py

@@ -20 +20 @@
     parser.add_option('-l', '--level',  action='store',dest='lvl',       type="string",  default="0",   help='level / start,stop,step (-i 2: p,mb)(-i 3,4: z,km) [0]')
     parser.add_option('-i', '--interp', action='store',dest='itp',       type="int",     default=None,  help='interpolation (2: p, 3: z-amr, 4:z-als)')
+    parser.add_option('--intas',        action='store',dest='intas',     type="string",  default=None,  help='specify "mcs" or "tes" for specific gcm pressure interpolation grid')
     parser.add_option('-N', '--no-api', action='store_true',dest='nocall',               default=False, help='do not recreate api file')
 

trunk/UTIL/PYTHON/pp.py

@@ -105 +105 @@
             interpolated_files=call_zrecast(interp_mode=opt.itp,\
                     input_name=zefiles,\
-                    fields=zevars)
+                    fields=zevars,\
+                    predifined=opt.intas)
 
             zefiles=interpolated_files
@@ -112 +113 @@
                interpolated_ref=call_zrecast(interp_mode=opt.itp,\
                     input_name=[opt.fref],\
-                    fields=zevars)
+                    fields=zevars,\
+                    predifined=opt.intas)
 
                reffile=interpolated_ref[0]

trunk/UTIL/PYTHON/zrecast_wrapper.py

@@ -4 +4 @@
 def call_zrecast (  interp_mode   = '4', \
                     input_name      = None, \
-                    fields  = 'all'):
+                    fields  = 'all', \
+                    predifined = None):
 
     import numpy as np
     from myplot import separatenames
     from os import system
-
+    pressure_axis_tes=[1658.152,1291.37,1005.72,783.2555,610.,475.0685,369.9837,288.1436,224.4065,174.7679,136.1094,106.0021,82.55452,64.29353,50.07185,38.99599,30.37011,23.65227,18.4204,14.34583,11.17254]
+    pressure_axis_mcs=[1878.9, 1658.2, 1463.3, 1291.4, 1139.6, 1005.7, 887.54, 783.26,691.22, 610, 538.32, 475.07, 419.25, 369.98, 326.51, 288.14, 254.29,224.41, 198.04, 174.77, 154.23, 136.11, 120.12, 106, 93.547, 82.555,72.854, 64.294, 56.739, 50.072, 44.188, 38.996, 34.414, 30.37, 26.802,23.652, 20.873, 18.42, 16.256, 14.346, 12.66, 11.173, 9.8597, 8.7012,7.6788, 6.7765, 5.9802, 5.2775, 4.6574, 4.1101, 3.6272, 3.201, 2.8249,2.4929, 2.2, 1.9415, 1.7134, 1.512, 1.3344, 1.1776, 1.0392, 0.9171,0.80934, 0.71424, 0.63031, 0.55625, 0.49089, 0.43321, 0.3823, 0.33738,0.29774, 0.26275, 0.23188, 0.20463, 0.18059, 0.15937, 0.14064, 0.12412,0.10953, 0.096661, 0.085303, 0.07528, 0.066434, 0.058628, 0.051739,0.04566, 0.040294, 0.03556, 0.031381, 0.027694, 0.02444, 0.021568,0.019034, 0.016797, 0.014824, 0.013082, 0.011545, 0.010188, 0.0089909,0.0079345, 0.0070021, 0.0061794, 0.0054533, 0.0048125, 0.004247]
     system("rm -f zrecast.auto.def")
     system("touch zrecast.auto.def")
@@ -26 +28 @@
         elif interp_mode == 2:
              append="_P"
-             system("echo 1 >> zrecast.auto.def")
-             system("echo yes >> zrecast.auto.def")
-             system("echo 370 0.1 >> zrecast.auto.def")  #I put that randomly! (a.c.)
-             system("echo 20 >> zrecast.auto.def")
+             if predifined is "TES":
+                system("echo 1 >> zrecast.auto.def")
+                system("echo no >> zrecast.auto.def")
+                system("echo "+len(pressure_axis_tes)+" >> zrecast.auto.def")
+                for zp in pressure_axis_tes:
+                   system("echo "+zp+" >> zrecast.auto.def")
+             elif predifined is "MCS":
+                system("echo no >> zrecast.auto.def")
+                system("echo "+len(pressure_axis_mcs)+" >> zrecast.auto.def")
+                for zp in pressure_axis_mcs:
+                   system("echo "+zp+" >> zrecast.auto.def")
+             else:
+                system("echo 1 >> zrecast.auto.def")
+                system("echo yes >> zrecast.auto.def")
+                system("echo 370 0.1 >> zrecast.auto.def")  #I put that randomly! (a.c.)
+                system("echo 20 >> zrecast.auto.def")
         else:
              print "zrecast interp option unsupported for now. Exiting."