Changeset 625 in lmdz_wrf

Timestamp:

Sep 13, 2015, 7:37:36 PM (9 years ago)

Author:

lfita

Message:

Adding radial averaged plots

Location:

trunk/tools

Files:

• drawing_tools.py (modified) (2 diffs)
• nc_var_tools.py (modified) (5 diffs)

trunk/tools/drawing_tools.py

@@ -57 +57 @@
 # plot_2D_shadow_line:
 # plot_lines: Function to plot a collection of lines
+# plot_ZQradii: Function to plot following radial averages only at exact grid poins
 
 # From nc_var_tools.py
@@ -5448 +5449 @@
 
     return
+
+def plot_ZQradii(Zmeans, graphtit, kfig, figname):
+    """ Function to plot following radial averages only at exact grid poins
+      Zmeans= radial means
+      radii= values of the taken radii
+      graphtit= title of the graph ('|', for spaces)
+      kfig= kind of figure
+      figname= name of the figure
+    """
+
+    fname = 'plot_ZQradii'
+
+    output_kind(kfig, figname, True)
+
+    return

trunk/tools/nc_var_tools.py

@@ -15155 +15155 @@
     return radpos
 
+def grid_combinations(x,y):
+    """ Function to provide all the possible grid points combination for a given pair of values
+      x,y= pair of grid points
+    >>>grid_combinations(1,2)
+    [[1, 2], [-1, 2], [-1, -2], [1, -2], [2, 1], [-2, 1], [-2, -1], [2, -1]]
+    """
+    fname = 'grid_combinations'
+
+    gridcomb = []
+    gridcomb.append([x,y])
+    gridcomb.append([-x,y])
+    gridcomb.append([-x,-y])
+    gridcomb.append([x,-y])
+
+    if x != y:
+        gridcomb.append([y,x])
+        gridcomb.append([-y,x])
+        gridcomb.append([-y,-x])
+        gridcomb.append([y,-x])
+
+    return gridcomb
+
+def radii_points(xpos,ypos,Lrad,Nang,dx,dy):
+    """ Function to provide the grid points for radial cross sections, by angle and radii and `squared' radii
+      xpos= x position of the center
+      ypos= y position of the center
+      Lrad= length of the maximum radi in grid points
+      Nang= number of equivalent angles
+      dx= x size of the domain of values
+      dy= y size of the domain of values
+    """
+    fname = 'radii_points'
+
+    radiipos = np.zeros((Nangle, radii, 2), dtype = int)
+    SQradiipos = {}
+
+# Getting the range of radii:
+    pairsSQradii = read_SQradii(Npt=Lrad)
+    SQradii = pairsSQradii.keys()
+    Nradii = len(SQradii)
+
+# Angle loop
+    for ia in range(Nang):
+        angle = 2.*np.pi*ia/Nang
+
+# Points at that given angle
+        pangle = radial_points(angle,radii)
+        posangle[:,0] = pangle[:,0] + xpos
+        posangle[:,1] = pangle[:,1] + ypos
+
+        for ip in range(radii):
+            iipos = int(posangle[ip,0])
+            jjpos = int(posangle[ip,1])
+
+# Creation of a matrix with all the grid points at each time-step
+            if iipos >= 0 and iipos < dx and jjpos >= 0 and jjpos < dy:
+                radiipos[ia,ip,0] = iipos
+                radiipos[ia,ip,1] = jjpos
+
+# SQ radii values
+
+# Radii loop (avoiding 0)
+    pairswithin = {}
+    for ir in range(1,Nradii):
+        pairs = pairsSQradii[SQradii[ir]]
+
+# pairs loop
+        pairsw = []
+        for ip in range(len(pairs)):
+            allpairs = grid_combinations(pairs[ip*2],pairs[ip*2+1])
+
+            for iap in range(len(allpairs)):
+                ppos = allpairs[iap]
+                iipos = trajvals[it,1] + ppos[0]
+                jjpos = trajvals[it,1] + ppos[0]
+
+# Creation of a matrix with all the grid points at each time-step
+                if iipos >= 0 and iipos < dx and jjpos >= 0 and jjpos < dy:
+                    pairsw.append([iipos,jjpos])
+
+        SQradiipos[SQradii[ir]] = pairsw
+
+    return radiipos, SQradiipos
+
 def compute_tevolboxtraj_radialsec(values, ncfile, varn):
     """ Function to compute tevolboxtraj_radialsec: temporal evolution at a given point along 
@@ -15245 +15329 @@
 # Trajectory values/slices
 ##
+# This is how we're going to proceed:
+#    Creation of 1 matrix and 1 dictionay with the grid points where values have to be taken
+#      trjradii(dimt,Nangle,Nrad,2): x,y pairs at each radii, angle and time step. 
+#        Ditance interpolation might be required!
+#      SQtrjradii(dimt,2): x,y pairs at each radii, number of pairs (Npair) and time
+#        step. Flip/rotation of points might be required!
+
+    trjradii = np.zeros((dimt, Nangle, radii, 2), dtype = int)
+    trjSQradii = {}
+
     iline = 0
     it = 0
-
-    xrangeslice = []
-    yrangeslice = []
-    xrangeslice2D = []
-    yrangeslice2D = []
-
-    cslicev = []
-    cslice2D = []
-    cslicevnoT = []
-
-    gtrajvals = np.zeros((Ttraj,3), dtype=int)
-    trajvals = np.zeros((Ttraj,3), dtype=np.float)
 
     trajpos = np.zeros(Ttraj,Nangle,raddi)
@@ -15283 +15365 @@
                 trajvals[it,0] = realdatetime1_CFcompilant(gdate, '19491201000000',  \
                   'hours')
+
                 tunits = 'hours since 1949-12-01 00:00:00'
             elif timekind == 'cf':
@@ -15292 +15375 @@
                 quit(-1)
 
-            trajvals[it,1] = lonv[gtrajvals[it,2],gtrajvals[it,1]]
-            trajvals[it,2] = latv[gtrajvals[it,2],gtrajvals[it,1]]
-
 #            print iline, it,'time:',trajvals[it,0],'lon lat:', trajvals[it,1],       \
 #              trajvals[it,2]
@@ -15305 +15385 @@
             posangle = np.zeros((radii,2), dtype=np.float)
 
-# Angle loop
-            for ia in range(Nangle):
-                angle = 2.*np.pi*ia/Nangle
-
-# Points at that given angle
-                pangle = radial_points(angle,radii)
-                posangle[:,0] = pangle[:,0] + trajvals[it,1]
-                posangle[:,1] = pangle[:,1] + trajvals[it,2]
-
-                for ip in range(radii):
-                    iipos = int(posangle[ip,0])
-                    jjpos = int(posangle[ip,1])
-
-# Creation of a matrix with all the grid points at each time-step
-                    if iipos >= 0 and iipos < dimx and jjpos >= 0 and jjpos < dimy:
-                         print 'Hola'
-                         quit()
-#                        varvalst[ia,ip] = 
-
-
-                if gtrajvals[it,2]-radii < 0: 
-                    yinit = 0
-                    yinit2D = radii - gtrajvals[it,2]
-                else:
-                    yinit = gtrajvals[it,2]-radii
-                    yinit2D = 0
-
-                if gtrajvals[it,2]+radii + 1 > dimy + 1: 
-                    yend = dimy + 1
-                    yend2D = dimy + 1 - gtrajvals[it,2] + radii
-                else:
-                    yend = gtrajvals[it,2]+radii + 1
-                    yend2D = radi
-
-                if gtrajvals[it,1]-box2 < 0: 
-                    xinit = 0
-                    xinit2D = box2 - gtrajvals[it,1]
-                else:
-                    xinit = gtrajvals[it,1]-box2
-                    xinit2D = 0
-
-                if gtrajvals[it,1]+box2 + 1 > dimx + 1: 
-                    xend = dimx + 1
-                    xend2D = dimx + 1 - gtrajvals[it,1] - box2
-                else:
-                    xend = gtrajvals[it,1]+box2 + 1
-                    xend2D = boxs
-
-                yrangeslice.append([yinit, yend])
-                xrangeslice.append([xinit, xend])
-                yrangeslice2D.append([yinit2D, yend2D])
-                xrangeslice2D.append([xinit2D, xend2D])
-            else:
-                yrangeslice.append([gtrajvals[it,2]-box2, gtrajvals[it,2]+box2 + 1])
-                xrangeslice.append([gtrajvals[it,1]-box2, gtrajvals[it,1]+box2 + 1])
-                yrangeslice2D.append([gtrajvals[it,2]-box2, gtrajvals[it,2]+box2 + 1])
-                xrangeslice2D.append([gtrajvals[it,1]-box2, gtrajvals[it,1]+box2 + 1])
-
-# circle values
-            circdist[it,:,:] = radius_dist(dimy,dimx,gtrajvals[it,2],gtrajvals[it,1])
-
-            if gtrajvals[it,2]-Nrad < 0 or gtrajvals[it,2]+Nrad + 1 > dimy + 1 \
-              or gtrajvals[it,1]-Nrad < 0 or gtrajvals[it,1]+Nrad + 1 > dimx + 1:
-
-                if gtrajvals[it,2]-Nrad < 0: 
-                    yinit = 0
-                    yinit2D = Nrad - gtrajvals[it,2]
-                else:
-                    yinit = gtrajvals[it,2]-Nrad
-                    yinit2D = 0
-
-                if gtrajvals[it,2]+Nrad + 1 > dimy + 1: 
-                    yend = dimy + 1
-                    yend2D = dimy + 1 - gtrajvals[it,2] + Nrad
-                else:
-                    yend = gtrajvals[it,2]+Nrad + 1
-                    yend2D = 2*Nrad+1
-
-                if gtrajvals[it,1]-Nrad < 0: 
-                    xinit = 0
-                    xinit2D = Nrad - gtrajvals[it,1]
-                else:
-                    xinit = gtrajvals[it,1]-Nrad
-                    xinit2D = 0
-
-                if gtrajvals[it,1]+Nrad + 1 > dimx + 1: 
-                    xend = dimx + 1
-                    xend2D = dimx + 1 - gtrajvals[it,1] - Nrad
-                else:
-                    xend = gtrajvals[it,1]+Nrad + 1
-                    xend2D = 2*Nrad+1
-
-                cyrangeslice.append([yinit, yend])
-                cxrangeslice.append([xinit, xend])
-                cyrangeslice2D.append([yinit2D, yend2D])
-                cxrangeslice2D.append([xinit2D, xend2D])
-            else:
-                cyrangeslice.append([gtrajvals[it,2]-Nrad, gtrajvals[it,2]+Nrad + 1])
-                cxrangeslice.append([gtrajvals[it,1]-Nrad, gtrajvals[it,1]+Nrad + 1])
-                cyrangeslice2D.append([gtrajvals[it,2]-Nrad, gtrajvals[it,2]+Nrad+1])
-                cxrangeslice2D.append([gtrajvals[it,1]-Nrad, gtrajvals[it,1]+Nrad+1])
+            trjradii[it,:,:,:], trjSQradii[it] = radii_points(trajvals[it,0],        \
+              trajvals[it,1],radii,Nangle,dimx,dimy)
+
             it = it + 1
             iline = iline + 1
+            quit()
 
     trajobj.close()