Changeset 814 in lmdz_wrf

Timestamp:

Jun 12, 2016, 2:12:37 PM (9 years ago)

Author:

lfita

Message:

Adding:

• `lonlatvarsFile': Functino to CF-define the longitudes and latitudes variables within a file
• `subbasins': Function to retrieve the subbasin (all the sub-flows untila a given lon,lat)

File:

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

trunk/tools/nc_var_tools.py

@@ -14852 +14852 @@
 #Partialmap_Entiremap('longitude,latitude,std,5000.,lonlat_dxdyFix', '/home/lluis/etudes/DYNAMICO/ORCHIDEE/interpolation/data/carteveg5km.nc', 'vegetation_map')
 #Partialmap_Entiremap('longitude,latitude,std,5000.,lonlat', '/home/lluis/etudes/DYNAMICO/ORCHIDEE/interpolation/data/carteveg5km.nc', 'vegetation_map')
+
+def lonlatvarsFile(lonvals,latvals,dnx,dny,oLlnc):
+    """ Functino to CF-define the longitudes and latitudes variables within a file
+      lonvals= 2D matrix with longitudes
+      latvals= 2D matrix with latitudes
+      dn[x/y]= names of the x and y dimensions
+      oLlnc= netCDF file object where to write the lon,lats
+    """
+    fname = 'lonlatvarsFile'
+
+    if len(lonvals.shape) != 2:
+        print errormsg
+        print '  ' + fname +': wrong shape:',lonvals.shape," of 'lons' variable, " + \
+          'must be 2D !!'
+        quit(-1)
+    if len(latvals.shape) != 2:
+        print errormsg
+        print '  ' + fname +': wrong shape:',latvals.shape," of 'lats' variable, " + \
+          'must be 2D !!'
+        quit(-1)
+    dx = lonvals.shape[1]
+    dy = lonvals.shape[0]
+
+    if not gen.searchInlist(oLlnc.dimensions, dnx):
+        print errormsg
+        print '  ' + fname + ': netCDF object does not have x dimension: ', dnx,"' !!"
+        quit(-1)
+    if not gen.searchInlist(oLlnc.dimensions, dny):
+        print errormsg
+        print '  ' + fname + ': netCDF object does not have y dimension: ', dny,"' !!"
+        quit(-1)
+ 
+# Creation of dimension variables
+    newvar = oLlnc.createVariable('lon', 'f8', (dny,dnx))
+    basicvardef(newvar, 'longitude', 'Longitude','degrees_East')
+    newvar[:] = lonvals
+    newvar.setncattr('axis', 'X')
+    newvar.setncattr('_CoordinateAxisType', 'Lon')
+
+    newvar = oLlnc.createVariable('lat', 'f8', (dny,dnx))
+    basicvardef(newvar, 'latitude', 'Latitude','degrees_North')
+    signlat = latvals[1,0] - latvals[0,0]
+    newvar[:] = latvals
+
+    newvar.setncattr('axis', 'Y')
+    newvar.setncattr('_CoordinateAxisType', 'Lat')   
+    
+    return
+
+def subbasin(values, ncfile):
+    """ Function to retrieve the subbasin (all the sub-flows untila a given lon,lat)
+      values=[subname],[sublon],[sublat]
+        subname: name for the subbasin ('!' for spaces)
+        sub[lon/lat]: longitude and latitude of the subbasin point 
+    """
+    fname = 'subbasin'
+    
+    if values == 'h':
+        print fname + '_____________________________________________________________'
+        print subbasin.__doc__
+        quit()
+
+    arguments = '[subnmame],[sublon],[sublat]'
+    gen.check_arguments(fname, values, arguments, ',')
+
+    ofile = 'subbasin.nc'
+
+    onc = NetCDFFile(ncfile, 'r')
+
+    subname = values.split(',')[0].replace('!',' ')
+    sublon = np.float(values.split(',')[1])
+    sublat = np.float(values.split(',')[2])
+
+    if ncfile == 'routing.nc':
+        lonname = 'nav_lon'
+        latname = 'nav_lat'
+        tripname = 'trip'
+    else:
+        print errormsg
+        print '  ' + fname + ": file '" + ncfile + "' not ready!!"
+        quit(-1)
+
+    olons = onc.variables[lonname]
+    olats = onc.variables[latname]
+    ooutflows = onc.variables[tripname]
+# Flipping the routing file
+    lons = olons[:]
+    lats = olats[:]
+    outflows = ooutflows[:]
+
+    onc.close()
+
+    difflonlat = np.sqrt((lons - sublon)**2 + (lats - sublat)**2)
+    mindifflonlat = np.min(difflonlat)
+    isub = gen.index_mat(difflonlat, mindifflonlat)
+    print '    ' + fname + "initial point of subbasin '" + subname + "' in the " +   \
+      'matrix:', isub
+
+    masksubbasin, subflows, subflowspt = gen.subbasin_point(outflows, isub)
+
+    onewnc = NetCDFFile(ofile, 'w')
+
+    dx = outflows.shape[1]
+    dy = outflows.shape[0]
+
+# Number of sub-flows within the sub-basin
+    Nsubflows = len(subflows.keys())
+
+# Maximum number of grid points of a sub-flow
+    Nmaxgridsubflow = -1
+    Ngridpoint = 0
+    for subflow in subflows.keys():
+        Ngridsubflow = len(subflows[subflow])
+        Ngridpoint = Ngridpoint + Ngridsubflow
+        if Ngridsubflow > Nmaxgridsubflow: Nmaxgridsubflow = Ngridsubflow
+
+# Creation of dimensions
+    Lstr = 50
+    newdim = onewnc.createDimension('x', dx)
+    newdim = onewnc.createDimension('y', dy)
+    newdim = onewnc.createDimension('Nsubflow', Nsubflows)
+    newdim = onewnc.createDimension('Lstring', Lstr)
+    newdim = onewnc.createDimension('Ngridsubflow', Nmaxgridsubflow)
+    
+    if ncfile == 'routing.nc':
+        lonlatvarsFile(lons, lats[::-1,:], 'x', 'y', onewnc)
+    else:
+        lonlatvarsFile(lons, lats, 'x', 'y', onewnc)
+
+# Creation of variables
+    newvar = onewnc.createVariable('subbasin', 'c', ('Lstring'))
+    for ic in range(len(subname)):
+        newvar[ic] = subname[ic:ic+1]
+    newattr = set_attributek(newvar,'Ngridpoint',Ngridpoint,'I')
+    newattr = set_attributek(newvar,'Nsubflow',Nsubflows,'I')
+    newattr = set_attributek(newvar,'lon',sublon,'R')
+    newattr = set_attributek(newvar,'lat',sublat,'R')
+    onewnc.sync()
+
+    newvar = onewnc.createVariable('subbasinmask','i',('y','x'))
+    vardef = basicvardef(newvar, 'subbasinmask', 'mask of the subbasin','1')
+    if ncfile == 'routing.nc':
+        newvar[:] = masksubbasin[::-1,:]
+    else:
+        newvar[:] = masksubbasin
+    newattr = set_attributek(newvar,'coordinates','lon lat','S')
+
+    onewnc.sync()
+
+    nsubflows = subflows.keys()
+    nsubflows.sort
+    newvar = onewnc.createVariable('subflow', 'c', ('Nsubflow', 'Lstring'))
+    vardef = basicvardef(newvar, 'subflow', 'reference of the sub-flow','-')
+    newvals = writing_str_nc(newvar, nsubflows, Lstr)
+    onewnc.sync()
+
+    newvar = onewnc.createVariable('Nsubflow', 'i', ('Nsubflow'))
+    basicvardef(newvar, 'Nsubflow', 'Number of grid points of the sub-flow', '-')
+    newvarlon = onewnc.createVariable('lonsubflow', 'f4', ('Nsubflow','Ngridsubflow'))
+    vardef = basicvardef(newvarlon, 'lonsubflow', 'Longitude of the grid point of '+ \
+      'the sub-flow','degrees_East')
+    newvarlat = onewnc.createVariable('latsubflow', 'f4', ('Nsubflow','Ngridsubflow'))
+    vardef = basicvardef(newvarlat, 'latsubflow', 'Latitude of the grid point of ' + \
+      'the sub-flow','degrees_North')
+    isub = 0
+    for subf in nsubflows:
+        grids = subflows[subf]
+        Ngrid = len(grids)
+        newvar[isub] = Ngrid
+        print '  ' + fname + 'Lluis; isub:', isub,'subf:',subf,'grids:',grids,'Ngrid',Ngrid
+        for igrid in range(Ngrid):
+            jigrid = subflowspt[grids[igrid]]
+            lonv = lons[jigrid[0],jigrid[1]]
+            latv = lats[jigrid[0],jigrid[1]]
+            newvarlon[igrid,isub] = lonv
+            newvarlat[igrid,isub] = latv
+        isub = isub + 1
+    onewnc.sync()
+    onewnc.close()
+
+    print fname + ": succesfull written of '" + ofile + "' !!"
+
+    return
+
 #quit()
 
@@ -14872 +15056 @@
        'missing': as a missing value
   addvattr, add a new attribute to any given variable: addvattr -S [attrname]|[attrvalue]
-  addvattrk, add a new attribute to any given variable: addvattrk -S [attrname]|[attrvalue]|[kind(S (!, white spaces),I,R,D)]
+  addvattrk, add a new attribute to any given variable: addvattr
+k -S [attrname]|[attrvalue]|[kind(S (!, white spaces),I,R,D)]
   checkNaNs, checks for NaN values over all variables in a file
   checkAllValues, check for variables with along all their dimensions with the same value in a file