Changeset 628 in lmdz_wrf for trunk

Timestamp:

Sep 14, 2015, 6:56:14 PM (9 years ago)

Author:

lfita

Message:

Working on the radi_traj

File:

• trunk/tools/nc_var_tools.py (modified) (10 diffs)

trunk/tools/nc_var_tools.py

@@ -13 +13 @@
 fillValueC = '-'
 fillValueI = -99999
+
+def values_line(line, splitv, chars):
+    """ Function to retrieve values from a line of ASCII text removing a series of characters
+      line= line of ASCII text to work with
+      splitv= character to use to split with the line
+      chars= list of characters/strings to remove
+    >>> values_line('0.0 0.0 1 [ 0 0 ]', ' ', ['[', ']'])
+    ['0.0', '0.0', '1', '0', '0']
+    """
+    fname = 'values_line'
+
+    vals = line.replace('\n', '')
+    for ic in chars:
+        vals = vals.replace(ic, '')
+
+    vals0 = vals.split(splitv)
+    vals = []
+    for val in vals0:
+        if not len(val) == 0: vals.append(val)
+
+    return vals
 
 def filexist(filen, emsg, fmsg):
@@ -15155 +15176 @@
     return radpos
 
+def read_SQradii(filename='SQradii.dat', Npt=-1, maxradii=None):
+    """ Function to read the outcome of the function `squared_radial'
+      filename= Name of the outcome file of the function 
+      Npt= Number of points on each grid
+      maxradii= allowed maximum radii
+      'SQradii.dat': file with the radii and paris of integers after ran squared_radial(1000)
+    """
+    fname = 'read_SQradii'
+
+    if not os.path.isfile(filename):
+        print errormsg
+        print '  ' + fname + ": file '" + filename + "' not found!!"
+        quit(-1)
+
+    infile = open(filename, 'r')
+
+    gridradii = {}
+    Nradii = {}
+    radii = []
+    for line in infile:
+        values = values_line(line, ' ', ['[', ']'])
+        rad = np.float(values[0])
+# Limited by a given radii
+        if maxradii is not None and rad > maxradii: break
+        Nrad = int(values[2])
+        couples = None
+        Nradgood = 0
+        if Nrad > 1:
+           Ncouples = 0
+           couples = []
+           for ic in range(Nrad):
+               ix = int(values[3 + ic*2])
+               iy = int(values[3 + ic*2 + 1])
+               if ix <= Npt and iy <= Npt:
+                   couples.append(np.array([ix,iy], dtype=int))
+                   Nradgood = Nradgood + 1
+               elif ix > Npt and iy > Npt: break
+        else:
+             ix = int(values[3])
+             iy = int(values[4])
+             if ix <= Npt and iy <= Npt:
+                 couples = np.array([ix,iy], dtype=int)
+                 Nradgood = Nradgood + 1
+             elif ix > Npt and iy > Npt: break
+
+        if Nradgood > 0:
+            Nradii[rad] = Nradgood
+            gridradii[rad] = np.array(couples)
+
+#    for rad in sorted(gridradii.keys()):
+#        print rad, gridradii[rad]
+
+    return gridradii
+
 def grid_combinations(x,y):
     """ Function to provide all the possible grid points combination for a given pair of values
@@ -15188 +15263 @@
     fname = 'radii_points'
 
-    radiipos = np.zeros((Nangle, radii, 2), dtype = int)
+    radiipos = np.zeros((Nang, Lrad, 2), dtype = int)
     SQradiipos = {}
 
@@ -15199 +15274 @@
     for ia in range(Nang):
         angle = 2.*np.pi*ia/Nang
+        posangle = np.zeros((Lrad,2), dtype=np.float)
 
 # Points at that given angle
-        pangle = radial_points(angle,radii)
+        pangle = radial_points(angle,Lrad)
         posangle[:,0] = pangle[:,0] + xpos
         posangle[:,1] = pangle[:,1] + ypos
 
-        for ip in range(radii):
+        for ip in range(Lrad):
             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 < dx and jjpos >= 0 and jjpos < dy:
+            if iipos >= 0 and iipos < dx-1 and jjpos >= 0 and jjpos < dy-1:
                 radiipos[ia,ip,0] = iipos
                 radiipos[ia,ip,1] = jjpos
@@ -15222 +15298 @@
 
 # pairs loop
+        Npairs = len(pairs.flatten())/2
         pairsw = []
-        for ip in range(len(pairs)):
-            allpairs = grid_combinations(pairs[ip*2],pairs[ip*2+1])
-
+        print 'Lluis ',SQradii[ir],'pairs:',pairs, type([pairs])
+        for ip in range(Npairs):
+            if Npairs == 0: 
+                print '  Lluis NO pair:', SQradii[ir],':',pairs
+                break
+            elif Npairs == 1: 
+                print '  Lluis single pair:',pairs[ip*2],pairs[ip*2+1]
+                allpairs = grid_combinations(pairs[ip*2],pairs[ip*2+1])
+            else:
+                print ip,'  Lluis 1 pair:',pairs[ip,0],pairs[ip,1]
+                allpairs = grid_combinations(pairs[ip,0],pairs[ip,1])
+                
             for iap in range(len(allpairs)):
                 ppos = allpairs[iap]
-                iipos = trajvals[it,1] + ppos[0]
-                jjpos = trajvals[it,1] + ppos[0]
+                iipos = xpos + ppos[0]
+                jjpos = ypos + ppos[1]
 
 # Creation of a matrix with all the grid points at each time-step
@@ -15317 +15403 @@
         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)
+# Removing this checking
+#    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
@@ -15342 +15429 @@
     it = 0
 
-    trajpos = np.zeros(Ttraj,Nangle,raddi)
+    trajpos = np.zeros((Ttraj,Nangle,radii), dtype=np.float)
+    gtrajvals = np.zeros((Ttraj,3), dtype=int)
+    trajvals = np.zeros((Ttraj,3), dtype=np.float)
 
 # Grid position at each time-step and angle
-    varvalst = np.ones((Traj, Nangle, radii, 2), dtype=np.float)*fillValue
-    varvalstSQ = np.zeros((Traj, radii, radii), dtype=bool)
+#    varvalst = np.ones((Ttraj, Nangle, radii, 2), dtype=np.float)*fillValue
+#    varvalstSQ = np.zeros((Ttraj, radii, radii), dtype=bool)
 
     it = 0
@@ -16036 +16125 @@
 
 compute_tevolboxtraj_radialsec('/bdd/PCER/workspace/lfita/etudes/FF/medic051215_2km/run/trajectory.dat@6,XLONG,XLAT,ZNU,Times,wrf,4,10',       \
-  '001/full_concatenated.nc', 'U10')
+  '/bdd/PCER/workspace/lfita/etudes/FF/medic051215_2km/run/wrfout/wrfout_d02_2005-12-12_00:00:00', 'U10')
 quit()
 
@@ -16443 +16532 @@
     return fvalues
 
-def values_line(line, splitv, chars):
-    """ Function to retrieve values from a line of ASCII text removing a series of characters
-      line= line of ASCII text to work with
-      splitv= character to use to split with the line
-      chars= list of characters/strings to remove
-    >>> values_line('0.0 0.0 1 [ 0 0 ]', ' ', ['[', ']'])
-    ['0.0', '0.0', '1', '0', '0']
-    """
-    fname = 'values_line'
-
-    vals = line.replace('\n', '')
-    for ic in chars:
-        vals = vals.replace(ic, '')
-
-    vals0 = vals.split(splitv)
-    vals = []
-    for val in vals0:
-        if not len(val) == 0: vals.append(val)
-
-    return vals
-
 def reduce_last_spaces(string):
     """ Function to reduce the last right spaces from a string
