Changeset 2556 in lmdz_wrf

Timestamp:

May 25, 2019, 7:04:18 PM (6 years ago)

Author:

lfita

Message:

Adding:

• `cut_between_ypolygon': Function to cut a polygon between 2 given value of the y-axis

File:

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

trunk/tools/geometry_tools.py

@@ -23 +23 @@
 
 ####### Contents:
+# cut_between_[x/y]polygon: Function to cut a polygon between 2 given value of the y-axis
 # cut_[x/y]polygon: Function to cut a polygon from a given value of the [x/y]-axis
 # deg_deci: Function to pass from degrees [deg, minute, sec] to decimal angles [rad]
@@ -465 +466 @@
                     dy = polygon[eep,0] - polygon[ip,0]
                     dd = yval - polygon[ip,0]
-                    print polygon[ip,:], ':', polygon[eep,:]
                     ept = [yval, polygon[ip,1]+dx*dd/dy]
     else:
@@ -514 +514 @@
         cutline[ip,:] = ipt + np.array([dy*ip,dx*ip])
     cutline[Nadd-1,:] = ept
-    cutpolygon[iip:iip+Nadd,:] = cutline
     if keep == 'bottom':
+        cutpolygon[iip:iip+Nadd,:] = cutline
         if type(polygon) == type(gen.mamat):
             cutpolygon[iip+Nadd:Npts,:] = polygon[ecut+1:N,:]
         else:
             cutpolygon[iip+Nadd:Npts,:] = polygon[ecut+1:N,:]
+    else:
+        cutpolygon[iip:iip+Nadd,:] = cutline[::-1,:]
 
     rmpolygon = []
@@ -541 +543 @@
     return Npts, cutpolygon
 
-def cut_xpolygon(polygon, xval, keep='left', Nadd=20):
-    """ Function to cut a polygon from a given value of the x-axis
+def cut_between_ypolygon(polygon, yval1, yval2, Nadd=20):
+    """ Function to cut a polygon between 2 given value of the y-axis
       polygon: polygon to cut
-      xval: value to use to cut the polygon
-      keep: part to keep from the cut ('left', default)
-         'left': left from the cut
-         'right': right from the cut
+      yval1: first value to use to cut the polygon
+      yval2: first value to use to cut the polygon
       Nadd: additional points to add to draw the line (20, default)
     """
-    fname = 'cut_xpolygon'
+    fname = 'cut_betwen_ypolygon'
 
     N = polygon.shape[0]
-    availkeeps = ['left', 'right']
-
-    if not gen.searchInlist(availkeeps, keep):
-        print errormsg
-        print '  ' + fname + ": wring keep '" + keep + "' value !!"
-        print '    available ones:', availkeeps
-        quit(-1)
 
     ipt = None
     ept = None
 
-    if type(polygon) == type(gen.mamat):
-        # Assuming clockwise polygons
-        for ip in range(N):
-            if not polygon.mask[ip,0]:
+    dx = np.zeros((2), dtype=np.float)
+    dy = np.zeros((2), dtype=np.float)
+    icut = np.zeros((2), dtype=int)
+    ecut = np.zeros((2), dtype=int)
+    ipt = np.zeros((2,2), dtype=np.float)
+    ept = np.zeros((2,2), dtype=np.float)
+
+    if yval1 > yval2:
+        print errormsg
+        print '  ' + fname + ': wrong between cut values !!'
+        print '     it is expected yval1 < yval2'
+        print '     values provided yval1: (', yval1, ')> yval2 (', yval2, ')'
+        quit(-1)
+
+    yvals = [yval1, yval2]
+
+    for ic in range(2):
+        yval = yvals[ic]
+        if type(polygon) == type(gen.mamat):
+            # Assuming clockwise polygons
+            for ip in range(N-1):
+                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[ic] = ip
+                        dx[ic] = polygon[eep,1] - polygon[ip,1]
+                        dy[ic] = polygon[eep,0] - polygon[ip,0]
+                        dd = yval - polygon[ip,0]
+                        ipt[ic,:] = [yval, polygon[ip,1]+dx[ic]*dd/dy[ic]]
+
+                    if polygon[ip,0] >= yval and polygon[eep,0] <= yval:
+                        ecut[ic] = ip
+                        dx[ic] = polygon[eep,1] - polygon[ip,1]
+                        dy[ic] = polygon[eep,0] - polygon[ip,0]
+                        dd = yval - polygon[ip,0]
+                        ept[ic,:] = [yval, polygon[ip,1]+dx[ic]*dd/dy[ic]]
+        else:
+            # Assuming clockwise polygons
+            for ip in range(N-1):
                 eep = ip + 1
                 if eep == N: eep = 0
       
-                if polygon[ip,1] <= xval and polygon[eep,1] >= xval:
-                    icut = ip
-                    dx = polygon[eep,1] - polygon[ip,1]
-                    dy = polygon[eep,0] - polygon[ip,0]
-                    dd = xval - polygon[ip,1]
-                    ipt = [ polygon[ip,0]+dy*dd/dx, xval]
-
-                if polygon[ip,1] >= xval and polygon[eep,1] <= xval:
-                    ecut = ip
-                    dx = polygon[eep,1] - polygon[ip,1]
-                    dy = polygon[eep,0] - polygon[ip,0]
-                    dd = xval - polygon[eep,1]
-                    ept = [polygon[eep,0]+dy*dd/dx, xval]
-    else:
-        # Assuming clockwise polygons
-        for ip in range(N):
-            eep = ip + 1
-            if eep == N: eep = 0
-      
-            if polygon[ip,1] <= xval and polygon[eep,1] >= xval:
-                icut = ip
-                dx = polygon[eep,1] - polygon[ip,1]
-                dy = polygon[eep,0] - polygon[ip,0]
-                dd = xval - polygon[ip,1]
-                ipt = [polygon[ip,0]+dy*dd/dx, xval]
-
-            if polygon[ip,1] >= xval and polygon[eep,1] <= xval:
-                ecut = ip
-                dx = polygon[eep,1] - polygon[ip,1]
-                dy = polygon[eep,0] - polygon[ip,0]
-                dd = xval - polygon[ip,1]
-                ept = [polygon[ip,0]+dy*dd/dx, xval]
-
-    if ipt is None or ept is None:
-        print errormsg
-        print '  ' + fname + ': no cutting for polygon at y=', xval, '!!'
-        irmv = 0
-        rmpolygon = []
-        if keep == 'left':
-            for ip in range(N):
-                if polygon[ip,1] > xval:
-                    rmpolygon.append([gen.fillValueF, gen.fillValueF])
-                else:
-                    rmpolygon.append(polygon[ip,:])
-        else:
-            for ip in range(N):
-                if polygon[ip,0] < xval:
-                    rmpolygon.append([gen.fillValueF, gen.fillValueF])
-                else:
-                    rmpolygon.append(polygon[ip,:])
-
-        Npts = len(rmpolygon)
-        cutpolygon = np.array(rmpolygon)
-    else:
-        if keep == 'left':
-            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,:]
-            iip = icut+1
-        else:
-            Npts = ecut - icut + Nadd-1
-
-            cutpolygon = np.zeros((Npts,2), dtype=np.float)
-            if type(polygon) == type(gen.mamat):
-                cutpolygon[0:ecut-icut-1,:] = polygon[icut+1:ecut,:].filled(gen.fillValueF)
-            else:
-                cutpolygon[0:ecut-icut-1,:] = polygon[icut+1:ecut,:]
-            iip = ecut-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[iip:iip+Nadd,:] = cutline
-        if keep == 'left':
-            cutpolygon[iip+Nadd:Npts,:] = polygon[ecut+1:N,:]
+                if polygon[ip,0] <= yval and polygon[eep,0] >= yval:
+                    icut[ic] = ip
+                    dx[ic] = polygon[eep,1] - polygon[ip,1]
+                    dy[ic] = polygon[eep,0] - polygon[ip,0]
+                    dd = yval - polygon[ip,0]
+                    ipt[ic,:] = [yval, polygon[ip,1]+dx[ic]*dd/dy[ic]]
+
+                if polygon[ip,0] >= yval and polygon[eep,0] <= yval:
+                    ecut[ic] = ip
+                    dx[ic] = polygon[eep,1] - polygon[ip,1]
+                    dy[ic] = polygon[eep,0] - polygon[ip,0]
+                    dd = yval - polygon[ip,0]
+                    ept[ic,:] = [yval, polygon[ip,1]+dx[ic]*dd/dy[ic]]
+
+        if ipt is None or ept is None:
+            print errormsg
+            print '  ' + fname + ': no cutting for polygon at y=', yval, '!!'
+
+    Npts = icut[1] - icut[0] + Nadd + ecut[0] - ecut[1]
+    cutpolygon = np.zeros((Npts,2), dtype=np.float)
+    cutpolygon[0,:] = ipt[0,:]
+    cutpolygon[1:icut[1]-icut[0]+1,:] = polygon[icut[0]+1:icut[1]+1,:]
+    iip = icut[1]-icut[0]
+
+    # cutting lines
+    Nadd2 = int(Nadd/2)
+    cutlines = np.zeros((2,Nadd2,2), dtype=np.float)
+
+    for ic in range(2):
+        dx = (ept[ic,1] - ipt[ic,1])/(Nadd2-2)
+        dy = (ept[ic,0] - ipt[ic,0])/(Nadd2-2)
+        cutlines[ic,0,:] = ipt[ic,:]
+        for ip in range(1,Nadd2-1):
+            cutlines[ic,ip,:] = ipt[ic,:] + np.array([dy*ip,dx*ip])
+        cutlines[ic,Nadd2-1,:] = ept[ic,:]
+
+    cutpolygon[iip:iip+Nadd2,:] = cutlines[1,:,:]
+    iip = iip + Nadd2
+    cutpolygon[iip:iip+(ecut[0]-ecut[1]),:] = polygon[ecut[1]+1:ecut[0]+1,:]
+    iip = iip + ecut[0]-ecut[1]
+    cutpolygon[iip:iip+Nadd2,:] = cutlines[0,::-1,:]
 
     cutpolygon = ma.masked_equal(cutpolygon, gen.fillValueF)
@@ -1682 +1669 @@
     Height = np.max(buoy1v[:,0])
 
-    print 'Lluis Height', Height
-
     Ncut, halfdown = cut_ypolygon(buoy1v, yval=Height/2., keep='bottom')
     Ncut, halfup = cut_ypolygon(buoy1v, yval=Height/2., keep='above')