Changeset 1398 in lmdz_wrf

Timestamp:

Dec 20, 2016, 5:43:18 PM (8 years ago)

Author:

lfita

Message:

Adding:

`nterp_curve': Function to interpolate (as weighted mean with 4-closest ones) a variable following a curve (only working now for 'gridline')

File:

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

trunk/tools/nc_var_tools.py

@@ -79 +79 @@
 # increaseDimvar: Function to increase with 1 dimension an existing variable within a netcdf file. Values of the variable will be repeated along the new dimension
 # insert_variable: Function to insert a variable in an existing file
+# interp_curve: Function to interpolate (as weighted mean with 4-closest ones) a variable following a curve
 # isgattrs: Give a single global attribute of a file and its type
 # isvattrs: Give a single attribute of a variable
@@ -18272 +18273 @@
       varn= name of the variable to fill
     """
-    fname = 'addVar'
+    fname = 'setvar_nc'
 
     if values == 'h':
@@ -18385 +18386 @@
     return
 
+def interp_curve(ov, xdn, ydn, curvpos, wgts, Nwgt, ijcurv):
+    """ Function to interpolate (as weighted mean with 4-closest ones) a variable following a curve
+      using a 3x3 matrix of weights
+      ov= netcdf object to interpolate
+      xdn= x-dimension name of the curve
+      ydn= y-dimension name of the curve
+      curvpos= position of the curve within the grid
+      wgts= matrix with the weights to use to interpolate
+      Nwgt= number of useful weights (<=4)
+      ijcurv= j,i indices of the distance-sorted grid points from 3x3 
+    """
+    fname = 'interp_curve'
+    dimvar = ov.dimensions
+
+    Ncurv = curvpos.shape[0]
+
+    dictslice = {}
+    for dn in ov.dimensions:
+        if dn == xdn or dn == ydn: ds = 0
+        else: ds = -1
+        dictslice[dn] = ds
+
+    varslice, dimrefvals = SliceVarDict(ov, dictslice)
+    varref = ov[varslice]
+    dimref = list(varref.shape)
+    curvevar = np.zeros(tuple(dimref+[Ncurv]), dtype=np.float)
+
+    for icv in range(Ncurv):
+        # Getting all `useful' weights
+        ix = int(curvpos[icv,1])
+        iy = int(curvpos[icv,0])
+        sortedwgt = list(wgts[icv,:,:].flatten())
+        sortedwgt.sort(reverse=True)
+
+#        print fname + '; Lluis icv:', icv ,' _____', Nwgt[icv]
+        for isp in range(Nwgt[icv]):
+            dictslice[ydn] = iy + int(ijcurv[icv,0,isp])
+            dictslice[xdn] = ix + int(ijcurv[icv,1,isp])
+            varslice, dimvarvals = SliceVarDict(ov, dictslice)
+            iwgt = sortedwgt[isp]
+            curvevar[...,icv] = curvevar[...,icv] + ov[tuple(varslice)]*iwgt
+#            print '    isp:', isp, 'iwgt:', iwgt, '<>', iy + int(ijcurv[icv,0,isp]), \
+#              ix + int(ijcurv[icv,1,isp]), 'ival:', ov[tuple(varslice)][0,0], 'curvevar:', curvevar[0,0,icv]
+
+    return curvevar
+
+def curve_section(values,ncfile,varn):
+    """ Function to provide a section of a file following a given curve
+      values= [kindcurve]: kind of curve to follow
+        'gridline',[xdimn],[ydimn],[iBL],[jBL],[iUR],[jUR],[Ncurve]: line from starting `bottom-left' grid point 
+          [iBL],[jBL] to `right-up' grid point [iUR],[jUR] and [Ncurve] number of line-points with 
+          [xdimn],[ydimn] name of the x and y dimensions
+        'lonlatline',[londimn]|[lonvarname],[latdimn]|[latvarname],[iSW],[jSW],[iNE],[jNE],[Ncurve]: line from starting 
+          `SW' lonlat point [lonSW],[latSW] to `NE' lonlat point [lonNE],[latNE] and [Ncurve] line-points 
+          with '[lon/latdimn]|[lon/latvarn]' name of the dimension and variable-dimension of longitude
+          and latitude
+        'gridfilecurve',[xdimn],[ydimn],[curvefile]: line following grid-values in the file 
+          [curvefile](x y; two columns space separated, '#' commment) with [xdimn],[ydimn] name of the x and y dimensions
+        'lonlatfilecurve',[londimn]|[lonvarname],[latdimn]|[latvarname],[curvefile]: line following 
+          lonlat-values in the file [curvefile](lon lat; two columns space separated, '#' commment) with 
+          '[lon/latdimn]|[lon/latvarn]' name of the dimension and variable-dimension of longitude
+          and latitude
+      ncfile= name of the file with the variable to fill
+      varn= name of the variable to get ('all', for all variables)
+    """
+    fname = 'curve_section'
+
+    if values == 'h':
+        print fname + '_____________________________________________________________'
+        print curve_section.__doc__
+        quit()
+
+    availcurves = ['gridline', 'lonlatline', 'gridfilecurve', 'lonlatfilecurve']
+
+    if values.split(',')[0] == 'gridline':
+        expectargs = 'gridline,[xdimn],[ydimn],[iBL],[jBL],[iUR],[jUR],[Ncurve]'
+    elif values.split(',')[0] == 'lonlatline':
+        expectargs = 'lonlatline,[londimn]|[lonvarname],[latdimn]|[latvarname],' +   \
+          '[iSW],[jSW],[iNE],[jNE],[Ncurve]'
+    elif values.split(',')[0] == 'gridfilecurve':
+        expectargs = 'gridfilecurve,[xdimn],[ydimn],[curvefile]'
+    elif values.split(',')[0] == 'lonlatfilecurve':
+        expectargs = 'lonlatfilecurve,[londimn]|[lonvarname],[latdimn]|' +           \
+          '[latvarname],[curvefile]'
+    else:
+        print errormsg
+        print '  ' + fname + ": kind of curve: '" + values.split(',')[0] +           \
+          "' not ready !!"
+        print '    available ones:', availcurves
+        quit(-1)
+
+    gen.check_arguments(fname,values,expectargs,',')
+
+    ofile = 'curve_section.nc'
+
+    onc = NetCDFFile(ncfile, 'r')
+
+    if varn == 'all':
+        varns = onc.variables.keys()
+    else:
+        varns = [varn]
+
+    # Getting curve
+    if values.split(',')[0] == 'gridline':
+        # 'gridline,[xdimn],[ydimn],[iBL],[jBL],[iUR],[jUR],[Ncurve]'
+        xdimn = values.split(',')[1]
+        ydimn = values.split(',')[2]
+        iBL = np.float(values.split(',')[3])
+        jBL = np.float(values.split(',')[4])
+        iUR = np.float(values.split(',')[5])
+        jUR = np.float(values.split(',')[6])
+        Ncurve = int(values.split(',')[7])
+        Dxcurve = (iUR-iBL)/(Ncurve-1.)
+        Dycurve = (jUR-jBL)/(Ncurve-1.)
+
+        curve = np.zeros((Ncurve,2), dtype=np.float)
+        for icv in range(Ncurve):
+            curve[icv,0] = jBL + Dycurve*icv
+            curve[icv,1] = iBL + Dxcurve*icv
+
+    elif values.split(',')[0] == 'lonlatline':
+        # 'lonlatline,[londimn]|[lonvarname],[latdimn]|[latvarname],[iSW],[jSW],[iNE],[jNE],[Ncurve]'
+        Lcurve = 'lon/lat'
+
+        xdimn = values.split(',')[1].split('|')[0]
+        xvardimn = values.split(',')[1].split('|')[1]
+        ydimn = values.split(',')[2].split('|')[0]
+        yvardimn = values.split(',')[2].split('|')[1]
+        lonSW = np.float(values.split(',')[3])
+        latSW = np.float(values.split(',')[4])
+        lonNE = np.float(values.split(',')[5])
+        latNE = np.float(values.split(',')[6])
+        Ncurve = int(values.split(',')[7])
+        Dxcurve = (lonNE-lonSW)/(Ncurve-1.)
+        Dycurve = (latNE-latSW)/(Ncurve-1.)
+
+        curve = np.zeros((Ncurve,2), dtype=np.float)
+        for icv in range(Ncurve):
+            curve[icv,0] = lonSW + Dycurve*icv
+            curve[icv,1] = latSW + Dxcurve*icv
+
+    elif values.split(',')[0] == 'gridfilecurve':
+        # 'gridfilecurve,[xdimn],[ydimn],[curvefile]'
+        xdimn = values.split(',')[1]
+        ydimn = values.split(',')[2]
+        curvefile = values.split(',')[3]
+
+        Scurve = []
+        Ncurve = 0
+        ocurve = open(curvefile, 'r')
+        for line in ocurve:
+            if line[0:1] != '#':
+                Scurve.append(line.split(' '))
+                Ncurve = Ncurve + 1
+        ocurve.close()
+
+        curve = np.zeros((Ncurve,2), dtype=np.float)
+        for icv in range(Ncurve):
+            curv = Scurve[icv]
+            curve[icv,0] = np.float(curv[1])
+            curve[icv,1] = np.float(curv[0])
+
+    elif values.split(',')[0] == 'lonlatfilecurve':
+        # 'lonlatfilecurve,[londimn]|[lonvarname],[latdimn]|[latvarname],[curvefile]'
+        xdimn = values.split(',')[1].split('|')[0]
+        xvardimn = values.split(',')[1].split('|')[1]
+        ydimn = values.split(',')[2].split('|')[0]
+        yvardimn = values.split(',')[2].split('|')[1]
+        curvefile = values.split(',')[3]
+
+        Scurve = []
+        Ncurve = 0
+        ocurve = open(curvefile, 'r')
+        for line in ocurve:
+            if line[0:1] != '#':
+                Scurve.append(line.split(' '))
+                Ncurve = Ncurve + 1
+        ocurve.close()
+
+        curve = np.zeros((Ncurve,2), dtype=np.float)
+        for icv in range(Ncurve):
+            curv = Scurve[icv]
+            curve[icv,0] = np.float(curv[1])
+            curve[icv,1] = np.float(curv[0])
+
+    # Getting dimension-variables
+    onc = NetCDFFile(ncfile, 'r')
+    if not gen.searchInlist(onc.dimensions, xdimn):
+        print errormsg
+        print '  ' + fname + ": file '" + ncfile + "' does not have dimension '" +   \
+          xdimn + "' !!"
+        print '    available dimensions:', onc.dimensions
+        onc.close() 
+        quit(-1)
+
+    if not gen.searchInlist(onc.dimensions, ydimn):
+        print errormsg
+        print '  ' + fname + ": file '" + ncfile + "' does not have dimension '" +   \
+          ydimn + "' !!"
+        print '    available dimensions:', onc.dimensions
+        onc.close() 
+        quit(-1)
+    dimx = len(onc.dimensions[xdimn])
+    dimy = len(onc.dimensions[ydimn])
+
+    if values.split(',')[0]== 'lonlatline' or values.split(',')[0]=='lonlatfilecurve':
+        if not onc.variables.has_key(xvardimn):
+            print errormsg
+            print '  ' + fname + ": file '" + ncfile + "' does not have variable-" + \
+              "dimension '" + xvardimn + "' !!"
+            print '    available variables:', onc.variables.keys()
+            onc.close() 
+            quit(-1)
+        if not onc.variables.has_key(yvardimn):
+            print errormsg
+            print '  ' + fname + ": file '" + ncfile + "' does not have variable-" + \
+              "dimension '" + yvardimn + "' !!"
+            print '    available variables:', onc.variables.keys()
+            onc.close() 
+            quit(-1)
+        ovar = onc.variables[xvardimn]
+        xvals0 = ovar[:]
+        Lcurvex = xvardimn
+        cxunit = ovar.units
+
+        ovar = onc.variables[yvardimn]
+        yvals0 = ovar[:]
+        Lcurvey = yvardimn
+        cyunit = ovar.units
+
+        xvals, yvals = gen.lonlat2D(xvals0, yvals0)
+    else:
+        x1vals = np.arange((dimx), dtype=np.float)
+        y1vals = np.arange((dimy), dtype=np.float)
+        xvals, yvals = gen.lonlat2D(x1vals, y1vals)
+        Lcurvex = 'grid'
+        Lcurvey = 'grid'
+        cxunit = '1'
+        cyunit = '1'
+
+    # Getting curve-values
+    loccurve, curvewgts, Nwgts, curveloc3x3 = gen.curvelocalize_2D(curve,xvals,yvals)
+
+    # Creation of file
+    newnc = NetCDFFile(ofile, 'w')
+
+    # Dimensions
+    newdim = newnc.createDimension('curve', Ncurve)
+    newdim = newnc.createDimension('box', 3)
+    newdim = newnc.createDimension('pt', 2)
+    newdim = newnc.createDimension('boxbox', 9)
+
+    # Var-dimensions
+    newvar = newnc.createVariable('xcurve', 'f8', ('curve'))
+    newvar[:] = curve[:,1]
+    basicvardef(newvar,'xcurve', 'x-' +Lcurvex+ ' localization of the curve', cxunit)
+
+    newvar = newnc.createVariable('ycurve', 'f8', ('curve'))
+    newvar[:] = curve[:,0]
+    basicvardef(newvar,'ycurve', 'y-' +Lcurvey+ ' localization of the curve', cyunit)
+
+    # variables
+    newvar = newnc.createVariable('weights', 'f', ('box','box','curve'))
+    for icv in range(Ncurve):
+        newvar[:,:,icv] = curvewgts[icv,:,:]
+    basicvardef(newvar,'weights', 'weights on the 3x3 around box following curve','1')
+
+    # variables
+    newvar = newnc.createVariable('Nweights', 'f', ('curve'))
+    newvar[:] = Nwgts[:]
+    basicvardef(newvar,'Nweights', 'number of usable weights of the 3x3 around box '+\
+      'following curve', '-')
+
+    # variables
+    newvar = newnc.createVariable('coordweights', 'f', ('pt','boxbox','curve'))
+    for icv in range(Ncurve):
+        newvar[:,:,icv] = curveloc3x3[icv,:,:]
+    basicvardef(newvar,'coordweights', 'coordinates of the usable weights of the ' + \
+      '3x3 around box following curve', '-')
+
+    for vn in varns:
+        if not onc.variables.has_key(vn):
+            print errormsg
+            print '  ' + fname + ": file '" + ncfile + "' does not have variable '" +\
+              vn + "' !!"
+            print '    available ones:', onc.variables.keys()
+            onc.close()
+            quit(-1)
+
+        ovar = onc.variables[vn]
+        dimvar = list(ovar.dimensions)
+        if gen.searchInlist(dimvar, xdimn): dimvar.remove(xdimn)
+        if gen.searchInlist(dimvar, ydimn): dimvar.remove(ydimn)
+
+        varcurve = interp_curve(ovar, xdimn, ydimn, loccurve, curvewgts, Nwgts,      \
+          curveloc3x3)
+
+        # Writing values of variable along the curve
+        add_dims(onc,newnc,dimvar)
+
+        newvar = newnc.createVariable(vn+'curve', 'f4', tuple(dimvar+['curve']))
+        newvar[:] = varcurve[:]
+        add_varattrs(onc,newnc,[vn],[vn+'curve'])
+        set_attribute(newvar,'along','xcurve,ycurve')
+ 
+    # Global attributes
+    add_globattrs(onc,newnc,'all')
+    newnc.setncattr('author', 'L. Fita')
+    newattr = set_attributek(newnc, 'institution', unicode('Laboratoire de M' +      \
+      unichr(233) + 't' + unichr(233) + 'orologie Dynamique'), 'U')
+    newnc.setncattr('university', 'Pierre Marie Curie - Jussieu')
+    newnc.setncattr('center', 'Centre National de Recherches Scientifiques')
+    newnc.setncattr('city', 'Paris')
+    newnc.setncattr('country', 'France')
+    newnc.setncattr('source', 'http://www.xn--llusfb-5va.cat/python/PyNCplot/')
+    newnc.setncattr('script', 'nc_var_tools.py')
+    newnc.setncattr('function', fname)
+
+    print fname + ": successfull writing of '" + ofile + "' !!"
+
+    return
+
+values = 'gridline,x,y,8.,8.,16.,16.,32'
+curve_section(values,'/home/lluis/PY/test.nc','var')
+
 #quit()