Changeset 705 in lmdz_wrf for trunk

Timestamp:

Apr 6, 2016, 4:32:37 PM (9 years ago)

Author:

lfita

Message:

Adding draw_river_desc' to plot rivers' description after river_desc.nc' from ORCDHIEE

Location:

trunk/tools

Files:

• drawing.py (modified) (7 diffs)
• drawing_tools.py (modified) (3 diffs)

trunk/tools/drawing.py

@@ -28 +28 @@
 ## e.g. # drawing.py -o draw_vectors -f wrfout_d01_2001-11-11_00:00:00 -S 'T|Time|Times|2,Y|south_north|XLAT|-1,X|west_east|XLONG|-1:3@3,wind@rainbow,9:10m wind,ms-1:cyl,l:WRF 10 m winds:pdf:winds' -v U10,V10
 ## e.g. # drawing.py -o draw_basins -f routing.py -S 'Y|y|nav_lat|-1,X|x|nav_lon|-1:1@1,rainbow,9:basins,-:cyl,l:ORCDHIEE river-basins:pdf:basins_named' -v nav_lon,nav_lat,trip,basins
+## 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
 
 main = 'drawing.py'
@@ -40 +41 @@
   'draw_lines', 'draw_lines_time', 'draw_Neighbourghood_evol',                       \
   'draw_points', 'draw_points_lonlat',                                               \
-  'draw_ptZvals', 'draw_timeSeries',                                                 \
+  'draw_ptZvals', 'draw_river_desc', 'draw_timeSeries',                              \
   'draw_topo_geogrid',                                                               \
   'draw_topo_geogrid_boxes', 'draw_trajectories', 'draw_vals_trajectories',          \
@@ -3329 +3330 @@
     return
 
