Changeset 2597 in lmdz_wrf

Timestamp:

Jun 7, 2019, 6:44:13 AM (6 years ago)

Author:

lfita

Message:

Seem to completely work 'cut_[y/x]polygon' for multiple cuts

File:

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

trunk/tools/geometry_tools.py

@@ -513 +513 @@
                     ipt.append([yval, polygon[ip,1]+dx*dd/dy])
 
-                if val_consec_between(polygon[ip,0], polygon[eep,0], yval):
+                if val_consec_between(polygon[eep,0], polygon[ip,0], yval):
                     ecut.append(ip)
                     dx = polygon[eep,1] - polygon[ip,1]
@@ -533 +533 @@
                 ipt.append([yval, polygon[ip,1]+dx*dd/dy])
 
-            if val_consec_between(polygon[ip,0], polygon[eep,0], yval):
+            if val_consec_between(polygon[eep,0], polygon[ip,0], yval):
                 ecut.append(ip)
                 dx = polygon[eep,1] - polygon[ip,1]
@@ -541 +541 @@
                 Ncuts = Ncuts + 1
 
+    # Looking for repeated
+    newicut = icut + []
+    newecut = ecut + []
+    newipt = ipt + []
+    newept = ept + []
+    newNcuts = Ncuts
+    for ic in range(newNcuts-1):
+        for ic2 in range(ic+1,newNcuts):
+            if newipt[ic] == newipt[ic2]:
+                Ncuts = Ncuts-1
+                icut.pop(ic2)
+                ecut.pop(ic2)
+                ipt.pop(ic2)
+                ept.pop(ic2)
+                newNcuts = Ncuts + 0
+
     if ipt is None or ept is None or Ncuts == 0:
         print errormsg
@@ -546 +562 @@
     else:
         print '  ' + fname + ': found ', Ncuts, ' Ncuts'
-        print '    yval=', yval, 'cut, ip; ipt ep; ept ________'
-        for ic in range(Ncuts):
-            print '      ', ic, icut[ic], ';', ipt[ic], ecut[ic], ';', ept[ic] 
-
+        if Ncuts > 1 and keep == 'below':
+            # Re-shifting cuts by closest distance.
+            xis = []
+            xes = []
+            for ic in range(Ncuts):
+                xp = ipt[ic]
+                xis.append(xp[1])
+                xp = ept[ic]
+                xes.append(xp[1])
+            xs = xis + xes
+            xs.sort()
+            newicut = icut + []
+            newecut = ecut + []
+            newipt = ipt + []
+            newept = ept + []
+            icut = []
+            ecut = []
+            ipt = []
+            ept = []
+            for xv in xs:
+                ic = xis.count(xv)
+                if ic != 0: 
+                    icc = xis.index(xv)
+                    if len(icut) > len(ecut): 
+                        ecut.append(newicut[icc])
+                        ept.append(newipt[icc])
+                    else: 
+                        icut.append(newicut[icc])
+                        ipt.append(newipt[icc])
+                else:
+                    icc = xes.index(xv)
+                    if len(icut) > len(ecut): 
+                        ecut.append(newecut[icc])
+                        ept.append(newept[icc])
+                    else: 
+                        icut.append(newecut[icc])
+                        ipt.append(newept[icc])
+
+#            # 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 '    yval=', yval, '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-1)
+        Naddc = (Nadd-Ncuts)/(Ncuts)
+        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))
+        Nadds.append(Nadd-Naddc*(Ncuts-1))
     else:
         Nadds.append(Nadd)
 
-    iip = 0
+    # Total points cut polygon
+    Ntotpts = 0
+    Ncpts = []
+    for ic in range(Ncuts):
+        if keep == 'below':
+            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
+
+        Ntotpts = Ntotpts + dpts
+        Ncpts.append(ecut[ic] - icut[ic])
+
+    cutpolygon = np.ones((Ntotpts+Ncuts,2), dtype=np.float)*gen.fillValue
+
     iipc = 0
     for ic in range(Ncuts):
