Changeset 2606 in lmdz_wrf for trunk/tools

Timestamp:

Jun 16, 2019, 11:43:31 PM (6 years ago)

Author:

lfita

Message:

Working on the generic version of 'cut_between_xpolygon'

File:

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

trunk/tools/geometry_tools.py

@@ -1204 +1204 @@
     N = polygon.shape[0]
 
+    if xval1 > xval2:
+        print errormsg
+        print '  ' + fname + ': wrong between cut values !!'
+        print '     it is expected xval1 < xval2'
+        print '     values provided xval1: (', xval1, ')> xval2 (', xval2, ')'
+        quit(-1)
+
+    xvals = [xval1, xval2]
+
+    ipt = None
+    ept = None
+
+    cuts = {}
+    if type(polygon) == type(gen.mamat) and type(polygon.mask) !=                    \
+      type(gen.mamat.mask[1]):
+        for ic in range(2):
+            xval = xvals[ic]
+            # There might be more than 1 cut ...
+            icut = []
+            ecut = []
+            ipt = []
+            ept = []
+            Ncuts = 0
+            # Assuming clockwise polygons
+            for ip in range(N-1):
+                if not polygon.mask[ip,0]:
+                    eep = ip + 1
+                    if eep == N: eep = 0
+      
+                    if val_consec_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.append([polygon[ip,0]+dy*dd/dx, xval])
+
+                    if val_consec_between(polygon[eep,1], polygon[ip,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.append([polygon[ip,0]+dy*dd/dx, xval])
+                        Ncuts = Ncuts + 1
+
+            # Looking for repeated
+            newicut = icut + []
+            newecut = ecut + []
+            newipt = ipt + []
+            newept = ept + []
+            newNcuts = Ncuts
+            for icp in range(newNcuts-1):
+                for ic2 in range(icp+1,newNcuts):
+                    if newipt[icp] == newipt[ic2]:
+                        Ncuts = Ncuts-1
+                        icut.pop(ic2)
+                        ecut.pop(ic2)
+                        ipt.pop(ic2)
+                        ept.pop(ic2)
+                        newNcuts = Ncuts + 0
+
+            cuts[ic] = [icut, ecut, ipt, ept, Ncuts]
+    else:
+        for ic in range(2):
+            xval = xvals[ic]
+            # There might be more than 1 cut ...
+            icut = []
+            ecut = []
+            ipt = []
+            ept = []
+            Ncuts = 0
+            # Assuming clockwise polygons
+            for ip in range(N-1):
+                eep = ip + 1
+                if eep == N: eep = 0
+      
+                if val_consec_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.append([polygon[ip,0]+dy*dd/dx, xval])
+
+                if val_consec_between(polygon[eep,1], polygon[ip,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.append([polygon[ip,0]+dy*dd/dx, xval])
+                    Ncuts = Ncuts + 1
+            # Looking for repeated
+            newicut = icut + []
+            newecut = ecut + []
+            newipt = ipt + []
+            newept = ept + []
+            newNcuts = Ncuts
+            for icp in range(newNcuts-1):
+                for ic2 in range(icp+1,newNcuts):
+                    if newipt[icp] == newipt[ic2]:
+                        Ncuts = Ncuts-1
+                        icut.pop(ic2)
+                        ecut.pop(ic2)
+                        ipt.pop(ic2)
+                        ept.pop(ic2)
+                        newNcuts = Ncuts + 0
+
+            cuts[ic] = [icut, ecut, ipt, ept, Ncuts]
+
+    for iic in range(0):
+        cutvs = cuts[iic]
+        icut = cutvs[0]
+        ecut = cutvs[1]
+        ipt = cutvs[2]
+        ept = cutvs[3]
+        Ncuts = cutvs[4]
+        if Ncuts > 1:
+            # 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 + []
+            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])
+
+        cuts[iic] = [icut, ecut, ipt, ept, Ncuts]
+
+    Naddlines = {}
+    for icc in range(2):
+        cutv = cuts[icc]
+        Ncuts = cutv[4]
+        # Length of joining lines
+        Nadds = []
+        if 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(Nadd-Naddc*(Ncuts-1))
+        else:
+            Nadds.append(Nadd)
+
+        Naddlines[icc] = Nadds
+
+    # sides
+    sides = {}
+    for iic in range(2):
+        cutvs = cuts[iic]
+        icut = cutvs[0]
+        ecut = cutvs[1]
+        ipt = cutvs[2]
+        ept = cutvs[3]
+        Ncuts = cutvs[4]
+        Nadds = Naddlines[iic]
+        cutpolygon = []
+        # left side
+        if iic == 0:
+            for ic in range(Ncuts-1):
+                cutpolygon.append(ipt[ic])
+
+                if ic == 0: 
+                    for ip in range(icut[ic],N-1): cutpolygon.append(polygon[ip,:])
+                for ip in range(icut[ic],ecut[ic]): cutpolygon.append(polygon[ip,:])
+                # line
+                print 'Lluis ipt+1', ipt[ic+1], 'ept', ept[ic], 'Nadds:', Nadds[ic]
+                dx = (ipt[ic+1][1] - ept[ic][1])/(Nadds[ic]-1)
+                dy = (ipt[ic+1][0] - ept[ic][0])/(Nadds[ic]-1)
+                for ip in range(1,Nadds[ic]-1):
+                    cutpolygon.append([ept[ic][0]+dy*ip, ept[ic][0]+dy*ip])
+            for ip in range(icut[ic],ecut[ic]): cutpolygon.append(polygon[ip,:])
+        # right side
+        else:
+            for ic in range(Ncuts-1):
+                cutpolygon.append(ipt[ic])
+
+                for ip in range(icut[ic],ecut[ic]): cutpolygon.append(polygon[ip,:])
+                # line
+                dx = (ipt[ic+1][1] - ept[ic][1])/(Nadds[ic]-1)
+                dy = (ipt[ic+1][0] - ept[ic][0])/(Nadds[ic]-1)
+                for ip in range(1,Nadds[ic]-1):
+                    cutpolygon.append([ept[ic,0]+dy*ip, ept[ic,0]+dy*ip])
+            for ip in range(icut[ic],ecut[ic]): cutpolygon.append(polygon[ip,:])
+            for ip in range(ecut[ic],N-1): cutpolygon.append(polygon[ip,:])
+        sides[iic] = cutpolygon
+        
+    # joining sides by e1[Ncuts1-1] --> i2[0]; e2[Ncuts2-1] --> i1[0]
+    cutv1 = cuts[0]
+    Ncuts1 = cutv1[4]
+    ec1 = cutv1[2][Ncuts1-1]
+    ic1 = cutv1[1][0]
+    ept1 = cutv1[3][Ncuts1-1]
+    ipt1 = cutv1[4][0]
+
+    cutv2 = cuts[0]
+    Ncuts2 = cutv2[4]
+    ec2 = cutv2[2][Ncuts1-1]
+    ic2 = cutv2[1][0]
+    ept2 = cutv2[3][Ncuts1-1]
+    ipt2 = cutv2[4][0]
+
+    finalcutpolygon = sides[0]
+    for ip in range(ec1,ic2): finalcutpolygon.append(polygon[ip,:])
+    finalcutpolygon = finalcutpolygon + sides[1]
+    for ip in range(ec2,ic1): finalcutpolygon.append(polygon[ip,:])
+
+
+    finalcutpolygon = np.array(finalcutpolygon)
+    finalcutpolygon = ma.masked_equal(finalcutpolygon, gen.fillValueF)
+
+    Npts = finalcutpolygon.shape[0]
+
+    return Npts, finalcutpolygon
+
+def cut_between_xpolygon_old(polygon, xval1, xval2, Nadd=20):
+    """ Function to cut a polygon between 2 given value of the x-axis
+      polygon: polygon to cut
+      xval1: first value to use to cut the polygon
+      xval2: first value to use to cut the polygon
+      Nadd: additional points to add to draw the line (20, default)
+    """
+    fname = 'cut_betwen_xpolygon'
+
+    N = polygon.shape[0]
+
     ipt = None
     ept = None