Changeset 624 in lmdz_wrf

Timestamp:

Sep 8, 2015, 12:26:44 PM (10 years ago)

Author:

lfita

Message:

Adding calculation of integer squared radii

Location:

trunk/tools

Files:

• SQradii.dat (added)
• nc_var_tools.py (modified) (3 diffs)

trunk/tools/nc_var_tools.py

@@ -874 +874 @@
     return valpos
 
+def addfileInfile(origfile,destfile,addfile,addsign):
+    """ Function to add the content of a file [addfile] to another one [origfile] at 
+      a given position [addsign] creating a new file called [destfile]
+      origfile= original file
+      destfile= destination file
+      addfile= content of the file to add
+      addsign= sign where to add the file
+    """
+    fname = 'addfileInfile'
+
+    if not os.path.isfile(origfile):
+        print errormsg
+        print '  ' + fname + ": original file '" + origfile + "' does not exist !!"
+        quit()
+
+    if not os.path.isfile(addfile):
+        print errormsg
+        print '  ' + fname + ": adding file '" + addfile + "' does not exist !!"
+        quit()
+    
+    ofile = open(origfile, 'r')
+
+# Inspecting flag
+##
+    Nfound = 0
+    for line in ofile:
+        if line == addsign + '\n': Nfound = Nfound + 1
+
+    if Nfound == 0:
+        print errormsg
+        print '  ' + fname + ": adding sign '" + addsign + "' not found !!"
+        quit()
+    print '  '+ fname + ': found', Nfound," times sign '" + addsign + "' "
+
+    ofile.seek(0)
+    dfile = open(destfile, 'w')
+
+    for line in ofile:
+        if line != addsign + '\n':
+            dfile.write(line)
+        else:
+            afile = open(addfile,'r')
+            for aline in afile:
+                print aline
+                dfile.write(aline)
+
+            afile.close()
+
+    ofile.close()
+    dfile.close()
+
+    print main + " successful writting of file '" + destfile + "' !!"
+    return
+
+
 ####### ###### ##### #### ### ## #
 
@@ -16388 +16443 @@
     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
@@ -16414 +16490 @@
 
     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):
+    """ Function to provide the series of radii as composite of pairs (x,y) of gid cells
+      Npt= largest amount of grid points on x and y directions
+    >>> squared_radial(5)
+    [[0.0, 0, 0], [1.0, 1, 0], [1.4142135623730951, 1, 1], [2.0, 2, 0], [2.2360679774997898, 2, 1], [2.8284271247461903, 2, 2], [3.0, 3, 0], [3.1622776601683795, 3, 1], [3.6055512754639891, 3, 2], [4.0, 4, 0], [4.1231056256176606, 4, 1], [4.2426406871192848, 3, 3], [4.4721359549995796, 4, 2], [5.0, array([4, 3]), array([5, 0])], [5.0990195135927845, 5, 1], [5.3851648071345037, 5, 2], [5.6568542494923806, 4, 4], [5.8309518948453007, 5, 3], [6.4031242374328485, 5, 4], [7.0710678118654755, 5, 5]]
+    """
+    fname = 'squared_radial'
+
+    radii = []
+    gridradii = {}
+    singradii = []
+    Nradii = {}
+
+# Computing all radii
+##
+    for i in range(Npt+1):
+        if np.mod(i,100) == 0: print 'done: ',i
+        for j in range(i+1):
+#            ir = np.float(i)
+#            jr = np.float(j)
+            ir = np.float64(i)
+            jr = np.float64(j)
+            rad = np.sqrt(ir*ir + jr*jr)
+            if not searchInlist(singradii, rad):
+                singradii.append(rad)
+                gridradii[rad] = np.array([i,j], dtype=int)
+                Nradii[rad] = 1
+            else:
+#                print warnmsg
+#                print '  ' + fname + ': repeated radii',rad,'!!'
+#                print '    Previous values:', gridradii[rad]
+
+                oldvalues = gridradii[rad]
+                Nradii[rad] = Nradii[rad] + 1
+                newvalues = np.zeros((Nradii[rad],2), dtype=int)
+                if Nradii[rad] == 2:
+                    newvalues[0,:] = oldvalues
+                    Nstart = 1
+                else:
+                    Nstart = 0
+                for Nr in range(Nstart,Nradii[rad]-1):
+                    newvalues[Nr,:] = oldvalues[Nr,:]
+                newvalues[Nradii[rad]-1,:] = np.array([i,j], dtype=int)
+                gridradii[rad] = newvalues
+
+# Sorting values
+##
+## Only required to downwrite the output file
+#    radiif = open('SQradii.dat', 'w')
+#    radii = []
+#    prevrad = 0.
+#    for rad in sorted(gridradii.keys()):
+#        dradii = rad - prevrad
+## Check fo repeated values
+#        radii.append([rad, gridradii[rad][0], gridradii[rad][1]])
+#        radiif.write(str(rad) + ' ' + str(dradii) + ' ' + str(Nradii[rad]) + ' [ '+  \
+#          numVector_String(gridradii[rad].flatten(), ' ') + ' ] \n')
+#        prevrad = rad.copy()
+
+#    radiif.close()
+    return gridradii
 
 def getvalues_lonlat(values, ncfile):