Changeset 2627 in lmdz_wrf

Timestamp:

Jun 23, 2019, 8:33:06 PM (5 years ago)

Author:

lfita

Message:

First working version of `pile_polygons'

File:

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

trunk/tools/geometry_tools.py

@@ -1595 +1595 @@
         linti = list(inti)
         for ip in range(iNpts):
+            iip1 = ip+1
+            if ip == iNpts-1: iip1 = 0
+            print ip, ipoly[ip,:], ':', ipoly[iip1,:]
             Nc = linti.count(ip)
-            dists = np.zeros((Nc), dtype=np.float)
+            ldists = []
+            ddists = {}
             if Nc > 0:
                 iic = gen.multi_index_vec(inti,ip)
-                # Sortering from distance respect the vertex ip
+                mindist = 1000000.
+                # Sorting from distance respect the vertex ip
                 for ic in range(Nc):
-                    dists[ic] = dist_points(ipoly[ip,:], pts[inti[ic],:])
+                    ddists[iic[ic]] = dist_points(ipoly[ip,:], pts[iic[ic],:])
+                    ldists.append(ddists[iic[ic]])
+                    print '  ', ic, ';', iic[ic], '=', pts[iic[ic],:], ddists[iic[ic]]
+                ldists.sort()
+                print '    ldists', ldists
+                newipoly.append(ipoly[ip,:])
+                for ic in range(Nc):
+                    iic = gen.dictionary_key(ddists,ldists[ic])
+                    newipoly.append(pts[iic,:])
+                    print '  ', ic, '|', iic, ';', pts[iic,:]
+                    if ic < Nc-1: newipoly.append([gen.fillValueF, gen.fillValueF])
+                newipoly.append(ipoly[iip1,:])
+
             else:
                 newipoly.append(ipoly[ip,:])
-
-
-
-        Npi=0
-        for ip in range(Nint):
-            # Determine side of the point of the polygon to cut
-            #   Here are assumed that all polygons are clockwise created, thus any
-            #   cutting point will lay to the left side of the face which cuts
-            #   but we need to deterime if it is the starting point or the ending point
-            #   the one which lays to the left of the over pile polygon
-            iip = inti[ip]
-            iip1 = iip + 1
-            if iip1 > iNpts-1: iip1=0
-            pip = intp[ip]
-            pip1 = pip + 1
-            if pip1 > Npts-1: pip1=0
-
-            print ip, ': Lluis iip, iip1', iip, iip1, 'pip, pip1', pip, pip1
-
-            iface = ipoly[iip,:]-poly[pip,:]
-            pface = poly[pip1,:]-poly[pip,:]
-            ifangle=angle_vectors2D(pface,iface)
-
-            iface = ipoly[iip1,:]-poly[pip,:]
-            pface = poly[pip1,:]-poly[pip,:]
-            efangle=angle_vectors2D(pface,iface)
-
-            if ifangle < 0.: 
-                # Initial point of the cuted face to the left of the overlying polygon
-                for iv in range(Npi,Npi+iip): newipoly.append(ipoly[iv,:])
-                Npi = Npi + iip
-                newipoly.append(pts[ip])
-                newipoly.append([gen.fillValueF, gen.fillValueF])
-
-            elif efangle < 0.:
-                # Ending point of the cuted face to the left of the overlying polygon
-                newipoly.append([gen.fillValueF, gen.fillValueF])
-                newipoly.append(pts[ip])
-                newipoly.append(ipoly[iip1,:])
-                Npi = iip1
-            else:
-                print errormsg
-                print '  ' + fname + ': No left point found !!'
-                print '    crossing point number', ip, 'below polygon vertex:',      \
-                  inti[ip], 'above', intp[ip]
-                print '    for pile face', poly[pip,:], '->', poly[pip1,:],          \
-                  'and below face', ipoly[iip,:], '->', ipoly[iip1,:]
     
         pilepolygons[polyns[Npolys-1]] = newipoly