Changeset 2580 in lmdz_wrf

Timestamp:

Jun 2, 2019, 5:47:28 PM (6 years ago)

Author:

lfita

Message:

Adding:

• `val_between': Function to provide if a given value is between two consecutive ones

File:

• trunk/tools/geometry_tools.py (modified) (7 diffs)

trunk/tools/geometry_tools.py

@@ -43 +43 @@
 # spheric_line: Function to transform a series of locations in lon, lat coordinates 
 #   to x,y,z over an 3D spaceFunction to provide coordinates of a line  on a 3D space
+# val_between: Function to provide if a given value is between two consecutive ones
 # write_join_poly: Function to write an ASCII file with the combination of polygons
 
@@ -435 +436 @@
     return polys
 
+def val_between(valA, valB, val):
+    """ Function to provide if a given value is between two consecutive ones
+      valA: first value
+      valB: second value
+      val: value to determine if it is between
+      >>> val_between(0.5,1.5,0.8)
+      True
+      >>> val_between(0.5,1.5.,-0.8)
+      False
+      >>> val_between(0.5,1.5,0.5)
+      True
+    """
+    fname = 'val_between'
+
+    btw = False
+    if (valA <= val and valB > val) or (valA < val and valB >= val): btw =True
+
+    return btw
+
 def cut_ypolygon(polygon, yval, keep='bottom', Nadd=20):
     """ Function to cut a polygon from a given value of the y-axis
@@ -465 +485 @@
     ept = []
 
-    if type(polygon) == type(gen.mamat):
+    if type(polygon) == type(gen.mamat) and type(polygon.mask) !=                    \
+      type(gen.mamat.mask[1]):
         # Assuming clockwise polygons
         for ip in range(N-1):
@@ -472 +493 @@
                 if eep == N: eep = 0
       
-                if polygon[ip,0] <= yval and polygon[eep,0] >= yval:
+                if val_between(polygon[ip,0], polygon[eep,0], yval):
                     icut.append(ip)
                     dx = polygon[eep,1] - polygon[ip,1]
@@ -485 +506 @@
                     dd = yval - polygon[ip,0]
                     ept.append([yval, polygon[ip,1]+dx*dd/dy])
-                    Ncuts = Nctus + 1
+                    Ncuts = Ncuts + 1
     else:
         # Assuming clockwise polygons
@@ -505 +526 @@
                 dd = yval - polygon[ip,0]
                 ept.append([yval, polygon[ip,1]+dx*dd/dy])
-
-    if ipt is None or ept is None:
+                Ncuts = Ncuts + 1
+
+    if ipt is None or ept is None or Ncuts == 0:
         print errormsg
         print '  ' + fname + ': no cutting for polygon at y=', yval, '!!'
+    else:
+        print '  ' + fname + ': found ', Ncuts, ' Ncuts'
+        print '    yval=', yval, 'cut, ip; ipt ep; ept ________'
+        for ic in range(Ncuts):
+            print ic, icut[ic], ';', ipt[ic], ecut[ic], ';', ept[ic] 
 
     Nadds = []
-    Naddc = Nadd/(Ncuts-1)
-    for ic in range(Ncuts-1):
-        Nadds.append(Naddc)
-
-    Nadds.append(N-Naddc*(Ncuts-1))
+    if Ncuts > 1:
+        Naddc = Nadd/(Ncuts-1)
+        for ic in range(Ncuts-1):
+            Nadds.append(Naddc)
+
+        Nadds.append(N-Naddc*(Ncuts-1))
+    else:
+        Nadds.append(Nadd)
 
     iip = 0
@@ -550 +580 @@
 
     rmpolygon = []
+    Npts = cutpolygon.shape[0]
     if keep == 'bottom':
         for ip in range(Npts):