Changeset 785 in lmdz_wrf

Timestamp:

May 28, 2016, 5:07:13 PM (9 years ago)

Author:

lfita

Message:

Adding `plot_river_pattern' to plot a given river and its subbasins with the stations using:

• routing.nc: discharge grid tracjectories
• river_desc.nc: river subbasins and description
• GRDC_Monthly_June14.nc: stations location

Location:

trunk/tools

Files:

• drawing.py (modified) (1 diff)
• drawing_tools.py (modified) (2 diffs)

trunk/tools/drawing.py

@@ -3333 +3333 @@
 
 def draw_basins(ncfile, values, varns):
-    """ Function to plot wind basins
+    """ Function to plot river basins with their discharge vector and basins id (from 'routing.nc')
       values= [lonlatbox]:[mapres]:[cbarname]:[xtrmbasin]:[mapdraw]:[veclength]:[freq]:
         [ifreq]:[plotcountry]:[basinidn]:[gtit]:[kindfig]:[figuren]

trunk/tools/drawing_tools.py

@@ -11 +11 @@
 from netCDF4 import Dataset as NetCDFFile
 import nc_var_tools as ncvar
+import generic_tools as gen
 
 errormsg = 'ERROR -- error -- ERROR -- error'
@@ -6397 +6398 @@
     return
 
+def plot_river_pattern(xvals, yvals, fvals, veccolor, descid, desclon, desclat,      \
+  stn, stlon, stlat, drawcountry, drawbasinid, colordescid, coloriver, colornoriver, \
+  colorst, graphtit, mapv, kfig, figname):
+    """ Function to plot vectors
+      xvals= values for the x-axis
+      yvals= values for the y-axis
+      fvals= values for the flow (1-8: N, NE, E, ..., NW; 97: sub-basin; 98: small to sea; 99: large to sea)
+      veccolor= basin id
+      descid= description id
+, desclon, desclat,
+      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
+      drawcountry= whether country lines should be plotted or not
+      drawbasinid= whether basins id should be plotted or not
+      graphtit= title of the graph ('|', for spaces)
+      kfig= kind of figure
+      figname= name of the figure
+    """
+    fname = 'plot_river_pattern'
+
+    dx=xvals.shape[1]
+    dy=xvals.shape[0]
+    xmin = np.min(xvals)
+    xmax = np.max(xvals)
+    ymin = np.min(yvals)
+    ymax = np.max(yvals)
+
+# flow direction
+
+    plt.rc('text', usetex=True)
+
+    ddx = np.abs(xvals[dy/2+1,dx/2] - xvals[dy/2,dx/2])
+    ddy = np.abs(yvals[dy/2,dx/2+1] - yvals[dy/2,dx/2])
+    fontcharac = {'family': 'serif', 'weight': 'normal', 'size': 3}
+
+    xlabpos = []
+    ylabpos = []
+    labels = []
+    labcol = []
+    flow = []
+    flowvals = []
+    xtrack = []
+    ytrack = []
+    colortrack = []
+    trackid = []
+    stnlL = []
+    stlonlL = []
+    stlatlL = []
+
+# Setting up colors for each label
+#   From: http://stackoverflow.com/questions/15235630/matplotlib-pick-up-one-color-associated-to-one-value-in-a-colorbar
+    my_cmap = plt.cm.get_cmap(colordescid)
+    vcolmin = np.min(descid)
+    vcolmax = np.max(descid)
+    print 'Lluis vcolmin:',vcolmin,'vcolmax:',vcolmax
+
+    norm = mpl.colors.Normalize(vcolmin, vcolmax)
+
+# Vector angles
+    lengthtrac = np.float(xvals[0,1] - xvals[0,0])
+    lengthtrac2 = lengthtrac*np.sqrt(2.)
+
+    for j in range(0,dy):
+        for i in range(0,dx):
+            if veccolor[j,i] != '--':
+#                ibeg = xvals[j,i]-lengthtrac/2.
+#                jbeg = yvals[j,i]-lengthtrac/2.
+                ibeg = xvals[j,i]
+                jbeg = yvals[j,i]
+                labels.append(int(veccolor[j,i]))
+                flowvals.append(fvals[j,i])
+                angle = (fvals[j,i] - 1)*np.pi/4
+                if gen.searchInlist([2,4,6,8], fvals[j,i]):
+                    iend = ibeg + lengthtrac2*np.sin(angle)
+                    jend = jbeg + lengthtrac2*np.cos(angle)
+                elif gen.searchInlist([1,3,5,7], fvals[j,i]):
+                    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
+
+                xlabpos.append(ibeg)
+                ylabpos.append(jbeg)
+                xtrack.append(ibeg)
+                xtrack.append(iend)
+                xtrack.append(None)
+                ytrack.append(jbeg)
+                ytrack.append(jend)
+                ytrack.append(None)
+                if len(desclon.shape) == 2:
+                    difflonlat = np.sqrt((desclon-xvals[j,i])**2 + (desclat-yvals[j,i])**2)
+                    mindiffLl = np.min(difflonlat)
+                    ilatlon = gen.index_mat(difflonlat, mindiffLl)
+                else:
+                    ilatlon = np.zeros((2), dtype=int)
+                    difflon = np.abs(desclon - xvals[j,i])
+                    mindiffl = np.min(difflon)
+                    ilatlon[1] = gen.index_vec(difflon,mindiffl)
+                    difflat = np.abs(desclat - yvals[j,i])
+                    mindiffL = np.min(difflat)
+                    ilatlon[0] = gen.index_vec(difflat,mindiffL)
+
+                if descid.mask[ilatlon[0],ilatlon[1]]: 
+                    labcol.append(colornoriver)
+                    colortrack.append(colornoriver)
+                else:
+                    if veccolor[j,i] != 6:
+                        print 'Lluis: veccol:', veccolor[j,i], 'mindiffl', mindiffl,'mindiffL:',mindiffL
+                    else:
+                        print 'Lluis right! mindiffl', mindiffl,'mindiffL:',mindiffL
+                    labcol.append(coloriver)
+                    colortrack.append(my_cmap(norm(descid[ilatlon[0],ilatlon[1]])))
+
+                trackid.append(descid[ilatlon[0],ilatlon[1]])
+
+    totvals = len(flowvals)
+
+    if not mapv is None:
+        lon00 = np.where(xvals[:] < 0., 360. + xvals[:], xvals[:])
+        lat00 = yvals[:]
+
+        map_proj=mapv.split(',')[0]
+        map_res=mapv.split(',')[1]
+
+        nlon = np.min(xvals)
+        xlon = np.max(xvals)
+        nlat = np.min(yvals)
+        xlat = np.max(yvals)
+
+        lon2 = xvals[dy/2,dx/2]
+        lat2 = yvals[dy/2,dx/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()
+        drawcountry = True
+        if drawcountry: m.drawcountries()
+
+        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(xlabpos)):
+        plt.plot(xtrack[j:j+2], ytrack[j:j+2], color=colortrack[i])
+        j = j + 3
+
+# Sea-flow
+    for i in range(len(xlabpos)):
+        if flowvals[i] == 97:
+            plt.plot(xlabpos[i], ylabpos[i], 'x', color=colortrack[i])
+        elif flowvals[i] == 98:
+            plt.plot(xlabpos[i], ylabpos[i], '*', color=colortrack[i])
+        elif flowvals[i] == 99:
+            plt.plot(xlabpos[i], ylabpos[i], 'h', color=colortrack[i])
+
+    plt.xlim(xmin,xmax)
+    plt.ylim(ymin,ymax)
+
+    if drawbasinid:
+        for i in range(len(xlabpos)):
+            plt.text(xlabpos[i]+0.05*lengthtrac, ylabpos[i]+0.05*lengthtrac, labels[i],    \
+              color=labcol[i], fontdict=fontcharac)
+#            plt.text(xlabpos[i]+0.1*lengthtrac, ylabpos[i]+0.1*lengthtrac, trackid[i],   \
+#              color=colortrack[i], fontdict=fontcharac)
+#            plt.text(xlabpos[i]+0.1*lengthtrac, ylabpos[i]+0.1*lengthtrac, str(flowvals[i]),   \
+#              color=colortrack[i], fontdict=fontcharac)
+
+    for ist in range(len(stlon)):
+        if stlon[ist] > xmin and stlon[ist] < xmax and stlat[ist] > ymin and stlat[ist] < ymax:
+            stname = '*'
+            for ic in range(len(stn[ist,:])):
+                stname  = stname + stn[ist,ic]
+#            plt.annotate(stname, xy=(stlon[ist],stlat[ist]), color=colorst, )
+            plt.text(stlon[ist], stlat[ist], stname, color=colorst, fontsize=3)
+
+#    cbar = plt.colorbar()
+    plt.title(graphtit.replace('&','\&'))
+
+    output_kind(kfig, figname, True)
+
+    return
+
 def var_3desc(ovar):
     """ Function to provide std_name, long_name and units from an object variable