Changeset 2564 in lmdz_wrf

Timestamp:

May 29, 2019, 9:56:33 PM (6 years ago)

Author:

lfita

Message:

Adding:

• `cut_between_xpolygon': Function to cut a polygon between 2 given value of the x-axis

File:

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

trunk/tools/geometry_tools.py

@@ -23 +23 @@
 
 ####### Contents:
-# cut_between_[x/y]polygon: Function to cut a polygon between 2 given value of the y-axis
+# cut_between_[x/y]polygon: Function to cut a polygon between 2 given value of the [x/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]
@@ -734 +734 @@
             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,:]
+        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
+
+def cut_between_xpolygon(polygon, xval1, xval2, Nadd=20):
+    """ Function to cut a polygon between 2 given value of the x-axis
+      polygon: polygon to cut
+      xval1: first value to use to cut the polygon
+      xval2: first value to use to cut the polygon
+      Nadd: additional points to add to draw the line (20, default)
+    """
+    fname = 'cut_betwen_xpolygon'
+
+    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 xval1 > xval2:
+        print errormsg
+        print '  ' + fname + ': wrong between cut values !!'
+        print '     it is expected xval1 < xval2'
+        print '     values provided xval1: (', xval1, ')> xval2 (', xval2, ')'
+        quit(-1)
+
+    xvals = [xval1, xval2]
+
+    for ic in range(2):
+        xval = xvals[ic]
+        if type(polygon) == type(gen.mamat):
+            # Assuming clockwise polygons
+            for ip in range(N-1):
+                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[ic] = ip
+                        dx[ic] = polygon[eep,1] - polygon[ip,1]
+                        dy[ic] = polygon[eep,0] - polygon[ip,0]
+                        dd = xval - polygon[ip,1]
+                        ipt[ic,:] = [polygon[ip,0]+dy[ic]*dd/dx[ic], xval]
+
+                    if polygon[ip,1] >= yval and polygon[eep,1] <= xval:
+                        ecut[ic] = ip
+                        dx[ic] = polygon[eep,1] - polygon[ip,1]
+                        dy[ic] = polygon[eep,0] - polygon[ip,0]
+                        dd = xval - polygon[ip,1]
+                        ept[ic,:] = [polygon[ip,0]+dy[ic]*dd/dx[ic], xval]
+        else:
+            # Assuming clockwise polygons
+            for ip in range(N-1):
+                eep = ip + 1
+                if eep == N: eep = 0
+      
+                if polygon[ip,1] <= xval and polygon[eep,1] >= xval:
+                    icut[ic] = ip
+                    dx[ic] = polygon[eep,1] - polygon[ip,1]
+                    dy[ic] = polygon[eep,0] - polygon[ip,0]
+                    dd = xval - polygon[ip,1]
+                    ipt[ic,:] = [polygon[ip,0]+dy[ic]*dd/dx[ic], xval]
+
+                if polygon[ip,1] >= xval and polygon[eep,1] <= xval:
+                    ecut[ic] = ip
+                    dx[ic] = polygon[eep,1] - polygon[ip,1]
+                    dy[ic] = polygon[eep,0] - polygon[ip,0]
+                    dd = xval - polygon[ip,1]
+                    ept[ic,:] = [polygon[ip,0]+dy[ic]*dd/dx[ic], xval]
+
+        if ipt is None or ept is None:
+            print errormsg
+            print '  ' + fname + ': no cutting for polygon at x=', xval, '!!'
 
     Npts = icut[1] - icut[0] + Nadd + ecut[0] - ecut[1]