-
 def draw_basins(ncfile, values, varns):
     """ Function to plot wind basins
@@ -3457 +3457 @@
     if len(varn.split(',')) != 3:
         print errormsg
-            print '  ' + fname + ": color of vectors should be according to '" +     \
-              coln + "' but a third varibale is not provided !!"
-            quit(-1)
+        print '  ' + fname + ": color of vectors should be according to '" +         \
+          coln + "' but a third varibale is not provided !!"
+        quit(-1)
 
     ocolvec = of.variables[varNs[3]]
@@ -3484 +3484 @@
 
     return
+
+def draw_river_desc(ncfile, values, riverns):
+    """ Function to plot rivers' description from ORCHIDEE's routing scheme
+      values= [dimname]|[vardimname]|[value]:[basinvals]:[upstreamvals]:[mapvalues]:
+        [gtit]:[kindfig]:[legloc]:[figuren]
+        'X/Y'|[dimname]|[vardimname]|[value]: ',', list for each basic dimension '|' separated of:
+          [dimname]: name of the dimension in the file
+          [vardimname]: name of the variable with the values for the dimension in the file
+          [value]: which value of the given dimension (-1, all; [ibeg]@[iend], i-range beginning, end)
+          No value takes all the range of the dimension
+        [basinsvals]= [colorline]
+          [basincolor]: ',' list of colors of the line to use to mark the basins contours (single value also possible)
+        [upstreamvals]= [upstreamvarn],[colorbar]
+          [upstreamcolor]: colorbar to use to plot the basins upstream values
+        [mapvalues]= 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
+        gtit= title of the graph ('|', for spaces)
+        kindfig= kind of figure
+        legloc= 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
+        figuren= name of the figure
+      ncfile= file to use
+      riverns= ',' list of the name of the rivers to plot
+    """
+    import numpy.ma as ma
+    fname = 'draw_river_desc'
+
+    if values == 'h':
+        print fname + '_____________________________________________________________'
+        print draw_river_desc.__doc__
+        quit()
+
+    expectargs = '[X/Y/Z/T]|[dimname]|[vardimname]|[value]:[basinvals]:' +           \
+      '[upstreamvals]:[mapvalues]:[gtit]:[kindfig]:[legloc]:[figuren]'
+ 
+    drw.check_arguments(fname,values,expectargs,':')
+
+    dimvals = values.split(':')[0]
+    basinvals = values.split(':')[1]
+    upstreamvals = values.split(':')[2]
+    mapvals = values.split(':')[3]
+    gtit = values.split(':')[4]
+    kindfig = values.split(':')[5]
+    legloc = int(values.split(':')[6])
+    figuren = values.split(':')[7]
+
+    basincol = basinvals
+    if basincol.find(',') != 1:
+        basincolor = basincol.split(',')
+    else:
+        basincolor = [basincol]
+
+    upstreamcolor = upstreamvals
+
+    of = NetCDFFile(ncfile,'r')
+
+    dims = {}
+    for dimv in dimvals.split(','):
+        dns = dimv.split('|')
+        dims[dns[0]] = [dns[1], dns[2], dns[3]]
+
+    varNs = []
+    for dn in dims.keys():
+        if dn == 'X':
+            varNs.append(dims[dn][1])
+            dimx = len(of.dimensions[dims[dn][0]])
+        elif dn == 'Y':
+            varNs.append(dims[dn][1])
+            dimy = len(of.dimensions[dims[dn][0]])
+
+    if riverns.find(',') != -1:
+        riverns = riverns.split(',')
+    else:
+        riverns = [riverns]
+
+    rivers = []
+    riversubbasins = {}
+    riversupstream = {}
+    riversoutflow = {}
+    for rivern in riverns:
+        print rivern
+
+# subBasins
+        basinvar = rivern + '_coding'
+        if not drw.searchInlist(of.variables.keys(), basinvar):
+            print errormsg
+            print '  ' + fname + ": file does not have variable '" + basinvar + "' !!"
+            quit(-1)
+        rivers.append(rivern)
+        obasin = of.variables[basinvar]
+        riversubbasins[rivern] = obasin[:]
+        if rivern == riverns[0]:
+            finalmask = obasin[:].mask
+        else:
+            finalmask = finalmask * obasin[:].mask
+
+# upstream
+        upstreamvar = rivern + '_upstream'
+        if not drw.searchInlist(of.variables.keys(), upstreamvar):
+            print errormsg
+            print '  ' + fname + ": file does not have variable '" + upstreamvar + "' !!"
+            quit(-1)
+        oupstream = of.variables[upstreamvar]
+        riversupstream[rivern] = oupstream[:]
+        if rivern == riverns[0]:
+            uunits = oupstream.getncattr('units')
+
+# River metadata
+        fracvar = rivern + '_frac'
+        if not drw.searchInlist(of.variables.keys(), fracvar):
+            print errormsg
+            print '  ' + fname + ": file does not have variable '" + fracvar + "' !!"
+            quit(-1)
+        ofrac = of.variables[fracvar]
+        riversoutflow[rivern] = [ofrac.getncattr('Longitude_of_outflow_point'),      \
+          ofrac.getncattr('Latitude_of_outflow_point')]
+
+    ivar = 0
+    for varN in varNs:
+        varslice = []
+
+        ovarN = of.variables[varN]
+        vard = ovarN.dimensions
+        for vdn in vard:
+            found = False
+            for dd in dims.keys():
+                if dims[dd][0] == vdn:
+                    if dims[dd][2].find('@') != -1:
+                        rvals = dims[dd][2].split('@')
+                        varslice.append(slice(int(rvals[0]), int(rvals[1])))
+                    elif dims[dd][2] == '-1':
+                        varslice.append(slice(0,len(of.dimensions[dims[dd][0]])))
+                    else:
+                        varslice.append(int(dims[dd][2]))
+
+                    found = True
+                    break
+            if not found:
+                varslice.append(slice(0,len(of.dimensions[dims[dd][0]])))
+
+        if varN == dims['X'][1]:
+            lonvals0 = np.squeeze(ovarN[tuple(varslice)])
+        elif varN == dims['Y'][1]:
+            latvals0 = np.squeeze(ovarN[tuple(varslice)])
+
+        ivar = ivar + 1
+
+    if len(lonvals0.shape) == 1:
+        lonvals, latvals = np.meshgrid(lonvals0, latvals0)
+    else:
+        lonvals = lonvals0
+        latvals = latvals0
+
+# Masking only the lon,lat with rivers
+    malonvals = ma.masked_array(lonvals, mask=finalmask)
+    malatvals = ma.masked_array(latvals, mask=finalmask)
+
+    if mapvals == 'None':
+        mapvalues = None
+    else:
+        mapvalues = mapvals
+
+    drw.plot_river_desc(malonvals, malatvals, rivers, riversubbasins, riversupstream, riversoutflow,  \
+      basincolor, upstreamcolor, uunits, mapvalues, gtit, kindfig, legloc, figuren)
+
+    of.close()
+
 #quit()
 
@@ -3489 +3667 @@
 
 ngraphics = "'" + drw.numVector_String(namegraphics, "', '") + "'"
-
+  
 ### Options
 ##string_operation="operation to make: " + '\n' + " out, output values -S inidim1,[inidim2,...]:enddim1,[enddim2,...]"
@@ -3568 +3746 @@
 elif oper == 'draw_ptZvals':
     draw_ptZvals(opts.ncfile, opts.values, opts.varname)
+elif oper == 'draw_river_desc':
+    draw_river_desc(opts.ncfile, opts.values, opts.varname)
 elif oper == 'draw_timeSeries':
     draw_timeSeries(opts.ncfile, opts.values, opts.varname)

trunk/tools/drawing_tools.py

@@ -3 +3 @@
 # From http://stackoverflow.com/questions/13336823/matplotlib-python-error
 import numpy as np
-from matplotlib.pylab import *
 import matplotlib as mpl
 mpl.use('Agg')
+from matplotlib.pylab import *
 import matplotlib.pyplot as plt
 from mpl_toolkits.basemap import Basemap
@@ -1520 +1520 @@
         elif u == '1/m': lu='$m^{-1}$'
         elif u == 'm-1': lu='$m^{-1}$'
+        elif u == 'm^2': lu='$m^{2}$'
         elif u == 'm2/s': lu='$m2s^{-1}$'
         elif u == 'm2s-1': lu='$m2s^{-1}$'
@@ -6406 +6407 @@
 
     return stdn, lonn, un
+
+def plot_river_desc(lons, lats, rns, rss, rus, ros, bcolor, ucolor, uuts, mapv,      \
+  graphtit, kfig, lloc, figname):
+    """ Function to plot rivers from 'river_desc.nc' ORCDHIEE
+      lons= values for the x-axis
+      lats= values for the y-axis
+      rns= list of the name of the rivers to plot
+      rss= dictionary with the lon,lats of the subbasins of each river
+      rus= dictionary with the lon,lats of the upstream values of each river
+      ros= dictionary with the lon,lats of the outflow points of each river
+      bcolor= color of the lines for the subbasins
+      ucolor= bar color for the upstreams
+      uuts= units of the upstream
+      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
+      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_basins'
+
+    colvalues, linekinds, pointkinds, lwidths, psizes = ColorsLinesPointsStyles(     \
+      len(rns), bcolor, ['-'], ['.'], [1.], [1.],  None)
+
+    dx=lons.shape[1]
+    dy=lats.shape[0]
+
+    plt.rc('text', usetex=True)
+
+    nlon = np.min(lons)
+    xlon = np.max(lons)
+    nlat = np.min(lats)
+    xlat = np.max(lats)
+
+    dlon = xlon - nlon
+    dlat = xlat - nlat
+
+# Making bigger the area to map
+    nlon = np.min(lons)-dlon*0.1
+    xlon = np.max(lons)+dlon*0.1
+    nlat = np.min(lats)-dlat*0.1
+    xlat = np.max(lats)+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')
+
+    umin = 10.e20
+    umax = -umin
+    for rn in rns:
+        nrns = np.min(rus[rn])
+        xrns = np.max(rus[rn])
+
+        if umin > nrns: umin = nrns
+        if umax < xrns: umax = xrns
+
+    i = 0
+    for rn in rns:
+        plt.pcolormesh(lons, lats, rus[rn], vmin=umin, vmax=umax,                    \
+          cmap=plt.get_cmap(ucolor), norm=matplotlib.colors.LogNorm())
+#        plt.pcolormesh(lons, lats, rss[rn], cmap=plt.get_cmap(ucolor))
+        if rn == rns[0]: cbar = plt.colorbar()
+        riverarea = rss[rn]
+
+        riverarea.mask = ma.nomask       
+        rarea = np.where(riverarea == 999999999, -5, riverarea)
+        rarea = np.where(rarea > 0, 1, 0)
+        plt.contour(lons, lats, rarea, levels = [0.], colors=colvalues[i],           \
+          linewidths=1.5)
+
+        oflow = ros[rn]
+        plt.plot(oflow[0], oflow[1], 'h', color=colvalues[i], label=rn)
+#        plt.annotate(rn, xy=(0., i*0.05), xycoords='axes fraction', color=colvalues[i])
+        i = i+1
+
+
+# Attempts to plot sub-basins
+
+# Gradient subbasins (not working)
+#        gradriver = np.zeros((dy,dx), dtype=int)
+#        for j in range(dy-2):
+#            for i in range(dx-2):
+#                gradriver[j,i] = riverarea[j,i+1]-riverarea[j,i]+riverarea[j+1,i]-  \
+#                  riverarea[j,i]
+#        garea = np.where(gradriver != 0, 1, gradriver)
+#        print garea
+#        plt.contour(lons, lats, garea, levels = [0.], colors='red', linewidths=1.)
+
+# Not working
+#        subbasins = np.sort(np.unique(rarea))
+#        for subb in subbasins:
+#            if subb > 0.:
+#                sarea = np.where(rarea != subb, 0, rarea)
+#                plt.contour(lons, lats, sarea, levels = [1.], colors='gray', linewidths=1.)
+
+    cbar.set_label('upstream (' + units_lunits(uuts) + ')')
+    plt.legend(loc=lloc)
+
+#    plt.annotate(windn.replace('_','\_') +' (' + units_lunits(wuts) + ')',           \
+#      xy=(0.80,-0.15), xycoords='axes fraction', color='black')
+
+    plt.title(graphtit.replace('|',' ').replace('&','\&'))
+
+    output_kind(kfig, figname, True)
+
+    return
+