Changeset 2579 in lmdz_wrf

Timestamp:

Jun 2, 2019, 5:32:08 PM (5 years ago)

Author:

lfita

Message:

Adding to `cut_xpolygon' once polygon has no mask (all masked values as False, which gives a scalar)

File:

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

trunk/tools/geometry_tools.py

@@ -458 +458 @@
     ept = None
 
+    # There might be more than 1 cut...
+    Ncuts = 0
+    icut = []
+    ecut = []
+    ipt = []
+    ept = []
+
     if type(polygon) == type(gen.mamat):
         # Assuming clockwise polygons
@@ -466 +473 @@
       
                 if polygon[ip,0] <= yval and polygon[eep,0] >= yval:
-                    icut = ip
+                    icut.append(ip)
                     dx = polygon[eep,1] - polygon[ip,1]
                     dy = polygon[eep,0] - polygon[ip,0]
                     dd = yval - polygon[ip,0]
-                    ipt = [yval, polygon[ip,1]+dx*dd/dy]
+                    ipt.append([yval, polygon[ip,1]+dx*dd/dy])
 
                 if polygon[ip,0] >= yval and polygon[eep,0] <= yval:
-                    ecut = ip
+                    ecut.append(ip)
                     dx = polygon[eep,1] - polygon[ip,1]
                     dy = polygon[eep,0] - polygon[ip,0]
                     dd = yval - polygon[ip,0]
-                    ept = [yval, polygon[ip,1]+dx*dd/dy]
+                    ept.append([yval, polygon[ip,1]+dx*dd/dy])
+                    Ncuts = Nctus + 1
     else:
         # Assuming clockwise polygons
@@ -485 +493 @@
       
             if polygon[ip,0] <= yval and polygon[eep,0] >= yval:
-                icut = ip
+                icut.append(ip)
                 dx = polygon[eep,1] - polygon[ip,1]
                 dy = polygon[eep,0] - polygon[ip,0]
                 dd = yval - polygon[ip,0]
-                ipt = [yval, polygon[ip,1]+dx*dd/dy]
+                ipt.append([yval, polygon[ip,1]+dx*dd/dy])
 
             if polygon[ip,0] >= yval and polygon[eep,0] <= yval:
-                ecut = ip
+                ecut.append(ip)
                 dx = polygon[eep,1] - polygon[ip,1]
                 dy = polygon[eep,0] - polygon[ip,0]
                 dd = yval - polygon[ip,0]
-                ept = [yval, polygon[ip,1]+dx*dd/dy]
+                ept.append([yval, polygon[ip,1]+dx*dd/dy])
 
     if ipt is None or ept is None:
@@ -502 +510 @@
         print '  ' + fname + ': no cutting for polygon at y=', yval, '!!'
 
-    if keep == 'bottom':
-        Npts = icut + (N-ecut) + Nadd
-        cutpolygon = np.zeros((Npts,2), dtype=np.float)
-        if type(polygon) == type(gen.mamat):
-            cutpolygon[0:icut+1,:] = polygon[0:icut+1,:]
+    Nadds = []
+    Naddc = Nadd/(Ncuts-1)
+    for ic in range(Ncuts-1):
+        Nadds.append(Naddc)
+
+    Nadds.append(N-Naddc*(Ncuts-1))
+
+    iip = 0
+    iipc = 0
+    for ic in range(Ncuts):
+        if keep == 'bottom':
+            Npts = icut[ic] + (N-ecut[ic]) + Nadds[ic]
+            cutpolygon = np.zeros((Npts,2), dtype=np.float)
+            if type(polygon) == type(gen.mamat):
+                cutpolygon[iipc:icut[ic]+1,:] = polygon[iip:icut[ic]+1,:]
+            else:
+                cutpolygon[iipc:icut[ic]+1,:] = polygon[iip:icut[ic]+1,:]
+            iip = icut[ic]+1
         else:
-            cutpolygon[0:icut+1,:] = polygon[0:icut+1,:]
-        iip = icut+1
-    else:
-        Npts = ecut - icut + Nadd-1
-        cutpolygon = np.zeros((Npts,2), dtype=np.float)
-        cutpolygon[0,:] = ipt
-        cutpolygon[1:ecut-icut,:] = polygon[icut+1:ecut,:]
-        iip = ecut-icut-1
-
-    # cutting line
-    cutline = np.zeros((Nadd,2), dtype=np.float)
-    dx = (ept[1] - ipt[1])/(Nadd-2)
-    dy = (ept[0] - ipt[0])/(Nadd-2)
-    cutline[0,:] = ipt
-    for ip in range(1,Nadd-1):
-        cutline[ip,:] = ipt + np.array([dy*ip,dx*ip])
-    cutline[Nadd-1,:] = ept
-    if keep == 'bottom':
-        cutpolygon[iip:iip+Nadd,:] = cutline
-        if type(polygon) == type(gen.mamat):
-            cutpolygon[iip+Nadd:Npts,:] = polygon[ecut+1:N,:]
+            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
+
+        # 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)
+        cutline[0,:] = ipt[ic]
+        for ip in range(1,Nadds[ic]-1):
+            cutline[ip,:] = ipt[ic] + np.array([dy*ip,dx*ip])
+        cutline[Nadds[ic]-1,:] = ept[ic]
+        if keep == 'bottom':
+            cutpolygon[iip:iip+Nadds[ic],:] = cutline
+            cutpolygon[iip+Nadds[ic]:Npts,:] = polygon[ecut[ic]+1:N,:]
         else:
-            cutpolygon[iip+Nadd:Npts,:] = polygon[ecut+1:N,:]
-    else:
-        cutpolygon[iip:iip+Nadd,:] = cutline[::-1,:]
+            cutpolygon[iip:iip+Nadds[ic],:] = cutline[::-1,:]
 
     rmpolygon = []
@@ -577 +592 @@
     ept = None
 
-    if type(polygon) == type(gen.mamat):
+    if type(polygon) == type(gen.mamat) and type(polygon.mask) !=                    \
+      type(gen.mamat.mask[1]):
         # Assuming clockwise polygons
         for ip in range(N-1):
-            if not polygon.mask[ip,0]:
+            if not polygon.mask[ip,1]:
                 eep = ip + 1
                 if eep == N: eep = 0
@@ -619 +635 @@
     if ipt is None or ept is None:
         print errormsg
-        print '  ' + fname + ': no cutting for polygon at y=', yval, '!!'
+        print '  ' + fname + ': no cutting for polygon at x=', xval, '!!'
 
     if keep == 'left':