Adding some tools to manage the hydrological network file.

Tools/AddExtremes.py

+#
+#
+import numpy as np
+import os
+from netCDF4 import Dataset
+#
+filename="routing_MED.nc"
+#
+#################################################################################
+#
+def ncvarextrema(varin) :
+    expval = None
+    for (a,v) in varin.__dict__.items():
+        if ( a == "_FillValue" ) :
+            expval = v
+        elif (a == "missing_value" ) :
+            expval = v
+            
+    if expval == None :
+        minval = np.min(varin[:,:])
+        maxval = np.max(varin[:,:])
+    else :
+        x = varin[:,:]
+        minval = np.min(x[np.where(x < expval)])
+        maxval = np.max(x[np.where(x < expval)])
+
+    return minval, maxval
+#
+#
+#
+def addextrema(varin, vmin, vmax) :
+    exec('varin.%s = "%s"' % ('min',vmin))
+    exec('varin.%s = "%s"' % ('max',vmax))
+#
+#
+#
+nf=Dataset(filename, 'a', format='NETCDF4_CLASSIC')
+#
+vars=["basins", "topoind_h", "topoind", "fac", "disto"]
+
+for v in vars :
+    vmin, vmax = ncvarextrema(nf.variables[v])
+    addextrema(nf.variables[v], vmin, vmax)
+
+nf.close()

Tools/AddFacDisto.py

