Changeset 631 in lmdz_wrf for trunk

Timestamp:

Sep 16, 2015, 6:00:32 PM (9 years ago)

Author:

lfita

Message:

Working on the radii average

File:

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

trunk/tools/nc_var_tools.py

@@ -15332 +15332 @@
     return radiipos, SQradiipos
 
+def hor_weight_int(SWval, SEval, NWval, NEval, xpos, ypos):
+    """ Function to weighted by distance interpolate a horizontal value from it closest 4 neighbourgs 
+      field= a 4 points representing the environment of a given point
+      xpos, ypos: position to which one want to interpolate (relative to 0,0 at the corner SW)
+    >>> hor_weight_int(1., 1., 1., 2., 0.75, 0.75)
+    1.44888128193
+    """
+    fname = 'hor_weight_int'
+
+    vals = np.array([SWval, SEval, NWval, NEval])
+    print vals
+    if xpos > 1.: 
+        print errormsg
+        print '  ' + fname + ': Wrong position to interpolate! (',xpos,',',ypos,     \
+          ') must be within (1,1) !!'
+        quit(-1)
+
+    vertexs = np.zeros((4,2), dtype=np.float) 
+    vertexs[0,0] = 0.
+    vertexs[0,1] = 0.
+    vertexs[1,0] = 1.
+    vertexs[1,1] = 0.
+    vertexs[2,0] = 0.
+    vertexs[2,1] = 1.
+    vertexs[3,0] = 1.
+    vertexs[3,1] = 1.
+    print vertexs
+
+    dist = np.sqrt((vertexs[:,0]-xpos)**2 + (vertexs[:,1]-ypos)**2)
+    print dist
+    done = False
+    for id in range(4):
+        if dist[id] == 0.: 
+            intval = vals[id]
+            done = True
+
+    if not done: intval = np.sum(vals/dist)/np.sum(1./dist)
+
+    return intval
+
 def compute_tevolboxtraj_radialsec(values, ncfile, varn):
     """ Function to compute tevolboxtraj_radialsec: temporal evolution at a given point along 
@@ -15338 +15378 @@
         [trajfile]: ASCII file with the trajectory ('#' not readed)
           [time] [xpoint] [ypoint]
-        [Tbeg]: equivalent first time-step of the trajectory within the netCDF file
+        [Tbeg]: equivalent first time-step of the traqjectory within the netCDF file
         [lonname],[latname],[zname],[timename]: longitude, latitude, z and time variables names
         [timekind]: kind of time
@@ -15433 +15473 @@
     trjSQradii = {}
 
+    TOTradii = []
+
     iline = 0
     it = 0
@@ -15458 +15500 @@
             gtrajvals[it,1] = int(line.split(' ')[1])
             gtrajvals[it,2] = int(line.split(' ')[2])
-#            print it,'t:',gtrajvals[it,0],'x y:', gtrajvals[it,1], gtrajvals[it,2]
+            print it,'t:',gtrajvals[it,0],'x y:', gtrajvals[it,1], gtrajvals[it,2]
         
             if timekind == 'wrf':
@@ -15484 +15526 @@
             posangle = np.zeros((radii,2), dtype=np.float)
 
-            trjradii[it,:,:,:], trjSQradii[it] = radii_points(trajvals[it,0],        \
-              trajvals[it,1],radii,Nangle,dimx,dimy)
-            print trjradii[it,:,:,:]
-            print '-----------'
-            print trjSQradii[it,:,:,:]
-
+            trjradii[it,:,:,:], trjSQradii[it] = radii_points(gtrajvals[it,1],        \
+              gtrajvals[it,2],radii,Nangle,dimx,dimy)
+            if it == 0: 
+                TOTradii = list(trjSQradii[it].keys())
+            else:
+                for ir in trjSQradii[it].keys():
+                    if searchInlist(TOTradii,ir): TOTradii.append(ir)
             it = it + 1
             iline = iline + 1
 
-            quit()
 
     trajobj.close()
 # Creation of the netCDF file
 ##
+    NtotSQradii = len(TOTradii)
+    TOTSQradii = sorted(TOTradii)
+    print '  ' + fname + ':',NtotSQradii,' squared radii:',TOTSQradii
+
     if varn.find(',') != -1:
         varnS = 'multi-var'
@@ -15503 +15549 @@
         varnS = varn.replace(',','-')
 
-    ofile = 'tevolboxtraj_' + varnS + '.nc'
+    ofile = 'tevolradiitraj_' + varnS + '.nc'
     objofile = NetCDFFile(ofile, 'w')
 
+    boxs = radii*2+1
 # Dimensions
     newdim = objofile.createDimension('x', boxs)
     newdim = objofile.createDimension('y', boxs)
-    newdim = objofile.createDimension('xr', Nrad*2+1)
-    newdim = objofile.createDimension('yr', Nrad*2+1)
     newdim = objofile.createDimension('time', None)
+    newdim = objofile.createDimension('radii', radii)
+    newdim = objofile.createDimension('SQradii', NtotSQradii)
 
     stsn = ['min', 'max', 'mean', 'mean2', 'stdev', 'ac']
@@ -15521 +15568 @@
     cstatnames = []
     for i in range(Nsts):
-        statnames.append(stsn[i] + 'box')
-        cstatnames.append(stsn[i] + 'circle')
+        statnames.append(stsn[i] + 'radii')
+        cstatnames.append(stsn[i] + 'SQradii')
 
 # Getting values
 ##
     ivar = 0
+    box2=0
 
     for vn in varns:
@@ -15581 +15629 @@
  
             varvals = np.ones(tuple([Ttraj,dimz,boxs,boxs]), dtype=np.float)
-            lonvals = np.ones(tuple([Ttraj,dimz,boxs,boxs]), dtype=np.float)
+            lonvals = np.ones(tuple([Ttraj,boxs,boxs]), dtype=np.float)
             latvals = np.ones(tuple([Ttraj,boxs,boxs]), dtype=np.float)
-            rvarvals = np.ones(tuple([Ttraj,dimz,Nrad*2+1,Nrad*2+1]), dtype=np.float)
-            rlonvals = np.ones(tuple([Ttraj,dimz,Nrad*2+1,Nrad*2+1]), dtype=np.float)
-            rlatvals = np.ones(tuple([Ttraj,dimz,Nrad*2+1,Nrad*2+1]), dtype=np.float)
+            radiivarvals = np.ones(tuple([Ttraj,dimz,radii]), dtype=np.float)*fillValueF
+            SQradiivarvals = np.ones(tuple([Ttraj,dimz,NtotSQradii]), dtype=np.float)*fillValueF
 
 # One dimension plus for the values at the center of the trajectory
@@ -15593 +15640 @@
             for it in range(Ttraj):
                 it0 = Tbeg + it
+# Radii values
+# trjradii[it,:,:,:], trjSQradii[it]
+                for ir in range(radii):
+                    for ia in range(Nangle):
+                        xp = trjradii[it,ia,ir,0]
+                        yp = trjradii[it,ia,ir,1]
+                        xpI = int(xp*1.)
+                        ypI = int(yp*1.)
+                        print 'xp yp:',xp,yp,'xpI ypI:',xpI,ypI
+                        quit()
+                        mean = 0. 
+                        for iz in range(dimz):
+                            SWv = varobj[it0,iz,ypI,xpI]
+                            SEv = varobj[it0,iz,ypI,xpI+1]
+                            NWv = varobj[it0,iz,ypI+1,xpI]
+                            NEv = varobj[it0,iz,ypI+1,xpI+1]
+                            val = hor_weight_int(SEv, SWv, NEv, NWv, xp-xpI, yp-ypI)
+                            print 'Lluis hor_weight_int(',SEv,',', SWv,',', NEv,',', NWv,',', xp-xpI,',', yp-ypI,')'
+                            print val
+                            mean = mean + val 
+                        print 'mean:',mean
+                        quit()
 
                 slicev = []
@@ -15833 +15902 @@
             dimt = varobj.shape[0]
  
-            varvals = np.ones(tuple([Ttraj,boxs,boxs1]), dtype=np.float)
+            Nrad=1
+            varvals = np.ones(tuple([Ttraj,boxs,boxs]), dtype=np.float)
             lonvals = np.ones(tuple([Ttraj,boxs,boxs]), dtype=np.float)
             latvals = np.ones(tuple([Ttraj,boxs,boxs]), dtype=np.float)
@@ -16135 +16205 @@
 #        [radii]: length in grid points of the radius
 
-compute_tevolboxtraj_radialsec('/bdd/PCER/workspace/lfita/etudes/FF/medic051215_2km/run/trajectory.dat@6,XLONG,XLAT,ZNU,Times,wrf,4,10',       \
-  '/bdd/PCER/workspace/lfita/etudes/FF/medic051215_2km/run/wrfout/wrfout_d02_2005-12-12_00:00:00', 'U10')
+compute_tevolboxtraj_radialsec('/bdd/PCER/workspace/lfita/etudes/FF/medic051215_2km/run/trajectory_shorten.dat@0,XLONG,XLAT,ZNU,Times,wrf,4,2',       \
+  '/bdd/PCER/workspace/lfita/etudes/FF/medic051215_2km/run/wrfout/wrfout_d02_2005-12-13_00:00:00', 'QVAPOR')
 quit()