Changeset 813 in lmdz_wrf

Timestamp:

Jun 12, 2016, 12:55:17 AM (8 years ago)

Author:

lfita

Message:

Adding:

• `subbasin_point': Function to provide sub-basins given a grid point following a matrix of trips
• `contflow': Function to bring back the increment in j,i grid points according to a trip: (inflow directions)
• `dictionary_key': Function to provide the first key in a dictionay with a given value
• `dictionary_key_list': Function to provide the first key in a dictionay of lists with a given value
• `incomming_flow': Function to determine if a fgrid-flow inflows to the central grid point
• `vals_around': Function to provide the 3x3 values around a given j,i point

File:

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

trunk/tools/generic_tools.py

@@ -11 +11 @@
 
 fillValue = 1.e20
-fillValueF = 1.e20
+fillValueF = 1.e20 
 fillValueC = '-'
 fillValueI = -99999
@@ -19 +19 @@
 
 ####### Content
+# contflow: Function to bring back the increment in j,i grid points according to a trip: (inflow directions)
+# dictionary_key: Function to provide the first key in a dictionay with a given value
+# dictionary_key_list: Function to provide the first key in a dictionay of lists with a given value
 # 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
 # PolyArea: Function to compute the area of the polygon following 'Shoelace formula'
 # 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
+# vals_around: Function to provide the 3x3 values around a given j,i point
+# subbasin_point: Function to provide sub-basins given a grid point following a matrix of trips
 
 def values_line(line, splitv, chars):
@@ -6460 +6466 @@
 #print unitsDate('19490101000000', '19760217082932', 'minute')
 
