Changeset 832 in lmdz_wrf for trunk

Timestamp:

Jun 16, 2016, 1:13:57 PM (8 years ago)

Author:

lfita

Message:

Adding `subbasin_point_orig': as the original version and direction search
Some corrections

File:

• trunk/tools/generic_tools.py (modified) (12 diffs)

trunk/tools/generic_tools.py

@@ -7062 +7062 @@
 #quit()
 
-def subbasin_point(trips, pt, direction):
+def index_mat_way(mat,val,way):
+    """ Function to look for a value within a matrix following a way and direction
+      mat: matrix
+      val: value to search
+      way: way of search
+        'anticlockwise12': looking in an anti-clockwise direction starting at 12
+        'anticlockwise6': looking in an anti-clockwise direction starting at 6
+        'clockwise12': looking in an clockwise direction starting at 12
+        'clockwise6': looking in an clockwise direction starting at 6
+    >>> 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 direction 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 direction 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 direction 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 direction 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
+
+def subbasin_point_orig(trips, pt):
     """ Function to provide sub-basins given a grid point following a matrix of trips
       trips= matrix of trips values
         1: N,  2: NE,  3: E,  4: SE,  5: S,  6: SW:  7: W,  8: NW
       pt= point within the trip matrix
-      direction= searching direction river-up ('left': clockwise, 'right': anti-clockwise)
     """
     fname = 'subbasin_point'
@@ -7085 +7253 @@
     Nij = 1
 
-    subbasin[pt[0], pt[1]] = True 
+    subbasin[pt[0], pt[1]] = True
     ijsubbasin[Nij] = np.array(pt)
     ijfound[Nij] = False
@@ -7096 +7264 @@
         Njipt = dictionary_key(subfinished, False)
         if Njipt is not None:
+#            print Njipt, 'points subflow:', ijsubflow[Njipt]
             for ipt in ijsubflow[Njipt]:
                 if ijfound[ipt] == False:
                     jipt = ijsubbasin[ipt]
                     break
+#            print '  working from point:', ipt, 'ji pair:', jipt
             ijfound[ipt] = True
             Nij = ipt
         else:
 # Keep searching since there are grid-points not found!
+#            print '  ' + fname + ': Keep searching since there are grid-points not found!!'
+#            print ijfound
             Nij = dictionary_key(ijfound, False)
             if Nij is None:
@@ -7112 +7284 @@
             parentNjipt = dictionary_key_list(ijsubflow, Nij)
             Njipt = chainSnum_levnext(parentNjipt, ijsubflow.keys())
+#            print '  ' + fname + "new sub-flow '" + Njipt + "' !!"
             subfinished[Njipt] = False
             Nij = np.max(ijfound.keys())
 
+#        print 'Lluis Nij:', Nij
 # Looking for which point of the sub-flow retake the search
         if Nij == 1:
@@ -7121 +7295 @@
 
         ardtrips = vals_around(trips,jipt)
+#        print '  ' + fname + 'ardtrips _______'
+#        print ardtrips
 
         arrive = incomming_flow(ardtrips)
         Narrive = np.sum(arrive)
+#        print Nij, '  ' + fname + ' Narrive:', Narrive
         if Narrive == 0:
             ijfound[Nij] = True
             subfinished[Njipt] = True
         else:
+#            print '  ' + fname + 'arrive _______'
+#            print arrive
             followvals = np.zeros((3,3), dtype=bool)
             followvals = arrive
@@ -7133 +7312 @@
 
             for ifollow in range(Narrive):
+#                print 'ifollow:',ifollow,'/',Narrive
 
 # We only want to work with that ij, which have not yet been found
                 while np.sum(followvals) != 0:
-
-# Searching direction river-up. 'left': starting from the left, 'right' starting from the right
-                    if direction == 'left':
-                        jifollow = index_mat_way(followvals, True, 'clockwise6')
-                    elif direction == 'right':
-                        jifollow = index_mat_way(followvals, True, 'anticlockwise6')
-
-                    jiequiv = jifollow[ifollow] - [1,1]
+#                    print Nij,'  Looking for a non-located point in subbasin ________'
+#                    print subbasin[jipt[0]-1:jipt[0]+2,jipt[1]-1:jipt[1]+2]
+                    jifollow = index_mat(followvals, True)
+                    jiequiv = jifollow - [1,1]
                     if subbasin[jipt[0]+jiequiv[0], jipt[1]+jiequiv[1]] == False:
                         Nij = np.max(ijfound.keys()) + 1
                         jiptnew = jipt + jiequiv
-                        if ifollow != 0: 
+                        if ifollow != 0:
 # Avoiding repetition of sub-flow name
                             if len(ijsubflow.keys()) > 1:
                                 Njipt = chainSnum_levnext(parentNjipt, ijsubflow.keys())
                             else:
-                                # First subflow
                                 Njipt = '11'
                         else:
                             Njipt = parentNjipt
+#                        print '  ' + fname + "new sub-flow '" + Njipt + "' !!"
                         subfinished[Njipt] = False
                         subbasin[jipt[0]+jiequiv[0], jipt[1]+jiequiv[1]] = True
@@ -7172 +7348 @@
                             ijfound[Nij] = True
                             subfinished[Njipt] = True
+#                            print Nij,"  subflow '" + Njipt + "' finished!!"
                             break
 
@@ -7179 +7356 @@
                       'the number of trips', TOTtrips,'!!'
                     quit()
+
+#            if ijsubflow.has_key(Njipt): 
+#                print "Lluis points of subflow: '" + Njipt + "' _______=", ijsubflow[Njipt]
+#                for isub in ijsubflow[Njipt]:
+#                    print '  ' , isub , ':', ijsubbasin[isub], ijfound[isub]
+#            if Nij == 10: print 'Nij = 9:', ijfound[9]
+
+#        print subbasin
+#        if Nij > 4: quit()
+
+    return subbasin
+
+def subbasin_point(trips, pt, direction):
+    """ Function to provide sub-basins given a grid point following a matrix of trips
+      trips= matrix of trips values
+        1: N,  2: NE,  3: E,  4: SE,  5: S,  6: SW:  7: W,  8: NW
+      pt= point within the trip matrix
+      direction= searching direction river-up ('left': clockwise, 'right': anti-clockwise)
+    """
+    fname = 'subbasin_point'
+
+    subbasin = np.zeros((trips.shape), dtype=bool)
+    TOTtrips = len(trips.flatten())
+
+# Dictionary with the j,i points of all the subbasin
+    ijsubbasin = {}
+# Dictionary to tell if a given j,i sub-basin point has been localized (as number from ijsubbasin)
+    ijfound = {}
+# Dictionary to tell if a given sub-flow has been finished (when Narrive == 0)
+    subfinished = {}
+# Dictionary with the list of points which are contained in a sub-flow
+    ijsubflow = {}
+# Number of ji sub-basin points
+    Nij = 1
+
+    subbasin[pt[0], pt[1]] = True
+    ijsubbasin[Nij] = np.array(pt)
+    ijfound[Nij] = False
+    ijsubflow[str(Nij)] = [Nij]
+    subfinished[str(Nij)] = False
+
+# Loop around all sub-flows
+##
+    while np.sum(subfinished.values) != 0 and np.sum(ijfound.values) != 0:
+        Njipt = dictionary_key(subfinished, False)
+        if Njipt is not None:
+#            print Njipt, 'points subflow:', ijsubflow[Njipt]
+            for ipt in ijsubflow[Njipt]:
+                if ijfound[ipt] == False:
+                    jipt = ijsubbasin[ipt]
+                    break
+#            print '  working from point:', ipt, 'ji pair:', jipt
+            ijfound[ipt] = True
+            Nij = ipt
+        else:
+# Keep searching since there are grid-points not found!
+#            print '  ' + fname + ': Keep searching since there are grid-points not found!!'
+#            print ijfound
+            Nij = dictionary_key(ijfound, False)
+            if Nij is None:
+                return subbasin, ijsubflow, ijsubbasin
+
+            ijfound[Nij] = True
+            jipt = ijsubbasin[Nij]
+            parentNjipt = dictionary_key_list(ijsubflow, Nij)
+            Njipt = chainSnum_levnext(parentNjipt, ijsubflow.keys())
+#            print '  ' + fname + "new sub-flow '" + Njipt + "' !!"
+            subfinished[Njipt] = False
+            Nij = np.max(ijfound.keys())
+
+#        print 'Lluis Nij:', Nij
+# Looking for which point of the sub-flow retake the search
+        if Nij == 1:
+            ipt = 1
+            jipt = pt
+
+        ardtrips = vals_around(trips,jipt)
+#        print '  ' + fname + 'ardtrips _______'
+#        print ardtrips
+
+        arrive = incomming_flow(ardtrips)
+        Narrive = np.sum(arrive)
+#        print Nij, '  ' + fname + ' Narrive:', Narrive
+        if Narrive == 0:
+            ijfound[Nij] = True
+            subfinished[Njipt] = True
+        else:
+#            print '  ' + fname + 'arrive _______'
+#            print arrive
+            followvals = np.zeros((3,3), dtype=bool)
+            followvals = arrive
+            parentNjipt = Njipt
+
+            for ifollow in range(Narrive):
+#                print 'ifollow:',ifollow,'/',Narrive
+
+# We only want to work with that ij, which have not yet been found
+                while np.sum(followvals) != 0:
+#                    print Nij,'  Looking for a non-located point in subbasin ________'
+#                    print subbasin[jipt[0]-1:jipt[0]+2,jipt[1]-1:jipt[1]+2]
+
+# Searching direction river-up. 'left': starting from the left, 'right' starting from the right
+                    if direction == 'left':
+                        jifollow0 = index_mat_way(followvals, True, 'clockwise6')
+                    elif direction == 'right':
+                        jifollow0 = index_mat_way(followvals, True, 'anticlockwise6')
+                    if len(jifollow0) > 1 and type(jifollow0) == type([2]):
+                        jifollow = jifollow0[0]
+                    else:
+                        jifollow = jifollow0
+## orig ##                   jifollow = index_mat(followvals, True)
+                    jiequiv = jifollow - [1,1]
+                    if subbasin[jipt[0]+jiequiv[0], jipt[1]+jiequiv[1]] == False:
+                        Nij = np.max(ijfound.keys()) + 1
+                        jiptnew = jipt + jiequiv
+                        if ifollow != 0:
+# Avoiding repetition of sub-flow name
+                            if len(ijsubflow.keys()) > 1:
+                                Njipt = chainSnum_levnext(parentNjipt, ijsubflow.keys())
+                            else:
+                                Njipt = '11'
+                        else:
+                            Njipt = parentNjipt
+#                        print '  ' + fname + "new sub-flow '" + Njipt + "' !!"
+                        subfinished[Njipt] = False
+                        subbasin[jipt[0]+jiequiv[0], jipt[1]+jiequiv[1]] = True
+                        ijfound[Nij] = False
+                        ijsubbasin[Nij] = jiptnew
+                        if ijsubflow.has_key(Njipt):
+                            subflowpts = ijsubflow[Njipt]
+                            subflowpts.append(Nij)
+                            ijsubflow[Njipt] = subflowpts
+                        else:
+                            ijsubflow[Njipt] = [Nij]
+                        break
+                    else:
+                        followvals[jifollow[0],jifollow[1]] = False
+                        if np.sum(followvals) == 0:
+                            ijfound[Nij] = True
+                            subfinished[Njipt] = True
+#                            print Nij,"  subflow '" + Njipt + "' finished!!"
+                            break
+
+                if Nij > TOTtrips:
+                    print errormsg
+                    print '  ' + fname + ': sub-flow point', Nij, 'larger than ' +   \
+                      'the number of trips', TOTtrips,'!!'
+                    quit()
+
+#            if ijsubflow.has_key(Njipt): 
+#                print "Lluis points of subflow: '" + Njipt + "' _______=", ijsubflow[Njipt]
+#                for isub in ijsubflow[Njipt]:
+#                    print '  ' , isub , ':', ijsubbasin[isub], ijfound[isub]
+#            if Nij == 10: print 'Nij = 9:', ijfound[9]
+
+#        print subbasin
+#        if Nij > 4: quit()
 
     return subbasin
@@ -7206 +7540 @@
 #print rivers
 #print rivers[point[0], point[1]]
-
-#masksubbasin, subflows, subflowspt = subbasin_point(rivers, point)
+#direc = 'left'
+
+#masksubbasin, subflows, subflowspt = subbasin_point_dir(rivers, point, 'right')
 #print rivers
 #print masksubbasin
@@ -7214 +7549 @@
 #  for isub in subflows[subflow]:
 #      print '  ' , isub , ':', subflowspt[isub]
-#
+
 #print 'Total sub-basin:', np.sum(masksubbasin)
+#quit()
 
 def color_deg(colors, Nsteps, typedeg, values=None):
@@ -7248 +7584 @@
         cols[:,icol] = np.array(colors[icol])
     
-    if typedeg == 'std':
-        Nstepcols = np.int(Nsteps*1. / (Ncols-1.))
-        iicol = 0
-        for icol in range(Ncols-1):
-            dcolors = (cols[:,icol+1]- cols[:,icol])/(Nstepcols-1.)
+    if Ncols > Nsteps:
+        print warnmsg
+        print '  ' + fname + ': too few steps:', Nsteps, 'for the number of colors:',\
+          Ncols,'!!'
+        print '    directly assigning the given color'
+        Nfincol = Nsteps
+        if Nsteps > 1:
+            for istp in range(Nsteps):
+                colordeg[0,istp] = istp / (Nsteps-1.)
+                colordeg[1:4,istp] = cols[:,istp]
+        else:
+            colordeg[0,0] = 1.
+            colordeg[1:4,0] = cols[:,0]
+    else:
+        if typedeg == 'std':
+            Nstepcols = np.int(Nsteps*1. / (Ncols-1.))
+            iicol = 0
+            for icol in range(Ncols-1):
+                dcolors = (cols[:,icol+1]- cols[:,icol])/(Nstepcols-1.)
+# Toavoid repetition of the color
+                if icol == 0: 
+                    ibeg = 0
+                else:
+                    ibeg = 1
+                for i in range(ibeg,Nstepcols):
+                    colordeg[0,iicol] = iicol*1.
+                    colordeg[1:4,iicol] = cols[:,icol] + dcolors*i
+                    iicol = iicol + 1
+            Nfincol = iicol-1
+            colordeg[0,:] = colordeg[0,:]/(Nfincol*1.)
+
+        elif typedeg == 'vals':
+            if values is None:
+                print errormsg
+                print '  ' + fname + ": type '" + typedeg + "' require values !!"
+                quit(-1)
+            Nvals = len(values)
+            if Nvals != Ncols:
+                print errormsg
+                print '  ' + fname + ': different number of values:', Nvals,         \
+                  'and colors:', Ncols, '!!'
+                print '    provided values:', values
+                print '    provided colors:', cols
+                quit(-1)
+            TOTrange = values[Nvals-1] - values[0]
+            iicol = 0
+            for icol in range(Ncols-1):
+                pairange = values[icol+1] - values[icol]
+                drange = pairange / TOTrange
+                Nstepcols = np.int(Nsteps*drange)
+                dcolors = (cols[:,icol+1]-cols[:,icol])/(Nstepcols-1.)
 # To avoid repetition of the color
-            if icol == 0: 
-                ibeg = 0
-            else:
-                ibeg = 1
-            for i in range(ibeg,Nstepcols):
-                colordeg[0,iicol] = iicol*1.
-                colordeg[1:4,iicol] = cols[:,icol] + dcolors*i
-                iicol = iicol + 1
-        Nfincol = iicol-1
-        colordeg[0,:] = colordeg[0,:]/(Nfincol*1.)
-
-    elif typedeg == 'vals':
-        if values is None:
+                if icol == 0: 
+                    ibeg = 0
+                else:
+                    ibeg = 1
+                for i in range(ibeg,Nstepcols):
+                    colordeg[0,iicol] = values[icol]+i*pairange/(Nstepcols-1)
+                    colordeg[1:4,iicol] = cols[:,icol] + dcolors*i
+                    iicol = iicol + 1
+            Nfincol = iicol - 1
+        else:
             print errormsg
-            print '  ' + fname + ": type '" + typedeg + "' require values !!"
+            print '  ' + fname + ": degradation type '" + typedeg + "' not ready !!"
             quit(-1)
-        Nvals = len(values)
-        if Nvals != Ncols:
-            print errormsg
-            print '  ' + fname + ': different number of values:', Nvals,             \
-              'and colors:', Ncols, '!!'
-            print '    provided values:', values
-            print '    provided colors:', cols
-            quit(-1)
-        TOTrange = values[Nvals-1] - values[0]
-        iicol = 0
-        for icol in range(Ncols-1):
-            pairange = values[icol+1] - values[icol]
-            drange = pairange / TOTrange
-            Nstepcols = np.int(Nsteps*drange)
-            dcolors = (cols[:,icol+1]-cols[:,icol])/(Nstepcols-1.)
-# To avoid repetition of the color
-            if icol == 0: 
-                ibeg = 0
-            else:
-                ibeg = 1
-            for i in range(ibeg,Nstepcols):
-                colordeg[0,iicol] = values[icol]+i*pairange/(Nstepcols-1)
-                colordeg[1:4,iicol] = cols[:,icol] + dcolors*i
-                iicol = iicol + 1
-        Nfincol = iicol - 1
-    else:
-        print errormsg
-        print '  ' + fname + ": degradation type '" + typedeg + "' not ready !!"
-        quit(-1)
 
     degcolor = np.zeros((4,Nfincol), dtype=np.float)
@@ -7348 +7698 @@
     return hexcol
 
-def index_mat_way(mat,val,way):
-    """ Function to look for a value within a matrix following a way and direction
-      mat: matrix
-      val: value to search
-      way: way of search
-        'anticlockwise12': looking in an anti-clockwise direction starting at 12
-        'anticlockwise6': looking in an anti-clockwise direction starting at 6
-        'clockwise12': looking in an clockwise direction starting at 12
-        'clockwise6': looking in an clockwise direction starting at 6
-    >>> 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 direction 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 direction 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 direction 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 direction 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
-
 #quit()