Plot of graphs improved but still not perfect.

Diags/plotgraphs.py

@@ -4,6 +4,7 @@ import matplotlib.pyplot as plt
 from mpl_toolkits.basemap import Basemap, cm
 import numpy as np
 from netCDF4 import Dataset
+import sys
 
 NCfile="../MEDCORDEX_test_graph.nc"
 GraphicFile = "test.png"
@@ -39,6 +40,12 @@ lat_bnd = nc.variables["lat_bnd"][:,:,:]
 land = nc.variables["land"][:,:]
 routetogrid = nc.variables["routetogrid"][:,:,:]
 routetohtu = nc.variables["routetobasin"][:,:,:]
+#
+# Coding for routetobasin :
+# nbasmax + 1 = Returnflow to the grid
+# nbasmax + 2 = Coastal flow
+# nbasmax + 3 = River flow
+#
 cglon = nc.variables["CG_lon"][:,:,:]
 cglat = nc.variables["CG_lat"][:,:,:]
 
@@ -46,7 +53,7 @@ nc.close()
 
 nhtu, nj, ni = routetogrid.shape
 landindex = LandGrid(land)
-print "Nb land points = ", landindex.nbland, 'Nb HTU : ', nhtu
+print("Nb land points = ", landindex.nbland, 'Nb HTU : ', nhtu)
 
 #
 # Get some average size of the grid boxes.
@@ -75,6 +82,8 @@ m.drawmeridians(np.arange(-80.,80.,0.1),labels=[0,0,0,1]) # draw meridians
 cmap = 'RdBu'
 cmap = 'Pastel1'
 cmap = 'prism'
+lwidth=0.0005
+markersize=5
 #
 #
 #
@@ -90,24 +99,33 @@ for i in range(ni) :
 #
 for il in range(landindex.nbland) :
     j,i = landindex.land2ij(il)
-    print il,"== routetogrid ==",routetogrid[:,j,i]
-    print il,"== routetohtu  ==",routetohtu[:,j,i]
+    xt,yt=m(lon_full[j,i],lat_full[j,i])
+    plt.text(xt,yt,str(il))
+    #
     for ih in range(nhtu) :
         lo = cglon[ih,j,i]
         la = cglat[ih,j,i]
         xl,yl=m(lo, la)
         if int(routetohtu[ih,j,i]-1) < nhtu :
-            plt.plot(xl,yl, marker=11)
+            plt.plot(xl,yl, ms=markersize, marker=11)
+        elif int(routetohtu[ih,j,i]-1) == nhtu+1 :
+            plt.plot(xl,yl, ms=markersize, marker="x")
+        elif int(routetohtu[ih,j,i]-1) == nhtu+2 :
+            plt.plot(xl,yl, ms=markersize, marker="o")
+        elif int(routetohtu[ih,j,i]-1) == nhtu+3 :
+            plt.plot(xl,yl, ms=markersize, marker="p")
         else :
-            plt.plot(xl,yl, marker="o")
+            print("Unforseen coding of outlof HTU : ",int(routetohtu[ih,j,i]-1))
+            sys.exit()
+        #
         if int(routetohtu[ih,j,i]-1) < nhtu :
             ile =  int(routetogrid[ih,j,i]-1)
-            ie,je = landindex.land2ij(ile)
+            je,ie = landindex.land2ij(ile)
             ihe = int(routetohtu[ih,j,i]-1)
             loe = cglon[ihe,je,ie]
             lae = cglat[ihe,je,ie]
             xle,yle=m(loe, lae)
-            plt.plot([xl,xle],[yl,yle], color="b")
+            plt.arrow(xl, yl, xle-xl, yle-yl, color="b", width=lwidth, head_width=20*lwidth, length_includes_head=True)
 #
 #
 #

run.def

-[OverAll]
-#
-#
-EarthRadius = 6370000.
-#
-ModelGridFile = /home/polcher/WORK/Data/NewRouting/geo_em.d01.nc
-# Mallorca
-WEST_EAST = 2.3, 3.5
-SOUTH_NORTH = 39.00, 40.1
-HydroFile = /home/polcher/WORK/Data/NewRouting/routing_MED.nc
-#
-# FORTRAN interface parameters
-#
-Documentation = true
-#
-# Configuration for the graph to be generated
-#
-nbasmax = 35
-#
-# Output
-#
-GraphFile = MEDCORDEX_test_graph.nc
-#
-# Diagnostics
-# You need to provide an interval in longitude and Latitude.
-#
-DiagLon = 2.3, 3.5
-DiagLat = 39.0, 40.1