-#quit()
+def incomming_flow(dirs):
+    """ Function to determine if a fgrid-flow inflows to the central grid point
+      dirs = [3,3] matrix with the trip: (inflow directions)
+          1: N,  2: NE,  3: E,  4: SE,  5: S,  6: SW:  7: W,  8: NW
+        being:
+          0,0: NW point, 0,1: N point, 0,2: NE point, 
+          1,0: W point, 1,1: grid-point, 1,2: E point,
+          2,0: SW point, 2,1: S point, 2,2: SE point
+    >>> dirflow = np.array([5, 5, 5, 98, 7, 7, 1, 1, 1], dtype=int).reshape(3,3)
+    >>> print dirflow
+    [[5 5 5]
+     [98 7 7]
+     [1 1 1]]
+    >>> print incomming_flow(dirflow)
+    [[False  True False]
+     [False False  True]
+     [False  True False]]
+    """
+    fname = 'incomming_flow'
+
+    inflow = np.zeros((3,3), dtype=bool)
+
+    rightinflow = np.zeros((3,3), dtype=int)
+    rightinflow[0,0] = 4
+    rightinflow[0,1] = 5
+    rightinflow[0,2] = 6
+    rightinflow[1,0] = 3
+    rightinflow[1,1] = 0
+    rightinflow[1,2] = 7
+    rightinflow[2,0] = 2
+    rightinflow[2,1] = 1
+    rightinflow[2,2] = 8
+
+    inflow = np.where(dirs == rightinflow, True, False)
+
+    return inflow
+
+def contflow(dirflow):
+    """ Function to bring back the increment in j,i grid points according to a trip: (inflow directions)
+          1: N,  2: NE,  3: E,  4: SE,  5: S,  6: SW:  7: W,  8: NW
+    >>> contflow(2)
+    [1 1]
+    >>> contflow(7)
+    [0 -1]
+    >>> contflow(98)
+    [-1 -1]
+    """
+    fname = 'contflow'
+    inc = np.ones((2), dtype=int)*(-1)
+
+    if dirflow == 1: #N
+        inc[0] = 1
+        inc[1] = 0
+    elif dirflow == 2: #NE
+        inc[0] = 1
+        inc[1] = 1
+    elif dirflow == 3: #E
+        inc[0] = 0
+        inc[1] = 1
+    elif dirflow == 4: #SE
+        inc[0] = -1
+        inc[1] = 1
+    elif dirflow == 5: #S
+        inc[0] = -1
+        inc[1] = 0
+    elif dirflow == 6: #SW
+        inc[0] = -1
+        inc[1] = -1
+    elif dirflow == 7: #W
+        inc[0] = 0
+        inc[1] = -1
+    elif dirflow == 8: #NW
+        inc[0] = 1
+        inc[1] = -1
+    else: #
+        inc[0] = -1
+        inc[1] = -1
+
+    return inc
+
+def vals_around(vals,pt):
+    """ Function to provide the 3x3 values around a given j,i point
+      vals = matrix of values
+      pt = j,i point to serch around
+    >>> vals_around(np.arange(16).reshape(4,4), [3, 3])
+    [[10.0 11.0 None]
+     [14.0 15.0 None]
+     [None None None]]
+    """
+    fname = 'vals_around'
+
+    vals3x3 = np.array([None]*9).reshape(3,3)
+
+    dx = vals.shape[1]
+    dy = vals.shape[0]
+
+    nx = np.max([0,pt[1]-1])
+    xx = np.min([dx-1,pt[1]+1]) + 1
+    ny = np.max([0,pt[0]-1])
+    xy = np.min([dy-1,pt[0]+1]) + 1
+
+    vals3x3[1-(pt[0]-ny):1+(xy-pt[0]),1-(pt[1]-nx):1+(xx-pt[1])] = vals[ny:xy,nx:xx]
+
+    return vals3x3
+
+def dictionary_key(dictv, val):
+    """ Function to provide the first key in a dictionay with a given value
+      ditcv= dictionary
+      val= value to search
+    >>> dictionary_key({'i': 1, 'ii': 2, 'iii': 3, 'iv': 4}, 3)
+    'iii'
+    """
+    fname = 'dictionary_key'
+
+    keyval = None
+    for keyv in dictv.keys():
+        if dictv[keyv] == val:
+            keyval = keyv
+            break
+
+    return keyval
+
+def dictionary_key_list(dictv, val):
+    """ Function to provide the first key in a dictionay of lists with a given value
+      ditcv= dictionary
+      val= value to search
+    >>> dictionary_key_list({'i': [1, 0, -1], 'ii': [2, 8, 16], 'iii': [-3, 3, 9], 'iv': [4]}, 3)
+    'iii'
+    """
+    fname = 'dictionary_key'
+
+    keyval = None
+    for keyv in dictv.keys():
+        if searchInlist(dictv[keyv],val):
+            keyval = keyv
+            break
+
+    return keyval
+
+def subbasin_point(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
+    """
+    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]
+            Njipt = dictionary_key_list(ijsubflow, Nij)
+            if subfinished.has_key(Njipt + '1'): 
+                if subfinished.has_key(Njipt + '01'): 
+                    if subfinished.has_key(Njipt + '001'): 
+                        Njipt = Njipt + '0001'
+                    else:
+                        Njipt = Njipt + '001'
+                else:
+                    Njipt = Njipt + '01'
+            else:
+                Njipt = Njipt + '1'
+            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
+
+            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]
+                    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 subfinished.has_key(Njipt + str(ifollow)): 
+                                Njipt = Njipt + '0' + str(ifollow)
+                            else:
+                                Njipt = Njipt + str(ifollow)
+                            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
+
+# Caceres
+rivers = np.array([7, 1, 1, 1, 1, 5, 5, 1, \
+  7, 1, 6, 3, 5, 5, 6, 5, \
+  8, 7, 5, 3, 5, 6, 3, 5, \
+  1, 1, 3, 5, 5, 7, 5, 5, \
+  3, 3, 5, 6, 5, 5, 5, 7 ], dtype=int).reshape(5,8)
+
+point = [4,4]
+
+# Porto Do Alegre
+#rivers = np.array([5, 5, 1, 7, 1, 8, 8, 2, \
+#  5, 6, 5, 7, 7, 1, 3, 1, \
+#  6, 3, 5, 1, 3, 3, 3, 3, \
+#  7, 5, 5, 5, 5, 5, 5, 1, \
+#  5, 5, 7, 7, 7, 6, 7, 7, \
+#  5, 5, 7, 7, 7, 7, 8, 7, \
+#  6, 6, 7, 7, 7, 7, 7, 7, \
+#  5, 7, 7, 7, 8, 7, 7, 6, \
+#  6, 7, 7, 7, 7, 7, 7, 6], dtype=int).reshape(9,8)
+
+#point = [7,0]
+
+print rivers
+print rivers[point[0], point[1]]
+
+masksubbasin, subflows, subflowspt = subbasin_point(rivers, point)
+print rivers
+print masksubbasin
+for subflow in subflows:
+  print subflow, ':', subflows[subflow]
+  for isub in subflows[subflow]:
+      print '  ' , isub , ':', subflowspt[isub]
+
+print 'Total sub-basin:', np.sum(masksubbasin)
+
+quit()