Changeset 458 in lmdz_wrf for trunk/tools

Timestamp:

Jun 5, 2015, 3:10:26 PM (10 years ago)

Author:

lfita

Message:

Adding `read_ISD.py' to read observations in ASCII in ISD format (single no spaces line) to netCDF
Adding `wdismean' to compute mean distance weighted values
Adding `radial_points' to provide grid points following an angle and radii
Adding (NOT WORKING yet) `compute_tevolboxtraj_radialsec' to retrieve radial sections of values following a trajectory
l

Location:

trunk/tools

Files:

• nc_var_tools.py (modified) (1 diff)
• read_ISD.py (added)

trunk/tools/nc_var_tools.py

@@ -14864 +14864 @@
     return tB
 
+def wdismean(pos,val4):
+    """ Function to compute the mean value weighted to its 4 distances
+      pos= x,y position within the 'square'
+      val4= values at the four vertexs (2,2)
+    >>>position = [0.005,0.005]
+    >>>val4 = np.zeros((2,2), dtype=np.float)
+    >>>val4 = np.arange(4).reshape(2,2)
+    0.035707955234
+    """
+    fname = 'wdismean'
+
+    dist = np.zeros((2,2), dtype=np.float)    
+    dist[0,0] = np.sqrt(pos[0]*pos[0]+pos[1]*pos[1])
+    dist[0,1] = np.sqrt(pos[0]*pos[0]+(1.-pos[1])*(1.-pos[1]))
+    dist[1,0] = np.sqrt((1.-pos[0])*(1.-pos[0])+pos[1]*pos[1])
+    dist[1,1] = np.sqrt((1.-pos[0])*(1.-pos[0])+(1.-pos[1])*(1.-pos[1]))
+
+    if np.min(dist) == 0.:
+        pos0 = index_mat(dist, 0.)
+        dist[:,:] = 0.
+        posmean = val4[pos0[0], pos0[1]]
+    else:
+        totdist = np.sum(1./dist)
+        posmean = np.sum(val4.flatten()/dist.flatten())/totdist
+
+    return posmean
+
+def radial_points(angle,dist):
+    """ Function to provide a number of grid point positions for a given angle
+      angle= desired angle (rad)
+      dist= number of grid points
+    >>> radial_points(0.5*np.pi,5)
+    [[  6.12323400e-17   1.00000000e+00]
+     [  1.22464680e-16   2.00000000e+00]
+     [  1.83697020e-16   3.00000000e+00]
+     [  2.44929360e-16   4.00000000e+00]
+     [  3.06161700e-16   5.00000000e+00]]
+    """
+    fname = 'radial_points'
+    radpos = np.zeros((dist,2), dtype=np.float)
+
+    for ip in range(dist):
+        radpos[ip,0] = (ip+1.)*np.cos(angle)
+        radpos[ip,1] = (ip+1.)*np.sin(angle)
+
+    return radpos
+
+def compute_tevolboxtraj_radialsec(values, ncfile, varn):
+    """ Function to compute tevolboxtraj_radialsec: temporal evolution at a given point along 
+          a number of radii at a given frequency following a trajectory
+      values= [trajfile]@[Tbeg],[lonname],[latname],[zname],[timename],[timekind],[Nangle],[radii]
+        [trajfile]: ASCII file with the trajectory ('#' not readed)
+          [time] [xpoint] [ypoint]
+        [Tbeg]: equivalent first time-step of the trajectory within the netCDF file
+        [lonname],[latname],[zname],[timename]: longitude, latitude, z and time variables names
+        [timekind]: kind of time
+          'cf': cf-compilant
+          'wrf': WRF kind
+        [Nangle]: Number of angles
+        [radii]: length in grtid points of the angle
+      ncfile= netCDF file to use
+      varn= ',' list of variables' name ('all', for all variables)
+    """
+    import numpy.ma as ma
+    import subprocess as sub
+
+    fname='compute_tevolboxtraj_radialsec'
+
+    if values == 'h':
+        print fname + '_____________________________________________________________'
+        print compute_tevolboxtraj_radialsec.__doc__
+        quit()
+
+    arguments = '[trajfile]@[Tbeg],[lonname],[latname],[zname],[timename],' +        \
+      '[timekind],[Nangle],[radii]]'
+    check_arguments(fname,len(arguments.split(',')),arguments,',',8)
+
+    trajfile = values.split(',')[0].split('@')[0]
+    Tbeg = int(values.split(',')[0].split('@')[1])
+    lonn = values.split(',')[1]
+    latn = values.split(',')[2]
+    zn = values.split(',')[3]
+    timn = values.split(',')[4]
+    timekind = values.split(',')[5]
+    Nangle = int(values.split(',')[6])
+    radii = int(values.split(',')[7])
+
+    if not os.path.isfile(ncfile):
+        print errormsg
+        print '  ' + fname + " netCDF file: '" + ncfile + "' does not exist !!"
+        quit(-1)
+
+    if not os.path.isfile(trajfile):
+        print errormsg
+        print '  ' + fname + " trajectory file: '" + trajfile + "' does not exist !!"
+        quit(-1)
+
+    objfile = NetCDFFile(ncfile, 'r')
+    lonobj = objfile.variables[lonn]
+    latobj = objfile.variables[latn]
+    timobj = objfile.variables[timn]
+
+    if varn.find(',') != -1:
+        varns = varn.split(',')
+    else:
+        if varn == 'all':
+            varns = objfile.variables
+        else:
+            varns = [varn]
+
+    lonv, latv = lonlat2D(lonobj[:],latobj[:])
+
+    dimx = lonv.shape[1]
+    dimy = lonv.shape[0]
+
+    timv = timobj[:]
+
+# Selecting accordingly a trajectory
+##
+    Ttraj = file_nlines(trajfile,'#')
+    if timekind == 'wrf':
+        dimt = objfile.variables[timn].shape[0] 
+    else:
+        dimt = objfile.variables[timn].shape
+
+    if Tbeg + Ttraj > dimt:
+        print errormsg
+        print '  ' + fname + ': trajectory has ', Ttraj, ' time steps and starts ' + \
+        ' at',Tbeg,' and data',varobj.shape[0], '" !!!!!'
+        quit(-1)
+
+    print '    ' + fname + ': Number of time-steps in trajectory file: ',Ttraj
+
+    trajobj = open(trajfile,'r')
+
+# Trajectory values/slices
+##
+    iline = 0
+    it = 0
+
+    xrangeslice = []
+    yrangeslice = []
+    xrangeslice2D = []
+    yrangeslice2D = []
+
+    cslicev = []
+    cslice2D = []
+    cslicevnoT = []
+
+    gtrajvals = np.zeros((Ttraj,3), dtype=int)
+    trajvals = np.zeros((Ttraj,3), dtype=np.float)
+
+    trajpos = np.zeros(Ttraj,Nangle,raddi)
+
+    it = 0
+    iline = 0
+    for line in trajobj:
+
+## Skipping first line
+##        if not iline == 0:
+##        it = iline - 1
+##        it = iline
+
+# Slicing brings to reduce 1 time-step.... ???
+        if line[0:1] != '#':
+            gtrajvals[it,0] = Tbeg + iline
+            gtrajvals[it,1] = int(line.split(' ')[1])
+            gtrajvals[it,2] = int(line.split(' ')[2])
+#            print it,'t:',gtrajvals[it,0],'x y:', gtrajvals[it,1], gtrajvals[it,2]
+        
+            if timekind == 'wrf':
+                gdate = datetimeStr_conversion(timv[gtrajvals[it,0]],'WRFdatetime',  \
+                  'matYmdHMS')
+                trajvals[it,0] = realdatetime1_CFcompilant(gdate, '19491201000000',  \
+                  'hours')
+                tunits = 'hours since 1949-12-01 00:00:00'
+            elif timekind == 'cf':
+                trajvals[it,0] = timv[gtrajvals[it,0]]
+                tunits = timobj.getncattr('units')
+            else:
+                print errormsg
+                print '  ' + fname + ' time kind: "' + timekind + '" not ready !!'
+                quit(-1)
+
+            trajvals[it,1] = lonv[gtrajvals[it,2],gtrajvals[it,1]]
+            trajvals[it,2] = latv[gtrajvals[it,2],gtrajvals[it,1]]
+
+#            print iline, it,'time:',trajvals[it,0],'lon lat:', trajvals[it,1],       \
+#              trajvals[it,2]
+
+# Assuming time as the first dimension in a fortran like way (t=0 as 1)
+#            trajcut[0] = tstept - 1
+# Assuming time as the first dimension in a C/python like way (t=0)
+
+            varvalst = np.ones((Nangle, radii), dtype=np.float)*fillValue
+            posangle = np.zeros((radii,2), dtype=np.float)
+
+# Angle loop
+            for ia in range(Nangle):
+                angle = 2.*np.pi*ia/Nangle
+
+# Points at that given angle
+                pangle = radial_points(angle,radii)
+                posangle[:,0] = pangle[:,0] + trajvals[it,1]
+                posangle[:,1] = pangle[:,1] + trajvals[it,2]
+
+                for ip in range(radii):
+                    iipos = int(posangle[ip,0])
+                    jjpos = int(posangle[ip,1])
+
+# Creation of a matrix with all the grid points at each time-step
+                    if iipos >= 0 and iipos < dimx and jjpos >= 0 and jjpos < dimy:
+                         print 'Hola'
+                         quit()
+#                        varvalst[ia,ip] = 
+
+
+                if gtrajvals[it,2]-radii < 0: 
+                    yinit = 0
+                    yinit2D = radii - gtrajvals[it,2]
+                else:
+                    yinit = gtrajvals[it,2]-radii
+                    yinit2D = 0
+
+                if gtrajvals[it,2]+radii + 1 > dimy + 1: 
+                    yend = dimy + 1
+                    yend2D = dimy + 1 - gtrajvals[it,2] + radii
+                else:
+                    yend = gtrajvals[it,2]+radii + 1
+                    yend2D = radi
+
+                if gtrajvals[it,1]-box2 < 0: 
+                    xinit = 0
+                    xinit2D = box2 - gtrajvals[it,1]
+                else:
+                    xinit = gtrajvals[it,1]-box2
+                    xinit2D = 0
+
+                if gtrajvals[it,1]+box2 + 1 > dimx + 1: 
+                    xend = dimx + 1
+                    xend2D = dimx + 1 - gtrajvals[it,1] - box2
+                else:
+                    xend = gtrajvals[it,1]+box2 + 1
+                    xend2D = boxs
+
+                yrangeslice.append([yinit, yend])
+                xrangeslice.append([xinit, xend])
+                yrangeslice2D.append([yinit2D, yend2D])
+                xrangeslice2D.append([xinit2D, xend2D])
+            else:
+                yrangeslice.append([gtrajvals[it,2]-box2, gtrajvals[it,2]+box2 + 1])
+                xrangeslice.append([gtrajvals[it,1]-box2, gtrajvals[it,1]+box2 + 1])
+                yrangeslice2D.append([gtrajvals[it,2]-box2, gtrajvals[it,2]+box2 + 1])
+                xrangeslice2D.append([gtrajvals[it,1]-box2, gtrajvals[it,1]+box2 + 1])
+
+# circle values
+            circdist[it,:,:] = radius_dist(dimy,dimx,gtrajvals[it,2],gtrajvals[it,1])
+
+            if gtrajvals[it,2]-Nrad < 0 or gtrajvals[it,2]+Nrad + 1 > dimy + 1 \
+              or gtrajvals[it,1]-Nrad < 0 or gtrajvals[it,1]+Nrad + 1 > dimx + 1:
+
+                if gtrajvals[it,2]-Nrad < 0: 
+                    yinit = 0
+                    yinit2D = Nrad - gtrajvals[it,2]
+                else:
+                    yinit = gtrajvals[it,2]-Nrad
+                    yinit2D = 0
+
+                if gtrajvals[it,2]+Nrad + 1 > dimy + 1: 
+                    yend = dimy + 1
+                    yend2D = dimy + 1 - gtrajvals[it,2] + Nrad
+                else:
+                    yend = gtrajvals[it,2]+Nrad + 1
+                    yend2D = 2*Nrad+1
+
+                if gtrajvals[it,1]-Nrad < 0: 
+                    xinit = 0
+                    xinit2D = Nrad - gtrajvals[it,1]
+                else:
+                    xinit = gtrajvals[it,1]-Nrad
+                    xinit2D = 0
+
+                if gtrajvals[it,1]+Nrad + 1 > dimx + 1: 
+                    xend = dimx + 1
+                    xend2D = dimx + 1 - gtrajvals[it,1] - Nrad
+                else:
+                    xend = gtrajvals[it,1]+Nrad + 1
+                    xend2D = 2*Nrad+1
+
+                cyrangeslice.append([yinit, yend])
+                cxrangeslice.append([xinit, xend])
+                cyrangeslice2D.append([yinit2D, yend2D])
+                cxrangeslice2D.append([xinit2D, xend2D])
+            else:
+                cyrangeslice.append([gtrajvals[it,2]-Nrad, gtrajvals[it,2]+Nrad + 1])
+                cxrangeslice.append([gtrajvals[it,1]-Nrad, gtrajvals[it,1]+Nrad + 1])
+                cyrangeslice2D.append([gtrajvals[it,2]-Nrad, gtrajvals[it,2]+Nrad+1])
+                cxrangeslice2D.append([gtrajvals[it,1]-Nrad, gtrajvals[it,1]+Nrad+1])
+            it = it + 1
+            iline = iline + 1
+
+    trajobj.close()
+# Creation of the netCDF file
+##
+    if varn.find(',') != -1:
+        varnS = 'multi-var'
+    else:
+        varnS = varn.replace(',','-')
+
+    ofile = 'tevolboxtraj_' + varnS + '.nc'
+    objofile = NetCDFFile(ofile, 'w')
+
+# Dimensions
+    newdim = objofile.createDimension('x', boxs)
+    newdim = objofile.createDimension('y', boxs)
+    newdim = objofile.createDimension('xr', Nrad*2+1)
+    newdim = objofile.createDimension('yr', Nrad*2+1)
+    newdim = objofile.createDimension('time', None)
+
+    stsn = ['min', 'max', 'mean', 'mean2', 'stdev', 'ac']
+    vstlname = ['minimum value within', 'maximum value within', 'mean value within', \
+      'squared mean value within', 'standard deviation value within',                \
+      'accumulated value within']
+    Nsts = len(stsn)
+    statnames = []
+    cstatnames = []
+    for i in range(Nsts):
+        statnames.append(stsn[i] + 'box')
+        cstatnames.append(stsn[i] + 'circle')
+
+# Getting values
+##
+    ivar = 0
+
+    for vn in varns:
+        vnst = variables_values(vn)[0]
+        if not searchInlist(objfile.variables, vn):
+            print errormsg
+            print '  ' + fname + ": variable name '" + vn + "' is not in file " +    \
+              ncfile + '" !!!!!'
+            quit(-1)
+
+        varobj = objfile.variables[vn]
+        Nvardims = len(varobj.shape)
+
+        print '  ' + fname + ": getting '" + vn + "' ... .. ." 
+
+        slicev = []
+        slice2D = []
+        slicevnoT = []
+        cslicev = []
+        cslice2D = []
+        cslicevnoT = []
+
+        if Nvardims == 4:
+# Too optimistic, but at this stage...
+            if not objofile.dimensions.has_key('z'):
+                varzobj = objfile.variables[zn]
+                if len(varzobj.shape) == 1:
+                    dimz = varzobj.shape
+                    zvals = varzobj[:]
+                elif len(varzobj.shape) == 2:
+                    dimz = varzobj.shape[1]
+                    zvals = varzobj[0,:]
+                elif len(varzobj.shape) == 3:
+                    dimz = varzobj.shape[2]
+                    zvals = varzobj[0,0,:]
+
+                objofile.createDimension('z', dimz)
+                newvar = objofile.createVariable(zn, 'f4', ('z'),fill_value=fillValue)
+                if searchInlist(varzobj.ncattrs(),'standard_name'):
+                    vsname = varzobj.getncattr('standard_name')
+                else:
+                    vsname = variables_values(zn)[1]
+                if searchInlist(varzobj.ncattrs(),'long_name'):
+                    vlname = varzobj.getncattr('long_name')
+                else:
+                    vlname = variables_values(zn)[4].replace('|',' ')
+                if searchInlist(varzobj.ncattrs(),'units'):
+                    vunits = varzobj.getncattr('units')
+                else:
+                    vunits = variables_values(zn)[5].replace('|',' ')
+    
+                newattr = basicvardef(newvar, vsname, vlname, vunits)
+                newvar[:] = zvals
+ 
+            varvals = np.ones(tuple([Ttraj,dimz,boxs,boxs]), dtype=np.float)
+            lonvals = np.ones(tuple([Ttraj,dimz,boxs,boxs]), dtype=np.float)
+            latvals = np.ones(tuple([Ttraj,boxs,boxs]), dtype=np.float)
+            rvarvals = np.ones(tuple([Ttraj,dimz,Nrad*2+1,Nrad*2+1]), dtype=np.float)
+            rlonvals = np.ones(tuple([Ttraj,dimz,Nrad*2+1,Nrad*2+1]), dtype=np.float)
+            rlatvals = np.ones(tuple([Ttraj,dimz,Nrad*2+1,Nrad*2+1]), dtype=np.float)
+
+# One dimension plus for the values at the center of the trajectory
+            statvarvals = np.ones(tuple([Ttraj,dimz,Nsts+1]), dtype=np.float)
+            rstatvarvals = np.ones(tuple([Ttraj,dimz,Nsts+1]), dtype=np.float)
+
+            for it in range(Ttraj):
+                it0 = Tbeg + it
+
+                slicev = []
+                slice2D = []
+                slicevnoT = []
+                cslicev = []
+                cslice2D = []
+                cslice2Dhor = []
+                cslicevnoT = []
+                cslicevnoThor = []
+
+                slicev.append(gtrajvals[it,0])
+                if gtrajvals[it,2]-box2 < 0 or gtrajvals[it,2]+box2 + 1 > dimy + 1   \
+                  or gtrajvals[it,1]-box2 < 0 or gtrajvals[it,1]+box2 + 1 > dimx +1:
+# box values
+                    slicev.append(slice(0,dimz))
+                    slicev.append(slice(yrangeslice[it][0],yrangeslice[it][1]))
+                    slicev.append(slice(xrangeslice[it][0],xrangeslice[it][1]))
+
+                    slicevnoT.append(slice(0,dimz))
+                    slicevnoT.append(slice(yrangeslice[it][0],yrangeslice[it][1]))
+                    slicevnoT.append(slice(xrangeslice[it][0],xrangeslice[it][1]))
+
+                    slice2D.append(slice(0,dimz))
+                    slice2D.append(slice(0,yrangeslice[it][1]-yrangeslice[it][0]))
+                    slice2D.append(slice(0,xrangeslice[it][1]-xrangeslice[it][0]))
+
+                    rvarvalst = np.ones((dimz, Nrad*2+1, Nrad*2+1),dtype=np.float)*  \
+                      fillValue
+
+                    varvalst[tuple(slice2D)] = varobj[tuple(slicev)]
+                    varvals[it,:,:,:] = varvalst
+    
+# box stats values
+                    maskedvals = ma.masked_values (varvalst, fillValue)
+                    maskedvals2 = maskedvals*maskedvals
+                    for iz in range(dimz):
+                        statvarvals[it,iz,0] = varvalst[iz,box2,box2]
+                        statvarvals[it,iz,1] = np.min(varvalst[iz,:,:])
+                        statvarvals[it,iz,2] = np.max(varvalst[iz,:,:])
+                        statvarvals[it,iz,3] = np.mean(varvalst[iz,:,:])
+                        statvarvals[it,iz,4] = maskedvals2[iz,:,:].mean()
+                        statvarvals[it,iz,5] = np.sqrt(statvarvals[it,iz,4] -        \
+                          statvarvals[it,iz,3]*statvarvals[it,iz,3])
+                        statvarvals[it,iz,6] = np.sum(varvalst[iz,:,:])
+
+                else:
+                    slicev.append(slice(0,dimz))
+                    slicev.append(slice(gtrajvals[it,2]-box2, gtrajvals[it,2]+box2+1))
+                    slicev.append(slice(gtrajvals[it,1]-box2, gtrajvals[it,1]+box2+1))
+                    slicevnoT.append(slice(0,dimz))
+                    slicevnoT.append(slice(gtrajvals[it,2]-box2, gtrajvals[it,2]+    \
+                      box2+1))
+                    slicevnoT.append(slice(gtrajvals[it,1]-box2, gtrajvals[it,1]+    \
+                      box2+1))
+                    slice2D.append(slice(0,dimz))
+                    slice2D.append(slice(gtrajvals[it,2]-box2, gtrajvals[it,2] +     \
+                      box2 + 1))
+                    slice2D.append(slice(gtrajvals[it,1]-box2, gtrajvals[it,1] +     \
+                      box2 + 1))
+
+                    varvalst = varobj[tuple(slicev)]
+# box values
+
+                    varvals[it,:,:,:] = varvalst
+#                    print 'values at time t______'
+#                    print varvalst
+
+# box stats values
+                    for iz in range(dimz):
+                        statvarvals[it,:,0] = varvalst[:,box2,box2]
+                        statvarvals[it,iz,1] = np.min(varvalst[iz,:,:])
+                        statvarvals[it,iz,2] = np.max(varvalst[iz,:,:])
+                        statvarvals[it,iz,3] = np.mean(varvalst[iz,:,:])
+                        statvarvals[it,iz,4] = np.mean(varvalst[iz,:,:]*             \
+                          varvalst[iz,:,:])
+                        statvarvals[it,iz,5] = np.sqrt(statvarvals[it,iz,4] -        \
+                          statvarvals[it,iz,3]*statvarvals[it,iz,3])
+                        statvarvals[it,iz,6] = np.sum(varvalst[iz,:,:])
+
+# Circle values
+                cslicev.append(gtrajvals[it,0])
+                if gtrajvals[it,2]-Nrad < 0 or gtrajvals[it,2]+Nrad + 1 >= dimy       \
+                  or gtrajvals[it,1]-Nrad < 0 or gtrajvals[it,1]+Nrad + 1 >= dimx:
+
+                    maxx = np.min([cxrangeslice[it][1], dimx])
+                    maxy = np.min([cyrangeslice[it][1], dimy])
+                    cslicev.append(slice(0,dimz))
+                    cslicev.append(slice(cyrangeslice[it][0], maxy))
+                    cslicev.append(slice(cxrangeslice[it][0], maxx)) 
+
+                    cslicevnoT.append(slice(0,dimz))
+                    cslicevnoT.append(slice(cyrangeslice[it][0], cyrangeslice[it][1]))
+                    cslicevnoT.append(slice(cxrangeslice[it][0], cxrangeslice[it][1]))
+
+                    cslice2D.append(slice(0,dimz))
+                    cslice2D.append(slice(0,maxy-cyrangeslice[it][0]))
+                    cslice2D.append(slice(0,maxx-cxrangeslice[it][0]))
+                    cslice2Dhor.append(slice(0, maxy - cyrangeslice[it][0]))
+                    cslice2Dhor.append(slice(0, maxx -cxrangeslice[it][0]))
+
+                    rvarvalst = np.ones((dimz,Nrad*2+1,Nrad*2+1),dtype=np.float)*    \
+                      fillValue
+                    rvarvalst[tuple(cslice2D)] = varobj[tuple(cslicev)]
+                    for iz in range(dimz):
+                        tslice = [slice(it,it)]+cslice2Dhor
+                        zslice = [slice(iz,iz)]+cslice2Dhor
+                        rvarvalst[tuple(zslice)] = np.where(circdist[tuple(tslice)] > \
+                          np.float(Nrad), fillValue, rvarvalst[tuple(zslice)]) 
+
+                    rvarvals[it,:,:,:] = rvarvalst
+
+# circle stats values
+                    maskedvals = ma.masked_values (rvarvalst, fillValue)
+                    maskedvals2 = maskedvals*maskedvals
+                    for iz in range(dimz):
+                        rstatvarvals[it,iz,0] = varvalst[iz,box2,box2]
+                        rstatvarvals[it,iz,1] = np.min(varvalst[iz,:,:])
+                        rstatvarvals[it,iz,2] = np.max(varvalst[iz,:,:])
+                        rstatvarvals[it,iz,3] = np.mean(varvalst[iz,:,:])
+                        rstatvarvals[it,iz,4] = maskedvals2[iz,:,:].mean()
+                        rstatvarvals[it,iz,5] = np.sqrt(rstatvarvals[it,iz,4] -      \
+                          rstatvarvals[it,iz,3]*rstatvarvals[it,iz,3])
+                        rstatvarvals[it,iz,6] = np.sum(varvalst[iz,:,:])
+
+                else:
+                    cslicev.append(slice(0,dimz))
+                    cslicev.append(slice(gtrajvals[it,2]-Nrad,gtrajvals[it,2]+Nrad+1))
+                    cslicev.append(slice(gtrajvals[it,1]-Nrad,gtrajvals[it,1]+Nrad+1))
+                    cslicevnoT.append(slice(0,dimz))
+                    cslicevnoT.append(slice(gtrajvals[it,2]-Nrad, gtrajvals[it,2]+   \
+                      Nrad+1))
+                    cslicevnoT.append(slice(gtrajvals[it,1]-Nrad, gtrajvals[it,1]+   \
+                      Nrad+1))
+                    cslicevnoThor.append(slice(gtrajvals[it,2]-Nrad,                 \
+                      gtrajvals[it,2] + Nrad+1))
+                    cslicevnoThor.append(slice(gtrajvals[it,1]-Nrad,                 \
+                      gtrajvals[it,1] + Nrad+1))
+                    cslice2D.append(slice(0,dimz))
+                    cslice2D.append(slice(gtrajvals[it,2]-Nrad,gtrajvals[it,2] +     \
+                      Nrad+1))
+                    cslice2D.append(slice(gtrajvals[it,1]-Nrad,gtrajvals[it,1] +     \
+                      Nrad+1))
+
+                    rvarvalst = varobj[tuple(cslicev)]
+# circle values
+                    for iz in range(dimz):
+                        tslice = [it]+cslicevnoThor
+                        rvarvalst[iz,:,:] = np.where(circdist[tuple(tslice)] >       \
+                          np.float(Nrad), fillValue, rvarvalst[iz,:,:]) 
+
+                    rvarvals[it,:,:,:] = rvarvalst
+
+# circle stats values
+                    maskedvals = ma.masked_values (rvarvalst, fillValue)
+                    maskedvals2 = maskedvals*maskedvals
+                    for iz in range(dimz):
+                        rstatvarvals[it,iz,0] = rvarvalst[iz,Nrad,Nrad]
+                        rstatvarvals[it,iz,1] = maskedvals[iz,:,:].min()
+                        rstatvarvals[it,iz,2] = maskedvals[iz,:,:].max()
+                        rstatvarvals[it,iz,3] = maskedvals[iz,:,:].mean()
+                        rstatvarvals[it,iz,4] = maskedvals2[iz,:,:].mean()
+                        rstatvarvals[it,iz,5] = np.sqrt(rstatvarvals[it,iz,4] -      \
+                          rstatvarvals[it,iz,3]*rstatvarvals[it,iz,3])
+                        rstatvarvals[it,iz,6] = maskedvals[iz,:,:].sum()
+
+#            print 'statistics:',rstatvarvals[it,:] 
+
+# variable box values
+            newvar = objofile.createVariable(vnst + 'box', 'f4', ('time','z','y','x'),\
+              fill_value=fillValue)
+            if searchInlist(varobj.ncattrs(),'standard_name'):
+                vsname = varobj.getncattr('standard_name')
+            else:
+                vsname = variables_values(vn)[1]
+            if searchInlist(varobj.ncattrs(),'long_name'):
+                vlname = varobj.getncattr('long_name')
+            else:
+                vlname = variables_values(vn)[4].replace('|',' ')
+            if searchInlist(varobj.ncattrs(),'units'):
+                vunits = varobj.getncattr('units')
+            else:
+                vunits = variables_values(vn)[5].replace('|',' ')
+
+            newattr = basicvardef(newvar, vsname, vlname, vunits)
+            newattr = newvar.setncattr('projection','lon lat')
+            newvar[:] = varvals
+
+# center of the trajectory
+            newvar = objofile.createVariable('trj_' + vnst, 'f', ('time','z'),     \
+              fill_value=fillValue)
+            newattr = basicvardef(newvar, 'trj_' + vsname, 'value along the ' +    \
+              'trajectory of '+ vn, vunits)
+            newvar[:] = statvarvals[:,:,0]
+
+# variable box statistics
+            ist = 0
+            for statn in statnames:
+                newvar = objofile.createVariable(statn + '_' + vnst,'f',('time','z'),\
+                  fill_value=fillValue)
+                newattr = basicvardef(newvar, statn + '_' + vsname, vstlname[ist] +  \
+                  ' the box ('+str(boxs)+'x'+str(boxs)+') of ' + vnst, vunits)
+                newvar[:] = statvarvals[:,:,ist+1]
+#                newattr = set_attributek(newvar,'_FillValue',fillValue,'npfloat')
+                ist = ist + 1
+        
+# variable circle values
+            newvar = objofile.createVariable(vnst + 'circle','f4',('time','z','yr',    \
+              'xr'), fill_value=fillValue)
+            if searchInlist(varobj.ncattrs(),'standard_name'):
+                vsname = varobj.getncattr('standard_name')
+            else:
+                vsname = variables_values(vn)[1]
+            if searchInlist(varobj.ncattrs(),'long_name'):
+                vlname = varobj.getncattr('long_name')
+            else:
+                vlname = variables_values(vn)[4].replace('|',' ')
+            if searchInlist(varobj.ncattrs(),'units'):
+                vunits = varobj.getncattr('units')
+            else:
+                vunits = variables_values(vn)[5].replace('|',' ')
+
+            newattr = basicvardef(newvar, vsname, vlname, vunits)
+            newattr = newvar.setncattr('projection','lon lat')
+            newvar[:] = rvarvals
+
+# variable circle statistics
+            ist = 0
+            for statn in cstatnames:
+                newvar = objofile.createVariable(statn + '_' + vnst,'f',('time','z'),\
+                  fill_value=fillValue)
+                newattr = basicvardef(newvar, statn + '_' + vsname, vstlname[ist] +  \
+                  ' the circle of radius ('+ str(circler)+') of ' + vnst, vunits)
+                newvar[:] = rstatvarvals[:,:,ist+1]
+#                newattr = set_attributek(newvar,'_FillValue',fillValue,'npfloat')
+                ist = ist + 1
+
+        elif Nvardims == 3:
+# Too optimistic, but at this stage...
+            dimt = varobj.shape[0]
+ 
+            varvals = np.ones(tuple([Ttraj,boxs,boxs]), dtype=np.float)
+            lonvals = np.ones(tuple([Ttraj,boxs,boxs]), dtype=np.float)
+            latvals = np.ones(tuple([Ttraj,boxs,boxs]), dtype=np.float)
+            rvarvals = np.ones(tuple([Ttraj,Nrad*2+1,Nrad*2+1]), dtype=np.float)
+            rlonvals = np.ones(tuple([Ttraj,Nrad*2+1,Nrad*2+1]), dtype=np.float)
+            rlatvals = np.ones(tuple([Ttraj,Nrad*2+1,Nrad*2+1]), dtype=np.float)
+
+# One dimension plus for the values at the center of the trajectory
+            statvarvals = np.ones(tuple([Ttraj,Nsts+1]), dtype=np.float)
+            rstatvarvals = np.ones(tuple([Ttraj,Nsts+1]), dtype=np.float)
+
+            for it in range(Ttraj):
+                it0 = Tbeg + it
+
+                slicev = []
+                slice2D = []
+                slicevnoT = []
+                cslicev = []
+                cslice2D = []
+                cslicevnoT = []
+                cslicevC = []
+
+                slicev.append(gtrajvals[it,0])
+                if gtrajvals[it,2]-box2 < 0 or gtrajvals[it,2]+box2 + 1 > dimy + 1   \
+                  or gtrajvals[it,1]-box2 < 0 or gtrajvals[it,1]+box2 + 1 > dimx + 1:
+# box values
+                    slicev.append(slice(yrangeslice[it][0],yrangeslice[it][1]))
+                    slicev.append(slice(xrangeslice[it][0],xrangeslice[it][1]))
+
+                    slicevnoT.append(slice(yrangeslice[it][0],yrangeslice[it][1]))
+                    slicevnoT.append(slice(xrangeslice[it][0],xrangeslice[it][1]))
+
+                    slice2D.append(slice(0,yrangeslice[it][1]-yrangeslice[it][0]))
+                    slice2D.append(slice(0,xrangeslice[it][1]-xrangeslice[it][0]))
+
+                    rvarvalst = np.ones((Nrad*2+1, Nrad*2+1),dtype=np.float)*fillValue
+
+                    varvalst[tuple(slice2D)] = varobj[tuple(slicev)]
+                    varvals[it,:,:] = varvalst
+    
+# box stats values
+                    maskedvals = ma.masked_values (varvalst, fillValue)
+                    statvarvals[it,0] = varvalst[box2,box2]
+                    statvarvals[it,1] = maskedvals.min()
+                    statvarvals[it,2] = maskedvals.max()
+                    statvarvals[it,3] = maskedvals.mean()
+                    maskedvals2 = maskedvals*maskedvals
+                    statvarvals[it,4] = maskedvals2.mean()
+                    statvarvals[it,5] = np.sqrt(statvarvals[it,4] -                   \
+                      statvarvals[it,3]*statvarvals[it,3])
+                    statvarvals[it,6] = maskedvals.sum()
+
+                    varvalst[tuple(slice2D)] = lonv[tuple(slicevnoT)]
+                    lonvals[it,:,:] = varvalst
+                    varvalst[tuple(slice2D)] = latv[tuple(slicevnoT)]           
+                    latvals[it,:,:] = varvalst
+
+                else:
+                    slicev.append(slice(gtrajvals[it,2]-box2, gtrajvals[it,2]+box2+1))
+                    slicev.append(slice(gtrajvals[it,1]-box2, gtrajvals[it,1]+box2+1))
+                    slicevnoT.append(slice(gtrajvals[it,2]-box2, gtrajvals[it,2]+    \
+                      box2+1))
+                    slicevnoT.append(slice(gtrajvals[it,1]-box2, gtrajvals[it,1]+    \
+                      box2+1))
+
+                    slice2D.append(slice(gtrajvals[it,2]-box2, gtrajvals[it,2]+box2 +\
+                      1))
+                    slice2D.append(slice(gtrajvals[it,1]-box2, gtrajvals[it,1]+box2 +\
+                      1))
+
+                    varvalst = varobj[tuple(slicev)]
+# box values
+
+                    varvals[it,:,:] = varvalst
+#                    print 'values at time t______'
+#                    print varvalst
+
+# box stats values
+                    statvarvals[it,0] = varvalst[box2,box2]
+                    statvarvals[it,1] = np.min(varvalst)
+                    statvarvals[it,2] = np.max(varvalst)
+                    statvarvals[it,3] = np.mean(varvalst)
+                    statvarvals[it,4] = np.mean(varvalst*varvalst)
+                    statvarvals[it,5] = np.sqrt(statvarvals[it,4] -                  \
+                      statvarvals[it,3]*statvarvals[it,3])
+                    statvarvals[it,6] = np.sum(varvalst)
+
+                    varvalst = lonv[tuple(slicevnoT)]
+                    lonvals[it,:,:] = varvalst
+                    varvalst = latv[tuple(slicevnoT)]           
+                    latvals[it,:,:] = varvalst
+
+# Circle values
+                cslicev.append(gtrajvals[it,0])
+                if gtrajvals[it,2]-Nrad < 0 or gtrajvals[it,2]+Nrad + 1 >= dimy      \
+                  or gtrajvals[it,1]-Nrad < 0 or gtrajvals[it,1]+Nrad + 1 >= dimx:
+                    maxx = np.min([cxrangeslice[it][1], dimx])
+                    maxy = np.min([cyrangeslice[it][1], dimy])
+                    cslicev.append(slice(cyrangeslice[it][0],maxy))
+                    cslicev.append(slice(cxrangeslice[it][0],maxx))
+
+                    cslicevnoT.append(slice(cyrangeslice[it][0],cyrangeslice[it][1]))
+                    cslicevnoT.append(slice(cxrangeslice[it][0],cxrangeslice[it][1]))
+
+                    cslice2D.append(slice(0,maxy - cyrangeslice[it][0]))
+                    cslice2D.append(slice(0,maxx - cxrangeslice[it][0]))
+
+                    rvarvalst = np.ones((Nrad*2+1,Nrad*2+1),dtype=np.float)*fillValue
+                    rvarvalst[tuple(cslice2D)] = varobj[tuple(cslicev)]
+                    rvarvalst[tuple(cslice2D)] = np.where(circdist[tuple(cslicev)] >\
+                      np.float(Nrad), fillValue, rvarvalst[tuple(cslice2D)])
+
+                    rvarvals[it,:,:] = rvarvalst
+
+# circle stats values
+                    maskedvals = ma.masked_values (rvarvalst, fillValue)
+                    rstatvarvals[it,0] = rvarvalst[Nrad,Nrad]
+                    rstatvarvals[it,1] = maskedvals.min()
+                    rstatvarvals[it,2] = maskedvals.max()
+                    rstatvarvals[it,3] = maskedvals.mean()
+                    maskedvals2 = maskedvals*maskedvals
+                    rstatvarvals[it,4] = maskedvals2.mean()
+                    rstatvarvals[it,5] = np.sqrt(rstatvarvals[it,4] -                \
+                      rstatvarvals[it,3]*rstatvarvals[it,3])
+                    rstatvarvals[it,6] = maskedvals2.sum()
+
+                else:
+                    cslicev.append(slice(gtrajvals[it,2]-Nrad,gtrajvals[it,2]+Nrad+1))
+                    cslicev.append(slice(gtrajvals[it,1]-Nrad,gtrajvals[it,1]+Nrad+1))
+                    cslicevnoT.append(slice(gtrajvals[it,2]-Nrad, gtrajvals[it,2]+    \
+                      Nrad+1))
+                    cslicevnoT.append(slice(gtrajvals[it,1]-Nrad, gtrajvals[it,1]+    \
+                      Nrad+1))
+
+                    cslice2D.append(slice(gtrajvals[it,2]-Nrad,gtrajvals[it,2]+Nrad+1))
+                    cslice2D.append(slice(gtrajvals[it,1]-Nrad,gtrajvals[it,1]+Nrad+1))
+
+                    cslicevC.append(it)
+                    cslicevC.append(slice(gtrajvals[it,2]-Nrad,gtrajvals[it,2]+Nrad+1))
+                    cslicevC.append(slice(gtrajvals[it,1]-Nrad,gtrajvals[it,1]+Nrad+1))
+
+                    rvarvalst = varobj[tuple(cslicev)]
+                    cdist = circdist[tuple(cslicevC)]
+# circle values
+                    rvarvalst = np.where(cdist > np.float(Nrad),fillValue,rvarvalst)
+                    rvarvals[it,:,:] = rvarvalst
+
+# circle stats values
+                    maskedvals = ma.masked_values (rvarvalst, fillValue)
+                    rstatvarvals[it,0] = rvarvalst[Nrad,Nrad]
+                    rstatvarvals[it,1] = maskedvals.min()
+                    rstatvarvals[it,2] = maskedvals.max()
+                    rstatvarvals[it,3] = maskedvals.mean()
+                    maskedvals2 = maskedvals*maskedvals
+                    rstatvarvals[it,4] = maskedvals2.mean()
+                    rstatvarvals[it,5] = np.sqrt(rstatvarvals[it,4] -                \
+                      rstatvarvals[it,3]*rstatvarvals[it,3])
+                    rstatvarvals[it,6] = maskedvals.sum()
+
+#            print 'statistics:',rstatvarvals[it,:] 
+
+# variable box values
+            newvar = objofile.createVariable(vnst + 'box', 'f4', ('time', 'y', 'x'), \
+              fill_value=fillValue)
+            if searchInlist(varobj.ncattrs(),'standard_name'):
+                vsname = varobj.getncattr('standard_name')
+            else:
+                vsname = variables_values(vn)[1]
+            if searchInlist(varobj.ncattrs(),'long_name'):
+                vlname = varobj.getncattr('long_name')
+            else:
+                vlname = variables_values(vn)[4].replace('|',' ')
+            if searchInlist(varobj.ncattrs(),'units'):
+                vunits = varobj.getncattr('units')
+            else:
+                vunits = variables_values(vn)[5].replace('|',' ')
+
+            newattr = basicvardef(newvar, vsname, vlname, vunits)
+            newattr = newvar.setncattr('projection','lon lat')
+            newvar[:] = varvals
+
+# center of the trajectory
+            newvar = objofile.createVariable('trj_' + vnst, 'f', ('time'),           \
+              fill_value=fillValue)
+            newattr = basicvardef(newvar, 'trj_' + vsname, 'value along the ' +      \
+              'trajectory of '+ vnst, vunits)
+            newvar[:] = statvarvals[:,0]
+
+# variable box statistics
+            ist = 0
+            for statn in statnames:
+                newvar = objofile.createVariable(statn + '_' + vnst, 'f', ('time'),  \
+                  fill_value=fillValue)
+                newattr = basicvardef(newvar, statn + '_' + vsname, vstlname[ist] +  \
+                  ' the box ('+str(boxs)+'x'+str(boxs)+') of ' + vnst, vunits)
+                newvar[:] = statvarvals[:,ist+1]
+#                newattr = set_attributek(newvar,'_FillValue',fillValue,'npfloat')
+                ist = ist + 1
+        
+# variable circle values
+            newvar = objofile.createVariable(vnst + 'circle','f4',('time','yr','xr'),\
+              fill_value=fillValue)
+            if searchInlist(varobj.ncattrs(),'standard_name'):
+                vsname = varobj.getncattr('standard_name')
+            else:
+                vsname = variables_values(vn)[1]
+            if searchInlist(varobj.ncattrs(),'long_name'):
+                vlname = varobj.getncattr('long_name')
+            else:
+                vlname = variables_values(vn)[4].replace('|',' ')
+            if searchInlist(varobj.ncattrs(),'units'):
+                vunits = varobj.getncattr('units')
+            else:
+                vunits = variables_values(vn)[5].replace('|',' ')
+
+            newattr = basicvardef(newvar, vsname, vlname, vunits)
+            newattr = newvar.setncattr('projection','lon lat')
+            newvar[:] = rvarvals
+
+# variable circle statistics
+            ist = 0
+            for statn in cstatnames:
+                newvar = objofile.createVariable(statn + '_' + vnst, 'f', ('time'),  \
+                  fill_value=fillValue)
+                newattr = basicvardef(newvar, statn + '_' + vsname, vstlname[ist] +  \
+                  ' the circle of radius ('+ str(circler)+') of ' + vnst, vunits)
+                newvar[:] = rstatvarvals[:,ist+1]
+#                newattr = set_attributek(newvar,'_FillValue',fillValue,'npfloat')
+                ist = ist + 1
+        else:
+# Other variables
+            print warnmsg
+            print '  ' + fname + ": variable '" + vn + "' shape:",varobj.shape,' not'\
+              + ' ready!!'
+            print '    skipping variable'
+            if len(varns) == 1:
+                objofile.close()
+                sub.call(['rm', ofile])
+                print '    uniq variable! removing file and finishing program'
+                quit()
+
+        if not objofile.variables.has_key('trlon') and Nvardims == 3:
+# var dimensions
+            newvar = objofile.createVariable('trlon', 'f8', ('time'))
+            newattr = basicvardef(newvar,'trlon','trajectory longitude',             \
+              'degrees west_east')
+            newvar[:] = trajvals[:,1]
+
+            newvar = objofile.createVariable('trlat', 'f8', ('time'))
+            newattr = basicvardef(newvar,'trlat','trajectory latitude',              \
+              'degrees north_south')
+            newvar[:] = trajvals[:,2]
+
+            newvar = objofile.createVariable('time', 'f8', ('time'))
+            newattr = basicvardef(newvar, 'time', 'time', tunits)
+            newvar[:] = trajvals[:,0]
+
+            newvar = objofile.createVariable('lon', 'f8', ('time', 'y', 'x'),        \
+              fill_value=fillValue)
+            newattr = basicvardef(newvar, 'longitude', 'longitude',                  \
+              'degrees west_east')
+            newvar[:] = lonvals
+
+            newvar = objofile.createVariable('lat', 'f8', ('time', 'y', 'x'),        \
+              fill_value=fillValue)
+            newattr = basicvardef(newvar, 'latitude', 'latitude',                    \
+              'degrees north_south')
+            newvar[:] = latvals
+           
+        ivar = ivar + 1
+
+
+# global attributes
+    objofile.setncattr('author', 'L. Fita')
+    objofile.setncattr('institution', 'Laboratire de Meteorologie Dynamique')
+    objofile.setncattr('university', 'Pierre Marie Curie - Jussieu')
+    objofile.setncattr('center', 'Centre National de Recherches Scientifiques')
+    objofile.setncattr('city', 'Paris')
+    objofile.setncattr('country', 'France')
+    objofile.setncattr('script', 'nc_var_tools.py')
+    objofile.setncattr('function', 'compute_tevolboxtraj')
+    objofile.setncattr('version', '1.0')
+    objofile.setncattr('data_file',ncfile)
+
+    objfile.close()
+
+    objofile.sync()
+    objofile.close()
+
+    print '  ' + fname + ': successful creation of file "' + ofile + '" !!!'
+
+    return 
+
 
 #quit()