Changeset 2628 in lmdz_wrf

Timestamp:

Jun 23, 2019, 9:08:30 PM (5 years ago)

Author:

lfita

Message:

Working version, but not for more than 2 polygons

File:

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

trunk/tools/geometry_tools.py

@@ -547 +547 @@
     prevpt = polygon[0,:]
     newpolygon.append(prevpt)
-    for ip in range(1,polygon.shape[0]-1):
+    for ip in range(1,polygon.shape[0]):
         if polygon[ip,0] != prevpt[0] or polygon[ip,1] != prevpt[1]:
             prevpt = polygon[ip,:]
@@ -1567 +1567 @@
     >>> plgs = {'sqra': polya, 'sqrb': polyb}
     >>> pile_polygons(pns, plgs)
+    #  sqrb : 
+    [[-0.75 -0.5]
+     [-0.5 -0.5]
+     [-- --]
+     [0.5 -0.5]
+     [0.75 -0.5]
+     [0.75 -0.5]
+     [0.75 0.5]
+     [0.5 0.5]
+     [-- --]
+     [-0.5 0.5]
+     [-0.75 0.5]
+     [-0.75 0.5]
+     [-0.75 -0.5]]
+    #  sqra : [[-0.5  -0.75]
+ [ 0.5  -0.75]
+ [ 0.5   0.75]
+ [-0.5   0.75]
+ [-0.5  -0.75]]
+
+
     [[-0.75 -0.5 ]
      [ 0.75 -0.5 ]
@@ -1574 +1595 @@
     fname = 'pile_polygons'
     pilepolygons = dict(polygons)
-
     Npolys = len(polyns)
-    ipoly = polygons[polyns[Npolys-1]]
-    iNpts = ipoly.shape[0]
-
-    pilesortpolys = polyns[::-1]
-    pilesortpolys.remove(polyns[Npolys-1])
-    newipoly = []
-    for polyn in pilesortpolys:
-        poly = polygons[polyn]
+
+    for ipolyp in range(Npolys-2,-1,-1):
+        polyn = polyns[ipolyp]
+        poly = pilepolygons[polyn]
         Npts = poly.shape[0]
-        Nint, inti, intp, pts = fsci.module_scientific.crossingpoints_polys(         \
-          nvertexa=iNpts, nvertexb=Npts, nvertexab=iNpts*Npts, polya=ipoly,          \
-          polyb=poly)
-        # We're in C-mode !
-        inti = inti-1
-        intp = intp-1
-
-        # Re-constructing below polygon looking for respective crossings
-        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)
-            ldists = []
-            ddists = {}
-            if Nc > 0:
-                iic = gen.multi_index_vec(inti,ip)
-                mindist = 1000000.
-                # Sorting from distance respect the vertex ip
-                for ic in range(Nc):
-                    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,:])
+        for ipolyi in range(ipolyp+1,Npolys,1):
+            ipolyn = polyns[ipolyi]
+            #print '  Lluis ' + polyn + ' above ' + ipolyn
+            ipoly = pilepolygons[ipolyn]
+            iNpts = ipoly.shape[0]
+            newipoly = []
+
+            Nint, inti, intp, pts = fsci.module_scientific.crossingpoints_polys(     \
+              nvertexa=iNpts, nvertexb=Npts, nvertexab=iNpts*Npts, polya=ipoly,      \
+              polyb=poly)
+            # We're in C-mode !
+            inti = inti-1
+            intp = intp-1
+
+            # Re-constructing below polygon looking for respective crossings
+            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)
+                ldists = []
+                ddists = {}
+                if Nc > 0:
+                    iic = gen.multi_index_vec(inti,ip)
+                    mindist = 1000000.
+                    # Sorting from distance respect the vertex ip
+                    for ic in range(Nc):
+                        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,:])
     
-        pilepolygons[polyns[Npolys-1]] = newipoly
+            newipoly = np.array(newipoly)
+            pilepolygons[polyns[Npolys-1]] = rm_consecpt_polygon(newipoly)
     
-    for polyn in polyns:
-        poly = pilepolygons[polyn]
-        poly = ma.masked_equal(poly, gen.fillValueF)
-        pilepolygons[polyn] = poly
+        for polyn in polyns:
+            poly = pilepolygons[polyn]
+            poly = ma.masked_equal(poly, gen.fillValueF)
+            pilepolygons[polyn] = poly
 
     return pilepolygons