Changeset 667 in lmdz_wrf for trunk/tools/nc_var_tools.py

Timestamp:

Jan 14, 2016, 10:45:54 AM (9 years ago)

Author:

lfita

Message:

Changing coincidence with 'missing' value to be with 'checking' value in `VarVal_FillValue'

File:

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

trunk/tools/nc_var_tools.py

@@ -13 +13 @@
 fillValueC = '-'
 fillValueI = -99999
+fillValueI16 = -99999
+fillValueF64 = 1.e20
+fillValueI32 = -99999
 
 def values_line(line, splitv, chars):
@@ -9875 +9878 @@
           ...
           [dim1]_N
+          NOTE: for character values, use '!' for space
         ncfile= netCDF file to use
         varname= variable to use
@@ -9917 +9921 @@
         for line in objasciif:
             val = (line.replace('\n','').replace(' ','').replace('\t',''))
-            valfinal[0] = retype(val, vartype)
-            if it <= values.shape[0]-1:
-                values[it] = retype(valfinal[0], vartype)
-                it = it +1
+            if vartype != '|S1':
+                valfinal[0] = retype(val, vartype)
+                if it <= values.shape[0]-1:
+                    values[it] = retype(valfinal[0], vartype)
+                    it = it +1
+            else:
+                for i1 in range(len(val)):
+                    values[it] = val.replace('!',' ')
     elif Ndims == 2:
         iline=0
         for line in objasciif:
             vals = line.replace('\n','').replace('\t','').split(' ')
-            if len(vals) != values.shape[1]:
-                print errormsg
-                print '  ' + fname + ': given: ',len(vals),' but variable reaquires: ',values.shape[1],'!!!'
-                exit(-1)
-            for i1 in range(len(vals)):
-                values[iline,i1] = retype(vals[i1], vartype)
+            if vartype != '|S1':
+                if len(vals) != values.shape[1]:
+                    print errormsg
+                    print '  ' + fname + ': given: ',len(vals),' but variable ' +    \
+                      'reaquires: ',values.shape[1],'!!!'
+                    exit(-1)
+                for i1 in range(len(vals)):
+                    values[iline,i1] = retype(vals[i1], vartype)
+            else:
+                for i1 in range(len(vals)):
+                    values[iline,i1] = vals[i1].replace('!',' ')
             iline=iline+1
     else:
@@ -14469 +14482 @@
           if [dsize] = 'None', give a third value with the real size
         [attributes]: [std_name]@[long_name]@[units], standard name, long name and units of the variable
+          ('!' for spaces)
         [kind]: type of variable (standard netCDF4/C-like values, 'c', 'i', 'f', 'f8',...)
       ncfile= name of the file
@@ -14485 +14499 @@
 
     dimensions = values.split('|')[0].split(',')
-    attributes = values.split('|')[1]
+    attributes = values.split('|')[1].replace('!',' ')
     kind = values.split('|')[2]
 
@@ -14521 +14535 @@
         newvar = onc.createVariable(varn, 'c', tuple(dnames))
 #        newvar[:] = np.zeros(tuple(dsize), dtype=np.float)
-    if kind == 'f' or kind == 'f4':
+    elif kind == 'f' or kind == 'f4':
         newvar = onc.createVariable(varn, 'f4', tuple(dnames), fill_value=fillValue)
         newvar[:] = np.zeros(tuple(dsize), dtype=np.float)
@@ -17158 +17172 @@
               'inv': inverse of the values 1/[varname]
               'nothing': direct values
+              'prodinv': inverse of the product between 2 variables ([varname1]*[varname2]), 
+                 NOTE: [varname] = [varname1]:[varname2]
           * 'reglonlat',[lonname],[latname]: regular lon/lat projection to compute the area size
         [xdimname]: name of the dimension for the x-axis
@@ -17171 +17187 @@
         quit()
 
-    operations = ['inv','nothing']
+    operations = ['inv','nothing','prodinv']
 
     weightk = values.split(',')[0]
@@ -17229 +17245 @@
 
     if weightk == 'varnames':
-        if not onc.variables.has_key(varname):
-            print errormsg
-            print '  ' + fname + ": file '" + filen + "' does not have variable '" + \
-              varname + "' !!"
-            quit(-1)
-        iovarwgt = onc.variables[varname]
-
-        ivarwgtv = iovarwgt[:]
+        if oper == 'inv':
+            if not onc.variables.has_key(varname):
+                print errormsg
+                print '  ' + fname + ": file '" + filen + "' does not have " +       \
+                  "variable '" + varname + "' !!"
+                quit(-1)
+            iovarwgt = onc.variables[varname]
+            ivarwgtv = iovarwgt[:]
+        elif oper == 'nothing':
+            if not onc.variables.has_key(varname):
+                print errormsg
+                print '  ' + fname + ": file '" + filen + "' does not have " +       \
+                  "variable '" + varname + "' !!"
+                quit(-1)
+            iovarwgt = onc.variables[varname]
+            ivarwgtv = iovarwgt[:]
+        elif oper == 'prodinv':
+            varname1 = varname.split(':')[0]
+            varname2 = varname.split(':')[1]
+            if not onc.variables.has_key(varname1):
+                print errormsg
+                print '  ' + fname + ": file '" + filen + "' does not have " +       \
+                  "variable '" + varname1 + "' !!"
+                quit(-1)
+            if not onc.variables.has_key(varname2):
+                print errormsg
+                print '  ' + fname + ": file '" + filen + "' does not have " +       \
+                  "variable '" + varname2 + "' !!"
+                quit(-1)
+            iovarwgt = onc.variables[varname1]
+            ivarwgtv = iovarwgt[:]*onc.variables[varname2]
+
         if oper == 'inv':
             longvarname = 'using inverse of variable ' + varname + ' as space weights'
@@ -17265 +17305 @@
                             newvals[id1,id2,id3]= np.mean(ivarv[tuple(slicevalues)]/ \
                               ivarwgtv[tuple(slicewgt)])
+            outweightvals = 1./ivarwgtv[tuple(slicewgt)]
         elif oper == 'nothing':
             longvarname = 'using variable ' + varname + ' as space weights'
@@ -17293 +17334 @@
                             newvals[id1,id2,id3]=np.mean(ivarv[tuple(slicevalues)]*  \
                               ivarwgtv[tuple(slicewgt)])
+            outweightvals = ivarwgtv[tuple(slicewgt)]
+        if oper == 'prodinv':
+            longvarname = 'using inverse of product of variables ' + varname1 + '*'  \
+              + varname2 + ' as space weights'
+            if len(loopshape) == 1:
+                newvals = np.ones((loopshape[0]), dtype=np.float)*fillValueF
+                for id1 in range(loopshape[0]):
+                    slicevalues = SliceVar(iovar,dimsloop,[id1])
+                    slicewgt = SliceVar(iovarwgt,dimsloop,[id1])
+                    newvals[id1] = np.mean(ivarv[tuple(slicevalues)] /               \
+                      ivarwgtv[tuple(slicewgt)])
+            elif len(loopshape) == 2:
+                newvals = np.ones((loopshape[0],loopshape[1]),dtype=np.float)*       \
+                  fillValueF
+                for id1 in range(loopshape[0]):
+                    for id2 in range(loopshape[1]):
+                        slicevalues = SliceVar(iovar,dimsloop,[id1,id2])
+                        slicewgt = SliceVar(iovarwgt,dimsloop,[id1,id2])
+                        newvals[id1,id2] = np.mean(ivarv[tuple(slicevalues)] /       \
+                          ivarwgtv[tuple(slicewgt)])
+            elif len(loopshape) == 3:
+                newvals = np.ones((loopshape[0],loopshape[1],loopshape[2]),          \
+                  dtype=np.float)*fillValueF
+                for id1 in range(loopshape[0]):
+                    for id2 in range(loopshape[1]):
+                        for id3 in range(loopshape[1]):
+                            slicevalues = SliceVar(iovar,dimsloop,[id1,id2,id3])
+                            slicewgt = SliceVar(iovarwgt,dimsloop,[id1,id2,id3])
+                            newvals[id1,id2,id3]= np.mean(ivarv[tuple(slicevalues)]/ \
+                              ivarwgtv[tuple(slicewgt)])
+            outweightvals = 1./ivarwgtv[tuple(slicewgt)]
     elif weightk == 'reglonlat':
         longvarname = 'using cos(latitude) variable ' + latname + ' as space weights'
@@ -17313 +17385 @@
         lonv, latv = lonlat2D(ilonv, ilatv)
 
-        ivarwgtv = np.abs(np.cos(latv*np.pi/180.))
+        ivarwgtv = np.abs(np.cos(latv*np.pi/180.))/(lonv.shape[0]*lonv.shape[1])
         if len(loopshape) == 1:
             newvals = np.ones((loopshape[0]), dtype=np.float)*fillValueF
@@ -17334 +17406 @@
                         newvals[id1,id2,id3] = np.mean(ivarv[tuple(slicevalues)] *   \
                           ivarwgtv)
+        outweightvals = ivarwgtv
 
 # Writting output file
@@ -17351 +17424 @@
       longvarname, iovar.getncattr('units'))
     newvar[:] = newvals
+
+# Spatial weight
+##
+    newdim = onewnc.createDimension(ydimname, outweightvals.shape[0])
+    newdim = onewnc.createDimension(xdimname, outweightvals.shape[1])
+    newvar = onewnc.createVariable('spatialweight', 'f4',tuple([ydimname, xdimname]),\
+      fill_value=fillValueF)
+    basicvardef(newvar, 'spatialweight', 'space weight ' + longvarname, '-')
+    newvar[:] = outweightvals
 
 # Global attributes
@@ -17372 +17454 @@
 #SpatialWeightedMean('varnames,MAPFAC_M,inv,west_east,south_north','wrfout_d01_2001-11-11_00:00:00','HGT')
 #SpatialWeightedMean('reglonlat,XLONG,XLAT,west_east,south_north','wrfout_d01_2001-11-11_00:00:00','HGT')
+
+def fill_value(vartype, fillval0):
+    """ Function to provide the fiven fill_Value for a kind of variable
+      [vartype]= type of the variable
+      [fillval0]= value to take as fill value, 'std' for the standard value
+           Float = 1.e20
+           Character = '-'
+           Integer = -99999
+           Integer16 = -99999
+           Float64 = 1.e20
+           Integer32 = -99999
+    """
+    fname = 'fill_value'
+
+    if vartype == '|S':
+        if fillval0 == 'std':
+            fillval = fillvalueC
+        else:
+            fillval = fillval0
+    elif vartype == type(int(1)):
+        if fillval0 == 'std':
+            fillval = fillValueI
+        else:
+            fillval = int(fillval0)
+    elif vartype == type(np.int16(1)):
+        if fillval0 == 'std':
+            fillval = np.int16(fillValueI)
+        else:
+            fillval = np.int16(fillval0)
+    elif vartype == type(np.int32(1)):
+        if fillval0 == 'std':
+            fillval = np.int32(fillValueI)
+        else:
+            fillval = np.int32(fillval0)
+    elif vartype == type(np.float(1.)):
+        if fillval0 == 'std':
+            fillval = np.float(fillValueF)
+        else:
+            fillval = np.float(fillval0)
+    elif vartype == type(np.float32(1.)):
+        if fillval0 == 'std':
+            fillval = np.float32(fillValueF)
+        else:
+            fillval = np.float32(fillval0)
+    elif vartype == type(np.float64(1.)):
+        if fillval0 == 'std':
+            fillval = np.float64(fillValueF)
+        else:
+            fillval = np.float64(fillval0)
+    else:
+        print errormsg
+        print '  ' + fname + ': variable type', vartype, 'not ready !!'
+        quit(-1)
+
+    return fillval
+
+def VarVal_FillValue(values, filen, varn):
+    """ Function to transform a given value from a given variable to _FillValue in a netCDF file
+      values= [value],[fillVal]
+        [value]: value to pass to '_FillValue'
+        [fillVal]: value for '_FillValue', 'std' for the standard value
+           Float = 1.e20
+           Character = '-'
+           Integer = -99999
+           Float64 = np.float(Float)
+           Integer32 = np.int32(Integer)
+      filen= name of the netCDF file
+      varn= name of the variable
+    """
+    fname = 'VarVal_FillValue'
+
+    if values == 'h':
+        print fname + '_____________________________________________________________'
+        print VarVal_FillValue.__doc__
+        quit()
+
+    onc = NetCDFFile(filen, 'a')
+    if not onc.variables.has_key(varn):
+        print errormsg
+        print '  ' + fname + ": file '" + filen + "' does not have " +       \
+          "variable '" + varn + "' !!"
+        quit(-1)
+
+    value = values.split(',')[0]
+    fval0 = values.split(',')[1]
+
+    varo = onc.variables[varn]
+    varvals = varo[:]
+    vartype = varo.dtype
+    
+    if vartype == '|S':
+        valchk = value
+    elif vartype == type(int(1)):
+        valchk = int(value)
+    elif vartype == type(np.int32(1)):
+        valchk = np.int32(value)
+    elif vartype == type(np.float(1.)):
+        valchk = np.float(value)
+    elif vartype == type(np.float32(1.)):
+        valchk = np.float32(value)
+    elif vartype == type(np.float64(1.)):
+        valchk = np.float64(value)
+    else:
+        print errormsg
+        print '  ' + fname + ': variable type', vartype, 'not ready !!'
+        quit(-1)
+
+    fval = fill_value(vartype, fval0)
+    newvarvals = np.where(varvals==valchk, fval, varvals)
+
+    if searchInlist(varo.ncattrs(), 'missing_value'):
+        missval = varo.getncattr('missing_value')
+        difvals = fval*10.**(-np.log10(valchk)) - missval*10.**(-np.log10(missval))
+        if missval == valchk or np.abs(difvals) < 1.e-6:
+            print warnmsg
+            print '  ' + fname + ': same missing and checking Value!'
+            print '   renaming missing_value'
+            varo.renameAttribute('missing_value','_FillValue')
+        else:
+            set_attribute(varo, '_FillValue', fval)
+    else:
+        set_attribute(varo, '_FillValue', fval)
+
+    varo[:] = newvarvals
+
+    onc.sync()
+    onc.close()
+
+    print fname +': Successfull change to _FillValue !!'
+
+    return
+
+#VarVal_FillValue('0.,std', '/home/lluis/etudes/DYNAMICO/ORCHIDEE/interpolation/data/new_lai2D.nc', 'LAI')
+
+def Partialmap_Entiremap(values, filen, varn):
+    """ Function to transform from a partial global map (e.g.: only land points) to an entire one
+      values= [lonmame],[latname]
+        [lonname]: name of the longitude variable
+        [latname]: name of the latitude variable
+        [fillVal]: value for '_FillValue', 'std' for the standard value
+           Float = 1.e20
+           Character = '-'
+           Integer = -99999
+           Float64 = 1.e20
+           Integer32 = -99999
+      filen= name of the netCDF file
+      varn= name of the variable
+    """
+    import numpy.ma as ma
+    fname = 'Partialmap_Entiremap'
+
+    if values == 'h':
+        print fname + '_____________________________________________________________'
+        print Partialmap_Entiremap.__doc__
+        quit()
+
+    ofile = 'EntireGlobalMap.nc'
+
+    onc = NetCDFFile(filen, 'r')
+    if not onc.variables.has_key(varn):
+        print errormsg
+        print '  ' + fname + ": file '" + filen + "' does not have " +       \
+          "variable '" + varn + "' !!"
+        quit(-1)
+
+    lonname = values.split(',')[0]
+    latname = values.split(',')[1]
+    fval0 = values.split(',')[2]
+
+    if not onc.variables.has_key(lonname):
+        print errormsg
+        print '  ' + fname + ": file '" + filen + "' does not have " +       \
+          "longitude '" + lonname + "' !!"
+        quit(-1)
+    if not onc.variables.has_key(latname):
+        print errormsg
+        print '  ' + fname + ": file '" + filen + "' does not have " +       \
+          "latitude '" + latname + "' !!"
+        quit(-1)
+
+    olon = onc.variables[lonname]
+    olat = onc.variables[latname]
+
+    lonvs = olon[:]
+    latvs = olat[:]
+
+    nlat = np.min(latvs)
+    nlon = np.min(lonvs)
+
+    minlat = np.min(np.abs(latvs))
+    minlon = np.min(np.abs(lonvs))
+
+    print '  ' + fname + ': Closest latitude to Equator:', minlat
+    print '  ' + fname + ': Closest longitude to Greenwich Meridian:', minlon
+
+# Existing longitudes along/closest to the Equator
+    eqlons = []
+    for i in range(len(latvs)):
+        if latvs[i] == minlat: eqlons.append(lonvs[i])
+
+# Existing latitudes along/closest to the Greenwich Meridian
+    grlats = []
+    for i in range(len(lonvs)):
+        if lonvs[i] == minlon: grlats.append(latvs[i])
+
+    sorteqlons = np.sort(eqlons)
+    sortgrlats = np.sort(grlats)
+ 
+    Neqlons = len(sorteqlons)
+    Ngrlats = len(sortgrlats)
+
+    if Neqlons > 1:
+        diffeqlons = sorteqlons[1:Neqlons-1] - sorteqlons[0:Neqlons-2]     
+        mindeqlon = np.min(diffeqlons)
+        print 'N sorteqlons:',Neqlons,'min deqlon:', mindeqlon
+    else:
+        mindeqlon = None
+
+    if Ngrlats > 1:
+        mindgrlat = np.min(diffgrlats)
+        diffgrlats = sortgrlats[1:Ngrlats-1] - sortgrlats[0:Ngrlats-2]
+        print 'N sortgrlats:',Ngrlats,'min dgrlat:', mindgrlat
+        latmap = np.range(0,360.,360./mindeqlon)
+    else:
+        mindgrlat = None
+
+# Global sorting
+
+#    TOTlons = lonvs.shape[0]
+#    TOTlats = lonvs.shape[0]
+
+#    sortlons = np.sort(lonvs)
+#    sortlats = np.sort(latvs)
+
+#    difflons = sortlons[1:TOTlons-1] - sortlons[0:TOTlons-2]
+#    difflats = sortlats[1:TOTlons-1] - sortlats[0:TOTlons-2]
+
+#    madifflons = ma.masked_equal(difflons, 0.)
+#    madifflats = ma.masked_equal(difflats, 0.)
+
+#    print 'Nmadifflons:',len(madifflons.compressed()),'Nmadifflats:',len(madifflats.compressed())
+
+#    sortdifflons = np.sort(madifflons)
+#    sortdifflats = np.sort(madifflats)
+
+#    print 'sortdifflons:',sortdifflons[0:10]
+#    print 'sortdifflats:',sortdifflats[0:10]
+
+# Fixing in case it has not worked
+    if mindeqlon is not None and mindgrlat is None: mindgrlat = mindeqlon
+    if mindeqlon is None and mindgrlat is not None: mindeqlon = mindgrlat
+    if mindeqlon is None and mindgrlat is None:
+        print errormsg
+        print fname + ': Not enough values along the Equator and Greenwich!!'
+        quit(-1)
+
+    if nlon < 0.:
+        lonmap = np.arange(-180.+mindeqlon,180.,mindeqlon)
+    else:
+        lonmap = np.arange(0.,360.-mindeqlon,mindeqlon)
+
+    if nlat < 0.:
+        latmap = np.arange(-90.+mindgrlat,90.,mindgrlat)
+    else:
+        latmap = np.arange(0.,180.-mindgrlat,mindgrlat)
+
+    print '  ' + fname + '_________'
+    print '     final longitudes:',lonmap
+    print '     final latitudes:',latmap
+
+    Nmaplon = len(lonmap)
+    Nmaplat = len(latmap)
+    print '  ' + fname + ': desired map will be:',Nmaplon,'x',Nmaplat
+
+    quit()
+
+# Matrix map creation
+##
+    ovar = onc.variables[varn]
+    vartype = ovar.dtype
+    varlongname = ovar.long_name
+    varunits = ovar.units
+
+    fval = fill_value(vartype, fval0)
+
+    mapval = np.ones((Nmaplat,Nmaplon), dtype=np.float)*fval
+
+    Ntotvals = lonvs.shape[0]
+    for iv in range(Ntotvals):
+        if np.mod(iv*100./Ntotvals, 10.) == 0: print '    done:',iv*100./Ntotvals,'%'
+        difflon = np.abs(lonmap - lonvs[iv])
+        difflat = np.abs(latmap - latvs[iv])
+
+#        print np.min(difflon), ':', np.min(difflat)
+
+        ilon = index_mat(difflon, np.min(difflon))
+        ilat = index_mat(difflat, np.min(difflat))
+        if np.min(difflon) != 0. or np.min(difflat) !=0.:
+            print 'No zero diffs!!:', lonvs[iv],latvs[iv],':',np.min(difflon),np.min(difflat),'i j',ilon,ilat
+            print lonmap[ilon], latmap[ilat]
+            quit()
+
+        if ilon > 0 and ilat > 0:
+            mapval[ilat,ilon] = ovar[iv]
+        else:
+            print errormsg
+            print '  ' + fname + ': point:',lonvs[iv],',',latvs[iv],' not relocated !!'
+            quit(-1)
+
+# Final map creation
+##
+    newnc = NetCDFFile(ofile, 'w')
+
+# dimensions
+    newdim = newnc.createDimension('lon',Nmaplon)
+    newdim = newnc.createDimension('lat',Nmaplat)
+
+# Variables
+    newvar = newnc.createVariable('lon','f8',('lon'))
+    basicvardef(newvar, 'lon', 'Longitudes','degrees East')
+    nevar[:] = lonmap
+
+    newvar = newnc.createVariable('lat','f8',('lat'))
+    basicvardef(newvar, 'lat', 'Latitudes','degrees North')
+    nevar[:] = latmap
+
+# map variable
+    newvar = newnc.createVariable(varn,'f4',('lat','lon'), fill_value=fval)
+    basicvardef(newvar, varn, 'Longitudes','degrees East')
+    nevar[:] = lonmap
+
+    newnc.sync()
+# Global attributes
+##
+    newnc.setncattr('script',  fname)
+    newnc.setncattr('version',  '1.0')
+    newnc.setncattr('author',  'L. Fita')
+    newnc.setncattr('institution',  'Laboratoire de Meteorology Dynamique')
+    newnc.setncattr('university',  'Pierre et Marie Curie')
+    newnc.setncattr('country',  'France')
+    for attrs in onc.ncattrs():
+        attrv = onc.getncattr(attrs)
+        attr = set_attribute(newnc, attrs, attrv)
+
+    newnc.sync()
+    onc.close()
+    newnc.close()
+
+    print fname + ": Successfull written of file: '" + ofile + "' !!"
+
+    return 
+
+#Partialmap_Entiremap('longitude,latitude,std', '/home/lluis/etudes/DYNAMICO/ORCHIDEE/interpolation/data/carteveg5km.nc', 'vegetation_map')
 #quit()