Changeset 2560 in lmdz_wrf

Timestamp:

May 28, 2019, 2:21:06 AM (6 years ago)

Author:

lfita

Message:

Starting to add `cut_xpolygon'

File:

• trunk/tools/geometry_tools.py (modified) (1 diff)

trunk/tools/geometry_tools.py

@@ -533 +533 @@
         for ip in range(Npts):
             if cutpolygon[ip,0] < yval:
+                rmpolygon.append([gen.fillValueF, gen.fillValueF])
+            else:
+                rmpolygon.append(cutpolygon[ip,:])
+    Npts = len(rmpolygon)
+    cutpolygon = np.array(rmpolygon)
+
+    cutpolygon = ma.masked_equal(cutpolygon, gen.fillValueF)
+
+    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
+      polygon: polygon to cut
+      yval: value to use to cut the polygon
+      keep: part to keep from the value ('left', default)
+         'left': left of the value
+         'right': right of the value
+      Nadd: additional points to add to draw the line (20, default)
+    """
+    fname = 'cut_xpolygon'
+
+    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-1):
+            if not polygon.mask[ip,0]:
+                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]
+    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[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=', yval, '!!'
+
+    print 'Lluis icut', icut, 'ipt:', ipt, 'ecut', ecut, 'ept:', ept, 'N', N
+    print 'Lluis 00', polygon[0:7,:]
+    print 'Lluis NN', polygon[ecut+1:N,:]
+
+    if keep == 'left':
+        Npts = icut + (N-ecut) + Nadd + 1
+        cutpolygon = np.zeros((Npts,2), dtype=np.float)
+        cutpolygon[0,:] = ept
+        cutpolygon[1:N-ecut,:] = polygon[ecut+1:N,:]
+        iip = N-ecut
+        cutpolygon[iip:iip+icut+1,:] = polygon[0:icut+1,:]
+        iip = iip + icut+1
+    else:
+        Npts = ecut - icut + Nadd-1
+        cutpolygon = np.zeros((Npts,2), dtype=np.float)
+        cutpolygon[0,:] = ipt
+        cutpolygon[1:ecut-icut,:] = 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
+    if keep == 'left':
+        cutpolygon[iip:iip+Nadd,:] = cutline
+#        cutpolygon[iip+Nadd:Npts,:] = polygon[ecut+1:N,:]
+    else:
+        cutpolygon[iip:iip+Nadd,:] = cutline[::-1,:]
+
+    rmpolygon = []
+    if keep == 'left':
+        for ip in range(Npts):
+            if cutpolygon[ip,1] > xval:
+                rmpolygon.append([gen.fillValueF, gen.fillValueF])
+            else:
+                rmpolygon.append(cutpolygon[ip,:])
+    else:
+        for ip in range(Npts):
+            if cutpolygon[ip,1] < xval:
                 rmpolygon.append([gen.fillValueF, gen.fillValueF])
             else: