Changeset 826 in lmdz_wrf for trunk/tools

Timestamp:

Jun 14, 2016, 7:24:36 PM (9 years ago)

Author:

lfita

Message:

Almost working `draw_subbasin'

Location:

trunk/tools

Files:

• drawing.py (modified) (4 diffs)
• drawing_tools.py (modified) (3 diffs)
• generic_tools.py (modified) (5 diffs)

trunk/tools/drawing.py

@@ -31 +31 @@
 ## e.g. # drawing.py -o draw_river_desc -f diver_desc.nc -S 'Y|lat|lat|-1,X|lon|lon|-1:red,green:Blues:cyl,l:ORCDHIEE rivers:pdf:0:or_rivers -v Amazon
 ## e.g. # drawing.py -o draw_vertical_levels -f wrfout_d01_2001-11-11_00:00:00 -S 'true:false:wrfout!vertical!levels!(standard!40):png:4' -v WRFz
+## e.g. $ drawing.py -o draw_subbasin -f Caceres_subbasin.nc -S 'Caceres:None:cyl,l:True:Caceres:pdf:0:Caceres_subbasin'
 
 main = 'drawing.py'
@@ -43 +44 @@
   'draw_lines', 'draw_lines_time', 'draw_Neighbourghood_evol',                       \
   'draw_points', 'draw_points_lonlat',                                               \
-  'draw_ptZvals', 'draw_river_desc', 'draw_timeSeries',                              \
+  'draw_ptZvals', 'draw_river_desc', 'draw_subbasin', 'draw_timeSeries',             \
   'draw_topo_geogrid',                                                               \
   'draw_topo_geogrid_boxes', 'draw_trajectories', 'draw_vals_trajectories',          \
@@ -3984 +3985 @@
     return
 
+def draw_subbasin(ncfile, values):
+    """ Function to plot subbasin from 'routnig.nc' ORCDHIEE
+      ncfile= file to use produced with nc_var.py#subbasin function 
+      values= [subasiname]:[rangecolors]:[mapv]:[drawsubid]:[gtit]:[figkind]:[legloc]:[figurename]
+        [subasiname]= name of the subbasin ('!' for spaces)
+        [rcolor]= '@', list of 'r|g|b' 1-based colors (as much as first level sub-flow). 'None' for automatic
+        [mapv]= map characteristics: [proj],[res]
+          see full documentation: http://matplotlib.org/basemap/
+          [proj]: projection
+            * 'cyl', cilindric
+            * 'lcc', lambert conformal
+          [res]: resolution:
+            * 'c', crude
+            * 'l', low
+            * 'i', intermediate
+            * 'h', high
+            * 'f', full
+        [drawsubid]= wehther sub-flow ids should be plot or not
+        [graphtit]= title of the graph ('|', for spaces)
+        [lloc]= location of the legend (0, automatic)
+          1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
+          5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
+          9: 'upper center', 10: 'center'      kfig= kind of figure
+        [figname]= name of the figure
+    """
+    fname = 'draw_subbasin'
+
+    if values == 'h':
+        print fname + '_____________________________________________________________'
+        print draw_subbasin.__doc__
+        quit()
+
+    expectargs = '[subasiname]:[rangecolors]:[mapv]:[drawsubid]:[gtit]:[figkind]:' + \
+      '[legloc]:[figurename]'
+ 
+    drw.check_arguments(fname,values,expectargs,':')
+
+    subbasiname = values.split(':')[0].replace('!',' ')
+    rangecolors = values.split(':')[1]
+    mapv = values.split(':')[2]
+    drawsubid = gen.Str_Bool(values.split(':')[3])
+    gtit = values.split(':')[4].replace('!',' ')
+    figkind = values.split(':')[5]
+    legloc = int(values.split(':')[6])
+    figurename = values.split(':')[7]
+
+    if not os.path.isfile(ncfile):
+        print errormsg
+        print '  ' + fname + ': file "' + ncfile + '" does not exist !!'
+        quit(-1)    
+
+    objf = NetCDFFile(ncfile, 'r')
+
+    searchvars = ['lon', 'lat', 'lonsubflow', 'latsubflow', 'outsubflow']
+    for searchvar in searchvars: 
+        if not gen.searchInlist(objf.variables,searchvar):
+            print errormsg
+            print '  ' + fname + ": WRF file '" + ncfile + "' does not have " +      \
+              "variable '" + searchvar + "' !!"
+            quit(-1)
+    
+# lon,lat
+    olon = objf.variables['lon']
+    olat = objf.variables['lat']
+    lon = olon[:]
+    lat = olat[:]
+
+# sub-flow names
+    osubnames = objf.variables['subflow']
+    subnames = drw.get_str_nc(osubnames, osubnames.shape[1])
+
+# sub-flow lat, lon
+    latlonsub = {}
+    outflowsub = {}
+    osublon = objf.variables['lonsubflow']
+    osublat = objf.variables['latsubflow']
+    oNsubflow = objf.variables['Nsubflow']
+    ooutsubflow = objf.variables['outsubflow']
+    Nsubflow = oNsubflow[:]
+    isub = 0
+    for Ssub in subnames:
+        sublatlon = []
+        suboutflow = []
+        for igrid in range(Nsubflow[isub]):
+            sublatlon.append([osublat[isub,igrid], osublon[isub,igrid]])
+            suboutflow.append(ooutsubflow[isub,igrid])
+        latlonsub[Ssub] = sublatlon
+        outflowsub[Ssub] = suboutflow
+        isub = isub + 1
+
+# colors
+    if rangecolors == 'None':
+        rangecols = None
+    else:
+        cols = rangecolors.split('@')
+        Ncols = len(cols)
+        rangecols = []
+        for icol in range(Ncols):
+            cval = cols[icol].split('|')
+            rangecols.append([np.float(cval[0]),np.float(cval[1]),np.float(cval[2])])
+
+    drw.plot_subbasin(subbasiname, lon, lat, subnames, latlonsub, outflowsub,        \
+      rangecols, mapv, drawsubid, gtit, figkind, legloc, figurename)
+
+    objf.close()
+
+    return
+
 #quit()
 
@@ -4070 +4179 @@
 elif oper == 'draw_river_desc':
     draw_river_desc(opts.ncfile, opts.values, opts.varname)
+elif oper == 'draw_subbasin':
+    draw_subbasin(opts.ncfile, opts.values)
 elif oper == 'draw_timeSeries':
     draw_timeSeries(opts.ncfile, opts.values, opts.varname)

trunk/tools/drawing_tools.py

@@ -278 +278 @@
 
     return varvalues, dimnslice
+
+def get_str_nc(varo, Lchar):
+    """ Function to get string values in a netCDF variable as a chains of 1char values
+    varo= netCDF variable object
+    Lchar = length of the string in the netCDF file
+    """
+    fname = 'get_str_nc'
+
+    Nvalues = varo.shape[0]
+
+    strings = []
+
+    for ival in range(Nvalues):
+        stringv = varo[ival,:]
+        string = str('')
+        for ich in range(Lchar):
+            charval = str(stringv[ich:ich+1]).replace('[','').replace(']','').replace("'","")
+            if charval == '--' or len(charval) == 0:
+                break
+            else:
+                if ich == 0:
+                    string = charval
+                else:
+                    string = string + charval
+
+        strings.append(string.strip())
+
+    return strings
 
 def interpolate_locs(locs,coords,kinterp):
@@ -6825 +6853 @@
     return
 
-
 def plot_vertical_lev(vertz, vertp, zlog, dzlog, plog, dplog, gtit, kfig, lloc):
     """ plotting vertical levels distribution
@@ -6913 +6940 @@
     return
 
+def plot_subbasin(subname, lons, lats, rsf, rLlf, rof, rcolors, mapv, drawsubid,     \
+  graphtit, kfig, lloc, figname):
+    """ Function to plot rivers from 'river_desc.nc' ORCDHIEE
+      subname= name of the subbasin
+      lons= values for the x-axis
+      lats= values for the y-axis
+      rsf= list of the name of the sub-flows to plot
+      rLlf= dictionary with the lat,lon of the subflows
+      rof= dictionary with the outflows of the subflows
+      rcolors= base of colors of the lines for the subflows (first level)
+      mapv= map characteristics: [proj],[res]
+        see full documentation: http://matplotlib.org/basemap/
+        [proj]: projection
+          * 'cyl', cilindric
+          * 'lcc', lambert conformal
+        [res]: resolution:
+          * 'c', crude
+          * 'l', low
+          * 'i', intermediate
+          * 'h', high
+          * 'f', full
+      drawsubid= wehther sub-flow ids should be plot or not
+      graphtit= title of the graph ('|', for spaces)
+      lloc= location of the legend (0, automatic)
+        1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
+        5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
+        9: 'upper center', 10: 'center'      kfig= kind of figure
+      figname= name of the figure
+    """
+    fname = 'plot_subbasin'
+
+    dx=lons.shape[1]
+    dy=lats.shape[0]
+
+    lengthtrac = np.float(lons[0,1] - lons[0,0])
+    lengthtrac2 = lengthtrac*np.sqrt(2.)
+
+    ddx = np.abs(lons[dy/2+1,dx/2] - lons[dy/2,dx/2])
+    ddy = np.abs(lats[dy/2,dx/2+1] - lats[dy/2,dx/2])
+    fontcharac = {'family': 'serif', 'weight': 'normal', 'size': 3}
+
+    print '  ' + fname + ' Lluis; rsf:', rsf, type(rsf)
+    subhr = gen.chainSnumHierarchy(rsf)
+
+# Getting colors
+    Nsubs = len(rsf)
+    maxNgrids = -1
+    
+    collev = {}
+    if rcolors is None:
+# 1st level colors from colorsauto
+        for isub in range(subhr.xlevs[1]):
+            collev[gen.num_chainSnum(isub+1)] = gen.colhex_dec(colorsauto[isub])
+    else:
+        if len(rcolors) != subhr.Llev1:
+            print errormsg
+            print '  ' + fname + ': number of first subbasins:', subhr.Llev1,        \
+              'and number of provided colors:', len(rcolors), 'differ !!'
+            print '    provided colors:', rcolors
+            quit(-1)
+        for isub in range(subhr.xlevs[1]):
+            collev[gen.num_chainSnum(isub+1)] = rcolors[isub]
+
+    print ' Lluis collev:', collev
+
+# All the colors should finish white...
+    endcol = [1., 1., 1.]
+    colorsub = {}
+    xlabpos = []
+    ylabpos = []
+    xtrack = []
+    ytrack = []
+    colortrack = []
+    outflows = []
+
+    for isub in rsf:
+        levs = gen.num_split(isub)
+        Nlevs = len(levs)
+        begcol = collev[levs[0]]
+        for ilev in range(2,Nlevs):
+# We do not want a last white value...
+            if Nlevs < subhr.Nlevs: 
+                maxNcol = subhr.xlevs[ilev]
+            else:
+                maxNcol = subhr.xlevs[ilev]+1
+            colors = gen.color_deg([begcol, endcol], maxNcol, 'std')
+            begcol = colors[1:4,gen.chainSnum_num(levs[ilev])-1]
+
+#        colorsub[isub] = gen.coldec_hex(begcolors)
+        colorsub[isub] = tuple(begcol)
+        Llfs = rLlf[isub]
+        outfs = rof[isub]
+        print '  Lluis; ' + fname + ': isub:', isub,' outfs:', outfs
+        for igrid in range(len(outfs)):
+            Llf = Llfs[igrid]
+            ibeg = Llf[1]
+            jbeg = Llf[0]
+            outflows.append(outfs[igrid])
+            angle = (outfs[igrid] - 1)*np.pi/4
+            if gen.searchInlist([2,4,6,8], outfs[igrid]):
+                iend = ibeg + lengthtrac2*np.sin(angle)
+                jend = jbeg + lengthtrac2*np.cos(angle)
+            elif gen.searchInlist([1,3,5,7], outfs[igrid]):
+                iend = ibeg + lengthtrac*np.sin(angle)
+                jend = jbeg + lengthtrac*np.cos(angle)
+            else:
+                ibeg = xvals[j,i]
+                jbeg = yvals[j,i]
+                iend = None
+                jend = None
+
+            xtrack.append(ibeg)
+            xtrack.append(iend)
+            xtrack.append(None)
+            ytrack.append(jbeg)
+            ytrack.append(jend)
+            ytrack.append(None)
+            colortrack.append(colorsub[isub])
+        xlabpos.append(ibeg)
+        ylabpos.append(jbeg)
+
+    plt.rc('text', usetex=True)
+
+    latlonss = rLlf.values()
+    nlon = np.min(latlonss[1])
+    xlon = np.max(latlonss[1])
+    nlat = np.min(latlonss[0])
+    xlat = np.max(latlonss[0])
+
+    dlon = xlon - nlon
+    dlat = xlat - nlat
+
+# Making bigger the area to map
+    nlon = nlon-dlon*0.1
+    xlon = xlon+dlon*0.1
+    nlat = nlat-dlat*0.1
+    xlat = xlat+dlat*0.1 
+
+    plt.xlim(nlon,xlon)
+    plt.ylim(nlat,xlat)
+
+    if not mapv is None:
+        lon00 = np.where(lons[:] < 0., 360. + lons[:], lons[:])
+        lat00 = lats[:]
+
+        map_proj=mapv.split(',')[0]
+        map_res=mapv.split(',')[1]
+
+        lon2 = (nlon + xlon)/2.
+        lat2 = (nlat + xlat)/2.
+
+        print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:',     \
+          xlon, ',', xlat
+
+        if map_proj == 'cyl':
+            m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat,         \
+              urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
+        elif map_proj == 'lcc':
+            m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \
+              llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
+        else:
+            print errormsg
+            print '  ' + fname + ": projection '" + map_proj + "' not ready!!"
+            print '    projections available: cyl, lcc'
+            quit(-1)
+
+        m.drawcoastlines()
+
+        meridians = pretty_int(nlon,xlon,5)
+        m.drawmeridians(meridians,labels=[True,False,False,True],color="black")
+
+        parallels = pretty_int(nlat,xlat,5)
+        m.drawparallels(parallels,labels=[False,True,True,False],color="black")
+
+        plt.xlabel('W-E')
+        plt.ylabel('S-N')
+
+    j = 0
+    for i in range(len(rsf)):
+        plt.plot(xtrack[j:j+2], ytrack[j:j+2], color=colortrack[i])
+        j = j + 3
+
+    print 'Lluis : xtrack shape', len(xtrack)/3, 'xlabpos:', len(xlabpos)
+
+# Sea-flow
+    for i in range(len(outflows)):
+        if outflows[i] == 97:
+            plt.plot(xlabpos[i], ylabpos[i], 'x', color=colorsub[i])
+        elif outflows[i] == 98:
+            plt.plot(xlabpos[i], ylabpos[i], '*', color=colorsub[i])
+        elif outflows[i] == 99:
+            plt.plot(xlabpos[i], ylabpos[i], 'h', color=colorsub[i])
+
+    if drawsubid:
+        for i in range(len(xlabpos)):
+            print ' Lluis:', rsf[i], ':', colorsub[rsf[i]]
+            plt.text(xlabpos[i]+0.05*lengthtrac, ylabpos[i]+0.05*lengthtrac, rsf[i], \
+              color=colorsub[rsf[i]], fontdict=fontcharac)
+
+    plt.title(graphtit)
+    output_kind(kfig, figname, True)
+
+    return
+

trunk/tools/generic_tools.py

@@ -19 +19 @@
 
 ####### Content
+# chainSnumHierarchy: Class to provide the structure of a `ChainStrNum' hierarchy
 # chainSnum_levnext: Function to provide the next value for a given level from a chainStrnum hirerarchy of numbers
 # chainSnum_num: Function to pass a `ChainStrNum' string to a number
+# coldec_hex: Function to pas a decimal ([r,g,b]; [0.,1.]) color to hexadecimal (#[RR][GG][BB], 00-64, 0A-FF)
+# colhex_dec: Function to pas a hexadecimal (#[RR][GG][BB]; 00-64, 0A-FF) color to decimal ([0.,1.])
+# color_deg: Function to generate a degradation of colors in rgb base
 # contflow: Function to bring back the increment in j,i grid points according to a trip: (inflow directions)
 # dictionary_key: Function to provide the first key in a dictionay with a given value
@@ -6661 +6665 @@
 def chainSnum_num(cSnum):
     """ Function to pass a `ChainStrNum' string to number
-        `ChainStrNum:' hirerarchy classification where each levels is a significant number within a string:
+        `ChainStrNum:' hierarchy classification where each levels is a significant number within a string:
             1 lev.: [n]
             2 lev.: [n][m]
@@ -6725 +6729 @@
       prevlev: parent level of the desired level
       ChainStrNums: list of strings as `ChainStrNum'
-        `ChainStrNum:' hyerarchic classification where each levels is a significant number within a string:
+        `ChainStrNum:' hierarchyc classification where each levels is a significant number within a string:
             1 lev.: [n]
             2 lev.: [n][m]
@@ -6779 +6783 @@
 
     return parentlev + newchainSnum
+
+class chainSnumHierarchy(object):
+    """ Function to provide the structure of a `ChainStrNum' hierarchy
+      ChainStrNums: list of strings as `ChainStrNum'
+        Nlevs: number of levels
+        levs: all the strings
+        xlevs: list with the maximum number of each level
+        Llev[N]: numbers of values in level N (up to 10 levels)
+        xlev[N]: maximun numbers of a level N (up to 10 levels)
+        Slev[N]: strings of level N (up to 10 levels)
+        Nlev[N]: numbers of level N (up to 10 levels)
+
+        # methods
+        levincr([ChainStrNum]): method to increase `ChainStrNum' level
+
+        `ChainStrNum:' hierarchyc classification where each levels is a significant number within a string:
+            1 lev.: [n]
+            2 lev.: [n][m]
+            3 lev.: [n][m][l]
+            ...
+          NOTE: when a given level has more than 9 values in a given potency of 10, then 'a' is attached in the front: 'a'[i]
+          NOTE: when a given level has more than 19 values in a given potency of 10, then 'b' is attached in the front: 'b'[i]
+          NOTE: when a given level has more than 29 values in a given potency of 10, then 'c' is attached in the front: 'c'[i]
+          (...)
+          NOTE: when a given level has more than 89 values in a given potency of 10, then 'i' is attached in the front: 'i'[i]
+          NOTE: 'o' is for that potencies of 10 without value
+    """
+    fname = 'chainSnumHierarchy'
+
+    def __init__(self,ChainStrNums):
+        self.Nchars = len(ChainStrNums)
+        self.Nlevs = 0
+        self.levs = []
+        self.xlevs = []
+        self.Llev1 = None
+        self.Llev2 = None
+        self.Llev3 = None
+        self.Llev4 = None
+        self.Llev5 = None
+        self.Llev6 = None
+        self.Llev7 = None
+        self.Llev8 = None
+        self.Llev9 = None
+        self.Llev10 = None
+        self.xlev1 = None
+        self.xlev2 = None
+        self.xlev3 = None
+        self.xlev4 = None
+        self.xlev5 = None
+        self.xlev6 = None
+        self.xlev7 = None
+        self.xlev8 = None
+        self.xlev9 = None
+        self.xlev10 = None
+        self.Slev1 = []
+        self.Slev2 = []
+        self.Slev3 = []
+        self.Slev4 = []
+        self.Slev5 = []
+        self.Slev6 = []
+        self.Slev7 = []
+        self.Slev8 = []
+        self.Slev9 = []
+        self.Slev10 = []
+        self.Nlev1 = []
+        self.Nlev2 = []
+        self.Nlev3 = []
+        self.Nlev4 = []
+        self.Nlev5 = []
+        self.Nlev6 = []
+        self.Nlev7 = []
+        self.Nlev8 = []
+        self.Nlev9 = []
+        self.Nlev10 = []
+
+        levsNvals = np.zeros((11,self.Nchars), dtype=int)
+        ilevvals = np.zeros((11), dtype=int)
+        xtrmlevs = []
+
+        Svals1 = []
+        Svals2 = []
+        Svals3 = []
+        Svals4 = []
+        Svals5 = []
+        Svals6 = []
+        Svals7 = []
+        Svals8 = []
+        Svals9 = []
+        Svals10 = []
+        for istr in ChainStrNums:
+            levs = num_split(istr)
+            Nlevs = len(levs)
+            if Nlevs > self.Nlevs: self.Nlevs = Nlevs
+            if Nlevs > 10:
+                print errormsg
+                print '  ' + fname + ': too levels:', Nlevs, " in StringNum '" +     \
+                  istr + "' !!"
+                print '    increase in the code the number of levels in the class'
+                quit(-1)
+            ilevvals[Nlevs] = ilevvals[Nlevs] + 1
+            if Nlevs == 1: Svals1.append(istr)
+            if Nlevs == 2: Svals2.append(istr)
+            if Nlevs == 3: Svals3.append(istr)
+            if Nlevs == 4: Svals4.append(istr)
+            if Nlevs == 5: Svals5.append(istr)
+            if Nlevs == 6: Svals6.append(istr)
+            if Nlevs == 7: Svals7.append(istr)
+            if Nlevs == 8: Svals8.append(istr)
+            if Nlevs == 9: Svals9.append(istr)
+            if Nlevs == 10: Svals10.append(istr)
+
+            levsNvals[Nlevs,ilevvals[Nlevs]] = chainSnum_num(levs[Nlevs-1])
+
+        valslist = ChainStrNums.sort()
+        self.levs = valslist
+        # 1st lev
+        self.Llev1 = ilevvals[1]
+        Svlist = Svals1
+        Nvlist = list(levsNvals[1,1:self.Llev1+1])
+        self.xlev1 = np.max(Nvlist)
+        xtrmlevs.append(self.xlev1)
+        Svlist.sort()
+        Nvlist.sort()
+        self.Slev1 = Svlist
+        self.Nlev1 = Nvlist
+        # 2nd lev
+        if self.Nlevs > 1:
+            self.Llev2 = ilevvals[2]
+            Svlist = Svals2
+            Nvlist = list(levsNvals[2,1:self.Llev2+1])
+            self.xlev2 = np.max(Nvlist)
+            xtrmlevs.append(self.xlev2)
+            Svlist.sort()
+            Nvlist.sort()
+            self.Slev2 = Svlist
+            self.Nlev2 = Nvlist
+        # 3rd lev
+        if self.Nlevs > 2:
+            self.Llev3 = ilevvals[3]
+            Svlist = Svals3
+            Nvlist = list(levsNvals[3,1:self.Llev3+1])
+            self.xlev3 = np.max(Nvlist)
+            xtrmlevs.append(self.xlev3)
+            Svlist.sort()
+            Nvlist.sort()
+            self.Slev3 = Svlist
+            self.Nlev3 = Nvlist
+        # 4th lev
+        if self.Nlevs > 3:
+            self.Llev4 = ilevvals[4]
+            Svlist = Svals4
+            Nvlist = list(levsNvals[4,1:self.Llev4+1])
+            self.xlev4 = np.max(Nvlist)
+            xtrmlevs.append(self.xlev4)
+            Svlist.sort()
+            Nvlist.sort()
+            self.Slev4 = Svlist
+            self.Nlev4 = Nvlist
+        # 5th lev
+        if self.Nlevs > 4:
+            self.Llev5 = ilevvals[5]
+            Svlist = Svals5
+            Nvlist = list(levsNvals[5,1:self.Llev5+1])
+            self.xlev5 = np.max(Nvlist)
+            xtrmlevs.append(self.xlev5)
+            Svlist.sort()
+            Nvlist.sort()
+            self.Slev5 = Svlist
+            self.Nlev5 = Nvlist
+        # 6th lev
+        if self.Nlevs > 5:
+            self.Llev6 = ilevvals[6]
+            Svlist = Svals6
+            Nvlist = list(levsNvals[6,1:self.Llev6+1])
+            self.xlev6 = np.max(Nvlist)
+            xtrmlevs.append(self.xlev6)
+            Svlist.sort()
+            Nvlist.sort()
+            self.Slev6 = Svlist
+            self.Nlev6 = Nvlist
+        # 7th lev
+        if self.Nlevs > 6:
+            self.Llev7 = ilevvals[7]
+            Svlist = Svals7
+            Nvlist = list(levsNvals[7,1:self.Llev7+1])
+            self.xlev7 = np.max(Nvlist)
+            xtrmlevs.append(self.xlev7)
+            Svlist.sort()
+            Nvlist.sort()
+            self.Slev7 = Svlist
+            self.Nlev7 = Nvlist
+        # 8th lev
+        if self.Nlevs > 7:
+            self.Llev8 = ilevvals[8]
+            Svlist = Svals8
+            Nvlist = list(levsNvals[8,1:self.Llev8+1])
+            self.xlev8 = np.max(Nvlist)
+            xtrmlevs.append(self.xlev8)
+            Svlist.sort()
+            Nvlist.sort()
+            self.Slev8 = Svlist
+            self.Nlev8 = Nvlist
+        # 9th lev
+        if self.Nlevs > 9:
+            self.Llev9 = ilevvals[9]
+            Svlist = Svals9
+            Nvlist = list(levsNvals[9,1:self.Llev9+1])
+            self.xlev9 = np.max(Nvlist)
+            xtrmlevs.append(self.xlev9)
+            Svlist.sort()
+            Nvlist.sort()
+            self.Slev9 = Svlist
+            self.Nlev9 = Nvlist
+        # 10th lev
+        if self.Nlevs > 9:
+            self.Llev10 = ilevvals[10]
+            Svlist = Svals10
+            Nvlist = list(levsNvals[10,1:self.Llev10+1])
+            self.xlev10 = np.max(Nvlist)
+            xtrmlevs.append(self.xlev10)
+            Svlist.sort()
+            Nvlist.sort()
+            self.Slev10 = Svlist
+            self.Nlev10 = Nvlist
+        self.xlevs = xtrmlevs
+
+    def levincr(self,plev):
+        """ Method to increase a given plevel level
+        levincr(self,plev): 
+        >>> ChainStrNums = ['1', '2', '11', '12', '13', '14', '15', '16', '17', '18', '19', '1a0', '111', '112', '113', '121', '122' ]
+        >>> CSN_hyr = chainSnumHierarchy(ChainStrNums)
+        >>> CSN_hyr.levincr('11')
+        114
+        """
+        return chainSnum_levnext(plev, ChainStrNums)
+
+#ChainStrNums = ['1', '2', '11', '12', '13', '14', '15', '16', '17', '18', '19', '1a0', '111', '112', '113', '121', '122' ]
+#CSN_hyr = chainSnumHierarchy(ChainStrNums)
+
+#print CSN_hyr.levincr.__doc__
+#quit()
 
 def subbasin_point(trips, pt):
@@ -6951 +7196 @@
 #
 #print 'Total sub-basin:', np.sum(masksubbasin)
+
+def color_deg(colors, Nsteps, typedeg, values=None):
+    """ Function to generate a degradation of colors in rgb base
+      colors= list of colors to use [r,g,b] r,g,b <= 1.
+      Nsteps= number of steps (NOTE: final number of colors might differ from this value)
+      typedeg= type of degradate
+        'std': simple linear degradation between equi-spaced colors 
+        'vals': change of colors at a given value (it requires [values])
+      values= list of values to use as reference (when applicable)
+      provides: [colordeg], [Ncolors]
+        [colordeg] = (4, [Ncolors]) matrix with
+          [0,:]= assigned value to the color (for 'std': counter)
+          [1:4,:]= r,g,b color
+    >>> color_deg([ [1.,0.,0.], [0.,1.,0.], [0.,0.,1.] ], 10, 'std')
+    [[ 0.     0.125  0.25   0.375  0.5    0.625  0.75   0.875  1.    ]
+     [ 1.     0.75   0.5    0.25   0.     0.     0.     0.     0.    ]
+     [ 0.     0.25   0.5    0.75   1.     0.75   0.5    0.25   0.    ]
+     [ 0.     0.     0.     0.     0.     0.25   0.5    0.75   1.    ]]
+    defcolors = color_deg([ [1.,0.,0.], [0.,1.,0.], [0.,0.,1.]], 12, 'vals', [-3., 0., 3.])
+    [[-3.  -2.4 -1.8 -1.2 -0.6  0.   0.6  1.2  1.8  2.4  3. ]
+     [ 1.   0.8  0.6  0.4  0.2  0.   0.   0.   0.   0.   0. ]
+     [ 0.   0.2  0.4  0.6  0.8  1.   0.8  0.6  0.4  0.2  0. ]
+     [ 0.   0.   0.   0.   0.   0.   0.2  0.4  0.6  0.8  1. ]]
+    """
+    fname = 'color_deg'
+
+    colordeg = np.zeros((4,Nsteps), dtype=np.float)
+    Ncols = len(colors)
+    cols = np.zeros((3,Ncols), dtype=np.float)
+    for icol in range(Ncols):
+        cols[:,icol] = np.array(colors[icol])
+    
+    if typedeg == 'std':
+        Nstepcols = np.int(Nsteps*1. / (Ncols-1.))
+        iicol = 0
+        for icol in range(Ncols-1):
+            dcolors = (cols[:,icol+1]- cols[:,icol])/(Nstepcols-1.)
+# To avoid repetition of the color
+            if icol == 0: 
+                ibeg = 0
+            else:
+                ibeg = 1
+            for i in range(ibeg,Nstepcols):
+                colordeg[0,iicol] = iicol*1.
+                colordeg[1:4,iicol] = cols[:,icol] + dcolors*i
+                iicol = iicol + 1
+        Nfincol = iicol-1
+        colordeg[0,:] = colordeg[0,:]/(Nfincol*1.)
+
+    elif typedeg == 'vals':
+        if values is None:
+            print errormsg
+            print '  ' + fname + ": type '" + typedeg + "' require values !!"
+            quit(-1)
+        Nvals = len(values)
+        if Nvals != Ncols:
+            print errormsg
+            print '  ' + fname + ': different number of values:', Nvals,             \
+              'and colors:', Ncols, '!!'
+            print '    provided values:', values
+            print '    provided colors:', cols
+            quit(-1)
+        TOTrange = values[Nvals-1] - values[0]
+        iicol = 0
+        for icol in range(Ncols-1):
+            pairange = values[icol+1] - values[icol]
+            drange = pairange / TOTrange
+            Nstepcols = np.int(Nsteps*drange)
+            dcolors = (cols[:,icol+1]-cols[:,icol])/(Nstepcols-1.)
+# To avoid repetition of the color
+            if icol == 0: 
+                ibeg = 0
+            else:
+                ibeg = 1
+            for i in range(ibeg,Nstepcols):
+                colordeg[0,iicol] = values[icol]+i*pairange/(Nstepcols-1)
+                colordeg[1:4,iicol] = cols[:,icol] + dcolors*i
+                iicol = iicol + 1
+        Nfincol = iicol - 1
+    else:
+        print errormsg
+        print '  ' + fname + ": degradation type '" + typedeg + "' not ready !!"
+        quit(-1)
+
+    degcolor = np.zeros((4,Nfincol), dtype=np.float)
+    degcolor = colordeg[:,0:Nfincol+1]
+
+    return degcolor
+
+def colhex_dec(col):
+    """ Function to pas a hexadecimal (#[RR][GG][BB]; 00-64, 0A-FF) color to decimal ([0.,1.])
+      col= '#'[RR][GG][BB] hexadecimal color string
+      >>>colhex_dec('#000000')
+      [0.0, 0.0, 0.0]
+      >>>colhex_dec('#00FF00')
+      [0.0, 1.0, 0.0]
+      >>>colhex_dec('#3264EE')
+      [0.19607843137254902, 0.39215686274509803, 0.9333333333333333]
+    """
+    fname = 'colhex_dec'
+    hexred = col[1:3]
+    hexgreen = col[3:5]
+    hexblue = col[5:7]
+
+    deccol = [int(hexred, 16)/255., int(hexgreen, 16)/255., int(hexblue, 16)/255.,]
+
+    return deccol
+
+def coldec_hex(col):
+    """ Function to pas a decimal ([r,g,b]; [0.,1.]) color to hexadecimal (#[RR][GG][BB], 00-64, 0A-FF)
+      col= '#'[RR][GG][BB] hexadecimal color string
+      >>>coldec_hex([0.0, 0.0, 0.0])
+      #000000
+      >>>coldec_hex([0.0, 1.0, 0.0])
+      #00FF00
+      >>>coldec_hex([0.19607843137254902, 0.39215686274509803, 0.9333333333333333])
+      #3264EE
+    """
+    fname = 'coldec_hex'
+
+    hexred = '{0:X}'.format(int(col[0]*255))
+    hexgreen = '{0:X}'.format(int(col[1]*255))
+    hexblue = '{0:X}'.format(int(col[2]*255))
+
+    if hexred == '0': hexred = '00'
+    if hexgreen == '0': hexgreen = '00'
+    if hexblue == '0': hexblue = '00'
+
+    hexcol = '#' + hexred + hexgreen + hexblue
+
+    return hexcol
+
+#quit()
+