Changeset 828 in lmdz_wrf for trunk

Timestamp:

Jun 16, 2016, 12:13:05 AM (8 years ago)

Author:

lfita

Message:

Adding:

• `radius_angle': Function to generate a matrix with the angle at a given point
• `subbasin_point': Function to provide sub-basins given a grid point following a matrix of trips

Location:

trunk/tools

Files:

• generic.py (modified) (1 diff)
• generic_tools.py (modified) (5 diffs)

trunk/tools/generic.py

@@ -126 +126 @@
         print gen.datetimeStr_conversion(vals[0], vals[1], vals[2])
 
-'days_period'
+#'days_period'
 
 elif oper == 'grid_combinations':

trunk/tools/generic_tools.py

@@ -30 +30 @@
 # ijlonlat: Function to provide the imin,jmin imax,jmax of a lon,lat box
 # incomming_flow: Function to determine if a fgrid-flow inflows to the central grid point
+# index_mat_way: Function to look for a value within a matrix following a sense
 # num_chainSnum: Function to pass a value to a `ChainStrNum' number
 # num_split: Function to split a string at each numeric value keeping the number as the last character of each cut
 # PolyArea: Function to compute the area of the polygon following 'Shoelace formula'
+# radius_angle: Function to generate a matrix with the angle at a given point
+# radius_dist: Function to generate a matrix with the distance at a given point
 # significant_decomposition: Function to decompose a given number by its signifcant potencies
 # unitsdsDate: Function to know how many units of time are from a given pair of dates
@@ -4634 +4637 @@
       [dx/y]: dimension of the matrix
       [ptx/y]: grid point coordinates of the point