+        dcpt = Ncpts[ic]
         if keep == 'below':
-            Npts = icut[ic] + (N-ecut[ic]) + Nadds[ic]
-            cutpolygon = np.zeros((Npts,2), dtype=np.float)
-            cutpolygon[iipc:iipc+icut[ic]+1,:] = polygon[iip:iip+icut[ic]+1,:]
-            iip = iip+icut[ic]+1
-            iipc = iipc+icut[ic]+1
+            if ic == 0:
+                cutpolygon[0:icut[ic],:] = polygon[0:icut[ic],:]
+                iipc = icut[ic]
+            else:
+                cutpolygon[iipc:iipc+dcpt-1,:] = polygon[icut[ic]+1:ecut[ic],:]
+                iipc = iipc + dcpt -1
         else:
-            Npts = ecut[ec] - icut[ic] + Nadds[ic]-1
-            cutpolygon = np.zeros((Npts,2), dtype=np.float)
             cutpolygon[iipc,:] = ipt[ic]
-            cutpolygon[iipc+1:ecut[ic]-icut[ic],:] = polygon[icut[ic]+1:ecut[ic],:]
-            iip = ecut[ic]-icut[ic]-1
-            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):
@@ -586 +679 @@
         cutline[Nadds[ic]-1,:] = ept[ic]
         if keep == 'below':
-            cutpolygon[iip:iip+Nadds[ic],:] = cutline
-            cutpolygon[iip+Nadds[ic]:Npts,:] = polygon[ecut[ic]+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[iip:iip+Nadds[ic],:] = cutline[::-1,:]
+            cutpolygon[iipc:iipc+Nadds[ic],:] = cutline[::-1,:]
+            iipc = iipc+Nadds[ic]
+        iipc = iipc + 1
 
     rmpolygon = []
@@ -635 +736 @@
     ept = None
 
+    # There might be more than 1 cut ...
     icut = []
     ecut = []
@@ -683 +785 @@
                 Ncuts = Ncuts + 1
 
+    # Looking for repeated
+    newicut = icut + []
+    newecut = ecut + []
+    newipt = ipt + []
+    newept = ept + []
+    newNcuts = Ncuts
+    for ic in range(newNcuts-1):
+        for ic2 in range(ic+1,newNcuts):
+            if newipt[ic] == newipt[ic2]:
+                Ncuts = Ncuts-1
+                icut.pop(ic2)
+                ecut.pop(ic2)
+                ipt.pop(ic2)
+                ept.pop(ic2)
+                newNcuts = Ncuts + 0
+
     if ipt is None or ept is None or Ncuts == 0:
         print errormsg
@@ -689 +807 @@
         ##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, ....
+            # Re-shifting cuts by closest heigth.
+            yis = []
+            yes = []
+            for ic in range(Ncuts):
+                yp = ipt[ic]
+                yis.append(yp[0])
+                yp = ept[ic]
+                yes.append(yp[0])
+            ys = yis + yes
+            ys.sort()
             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]
-
+            icut = []
+            ecut = []
+            ipt = []
+            ept = []
+            for yv in ys:
+                ic = yis.count(yv)
+                if ic != 0: 
+                    icc = yis.index(yv)
+                    if len(icut) > len(ecut): 
+                        ecut.append(newicut[icc])
+                        ept.append(newipt[icc])
+                    else: 
+                        icut.append(newicut[icc])
+                        ipt.append(newipt[icc])
+                else:
+                    icc = yes.index(yv)
+                    if len(icut) > len(ecut): 
+                        ecut.append(newecut[icc])
+                        ept.append(newept[icc])
+                    else: 
+                        icut.append(newecut[icc])
+                        ipt.append(newept[icc])
         ##print '    xval=', xval, 'cut, ip; ipt ep; ept ________'
         ##for ic in range(Ncuts):