A first attempt to plot the graph of the river flow.

Diags/plotgraphs.py

+import matplotlib
+matplotlib.use('Agg')
+import matplotlib.pyplot as plt
+from mpl_toolkits.basemap import Basemap, cm
+import numpy as np
+from netCDF4 import Dataset
+
+NCfile="../MEDCORDEX_test_graph.nc"
+GraphicFile = "test.png"
+title="Test"
+#
+#####################################################################
+#
+class LandGrid :
+    def __init__(self, land) :
+        nj,ni = land.shape
+        self.indP=[]
+        self.nbland = 0
+        for i in range(ni) :
+            for j in range(nj) :
+                if (land[j,i] > 0 ) :
+                    self.nbland = self.nbland+1
+                    self.indP.append([j,i])
+        return
+                    
+    def land2ij(self, il) :
+        j,i = self.indP[il]
+        return j,i
+#
+#####################################################################
+#
+#
+nc = Dataset(NCfile, 'r')
+
+lon_full = nc.variables["lon"][:,:]
+lat_full = nc.variables["lat"][:,:]
+lon_bnd = nc.variables["lon_bnd"][:,:,:]
+lat_bnd = nc.variables["lat_bnd"][:,:,:]
+land = nc.variables["land"][:,:]
+routetogrid = nc.variables["routetogrid"][:,:,:]
+routetohtu = nc.variables["routetobasin"][:,:,:]
+cglon = nc.variables["CG_lon"][:,:,:]
+cglat = nc.variables["CG_lat"][:,:,:]
+
+nc.close()
+
+nhtu, nj, ni = routetogrid.shape
+landindex = LandGrid(land)
+print "Nb land points = ", landindex.nbland, 'Nb HTU : ', nhtu
+
+#
+# Get some average size of the grid boxes.
+#
+lon_res = np.nanmean(np.abs(np.gradient(lon_bnd, axis=0)))
+lat_res = np.nanmean(np.abs(np.gradient(lat_bnd, axis=0)))
+#
+#
+#
+fig = plt.figure(1)
+ax = fig.add_subplot(111)
+ax.set_title(title)
+ax.set_xlabel('Longitude',labelpad=20)
+ax.set_ylabel('Latitude',labelpad=20)
+ax.axis([0, 10, 0, 10])
+#
+m = Basemap(projection='cyl', resolution='i', \
+            llcrnrlon=np.nanmin(lon_bnd)-lon_res, llcrnrlat=np.nanmin(lat_bnd)-lat_res, \
+            urcrnrlon=np.nanmax(lon_bnd)+lon_res, urcrnrlat=np.nanmax(lat_bnd)+lat_res)
+
+m.drawcoastlines(linewidth=0.5)
+m.drawcountries(linewidth=0.5)
+m.drawparallels(np.arange(-180.,180.,0.1),labels=[1,0,0,0]) # draw parallels
+m.drawmeridians(np.arange(-80.,80.,0.1),labels=[0,0,0,1]) # draw meridians
+#
+cmap = 'RdBu'
+cmap = 'Pastel1'
+cmap = 'prism'
+#
+#
+#
+for i in range(ni) :
+    for j in range(nj) :
+        if land[j,i] > 0 :
+            lo = np.append(lon_bnd[:,j,i],lon_bnd[0,j,i])
+            la = np.append(lat_bnd[:,j,i],lat_bnd[0,j,i])
+            xl,yl=m(lo, la)
+            plt.plot(xl,yl, color="m")
+#
+# Draw the graph
+#
+for il in range(landindex.nbland) :
+    j,i = landindex.land2ij(il)
+    print il,"== routetogrid ==",routetogrid[:,j,i]
+    print il,"== routetohtu  ==",routetohtu[:,j,i]
+    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)
+        else :
+            plt.plot(xl,yl, marker="o")
+        if int(routetohtu[ih,j,i]-1) < nhtu :
+            ile =  int(routetogrid[ih,j,i]-1)
+            ie,je = 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.savefig(GraphicFile, dpi=450)
+plt.close()

Interface.py

@@ -177,6 +177,7 @@ class HydroGraph :
         NCFillValue=1.0e20
         vtyp=np.float64
         cornerind=[0,2,4,6]
+        nbcorners = len(cornerind)
         #
         if part.rank == 0 :
             outnf=Dataset(filename, 'w', format='NETCDF4_CLASSIC')
@@ -185,7 +186,7 @@ class HydroGraph :
             outnf.createDimension('y', globalgrid.nj)
             outnf.createDimension('land', len(procgrid.area))
             outnf.createDimension('htu', self.nbasmax)
-            outnf.createDimension('bnd', 4)
+            outnf.createDimension('bnd', nbcorners)
         #
         # Coordinates
         #
@@ -197,11 +198,6 @@ class HydroGraph :
             lon.title="Longitude"
             lon.axis="X"
             lon[:,:] = longitude[:,:]