@@ -16490 +16558 @@
 
     return newstr
-
-def read_SQradii(filename='SQradii.dat', Npt=1000, maxradii=None):
-    """ Function to read the outcome of the function `squared_radial'
-      filename= Name of the outcome file of the function 
-      Npt= Number of points on each grid
-      maxradii= allowed maximum radii
-      'SQradii.dat': file with the radii and paris of integers after ran squared_radial(1000)
-    """
-    fname = 'read_SQradii'
-
-    if not os.path.isfile(filename):
-        print errormsg
-        print '  ' + fname + ": file '" + filename + "' not found!!"
-        quit(-1)
-
-    infile = open(filename, 'r')
-
-    gridradii = {}
-    Nradii = {}
-    radii = []
-    for line in infile:
-        print 'line:',line
-        values = values_line(line, ' ', ['[', ']'])
-        rad = np.float(values[0])
-# Limited by a given radii
-        if maxradii is not None and rad > maxradii: break
-        Nrad = int(values[2])
-        couples = None
-        Nradgood = 0
-        if Nrad > 1:
-           Ncouples = 0
-           couples = []
-           for ic in range(Nrad):
-               ix = int(values[3 + ic*2])
-               iy = int(values[3 + ic*2 + 1])
-               if ix <= Npt and iy <= Npt:
-                   couples.append(np.array([ix,iy], dtype=int))
-                   Nradgood = Nradgood + 1
-               elif ix > Npt and iy > Npt: break
-        else:
-             ix = int(values[3])
-             iy = int(values[4])
-             if ix <= Npt and iy <= Npt:
-                 couples = np.array([ix,iy], dtype=int)
-                 Nradgood = Nradgood + 1
-             elif ix > Npt and iy > Npt: break
-
-        if Nradgood > 0:
-            Nradii[rad] = Nradgood
-            gridradii[rad] = np.array(couples)
-
-    for rad in sorted(gridradii.keys()):
-        print rad, gridradii[rad]
-
-    return gridradii
 
 def squared_radial(Npt):