Changeset 2545 in lmdz_wrf

Timestamp:

May 19, 2019, 3:23:05 PM (6 years ago)

Author:

lfita

Message:

Adding:

• `cut_ypolygon': Function to cut a polygon from a given value of the y-axis

File:

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

trunk/tools/geometry_tools.py

@@ -17 +17 @@
 import generic_tools as gen
 import numpy.ma as ma
+import module_ForSci as sci
 
 errormsg = 'ERROR -- error -- ERROR -- error'
@@ -22 +23 @@
 
 ####### Contents:
+# cut_ypolygon: Function to cut a polygon from a given value of the y-axis
 # deg_deci: Function to pass from degrees [deg, minute, sec] to decimal angles [rad]
 # dist_points: Function to provide the distance between two points
@@ -420 +422 @@
 
     return polys
+
+def cut_ypolygon(polygon, yval, keep='bottom', Nadd=20):
+    """ Function to cut a polygon from a given value of the y-axis
+      polygon: polygon to cut
+      yval: value to use to cut the polygon
+      keep: part to keep from the height ('bottom', default)
+      Nadd: additional points to add to draw the line (20, default)
+    """
+    fname = 'cut_ypolygon'
+
+    N = polygon.shape[0]
+
+    ipt = None
+    ept = None
+
+    if type(polygon) == type(gen.mamat):
+        # Assuming clockwise polygons
+        for ip in range(N):
+            if not polygon.mask[ip,0]:
+                eep = ip + 1
+                if eep == N: eep = 0
+      
+                if polygon[ip,0] <= yval and polygon[eep,0] >= yval:
+                    icut = ip
+                    dx = polygon[eep,1] - polygon[ip,1]
+                    dy = polygon[eep,0] - polygon[ip,0]
+                    dd = yval - polygon[ip,0]
+                    ipt = [yval, polygon[ip,1]+dx*dd/dy]
+
+                if polygon[ip,0] >= yval and polygon[eep,0] <= yval:
+                    ecut = ip
+                    dx = polygon[eep,1] - polygon[ip,1]
+                    dy = polygon[eep,0] - polygon[ip,0]
+                    dd = yval - polygon[eep,0]
+                    ept = [yval, polygon[eep,1]+dx*dd/dy]
+    else:
+        # Assuming clockwise polygons
+        for ip in range(N):
+            eep = ip + 1
+            if eep == N: eep = 0
+      
+            if polygon[ip,0] <= yval and polygon[eep,0] >= yval:
+                icut = ip
+                dx = polygon[eep,1] - polygon[ip,1]
+                dy = polygon[eep,0] - polygon[ip,0]
+                dd = yval - polygon[ip,0]
+                ipt = [yval, polygon[ip,1]+dx*dd/dy]
+                print ip, 'Lluis pi:', polygon[ip,:], 'pi+1', polygon[eep,:]
+                print '    yval', yval, 'dx', dx, 'dy', dy
+                print '      ipt', ipt
+
+            if polygon[ip,0] >= yval and polygon[eep,0] <= yval:
+                ecut = ip
+                dx = polygon[eep,1] - polygon[ip,1]
+                dy = polygon[eep,0] - polygon[ip,0]
+                dd = yval - polygon[ip,0]
+                ept = [yval, polygon[ip,1]+dx*dd/dy]
+                print ip, 'Lluis pi:', polygon[ip,:], 'pi+1', polygon[eep,:]
+                print '    yval', yval, 'dx', dx, 'dy', dy
+                print '      ept', ept
+
+    if ipt is None or ept is None:
+        print errormsg
+        print '  ' + fname + ': no cutting for polygon at y=', yval, '!!'
+        irmv = 0
+        rmpolygon = []
+        for ip in range(N):
+            if polygon[ip,0] > yval:
+                rmpolygon.append([gen.fillValueF, gen.fillValueF])
+            else:
+                rmpolygon.append(polygon[ip,:])
+        Npts = len(rmpolygon)
+        cutpolygon = np.array(rmpolygon)
+    else:
+        Npts = icut + (N-ecut) + Nadd
+
+        cutpolygon = np.zeros((Npts,2), dtype=np.float)
+        if type(polygon) == type(gen.mamat):
+            cutpolygon[0:icut+1,:] = polygon[0:icut+1,:].filled(gen.fillValueF)
+        else:
+            cutpolygon[0:icut+1,:] = polygon[0:icut+1,:]
+
+        # cutting line
+        cutline = np.zeros((Nadd,2), dtype=np.float)
+        dx = (ept[1] - ipt[1])/(Nadd-2)
+        dy = (ept[0] - ipt[0])/(Nadd-2)
+        cutline[0,:] = ipt
+        for ip in range(1,Nadd-1):
+            cutline[ip,:] = ipt + np.array([dy*ip,dx*ip])
+        cutline[Nadd-1,:] = ept
+        cutpolygon[icut+1:icut+1+Nadd,:] = cutline
+        cutpolygon[icut+1+Nadd:Npts,:] = polygon[ecut+1:N,:]
+
+    cutpolygon = ma.masked_equal(cutpolygon, gen.fillValueF)
+
+    return Npts, cutpolygon
 
 ####### ###### ##### #### ### ## #