Changeset 1869 in lmdz_wrf

Timestamp:

Mar 26, 2018, 11:24:40 PM (7 years ago)

Author:

lfita

Message:

Working with `getdim_listonc'

File:

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

trunk/tools/nc_var_tools.py

@@ -77 +77 @@
 # get_1str_nc: Function to get 1 string value in a netCDF variable as a chain of 1char values
 # get_attribute: Function to get an attribute from a netCDF file
+# getdim_listonc: Function to get a dimension object from a list of netCDF file object
 # getvar_listonc: Function to get a variable object from a list of netCDF file object
 # get_namelist_vars: Function to get namelist-like  values ([varname] = [value])
@@ -22867 +22868 @@
 #quit()
 
+def getdim_listonc(dimn, fns, listoc):
+    """ Function to get a dimension object from a list of netCDF file object
+      dimn= name of the dimension to be found
+      fns= list with the name of the files
+      listoc= list of NetCDF file objects
+    """
+    fname = 'getvar_listonc'
+
+    allfs = {}
+    for ioc in range(len(listoc)):
+        oc = listoc[ioc]
+        # Seems that ObjNetcdf.filepath() requires netcdf >= 4.1.2....
+        #  Thus avoiding to use the functionality
+        #allfs[oc.filepath()] = ...
+        allfs[fns[ioc]] = [oc.dimensions, oc.variables.keys()]
+
+    # Looping for all files
+    found = False
+    for ioc in range(len(listoc)):
+        oc = listoc[ioc]
+        # Seems that ObjNetcdf.filepath() requires netcdf >= 4.1.2....
+        #  Thus avoiding to use the functionality
+        oncvars = allfs[fns[ioc]]
+        oncvars.sort()
+        if gen.searchInlist(oc.dimensions,dimn):
+            od = oc.dimensions[dimn]
+            found = True
+            break
+
+    if not found:
+        print errormsg
+        print '  ' + fname + "': any netcdf file '" + ','.join(fns) + "' contain " + \
+          "dimension '" + dimn + "' !!"
+        print '    available ones _______'
+        ioc = 0
+        for fn in fns:
+            print fn, ':', allfs[fn][0]
+            listoc[ioc].close()
+            ioc = ioc + 1
+        quit(-1)
+
+    return od
+
+def projection2D_information(PROJn, AXinf, projinf, dx, dy, Fns, oncs):
+    """ Function to get the basic information of the 2D-projection: new dimensions, variables, ...
+      Reproducing: class Projection from wrfncxnj.py
+        PROJn= name of the projection
+        AXinf= dictionary with the informatino of the axes 'X', 'Y', 'Z', 'T'
+        projinf= dictionary with all the information of the projection
+        Fns= list with the name of the files
+        oncs= list of netCDF objects from which retrieve the information
+      Returns:
+        newdim: dictionary new dimensions to be added {dimn1: dimlength1, ..., dimnN: dimlengthN} 
+        newvar: dictionary new variables to be added {varn1: [stdn1, longn1, units1], ..., varnM: [stdnM, longnM, unitsM]}
+        newvarv: dictionary values new variables to be added {varn1: [values,...], ..., varnM: [values,...]}
+        nonCFattr: list of attributes from PROJn non-CF and non to be added to the file
+    """
+    fname = 'projection2D_information'
+
+    availproj = ['Lambert_Conformal', 'Rotated_Pole/Rotate_Latituded_Longitude',     \
+      'Mercator']
+
+    # Explanation of certain attributes in projection file
+    projfileexpl = {'x_resolution': 'Resolution in file alon x axis',                \
+      'y_resolution': 'Resolution in file alon y axis',                              \
+      'proj_units': 'units of the projection (e.g.: km, m, degrees,...)',            \
+      'fileXvals': 'variable with the reference values for the projection ' +        \
+        'along X axis (e.g.: values for rlon)',                                      \
+      'fileYvals': 'variable with the reference values for the projection ' +        \
+        'along Y axis (e.g.: values for rlat)'}
+    CFprojweb = 'http://cfconventions.org/Data/cf-conventions/cf-conventions-1.7/' + \
+      'build/ch05s06.html'
+
+    isrotated = False
+    if PROJn.find('rotated') != -1 or PROJn.find('Rotated') != -1:
+        isrotated = True
+        projpiec = PROJn.split('_')
+        if PROJn.find('rotated') != -1:
+            projpiec.remove('rotated')
+        else:
+            projpiec.remove('Rotated')
+        projn = '_'.join(projpiec)
+
+    if PROJn == "Lambert_Conformal":
+        proj_attr = ["grid_mapping_name", "cone_type", "northern_parallel",          \
+          "southern_parallel", "longitude_of_central_meridian",                      \
+          "latitude_of_projection_origin", 'x_resolution', 'y_resolution',           \
+          'proj_units']
+        # Removing the non-CF attributes
+        nonCFattr = ['x_resolution', 'y_resolution', 'proj_units']
+        for atn in proj_attr:
+            if not projinf.has_key(atn):
+                print errormsg
+                print '  ' + fname + ": to process projection '" + PROJn +           \
+                  "' is required to provide argument '" + atn + "' in the " +        \
+                  "'[projectfile]' !!"
+                if projfileexpl.has_key(atn): 
+                    print '    ', atn, ':', projfileexpl[atn]
+                else:
+                    print '    ', atn, ': CF-attribute, please visit: ', CFprojweb
+                print '    values provided:', projinf
+                quit(-1)
+        newdim = {}
+        newvar = {}
+        newvarv = {}
+
+        newdim['x'] = dx
+        newvar['x'] = ['projection_x_coordinate', 'x coordinate of projection', 'km']
+        Xres = np.float(projinf['x_resolution'][0])
+        if projinf['proj_units'][0] == 'km':
+            scale = 1.
+        elif projinf['proj_units'][0] == 'm':
+            scale = 1000.
+        else:
+            print errormsg
+            print '  ' + fname + ": Resolution units '" + projinf['proj_units'][0] + \
+              "' not ready !!"
+            print '    available ones:', ['km', 'm']
+            quit(-1)
+        newvarv['x'] = (np.arange(1,dx+1)-dx/2)*Xres/scale
+        # updating values on axes
+        axv = AXinf['X']
+        axv['dimn'] = 'x'
+        axv['vdimn'] = 'lon'
+        AXinf['X'] = axv
+
+        newdim['y'] = dy
+        newvar['y'] = ['projection_y_coordinate', 'y coordinate of projection', 'km']
+        Yres = np.float(projinf['y_resolution'][0])
+        newvarv['y'] = (np.arange(1,dy+1)-dy/2)*Yres/scale
+        # updating values on axes
+        axv = AXinf['Y']
+        axv['dimn'] = 'y'
+        axv['vdimn'] = 'lat'
+        AXinf['Y'] = axv
+
+    elif PROJn == "Rotated_Pole" or PROJn == 'Rotated_Latitude_Longitude':
+        proj_attr = ["grid_mapping_name", "grid_north_pole_latitude",                \
+          "grid_north_pole_longitude", 'fileXvals', 'fileYvals']
+        # Removing the non-CF attributes
+        nonCFattr = ['fileXvals', 'fileYvals']
+        for atn in proj_attr:
+            if not projinf.has_key(atn):
+                print errormsg
+                print '  ' + fname + ": to process projection '" + PROJn +           \
+                  "' is required to provide argument '" + atn + "' in the " +        \
+                  "'[projectfile]' !!"
+                if projfileexpl.has_key(atn): 
+                    print '    ', atn, ':', projfileexpl[atn]
+                else:
+                    print '    ', atn, ': CF-attribute, please visit: ', CFprojweb
+                print '    values provided:', projinf
+                quit(-1)
+        newdim = {}
+        newvar = {}
+        newvarv = {}
+
+        newdim['rlon'] = dx
+        newvar['rlon'] = ['grid_longitude', 'longitude in rotated pole grid',        \
+          'degrees']
+
+        ovar = getvar_listonc(projinf['fileXvals'][0], Fns, oncs)
+        oXvar = ovar[0]
+        if len(oXvar.shape) == 3:
+            newvarv['rlon'] = oXvar[0,0,:]
+        elif len(oXvar.shape) == 2:
+            newvarv['rlon'] = oXvar[0,:]
+        elif len(oXvar.shape) == 1:
+            newvarv['rlon'] = oXvar[:]
+        else:
+            print errormsg
+            print '  ' + fname + ": rank ", oXvar.shape, " of values for X ref " +   \
+              "not ready !!"
+            quit(-1)
+        # updating values on axes
+        axv = AXinf['X']
+        axv['dimn'] = 'rlon'
+        axv['vdimn'] = 'lon'
+        AXinf['X'] = axv
+
+        newdim['rlat'] = dy
+        newvar['rlat'] = ['grid_latitude', 'latitude in rotated pole grid', 'degrees']
+        ovar = getvar_listonc(projinf['fileYvals'][0], Fns, oncs)
+        oYvar = ovar[0]
+        if len(oYvar.shape) == 3:
+            newvarv['rlat'] = oYvar[0,0,:]
+        elif len(oYvar.shape) == 2:
+            newvarv['rlat'] = oYvar[0,:]
+        elif len(oYvar.shape) == 1:
+            newvarv['rlat'] = oYvar[:]
+        else:
+            print errormsg
+            print '  ' + fname + ": rank ", oYvar.shape, " of values for Y ref " +   \
+              "not ready !!"
+            quit(-1)
+        # updating values on axes
+        axv = AXinf['Y']
+        axv['dimn'] = 'rlat'
+        axv['vdimn'] = 'lat'
+        AXinf['Y'] = axv
+
+    elif PROJn == "Mercator":
+        proj_attr = ["grid_mapping_name", "longitude_of_projection_origin",          \
+          "standard_parallel", 'x_resolution', 'y_resolution', 'proj_units']
+        # Removing the non-CF attributes
+        nonCFattr = ['x_resolution', 'y_resolution', 'proj_units']
+        for atn in proj_attr:
+            if not projinf.has_key(atn):
+                print errormsg
+                print '  ' + fname + ": to process projection '" + PROJn +           \
+                  "' is required to provide argument '" + atn + "' in the " +        \
+                  "'[projectfile]' !!"
+                if projfileexpl.has_key(atn): 
+                    print '    ', atn, ':', projfileexpl[atn]
+                else:
+                    print '    ', atn, ': CF-attribute, please visit: ', CFprojweb
+                print '    values provided:', projinf
+                quit(-1)
+
+        newdim = {}
+        newvar = {}
+        newvarv = {}
+
 def getvar_listonc(varn, fns, listoc):
     """ Function to get a variable object from a list of netCDF file object
@@ -23529 +23753 @@
                 else:
                     CFvardimvalues.append(dn)
-                    varslice.append(slice(0,len(onc.dimensions[dn])))
+                    oaxisd = getdim_listonc(dn, afNs, allonc)
+                    varslice.append(slice(0,len(oaxisd)))
                     if dn == dimn and CFdimvalues['length'] != -1:
-                        CFdimvalues['length'] = len(onc.dimensions[dn])
+                        CFdimvalues['length'] = len(oaxisd)
             axsiv = oaxisv[tuple(varslice)]
         else: