Changeset 2592 in lmdz_wrf

Timestamp:

Jun 4, 2019, 4:38:08 AM (6 years ago)

Author:

lfita

Message:

First working version of `cut_xpolygon' for multiple cuts

File:

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

trunk/tools/geometry_tools.py

@@ -687 +687 @@
         print '  ' + fname + ': no cutting for polygon at x=', xval, '!!'
     else:
-        print '  ' + fname + ': found ', Ncuts, ' Ncuts'
-        print '    xval=', xval, 'cut, ip; ipt ep; ept ________'
-        for ic in range(Ncuts):
-            print '      ', ic, icut[ic], ';', ipt[ic], ecut[ic], ';', ept[ic] 
+        ##print '  ' + fname + ': found ', Ncuts, ' Ncuts'
+        if Ncuts > 1 and keep == 'left':
+            # Re-shifting cuts. 1st icut --> last ecut; 1st ecut as 1st icut; 
+            #    2nd icut --> last-1 ecut, ....
+            newicut = icut + []
+            newecut = ecut + []
+            newipt = ipt + []
+            newept = ept + []
+            for ic in range(Ncuts-1):
+                ecut[ic] = newecut[Ncuts-ic-1]
+                ept[ic] = newept[Ncuts-ic-1]
+                icut[ic+1] = newecut[ic]
+                ipt[ic+1] = newept[ic]
+
+            ecut[Ncuts-1] = newicut[Ncuts-1]
+            ept[Ncuts-1] = newipt[Ncuts-1]
+
+        ##print '    xval=', 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)
+        Naddc = (Nadd-Ncuts)/(Ncuts)
         if Naddc < 3:
             print errormsg
@@ -704 +720 @@
             Nadds.append(Naddc)
 
-        Nadds.append(N-Naddc*(Ncuts-1))
+        Nadds.append(Nadd-Naddc*(Ncuts-1))
     else:
         Nadds.append(Nadd)
@@ -712 +728 @@
     Ncpts = []
     for ic in range(Ncuts):
-        ip = ipt[ic]
-        if ic == 0: 
-            dpts = ecut[ic] - icut[ic] + Nadds[ic] + icut[ic]
+        if keep == 'left':
+            if ic == 0: 
+                dpts = icut[ic] + Nadds[ic] + (N - ecut[ic])
+            else:
+                dpts = ecut[ic] - icut[ic] + Nadds[ic] - 1
+
+            # Adding end of the polygon in 'left' keeps
+            if ic == Ncuts - 1: dpts = dpts + N-ecut[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+Ncuts-1,2), dtype=np.float)*gen.fillValue
-
-    iip = 0
+        Ncpts.append(ecut[ic] - icut[ic])
+
+    cutpolygon = np.ones((Ntotpts+Ncuts,2), dtype=np.float)*gen.fillValue
+
     iipc = 0
     for ic in range(Ncuts):
-        Npts = Ncpts[ic]
+        dcpt = Ncpts[ic]
         if keep == 'left':
             if ic == 0:
                 cutpolygon[0:icut[ic],:] = polygon[0:icut[ic],:]
-                iip = icut[ic]
                 iipc = icut[ic]
-                dcpt = ecut[ic]-icut[ic]
-                print ic, ':', 0, icut[ic], 'iipc:', iipc, 'dcpt', dcpt
             else:
-                dcpt = Ncpts[ic]-Nadds[ic]
-            print ic, 'dcpt:', dcpt, '<Ncpts>', Ncpts[ic], '<--', icut[ic]-1,ecut[ic]+1, '=',  iipc+1, iipc+dcpt+1
-            cutpolygon[iipc+1:iipc+dcpt+1,:] = polygon[icut[ic]-1:ecut[ic]+1,:]
-            iipc = iipc + dcpt
+                cutpolygon[iipc:iipc+dcpt-1,:] = polygon[icut[ic]+1:ecut[ic],:]
+                iipc = iipc + dcpt -1
         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
+            cutpolygon[iipc:iipc+dcpt-1,:]=polygon[icut[ic]+1:ecut[ic],:]
+            iipc = iipc+dcpt-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)
+        dx = (ept[ic][1] - ipt[ic][1])/(Nadds[ic]-1)
+        dy = (ept[ic][0] - ipt[ic][0])/(Nadds[ic]-1)
         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]
+        cutline[Nadds[ic]-1,:] = ept[ic]
         if keep == 'left':
-            cutpolygon[iipc:iipc+Nadds[ic],:] = cutline
-#            cutpolygon[iip+Nadd:Npts,:] = polygon[ecut+1:N,:]
+            if ic == 0: cutpolygon[iipc:iipc+Nadds[ic],:] = cutline
+            else: cutpolygon[iipc:iipc+Nadds[ic],:] = cutline[::-1]
+            iipc = iipc+Nadds[ic]
+            if ic == 0:
+                cutpolygon[iipc:iipc+N-ecut[ic]-1,:] = polygon[ecut[ic]+1:N,:]
+                iipc = iipc + N-ecut[ic]-1
+                cutpolygon[iipc,:] = polygon[0,:]
+                iipc = iipc + 1
         else:
             cutpolygon[iipc:iipc+Nadds[ic],:] = cutline[::-1,:]
+            iipc = iipc+Nadds[ic]
+        iipc = iipc + 1
 
     rmpolygon = []