Changeset 2547 in lmdz_wrf

Timestamp:

May 19, 2019, 4:46:20 PM (6 years ago)

Author:

lfita

Message:

Fixing `cut_ypolygon'

File:

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

trunk/tools/geometry_tools.py

@@ -23 +23 @@
 
 ####### Contents:
-# cut_ypolygon: Function to cut a polygon from a given value of the y-axis
+# cut_[x/y]polygon: Function to cut a polygon from a given value of the [x/y]-axis
 # deg_deci: Function to pass from degrees [deg, minute, sec] to decimal angles [rad]
 # dist_points: Function to provide the distance between two points
@@ -448 +448 @@
     if type(polygon) == type(gen.mamat):
         # Assuming clockwise polygons
-        for ip in range(N):
+        for ip in range(N-1):
             if not polygon.mask[ip,0]:
                 eep = ip + 1
@@ -464 +464 @@
                     dx = polygon[eep,1] - polygon[ip,1]
                     dy = polygon[eep,0] - polygon[ip,0]
-                    dd = yval - polygon[eep,0]
-                    ept = [yval, polygon[eep,1]+dx*dd/dy]
+                    dd = yval - polygon[ip,0]
+                    print polygon[ip,:], ':', polygon[eep,:]
+                    print 'dx', dx, 'dy', dy, 'dd', dd
+                    ept = [yval, polygon[ip,1]+dx*dd/dy]
     else:
         # Assuming clockwise polygons
-        for ip in range(N):
+        for ip in range(N-1):
             eep = ip + 1
             if eep == N: eep = 0
@@ -484 +486 @@
                 dy = polygon[eep,0] - polygon[ip,0]
                 dd = yval - polygon[ip,0]
+                print polygon[ip,:], ':', polygon[eep,:]
+                print 'dx', dx, 'dy', dy, 'dd', dd
                 ept = [yval, polygon[ip,1]+dx*dd/dy]
+
+    print 'Lluis ipt:', ipt, 'ept:', ept
 
     if ipt is None or ept is None:
         print errormsg
         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,:].filled(gen.fillValueF)
+        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)
+        if type(polygon) == type(gen.mamat):
+            cutpolygon[0:ecut-icut-1,:] = polygon[icut+1:ecut,:].filled(gen.fillValueF)
+        else:
+            cutpolygon[0:ecut-icut-1,:] = 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
+    cutpolygon[iip:iip+Nadd,:] = cutline
+    if keep == 'bottom':
+        if type(polygon) == type(gen.mamat):
+            cutpolygon[iip+Nadd:Npts,:] = polygon[ecut+1:N,:].filled(gen.fillValueF)
+        else:
+            cutpolygon[iip+Nadd:Npts,:] = polygon[ecut+1:N,:]
+
+    rmpolygon = []
+    if keep == 'bottom':
+        for ip in range(N):
+            if polygon[ip,0] > yval:
+                rmpolygon.append([gen.fillValueF, gen.fillValueF])
+            else:
+                rmpolygon.append(polygon[ip,:])
+    else:
+        for ip in range(N):
+            if polygon[ip,0] < yval:
+                rmpolygon.append([gen.fillValueF, gen.fillValueF])
+            else:
+                rmpolygon.append(polygon[ip,:])
+    Npts = len(rmpolygon)
+    cutpolygon = np.array(rmpolygon)
+
+    cutpolygon = ma.masked_equal(cutpolygon, gen.fillValueF)
+
+    return Npts, cutpolygon
+
+def cut_xpolygon(polygon, xval, keep='left', Nadd=20):
+    """ Function to cut a polygon from a given value of the x-axis
+      polygon: polygon to cut
+      xval: value to use to cut the polygon
+      keep: part to keep from the cut ('left', default)
+         'left': left from the cut
+         'right': right from the cut
+      Nadd: additional points to add to draw the line (20, default)
+    """
+    fname = 'cut_xpolygon'
+
+    N = polygon.shape[0]
+    availkeeps = ['left', 'right']
+
+    if not gen.searchInlist(availkeeps, keep):
+        print errormsg
+        print '  ' + fname + ": wring keep '" + keep + "' value !!"
+        print '    available ones:', availkeeps
+        quit(-1)
+
+    ipt = None
+    ept = None
+
+    if type(polygon) == type(gen.mamat):
+        # Assuming clockwise polygons
+        for ip in range(N):
+            if not polygon.mask[ip,0]:
+                eep = ip + 1
+                if eep == N: eep = 0
+      
+                if polygon[ip,1] <= xval and polygon[eep,1] >= xval:
+                    icut = ip
+                    dx = polygon[eep,1] - polygon[ip,1]
+                    dy = polygon[eep,0] - polygon[ip,0]
+                    dd = xval - polygon[ip,1]
+                    ipt = [ polygon[ip,0]+dy*dd/dx, xval]
+
+                if polygon[ip,1] >= xval and polygon[eep,1] <= xval:
+                    ecut = ip
+                    dx = polygon[eep,1] - polygon[ip,1]
+                    dy = polygon[eep,0] - polygon[ip,0]
+                    dd = xval - polygon[eep,1]
+                    ept = [polygon[eep,0]+dy*dd/dx, xval]
+    else:
+        # Assuming clockwise polygons
+        for ip in range(N):
+            eep = ip + 1
+            if eep == N: eep = 0
+      
+            if polygon[ip,1] <= xval and polygon[eep,1] >= xval:
+                icut = ip
+                dx = polygon[eep,1] - polygon[ip,1]
+                dy = polygon[eep,0] - polygon[ip,0]
+                dd = xval - polygon[ip,1]
+                ipt = [polygon[ip,0]+dy*dd/dx, xval]
+
+            if polygon[ip,1] >= xval and polygon[eep,1] <= xval:
+                ecut = ip
+                dx = polygon[eep,1] - polygon[ip,1]
+                dy = polygon[eep,0] - polygon[ip,0]
+                dd = xval - polygon[ip,1]
+                ept = [polygon[ip,0]+dy*dd/dx, xval]
+
+    if ipt is None or ept is None:
+        print errormsg
+        print '  ' + fname + ': no cutting for polygon at y=', xval, '!!'
         irmv = 0
         rmpolygon = []
-        if keep == 'bottom':
+        if keep == 'left':
             for ip in range(N):
-                if polygon[ip,0] > yval:
+                if polygon[ip,1] > xval:
                     rmpolygon.append([gen.fillValueF, gen.fillValueF])
                 else:
@@ -499 +624 @@
         else:
             for ip in range(N):
-                if polygon[ip,0] < yval:
+                if polygon[ip,0] < xval:
                     rmpolygon.append([gen.fillValueF, gen.fillValueF])
                 else:
@@ -507 +632 @@
         cutpolygon = np.array(rmpolygon)
     else:
-        if keep == 'bottom':
+        if keep == 'left':
             Npts = icut + (N-ecut) + Nadd
 
@@ -535 +660 @@
         cutline[Nadd-1,:] = ept
         cutpolygon[iip:iip+Nadd,:] = cutline
-        if keep == 'bottom':
+        if keep == 'left':
             cutpolygon[iip+Nadd:Npts,:] = polygon[ecut+1:N,:]
 
@@ -1567 +1692 @@
       'sign2': [buoy[N2+N32+2:N,:],'-','k',1.5]}
 
+    # painting it
+    Height = height + 1.8*bradii
+#    print 'Lluis Height', Height
+#    Ncut, halfdown = cut_ypolygon(buoy, yval=Height/2., keep='bottom')
+#    Ncut, halfup = cut_ypolygon(buoy, yval=Height/2., keep='above')
+
+#    buoysecs.append('halfk')
+#    buoydic['halfk'] = [halfup, '-', 'k', 1.]
+#    buoysecs.append('halfy')
+#    buoydic['halfy'] = [halfdown, '-', '#FFFF00', 1.]
+
     return buoy, buoysecs, buoydic