-        # # Longitude bounds
-        # lonbnd=outnf.createVariable("lon_bnd", vtyp, ('bnd','y','x'), fill_value=NCFillValue)
-        # lonbnd.units="grid box corners degrees east"
-        # lonbnd.title="Longitude Corners"
-        # lonbnd[:,:] = np.array(procgrid.polyll)[:,cornerind,0]
         #
         # Latitude
         latitude = part.gather(procgrid.lat_full)
@@ -212,12 +208,28 @@ class HydroGraph :
             lat.title="Latitude"
             lat.axis="Y"
             lat[:] = latitude[:,:]
-        # # Latitude bounds
-        # latbnd=outnf.createVariable("lat_bnd", vtyp, ('bnd','land'), fill_value=NCFillValue)
-        # latbnd.units="grid box corners degrees north"
-        # latbnd.title="Latitude Corners"
-        # latbnd[:,:] = np.array(procgrid.polyll)[:,cornerind,1]
         #
+        # Bounds of grid box
+        #
+        llonpoly=np.zeros((nbcorners,procgrid.nbland))
+        llatpoly=np.zeros((nbcorners,procgrid.nbland))
+        for i in range(procgrid.nbland) :
+            llonpoly[:,i] = [procgrid.polyll[i][ic][0] for ic in cornerind]
+            llatpoly[:,i] = [procgrid.polyll[i][ic][1] for ic in cornerind]
+        lon_bnd = procgrid.landscatter(llonpoly)
+        lat_bnd = procgrid.landscatter(llatpoly)
+        if part.rank == 0 :
+            lonbnd=outnf.createVariable("lon_bnd", vtyp, ('bnd','y','x'), fill_value=NCFillValue)
+            lonbnd.units="grid box corners degrees east"
+            lonbnd.title="Longitude of Corners"
+            latbnd=outnf.createVariable("lat_bnd", vtyp, ('bnd','y','x'), fill_value=NCFillValue)
+            latbnd.units="grid box corners degrees north"
+            latbnd.title="Latitude of Corners"
+        else :
+            lonbnd= np.zeros((1,1,1))
+            latbnd= np.zeros((1,1,1))
+        lonbnd[:,:,:] = part.gather(lon_bnd[:,:,:])
+        latbnd[:,:,:] = part.gather(lat_bnd[:,:,:])
         #
         # Land sea mask
         #

SampleScripts/BuildHTUs.pbs

@@ -5,8 +5,8 @@
 #PBS -j oe
 #PBS -l nodes=1:ppn=42
 #PBS -l walltime=48:00:00
-#PBS -l mem=120gb
-#PBS -l vmem=120gb
+#PBS -l mem=160gb
+#PBS -l vmem=160gb
 #
 cd ${PBS_O_WORKDIR}/..
 

diagplots.py

@@ -5,6 +5,11 @@ import matplotlib.pyplot as plt
 from mpl_toolkits.basemap import Basemap, cm
 import numpy as np
 
+import getargs
+log_master, log_world = getargs.getLogger()
+INFO, DEBUG, ERROR = log_master.info, log_master.debug, log_world.error
+INFO_ALL, DEBUG_ALL = log_world.info, log_world.debug
+
 inc=np.zeros((8,2), dtype=np.int8)
 inc[0,0] = 1; inc[0,1] =0
 inc[1,0] = 1; inc[1,1] =1
@@ -204,7 +209,7 @@ def SUPERMESHPLOT(filename, loninterval, latinterval, hoverlap, hsuper, atmpolys
     lat = hoverlap.lat_bx
     basin_sz = hsuper.basin_sz
     basin_pts = hsuper.basin_pts
-    outlet = hsuper.basin_coor
+    outlet = hsuper.basin_outcoor
     #
     if np.nanmin(lon) <= max(loninterval) and np.nanmax(lon) > min(loninterval) and \
        np.nanmin(lat) <= max(latinterval) and np.nanmax(lat) > min(latinterval) :
@@ -265,6 +270,7 @@ def SUPERMESHPLOT(filename, loninterval, latinterval, hoverlap, hsuper, atmpolys
                             if ibo == ib :
                                 # Same grid box
                                 if hsuper.outflow_basin[ib,iz] < undef_int :
+                                    INFO("hsuper.outflow_basin[ib,iz] = "+str(hsuper.outflow_basin[ib,iz]))
                                     ERROR("We cannot be here as we flow into an HTU of the same grid but we do not know the HTU index")
                                 else :
                                     izo = hsuper.outflow_basin[ib,iz]-1

run.def

@@ -24,5 +24,5 @@ GraphFile = MEDCORDEX_test_graph.nc
 # Diagnostics
 # You need to provide an interval in longitude and Latitude.
 #
-DiagLon = 3.9, 3.9
-DiagLat = 40.0, 40.0
+DiagLon = 2.3, 3.5
+DiagLat = 39.0, 40.1