Changeset 1391 in lmdz_wrf for trunk/tools

Timestamp:

Dec 19, 2016, 4:09:08 PM (8 years ago)

Author:

lfita

Message:

Changing no-value by masked array on vals_around', pass changes to linearint_3x3weights'
Adding:
curvelocalize_2D: Function to provide the localization a curve in a 2D field of positions via the equivalent

i,j-quads within which the curve lays (-1, no value for the given quad [2x2 points around]) and provide
the weights of the quad to perform a distance-weighted mean at the given curve point

File:

• trunk/tools/generic_tools.py (modified) (11 diffs)

trunk/tools/generic_tools.py

@@ -55 +55 @@
 # color_deg: Function to generate a degradation of colors in rgb base
 # contflow: Function to bring back the increment in j,i grid points according to a trip: (inflow directions)
+# curvelocalize_2D: Function to provide the localization a curve in a 2D field of positions via the equivalent
+#   i,j-quads within which the curve lays (-1, no value for the given quad [2x2 points around]) and provide 
+#   the weights of the quad to perform a distance-weighted mean at the given curve point
 # dictionary_key: Function to provide the first key in a dictionay with a given value
 # dictionary_key_list: Function to provide the first key in a dictionary of lists with a given value
@@ -7592 +7595 @@
      [None None None]]
     """
+    import numpy.ma as ma
     fname = 'vals_around'
 
-    vals3x3 = np.array([None]*9).reshape(3,3)
+    # Issues with 'None'
+    #vals3x3 = np.array([None]*9).reshape(3,3)
+    vals3x3 = ma.masked_equal(np.ones((3,3), dtype=np.float), 1.)
 
     dx = vals.shape[1]
@@ -10803 +10809 @@
     [  5.00000000e-04   9.99500000e-01]
     """
-    fname = 'linear_wegiths'
+    fname = 'linear_weights'
 
     if vals[0] > vals[1]: 
@@ -10839 +10845 @@
         * Number of closest grid points (<=4)
         * 2x9 matrix with distance-sorted positions within the 3x3 matrix
-    xvals = np.zeros(9).reshape(3,3)
+    xvals = ma.masked_equal(np.zeros(9).reshape(3,3),0.)
     xvals[:,0] = 1.
     xvals[:,1] = 2.
     xvals[:,2] = 3.
-    xvals[2,2] = None
+    xvals.mask[2,2] = True
     yvals = xvals.transpose()
     >>> linearint_3x3weights(xvals, yvals, [1.5, 2.5])
@@ -10858 +10864 @@
            [ 1,  2,  1,  2,  0,  0,  1,  0, -1]]))
     """
-    fname = 'linear_wegiths'
+    fname = 'linear_3x3wegiths'
 
     intval = np.array(intval0, dtype=np.float)
-    print 'type xvals:',  xvals.dtype
-
-    # Assuming that values come from `vals_around' function (None, when outside matrix size)
-    xvalsma = ma.masked_equal(xvals,None)
-    yvalsma = ma.masked_equal(yvals,None)
 
     # Check data
-    nx = np.min(xvals)
-    xx = np.max(xvals)
-    ny = np.min(yvals)
-    xy = np.max(yvals)
+    nx = xvals.min()
+    xx = xvals.max()
+    ny = yvals.min()
+    xy = yvals.max()
+
+    #print 'Lluis intval:', intval
+    #print '  ' + fname + ' Lluis xvals ____'
+    #print xvals
+
+    #print '  ' + fname + ' Lluis yvals ____'
+    #print yvals
 
     if not intval[1] >= nx and intval[1] <= xx:
         print errormsg
-        print '  ' + fname + ': wrong value:', intval[1] , ' to interpolate for range:'
+        print '  ' + fname + ': wrong value:', intval[1] , ' to interpolate for ' +  \
+          'x-range:'
         for i in range(3):
             print '    ', xvals[:,i]
@@ -10881 +10890 @@
     if not intval[0] >= ny and intval[0] <= xy:
         print errormsg
-        print '  ' + fname + ': wrong value:', intval[0] , ' to interpolate for range:'
+        print '  ' + fname + ': wrong value:', intval[0] , ' to interpolate for ' +  \
+          'y-range:'
         for i in range(3):
             print '    ', yvals[:,i]
         quit(-1)
 
-    print 'Lluis intval:', intval
-    print '  ' + fname + ' Lluis xvalsma ____'
-    print xvalsma
-
-    print '  ' + fname + ' Lluis yvalsma ____'
-    print yvalsma
-
     dist = np.sqrt( (xvals-intval[1])**2 + (yvals-intval[0])**2)
     inc = np.sum(dist)
-    print 'inc:', inc, 'intval', intval
-    for i in range(3):
-        print dist[:,i]
+    #print 'inc:', inc, 'intval', intval, 'dist ___'
+    #for i in range(3):
+    #    print dist[:,i]
 
     ijvals = dist.copy()
@@ -10903 +10906 @@
 
     # shortest minium distance
-    sortdist = list(dist.flatten())
+    sortdist = list(dist.compressed())
     sortdist.sort()
     iijj = 0    
@@ -10933 +10936 @@
 
     wgt = np.where(dist <= maxdist, dist / inc, 0.)
+    #print 'Lluis weights _______', sortdist[0:4]
+    #for i in range(3):
+    #    print wgt[:,i]
 
     return wgt, Nvals, distpos
-
+ 
 def curvelocalize_2D(curve,xvals,yvals):
     """ Function to provide the localization a curve in a 2D field of positions via the equivalent i,j-quads 