+#
+# Code to add flow accumulation and distance to ocean to the old routing.nc file
+#
+# fac is probably the accumulation of upstream grid boxes which flow through a given point.
+#
+import numpy as np
+import os, sys
+from netCDF4 import Dataset
+#
+rose=np.array([[0,0],[-1,0],[-1,+1],[0,+1],[+1,+1],[+1,0],[+1,-1],[0,-1],[-1,-1]])
+debug=False
+FillValue=1.0e+20
+#
+infile="routing.nc"
+outfile="routing_GLOB.nc"
+#
+# Copy variable
+def copyvar (innc, outnc, var, MinMax=False):
+    #
+    if ( debug ) : print("Working ", var)
+    shape=innc.variables[var].shape
+    #
+    # Scalar
+    #
+    if (len(shape) == 0 ) :
+        dtyp=innc.variables[var].datatype
+        outvar = outnc.createVariable(var,dtyp,())
+        outvar[0] = innc.variables[var][:]
+    #
+    # Vector
+    #
+    elif (len(shape) == 1 ) :
+        dtyp=innc.variables[var].datatype
+        ddim=innc.variables[var].dimensions
+        outvar = outnc.createVariable(var,dtyp,ddim)
+        outvar[:] = innc.variables[var][:]
+    #
+    # 2D matrix
+    #
+    elif (len(shape) == 2 ) :
+        dtyp=innc.variables[var].datatype
+        ddim=innc.variables[var].dimensions
+        outvar = outnc.createVariable(var,dtyp,ddim,fill_value=FillValue)
+        x = np.copy(innc.variables[var][:,:])
+        x[np.greater_equal(x, FillValue, where=~np.isnan(x))] = np.nan
+        outvar[:,:] = x
+        #
+        if MinMax :
+            outvar.min = np.nanmin(x)
+            outvar.max = np.nanmax(x)
+    #
+    else :
+        print('A 3D var is not yet foreseen in copyvar')
+        sys.exit(2)
+    #
+    # Transfer attributes
+    #
+    copyatt(innc.variables[var], outvar)
+    #
+    # Corrections to attributes
+    #
+    outvar.axis="YX"
+    outvar.coordinates="nav_lat, nav_lon"
+#
+# Function to copy attributes
+#
+#################################################################################
+
+# Function to copy attributes
+#
+#################################################################################
+#
+# Function which copied the attributes from one file netCDF file to another.
+#
+def copyatt(varin, varout):
+    for (a,v) in varin.__dict__.items():
+        if ( a != "_FillValue" ) :
+            if isinstance(v, str) : 
+                vv=v.replace('\n','').replace('_',' ')
+            else :
+                vv=v
+            if ( debug ) : print("Attribute ", a, "--", vv)
+            exec('varout.%s = "%s"' % (a,vv))
+#
+# Function to copy global attributes
+#
+def copyglobatt(innc, outnc):
+    for (a,v) in innc.__dict__.items():
+        if ( a != "_FillValue" ) :
+            if isinstance(v, str) : 
+                vv=v.replace('\n','').replace('_',' ')
+            else :
+                vv=v
+            if ( debug ) : print("Global Attribute ", a, "--", vv, type(vv))
+            setattr(outnc, a, vv)
+#
+#
+###################################################################################
+#
+#
+#
+###################################################################################
+#
+inf=Dataset(infile,'r')
+out=Dataset(outfile, 'w', format='NETCDF4_CLASSIC')
+#
+jdim,idim = inf.variables["nav_lon"].shape
+#
+# Creat output dimensions
+#
+xdim = out.createDimension('x', idim)
+ydim = out.createDimension('y', jdim)
+#
+#
+#
+lon=np.copy(inf.variables["nav_lon"][:,:])
+copyvar(inf, out, "nav_lon")
+lat=np.copy(inf.variables["nav_lat"][:,:])
+copyvar(inf, out, "nav_lat")
+#
+trip=np.copy(inf.variables["trip"][:,:])
+copyvar(inf, out, "trip")
+#
+basins=np.copy(inf.variables["basins"][:,:])
+copyvar(inf, out, "basins", MinMax=True)
+#
+topoind=np.copy(inf.variables["topoind"][:,:])
+copyvar(inf, out, "topoind", MinMax=True)
+#
+hdiff=np.copy(inf.variables["hdiff"][:,:])
+copyvar(inf, out, "hdiff", MinMax=True)
+#
+riverl=np.copy(inf.variables["riverl"][:,:])
+copyvar(inf, out, "riverl", MinMax=True)
+#
+orog=np.copy(inf.variables["orog"][:,:])
+copyvar(inf, out, "orog", MinMax=True)
+#
+varorog=np.copy(inf.variables["varorog"][:,:])
+copyvar(inf, out, "varorog", MinMax=True)
+#
+#
+# Compute the new variables
+#
+trip[np.greater_equal(trip, FillValue, where=~np.isnan(trip))] = np.nan
+fac=np.zeros((jdim,idim))
+disto=np.zeros((jdim,idim))
+disto[:,:]=np.nan
+for i in range(idim):
+    for j in range(jdim):
+        iflow=i
+        jflow=j
+        distotmp=0
+        while trip[jflow,iflow] > 0 and ~np.isnan(trip[jflow,iflow]) and trip[jflow,iflow] < 97:
+            # Add grid to flow accumulation
+            fac[jflow,iflow]=fac[jflow,iflow]+1
+            distotmp=distotmp+riverl[jflow,iflow]
+            # Move on
+            inci=rose[int(trip[jflow,iflow])][1]
+            incj=rose[int(trip[jflow,iflow])][0]
+            iflow=(iflow+inci)%idim
+            jflow=jflow+incj
+        #
+        # Save distance to ocean
+        #
+        if ~np.isnan(trip[j,i]) :
+            disto[j,i] = distotmp
+#
+# Add new fields to file
+#
+facnc = out.createVariable('fac',np.float,('y','x'),fill_value=FillValue)
+facnc.long_name="Flow accumulation"
+facnc.associate="nav_lat nav_lon"
+facnc.coordinates="nav_lat, nav_lon"
+facnc.units=""
+facnc.min=np.nanmin(fac)
+facnc.max=np.nanmax(fac)
+facnc[:,:] = fac[:,:]
+distonc = out.createVariable('disto',np.float,('y','x'),fill_value=FillValue)
+distonc.long_name="Distance to the ocean"
+distonc.associate="nav_lat nav_lon"
+distonc.coordinates="nav_lat, nav_lon"
+distonc.axis="YX"
+distonc.units="m"
+distonc.min=np.nanmin(disto)
+distonc.max=np.nanmax(disto)
+distonc[:,:] = disto[:,:]
+#
+# Close files
+#
+inf.close()
+out.close()
+