-    >>> radius_dist(3,5,2,2)
-    [[ 1.41421356  1.          1.41421356  2.23606798  3.16227766]
-     [ 1.          0.          1.          2.          3.        ]
-     [ 1.41421356  1.          1.41421356  2.23606798  3.16227766]]
+    >>> radius_dist(5,4,2,2)
+    [[ 2.82842712  2.23606798  2.          2.23606798  2.82842712]
+     [ 2.23606798  1.41421356  1.          1.41421356  2.23606798]
+     [ 2.          1.          0.          1.          2.        ]
+     [ 2.23606798  1.41421356  1.          1.41421356  2.23606798]]
     """
 
@@ -4648 +4652 @@
         quit(-1)
 
-    xdist =  np.zeros((dx,dy), dtype=np.float)
-    ydist =  np.zeros((dx,dy), dtype=np.float)
-    dist =  np.zeros((dx,dy), dtype=np.float)
+    xdist =  np.zeros((dy,dx), dtype=np.float)
+    ydist =  np.zeros((dy,dx), dtype=np.float)
+    dist =  np.zeros((dy,dx), dtype=np.float)
 
     for ix in range(dx):
-        xdist[ix,:] = np.float(ix-ptx)
+        xdist[:,ix] = np.float(ix-ptx)
     for iy in range(dy):
-        ydist[:,iy] = np.float(iy-pty)
+        ydist[iy,:] = np.float(pty-iy)
 
     dist = np.sqrt(xdist*xdist + ydist*ydist)
 
     return dist
+
+def radius_angle(dx,dy,ptx,pty):
+    """ Function to generate a matrix with the angle at a given point
+    radius_angle(dx,dy,ptx,pty)
+      [dx/y]: dimension of the matrix
+      [ptx/y]: grid point coordinates of the point
+    >>> radius_angle(5,4,2,2)*180./np.pi
+    [[ 315.          333.43494882    0.           26.56505118   45.        ]
+     [ 296.56505118  315.            0.           45.           63.43494882]
+     [ 270.          270.            0.           90.           90.        ]
+     [ 243.43494882  225.          180.          135.          116.56505118]]
+    """
+    fname = 'radius_angle'
+
+    if ptx < 0 or ptx > dx-1 or pty < 0 or pty > dy-1:
+        print errormsg
+        print '  ' + fname + ': wrong point coordinates:',dx,',',dy,'for matrix;',dx \
+         ,'x',dy
+        quit(-1)
+
+    xdist =  np.zeros((dy,dx), dtype=np.float)
+    ydist =  np.zeros((dy,dx), dtype=np.float)
+    angle =  np.zeros((dy,dx), dtype=np.float)
+
+    for ix in range(dx):
+        xdist[:,ix] = np.float(ix-ptx)
+    for iy in range(dy):
+        ydist[iy,:] = np.float(pty-iy)
+
+    angle = np.arctan2(xdist,ydist)
+    angle = np.where(angle < 0., angle + 2*np.pi, angle)
+
+    return angle
 
 def lonlat2D(lon,lat):
@@ -7320 +7357 @@
     hexblue = '{0:X}'.format(int(col[2]*255))
 
-    if hexred == '0': hexred = '00'
+    if hexred == '0': hexred = '00' 
     if hexgreen == '0': hexgreen = '00'
     if hexblue == '0': hexblue = '00'
@@ -7328 +7365 @@
     return hexcol
 
+def index_mat_way(mat,val,way):
+    """ Function to look for a value within a matrix following a sense
+      mat: matrix
+      val: value to search
+      way: way of search
+        'anticlockwise12': looking in an anti-clockwise sens starting at 12
+        'clockwise12': looking in an clockwise sens starting at 12
+    >>> mat = np.arange(20).reshape(4,5)
+    >>> mat[2,1] = 3
+    >>> index_mat_way(mat,3,'clockwise12')
+    [array([0, 3]), array([2, 1])]
+    >>> index_mat_way(mat,3,'anticlockwise12') 
+    [array([2, 1]), array([0, 3])]
+    """
+    fname = 'index_mat_way'
+
+    ways = ['anticlockwise12', 'anticlockwise6', 'clockwise12', 'clockwise6']
+
+    Ndims = len(mat.shape)
+    dx = mat.shape[Ndims-1]
+    dy = mat.shape[Ndims-2]
+
+    diffmat = np.abs(mat - val)
+    mindiff = np.min(diffmat)
+    if np.sum(diffmat == mindiff) == 1:
+        ijindex = index_mat(diffmat,mindiff)
+        return ijindex
+
+    if Ndims > 2: dz = mat.shape[Ndims-3]
+    if Ndims > 3: dt = mat.shape[Ndims-4]
+
+    if way == 'anticlockwise12':
+# Sorting following a double criteria, first angle and then distance taken from the
+#   center in the anti-clockwise sens starting at 12
+        if Ndims > 2:
+            print errormsg
+            print '  ' + fame + ': wiht more than 2-dims:', Ndims, "'" + way +       \
+              "' is not possible !!"
+            quit(-1)
+ 
+        distrad = radius_dist(dx,dy,dx/2,dy/2)
+        distangle = radius_angle(dx,dy,dx/2,dy/2)
+        distangle = np.where(distangle >= np.pi, distangle - 2.*np.pi, distangle)
+        distangle = np.where(distangle > 0., 2.* np.pi - distangle, -distangle)
+
+        sortedij = {}
+        minangle = 10000.
+        notfound = np.zeros((dy,dx), dtype=bool)
+        maskdistrad = ma.array(distrad, mask=notfound)
+        maskdistangle = ma.array(distangle, mask=notfound)
+        for ip in range(dx*dy):
+            minangle = np.min(maskdistangle)
+            jiangle = multi_index_mat(maskdistangle, minangle)
+            mindist = 10000.
+            for ia in range(len(jiangle)):
+                iaji = jiangle[ia]
+                if maskdistrad[iaji[0], iaji[1]] < mindist: 
+                    mindist = maskdistrad[iaji[0], iaji[1]]
+                    idistangle = iaji
+            notfound[idistangle[0],idistangle[1]] = True
+            sortedij[ip] = idistangle
+
+    elif way == 'anticlockwise6':
+# Sorting following a double criteria, first angle and then distance taken from the
+#   center in the anti-clockwise sens starting at 6
+        if Ndims > 2:
+            print errormsg
+            print '  ' + fame + ': wiht more than 2-dims:', Ndims, "'" + way +       \
+              "' is not possible !!"
+            quit(-1)
+ 
+        distrad = radius_dist(dx,dy,dx/2,dy/2)
+        distangle = radius_angle(dx,dy,dx/2,dy/2)
+        distangle = np.where(distangle > np.pi, distangle - 2.*np.pi, distangle)
+        distangle = np.where(distangle >= 0., np.pi - distangle, np.pi-distangle)
+
+        sortedij = {}
+        minangle = 10000.
+        notfound = np.zeros((dy,dx), dtype=bool)
+        maskdistrad = ma.array(distrad, mask=notfound)
+        maskdistangle = ma.array(distangle, mask=notfound)
+        for ip in range(dx*dy):
+            minangle = np.min(maskdistangle)
+            jiangle = multi_index_mat(maskdistangle, minangle)
+            mindist = 10000.
+            for ia in range(len(jiangle)):
+                iaji = jiangle[ia]
+                if maskdistrad[iaji[0], iaji[1]] < mindist: 
+                    mindist = maskdistrad[iaji[0], iaji[1]]
+                    idistangle = iaji
+            notfound[idistangle[0],idistangle[1]] = True
+            sortedij[ip] = idistangle
+
+    elif way == 'clockwise12':
+# Sorting following a double criteria, first angle and then distance taken from the
+#   center in the cloc-kwise sens starting at 12
+        if Ndims > 2:
+            print errormsg
+            print '  ' + fame + ': wiht more than 2-dims:', Ndims, "'" + way +       \
+              "' is not possible !!"
+            quit(-1)
+ 
+        distrad = radius_dist(dx,dy,dx/2,dy/2)
+        distangle = radius_angle(dx,dy,dx/2,dy/2)
+
+        sortedij = {}
+        minangle = 10000.
+        notfound = np.zeros((dy,dx), dtype=bool)
+        maskdistrad = ma.array(distrad, mask=notfound)
+        maskdistangle = ma.array(distangle, mask=notfound)
+        for ip in range(dx*dy):
+            minangle = np.min(maskdistangle)
+            jiangle = multi_index_mat(maskdistangle, minangle)
+            mindist = 10000.
+            for ia in range(len(jiangle)):
+                iaji = jiangle[ia]
+                if maskdistrad[iaji[0], iaji[1]] < mindist: 
+                    mindist = maskdistrad[iaji[0], iaji[1]]
+                    idistangle = iaji
+            notfound[idistangle[0],idistangle[1]] = True
+            sortedij[ip] = idistangle
+
+    elif way == 'clockwise6':
+# Sorting following a double criteria, first angle and then distance taken from the
+#   center in the clockwise sens starting at 6
+        if Ndims > 2:
+            print errormsg
+            print '  ' + fame + ': wiht more than 2-dims:', Ndims, "'" + way +       \
+              "' is not possible !!"
+            quit(-1)
+ 
+        distrad = radius_dist(dx,dy,dx/2,dy/2)
+        distangle = radius_angle(dx,dy,dx/2,dy/2)
+        distangle = np.where(distangle >= np.pi, distangle - 2.*np.pi, distangle)
+        distangle = np.where(distangle > 0., np.pi + distangle, np.pi+distangle)
+
+        sortedij = {}
+        minangle = 10000.
+        notfound = np.zeros((dy,dx), dtype=bool)
+        maskdistrad = ma.array(distrad, mask=notfound)
+        maskdistangle = ma.array(distangle, mask=notfound)
+        for ip in range(dx*dy):
+            minangle = np.min(maskdistangle)
+            jiangle = multi_index_mat(maskdistangle, minangle)
+            mindist = 10000.
+            for ia in range(len(jiangle)):
+                iaji = jiangle[ia]
+                if maskdistrad[iaji[0], iaji[1]] < mindist: 
+                    mindist = maskdistrad[iaji[0], iaji[1]]
+                    idistangle = iaji
+            notfound[idistangle[0],idistangle[1]] = True
+            sortedij[ip] = idistangle
+
+    else:
+        print errormsg
+        print '  ' + fname + ": way of search '" + way + "' not ready !!"
+        print '    available ways:', ways
+
+    indices = []
+    for ij in range(dx*dy):
+        sij = sortedij[ij]
+        sval = mat[sij[0], sij[1]]
+        if sval == val:
+            indices.append(sij)
+
+    return indices
+
+mat = np.arange(20).reshape(4,5)
+mat[2,1] = 3
+print index_mat_way(mat,3,'anticlockwise6')
+
 #quit()