@@ -10956 +10962 @@
       yvals= 2D matrix of y-coordinate values
       RETURNS: 
-        curvloc= [2,Npts] matrix with the closest i,j grid point in the 2D space of the field positions
+        curvloc= [Npts,2] matrix with the closest i,j grid point in the 2D space of the field positions
         curvweights= [Npts,3,3] 3x3 matrix weights around the curve position
         Nwgt= Number of weights
-        ijaorund3x3= [2,9] distance (in xval, yval units) to the curve position
+        ijaorund3x3= [Npts,2,9] distance (in xval, yval units) to the curve position
     """
     fname = 'curvelocalize_2D'
 
-    Npts = curve.shape[1]
+    Npts = curve.shape[0]
 
     curvloc = np.ones((curve.shape), dtype=int)*(-1)
-    Nquad = np.ones((Npts), dtype=int)*(-1)
-    curvweights = np.zeros((Npts,3,3), dtype=float)
+    curveweights = np.zeros((Npts,3,3), dtype=float)
+    ijaround3x3 = np.zeros((Npts,2,9), dtype=int)
 
     for icv in range(Npts):
         # Closest grid-point
-        dist = np.sqrt((xvals - curve[1,icv])**2 + (yvals - curve[0,icv])**2)
+        dist = np.sqrt((xvals - curve[icv,1])**2 + (yvals - curve[icv,0])**2)
         mindist = np.min(dist)
         ijloc = index_mat(dist,mindist)
 
-        curvloc[0,icv] = ijloc[0]
-        curvloc[1,icv] = ijloc[1]
+        if np.any(ijloc == -1): 
+            print errormsg
+            print '  ', fname, ': curve position:', curve[icv,:], 'not found !!'
+            quit(-1)
+        curvloc[icv,0] = ijloc[0]
+        curvloc[icv,1] = ijloc[1]
 
         # Values around the nearest 2D-point to the curve point
         xaroundvls = vals_around(xvals,ijloc)
         yaroundvls = vals_around(yvals,ijloc)
-        print 'Lluis type xaroundvls:', type(xaroundvls)
-
-        curveweights[icv,:,:], Nwgt, ijaorund3x3= linearint_3x3weights(xaroundvls,   \
-          yaroundvls, curve[:,icv])
+
+        curveweights[icv,:,:], Nwgt, ijaround3x3[icv,:,:,] = linearint_3x3weights(    \
+          xaroundvls, yaroundvls, curve[icv,:])
         
     return curvloc, curveweights, Nwgt, ijaround3x3
@@ -11013 +11022 @@
 #Dycurve = (elatcurve-ilatcurve)/(Ncurve-1.)
 
-#curve = np.zeros((2,Ncurve), dtype=np.float)
+#curve = np.zeros((Ncurve,2), dtype=np.float)
 #for icv in range(Ncurve):
-#    curve[0,icv] = ilatcurve + Dycurve*icv
-#    curve[1,icv] = iloncurve + Dxcurve*icv
+#    curve[icv,0] = ilatcurve + Dycurve*icv
+#    curve[icv,1] = iloncurve + Dxcurve*icv
 
 #print 'Curve origin:', ilatcurve, ',', iloncurve
@@ -11024 +11033 @@
 
 #loccurve, wgts, Nwgts, loc3x3 = curvelocalize_2D(curve,lons,lats)
-#print loccurve
+#for ipt in range(Ncurve):
+#    print 'pos:', ipt, 'curve:', curve[ipt,:],'loccurve:', loccurve[ipt,:], 'wgts:', wgts[ipt,:,:], 'Nwgts:', Nwgts, 'loc3x3:', loc3x3[ipt,:]
 #quit()