Changeset 2583 in lmdz_wrf

Timestamp:

Jun 2, 2019, 6:23:46 PM (6 years ago)

Author:

lfita

Message:

Adding multiple cuts in `cut_xpolygon'

File:

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

trunk/tools/geometry_tools.py

@@ -622 +622 @@
     ept = None
 
+    icut = []
+    ecut = []
+    ipt = []
+    ept = []
+    Ncuts = 0
     if type(polygon) == type(gen.mamat) and type(polygon.mask) !=                    \
       type(gen.mamat.mask[1]):
@@ -630 +635 @@
                 if eep == N: eep = 0
       
-                if polygon[ip,1] <= xval and polygon[eep,1] >= xval:
-                    icut = ip
+                if val_between(polygon[ip,1], polygon[eep,1], xval):
+                    icut.append(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
+                    ipt.append([polygon[ip,0]+dy*dd/dx, xval])
+
+                if val_between(polygon[ip,1], polygon[eep,1], xval):
+                    ecut.append(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]
+                    ept.append([polygon[ip,0]+dy*dd/dx, xval])
+                    Npts = Npts + 1
     else:
         # Assuming clockwise polygons
@@ -649 +655 @@
             if eep == N: eep = 0
       
-            if polygon[ip,1] <= xval and polygon[eep,1] >= xval:
-                icut = ip
+            if val_between(polygon[ip,1], polygon[eep,1], xval):
+                icut.append(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
+                ipt.append([polygon[ip,0]+dy*dd/dx, xval])
+
+            if val_between(polygon[ip,1], polygon[eep,1], xval):
+                ecut.append(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:
+                ept.append([polygon[ip,0]+dy*dd/dx, xval])
+
+    if ipt is None or ept is None or Ncuts == 0:
         print errormsg
         print '  ' + fname + ': no cutting for polygon at x=', xval, '!!'
-
-    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,:]
+        print '  ' + fname + ': found ', Ncuts, ' Ncuts'
+        print '    yval=', xval, 'cut, ip; ipt ep; ept ________'
+        for ic in range(Ncuts):
+            print '      ', ic, icut[ic], ';', ipt[ic], ecut[ic], ';', ept[ic] 
+
+    # Length of joining lines
+    Nadds = []    
+    if Ncuts > 1:
+        Naddc = (Nadd-Ncuts)/(Ncuts-1)
+        if Naddc < 3:
+            print errormsg
+            print '  ' + fname + ': too few points for jioning lines !!'
+            print '    increase Nadd at least to:', Ncuts*3+Ncuts
+            quit(-1)
+        for ic in range(Ncuts-1):
+            Nadds.append(Naddc)
+
+        Nadds.append(N-Naddc*(Ncuts-1))
     else:
-        cutpolygon[iip:iip+Nadd,:] = cutline[::-1,:]
+        Nadds.append(Nadd)
+
+   # Total points cut polygon
+    Ntotpts = 0
+    Ncpts = []
+    for ic in range(Ncuts):
+        ip = ipt[ic]
+        if ic == 0: 
+            dpts = icut[ic] + ecut[ic] + Nadds[ic]
+        else:
+            dpts = ecut[ic] - icut[ic] + Nadds[ic] - 1
+
+        # Adding end of the polygon in 'left' keeps
+        if keep == 'left' and ic == Ncuts - 1: dpts = dpts + N-ecut[ic]
+        Ncpts.append(dpts)
+        Ntotpts = Ntotpts + dpts
+
+    cutpolygon = np.ones((Ntotpts,2), dtype=np.float)*gen.fillValue
+
+    iip = 0
+    iipc = 0
+    for ic in range(Ncuts):
+        Npts = Ncpts[ic]
+        if keep == 'left':
+            if ic == 0:
+                cutpolygon[0:icut[ic]] = polygon[0:icut[ic],:]
+                iip = icut[ic]
+                iipc = icut[ic]
+            dcpt = Ncpts[ic]-Nadds[ic]
+            cutpolygon[iipc+1:iipc+dcpt,:] = polygon[icut[ic]-1:ecut[ic]+1,:]
+            iipc = iipc + dcpt
+        else:
+            cutpolygon[iipc,:] = ipt[ic]
+            cutpolygon[iipc+1:iipc+ecut[ic]-icut[ic],:]=polygon[icut[ic]+1:ecut[ic],:]
+            iipc = iipc+ecut[ic]-icut[ic]-1
+
+        # cutting line
+        cutline = np.zeros((Nadds[ic],2), dtype=np.float)
+        dx = (ept[ic][1] - ipt[ic][1])/(Nadds[ic]-2)
+        dy = (ept[ic][0] - ipt[ic][0])/(Nadds[ic]-2)
+        cutline[0,:] = ipt[ic]
+        for ip in range(1,Nadds[ic]-1):
+            cutline[ip,:] = ipt[ic] + np.array([dy*ip,dx*ip])
+        cutline[Nadd-1,:] = ept[ic]
+        if keep == 'left':
+            cutpolygon[iipc:iipc+Nadds[ic],:] = cutline
+#            cutpolygon[iip+Nadd:Npts,:] = polygon[ecut+1:N,:]
+        else:
+            cutpolygon[iipc:iipc+Nadds[ic],:] = cutline[::-1,:]
 
     rmpolygon = []