Changeset 2605 in lmdz_wrf for trunk/tools

Timestamp:

Jun 16, 2019, 10:37:06 PM (6 years ago)

Author:

lfita

Message:

Working version of `cut_between_ypolygon'
Removing '--' values when 'paint_filled'

File:

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

trunk/tools/geometry_tools.py

@@ -952 +952 @@
     N = polygon.shape[0]
 
-    ipt = None
-    ept = None
-
-    dx = np.zeros((2), dtype=np.float)
-    dy = np.zeros((2), dtype=np.float)
-    icut = np.zeros((2), dtype=int)
-    ecut = np.zeros((2), dtype=int)
-    ipt = np.zeros((2,2), dtype=np.float)
-    ept = np.zeros((2,2), dtype=np.float)
-
     if yval1 > yval2:
         print errormsg
@@ -971 +961 @@
     yvals = [yval1, yval2]
 
-    for ic in range(2):
-        yval = yvals[ic]
-        if type(polygon) == type(gen.mamat):
+    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):
+            yval = yvals[ic]
+            # There might be more than 1 cut ...
+            icut = []
+            ecut = []
+            ipt = []
+            ept = []
+            Ncuts = 0
             # Assuming clockwise polygons
             for ip in range(N-1):
@@ -980 +981 @@
                     if eep == N: eep = 0
       
-                    if polygon[ip,0] <= yval and polygon[eep,0] >= yval:
-                        icut[ic] = ip
-                        dx[ic] = polygon[eep,1] - polygon[ip,1]
-                        dy[ic] = polygon[eep,0] - polygon[ip,0]
+                    if val_consec_between(polygon[ip,0], polygon[eep,0], yval):
+                        icut.append(ip)
+                        dx = polygon[eep,1] - polygon[ip,1]
+                        dy = polygon[eep,0] - polygon[ip,0]
                         dd = yval - polygon[ip,0]
-                        ipt[ic,:] = [yval, polygon[ip,1]+dx[ic]*dd/dy[ic]]
-
-                    if polygon[ip,0] >= yval and polygon[eep,0] <= yval:
-                        ecut[ic] = ip
-                        dx[ic] = polygon[eep,1] - polygon[ip,1]
-                        dy[ic] = polygon[eep,0] - polygon[ip,0]
+                        ipt.append([yval, polygon[ip,1]+dx*dd/dy])
+
+                    if val_consec_between(polygon[eep,0], polygon[ip,0], yval):
+                        ecut.append(ip)
+                        dx = polygon[eep,1] - polygon[ip,1]
+                        dy = polygon[eep,0] - polygon[ip,0]
                         dd = yval - polygon[ip,0]
-                        ept[ic,:] = [yval, polygon[ip,1]+dx[ic]*dd/dy[ic]]
-        else:
+                        ept.append([yval, polygon[ip,1]+dx*dd/dy])
+                        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):
+            yval = yvals[ic]
+            # There might be more than 1 cut ...
+            icut = []
+            ecut = []
+            ipt = []
+            ept = []
+            Ncuts = 0
             # Assuming clockwise polygons
             for ip in range(N-1):
@@ -999 +1027 @@
                 if eep == N: eep = 0
       
-                if polygon[ip,0] <= yval and polygon[eep,0] >= yval:
-                    icut[ic] = ip
-                    dx[ic] = polygon[eep,1] - polygon[ip,1]
-                    dy[ic] = polygon[eep,0] - polygon[ip,0]
+                if val_consec_between(polygon[ip,0], polygon[eep,0], yval):
+                    icut.append(ip)
+                    dx = polygon[eep,1] - polygon[ip,1]
+                    dy = polygon[eep,0] - polygon[ip,0]
                     dd = yval - polygon[ip,0]
-                    ipt[ic,:] = [yval, polygon[ip,1]+dx[ic]*dd/dy[ic]]
-
-                if polygon[ip,0] >= yval and polygon[eep,0] <= yval:
-                    ecut[ic] = ip
-                    dx[ic] = polygon[eep,1] - polygon[ip,1]
-                    dy[ic] = polygon[eep,0] - polygon[ip,0]
+                    ipt.append([yval, polygon[ip,1]+dx*dd/dy])
+
+                if val_consec_between(polygon[eep,0], polygon[ip,0], yval):
+                    ecut.append(ip)
+                    dx = polygon[eep,1] - polygon[ip,1]
+                    dy = polygon[eep,0] - polygon[ip,0]
                     dd = yval - polygon[ip,0]
-                    ept[ic,:] = [yval, polygon[ip,1]+dx[ic]*dd/dy[ic]]
-
-        if ipt is None or ept is None:
-            print errormsg
-            print '  ' + fname + ': no cutting for polygon at y=', yval, '!!'
-
-    Npts = icut[1] - icut[0] + Nadd + ecut[0] - ecut[1]
-    cutpolygon = np.zeros((Npts,2), dtype=np.float)
-    cutpolygon[0,:] = ipt[0,:]
-    cutpolygon[1:icut[1]-icut[0]+1,:] = polygon[icut[0]+1:icut[1]+1,:]
-    iip = icut[1]-icut[0]
-
-    # cutting lines
-    Nadd2 = int(Nadd/2)
-    cutlines = np.zeros((2,Nadd2,2), dtype=np.float)
-
-    for ic in range(2):
-        dx = (ept[ic,1] - ipt[ic,1])/(Nadd2-2)
-        dy = (ept[ic,0] - ipt[ic,0])/(Nadd2-2)
-        cutlines[ic,0,:] = ipt[ic,:]
+                    ept.append([yval, polygon[ip,1]+dx*dd/dy])
+                    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]
+
+    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)
+
+        # Total points cut polygon
+        Ntotpts = 0
+        Ncpts = []
+        for ic in range(Ncuts):
+            dpts = ecut[ic] - icut[ic] + Nadds[ic] - 1
+
+            Ntotpts = Ntotpts + dpts
+            Ncpts.append(ecut[ic] - icut[ic])
+
+        Naddlines[icc] = [Nadds, Ntotpts, Ncpts]
+
+    cutv1 = cuts[0]
+    addv1 = Naddlines[0]
+    Nadds1 = addv1[0]
+    Ncuts1 = cutv1[4]
+
+    cutv2 = cuts[1]
+    addv2 = Naddlines[1]
+    Nadds2 = addv2[0]
+    Ncuts2 = cutv2[4]
+
+    if Ncuts1 != Ncuts2: 
+        print errormsg
+        print '  ' + fname + ": different number of cuts !!"
+        print '    yval1:', yval1, 'Ncuts=', Ncuts1
+        print '    yval2:', yval2, 'Ncuts=', Ncuts2
+        print '      I am not prepare to deal with it'
+        quit(-1)
+    #else:
+    #    print '  ' + fname + ' _______'
+    #    print '    yval1:', yval1, 'Ncuts=', Ncuts1
+    #    print '    yval2:', yval2, 'Ncuts=', Ncuts2
+
+    icut1 = cutv1[0]
+    ecut1 = cutv1[1]
+    ipt1 = cutv1[2]
+    ept1 = cutv1[3]
+    icut2 = cutv2[0]
+    ecut2 = cutv2[1]
+    ipt2 = cutv2[2]
+    ept2 = cutv2[3]
+
+    # Looking for pairs of cuts. Grouping for smallest x distance between initial    
+    #   points of each cut
+    cutpolygons = []
+    for iic1 in range(Ncuts1):
+        iip = 0
+        cutpolygon = []
+        ic1 = icut1[iic1]
+        ec1 = ecut1[iic1]
+        ip1 = ipt1[iic1]
+        ep1 = ept1[iic1]
+
+        ipx2s = []
+        for ip in range(Ncuts2): 
+            ip2 = ipt2[ip]
+            ipx2s.append(ip2[1])
+        dxps = ipx2s - ip1[1]
+        dxps = np.where(dxps < 0., gen.fillValueF, dxps)
+        ndxps = np.min(dxps)
+        iic12 = gen.index_vec(dxps,ndxps)
+
+        ic2 = icut2[iic12]
+        ec2 = ecut2[iic12]
+        ip2 = ipt2[iic12]
+        ep2 = ept2[iic12]
+
+        #print 'Lluis iic1', iic1, 'ic1', ic1, 'ec1', ec1, 'ipt1', ip1, 'ept1', ep1, 'Nadds1', Nadds1
+        #print '    iic12', iic12, 'ic2', ic2, 'ec2', ec2, 'ipt2', ip2, 'ept2', ep2, 'Nadds2', Nadds2
+
+        cutpolygon.append(ip1)
+        for ip in range(ic1+1,ic2-1):
+            cutpolygon.append(polygon[ip,:])
+        iip = ic2-ic1
+        # cutting line 1
+        Nadd2 = Nadds1[iic1]
+        cutlines = np.zeros((Nadd2,2), dtype=np.float)
+        dx = (ep2[1] - ip2[1])/(Nadd2-2)
+        dy = (ep2[0] - ip2[0])/(Nadd2-2)
+        cutlines[0,:] = ip2
         for ip in range(1,Nadd2-1):
-            cutlines[ic,ip,:] = ipt[ic,:] + np.array([dy*ip,dx*ip])
-        cutlines[ic,Nadd2-1,:] = ept[ic,:]
-
-    cutpolygon[iip:iip+Nadd2,:] = cutlines[1,:,:]
-    iip = iip + Nadd2
-    cutpolygon[iip:iip+(ecut[0]-ecut[1]),:] = polygon[ecut[1]+1:ecut[0]+1,:]
-    iip = iip + ecut[0]-ecut[1]
-    cutpolygon[iip:iip+Nadd2,:] = cutlines[0,::-1,:]
-
-    cutpolygon = ma.masked_equal(cutpolygon, gen.fillValueF)
-
-    return Npts, cutpolygon
+            cutlines[ip,:] = ip2 + np.array([dy*ip,dx*ip])
+        cutlines[Nadd2-1,:] = ep2
+        for ip in range(Nadd2): cutpolygon.append(cutlines[ip,:])
+        iip = iip + Nadd2
+
+        for ip in range(ec2,ec1):
+            cutpolygon.append(polygon[ip,:])
+        iip = iip + ec1-ec2
+        # cutting line 2
+        Nadd2 = Nadds2[iic12]
+        cutlines = np.zeros((Nadd2,2), dtype=np.float)
+        dx = (ep1[1] - ip1[1])/(Nadd2-2)
+        dy = (ep1[0] - ip1[0])/(Nadd2-2)
+        cutlines[0,:] = ip1
+        for ip in range(1,Nadd2-1):
+            cutlines[ip,:] = ip1 + np.array([dy*ip,dx*ip])
+        cutlines[Nadd2-1,:] = ep1
+        for ip in range(Nadd2-1,0,-1):
+            cutpolygon.append(cutlines[ip,:])
+
+        cutpolygon.append(ip1)
+
+        cutpolygon.append([gen.fillValueF,gen.fillValueF])
+        if len(cutpolygons) == 0: cutpolygons = cutpolygon
+        else: cutpolygons = cutpolygons + cutpolygon 
+
+    cutpolygons = np.array(cutpolygons)
+    cutpolygons = ma.masked_equal(cutpolygons, gen.fillValueF)
+
+    Npts = cutpolygons.shape[0]
+
+    return Npts, cutpolygons
 
 def cut_between_xpolygon(polygon, xval1, xval2, Nadd=20):
@@ -2389 +2537 @@
     buoy = ma.masked_equal(buoy, gen.fillValueF)
 
-    buoysecs = ['buoy', 'sign1', 'sign2', 'halfk', 'halfy']
+    buoysecs = ['buoy', 'sign1', 'sign2', 'half1', 'half2']
     buoydic = {'buoy': [buoy[0:N2,:],'-','k',1.5],                                   \
       'sign1': [buoy[N2+1:N2+N32+1,:],'-','k',1.5],                                  \
@@ -3185 +3333 @@
         secvals=objdic[secn]
         pvals = secvals[0]
-        plt.fill(pvals[:,1], pvals[:,0], color=secvals[2])
+        fillsec = []
+        for ip in range(pvals.shape[0]):
+            if type(pvals[ip,0]) != type(gen.mamat[1]) and type(pvals[ip,1]) != type(gen.mamat[1]): 
+                fillsec.append(pvals[ip,:])
+        #plt.fill(pvals[:,1], pvals[:,0], color=secvals[2])
+        fillsec = np.array(fillsec)
+        plt.fill(fillsec[:,1], fillsec[:,0], color=secvals[2])
 
     return