Changeset 1861 in lmdz_wrf for trunk

Timestamp:

Mar 26, 2018, 3:53:59 PM (7 years ago)

Author:

lfita

Message:

Adding:

• `projection2D_information': Function to get the basic information of the 2D-projection: new dimensions, variables, ...

File:

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

trunk/tools/nc_var_tools.py

@@ -114 +114 @@
 # Partialmap_EntiremapForExact: Function to transform from a partial global map (e.g.: only land points) to an entire one using Fortran code with exact location
 # pinterp: Function to vertically interpolate using subroutines from the p_interp.F90 NCAR program
+# projection2D_information: Function to get the basic information of the 2D-projection: new dimensions, variables, ...
 # put_variable_slice: Function to add a variable from a netcdf object to a new one following a slice
 # reconstruct_matrix_from_vector: Function to reconstruct a 2D matrix from a pair of syncronized vectors with the positions on x and y coordinates
@@ -475 +476 @@
                 var[i,:,:] = gen.valmodoper(varvals, values)
         elif Nvars == 4:
-            for i in range(varshape[0]):
+            for i in range(varshape[0]):
                 for j in range(varshape[1]):
                     varvals[:] = var[i,j,:,:]
                     var[i,j,:,:] = gen.valmodoper(varvals, values)
         elif Nvars == 5:
-            for i in range(varshape[0]):
+            for i in range(varshape[0]):
                 for j in range(varshape[1]):
                     for k in range(varshape[2]):
@@ -486 +487 @@
                         var[i,j,k,:,:] = gen.valmodoper(varvals, values)
         elif Nvars == 6:
-            for i in range(varshape[0]):
+            for i in range(varshape[0]):
                 for j in range(varshape[1]):
                     for k in range(varshape[2]):
@@ -22864 +22865 @@
 #quit()
 
+def projection2D_information(PROJn, AXinf, projinf, dx, dy, onc):
+    """ 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
+        onc= netCDF object 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', '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":
+        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']
+        if not gen.searchInlist(onc.variables.keys(),projinf['fileXrefvals'][0]):
+            print errormsg
+            print '  ' + fname + ": file does not have variable with X ref values '"+\
+              projinf['fileXrefvals'][0] + "' !!"
+            print '    available ones:', onc.variables.keys()
+            quit(-1)
+        oXvar = onc.variables[projinf['fileXvals'][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']
+        if not gen.searchInlist(onc.variables.keys(),projinf['fileYrefvals'][0]):
+            print errormsg
+            print '  ' + fname + ": file does not have variable with Y ref values '"+\
+              projinf['fileYrefvals'][0] + "' !!"
+            print '    available ones:', onc.variables.keys()
+            quit(-1)
+        oYvar = onc.variables[projinf['fileYvals'][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 = {}
+
+        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
+
+    else:
+        print errormsg
+        print '  ' + fname + ": projection '" + + "' not ready !!"
+        print '    available ones:', availproj
+        quit(-1)
+
+    return newdim, newvar, newvarv, nonCFattr
+
 def CFmorzization(values, ncfile, variable):
     """ Function to provide a CF-compilation version of a variable within a file
@@ -23190 +23416 @@
         CFaxisvardimvals[axn] = CFvardimvalues
 
-    # Renaming dimensions for 2D lon,lat variable-dimensions
+    # Getting further projection information
     if len(filedimvals['X'].shape) == 2 or len(filedimvals['Y'].shape) == 2:
-        print infmsg
-        print '    because Xvals have rank 2:', filedimvals['X'].shape,              \
-          'renaming lon,lat dimensions to x,y'
+        newDim, newVar, newVARv, nonCFproj = projection2D_information(gattr['grid'], \
+          CFdimvals, pattr, filedimvals['X'].shape[1], filedimvals['X'].shape[0], onc)
+
+    if len(newDim.keys()) > 0:
+        print '    Dealing with 2D lon,lat dimensions ...'
         addCFdimvals = {}
-        axv = CFdimvals['X']
-        axv['dimn'] = 'x'
-        axv['vdimn'] = 'lon'
-        CFdimvals['X'] = axv
-        addCFdimvals['x'] = {'dimn': 'x', 'vdimn': 'x',
-          'stdn': 'projection_x_coordinate','longname': 'x coordinate of projection',\
-          'units': pattr['proj_units'][0], 'maxrank': 1,                             \
-          'length': filedimvals['X'].shape[1]}
-        Xres = np.float(pattr['x_resolution'][0])
-        filedimvals['x'] = np.arange((filedimvals['X'].shape[1]),dtype=np.float)*Xres
-        axv = CFdimvals['Y']
-        axv['dimn'] = 'y'
-        axv['vdimn'] = 'lat'
-        CFdimvals['Y'] = axv
-        addCFdimvals['y'] = {'dimn': 'y', 'vdimn': 'y',
-          'stdn': 'projection_y_coordinate','longname': 'y coordinate of projection',\
-          'units': pattr['proj_units'][0], 'maxrank': 1,                             \
-          'length': filedimvals['X'].shape[0]}
-        Yres = np.float(pattr['y_resolution'][0])
-        filedimvals['y'] = np.arange((filedimvals['X'].shape[0]),dtype=np.float)*Yres
-        #axiv = axisinf['X']
-        #filecfdimn[axiv[0]] = 'x'
-        #axiv = axisinf['Y']
-        #filecfdimn[axiv[0]] = 'y'
+        for newdn in newDim.keys():
+            addCFdimvals[newdn] = {'dimn': newdn, 'vdimn': newdn,                    \
+              'stdn': newVar[newdn][0], 'longname': newVar[newdn][1],                \
+              'units': newVar[newdn][2], 'maxrank': 1,                               \
+              'length': filedimvals['X'].shape[1]}
+            filedimvals[newdn] = newVARv[newdn]
 
     print 'Values for variables-axes from file ________'
@@ -23456 +23666 @@
         newvar = onewnc.createVariable(gattr['grid'], 'i')
         for attrn in ppn:
-            pvattr = pattr[attrn]
-            set_attributek(newvar, attrn, pvattr[0], pvattr[1])
-            onewnc.sync()
+            if not gen.searchInlist(nonCFproj, attrn):
+                pvattr = pattr[attrn]
+                set_attributek(newvar, attrn, pvattr[0], pvattr[1])
+                onewnc.sync()
 
         # Global attributes