From 8e7b5a108529f8520f109b60b635e6678a485718 Mon Sep 17 00:00:00 2001
From: POLCHER Jan <jan.polcher@lmd.jussieu.fr>
Date: Fri, 22 May 2020 14:59:51 +0200
Subject: [PATCH] Adding some tools to manage the hydrological network file.
---
Tools/AddExtremes.py | 45 ++++++++++
Tools/AddFacDisto.py | 193 +++++++++++++++++++++++++++++++++++++++++++
2 files changed, 238 insertions(+)
create mode 100644 Tools/AddExtremes.py
create mode 100644 Tools/AddFacDisto.py
diff --git a/Tools/AddExtremes.py b/Tools/AddExtremes.py
new file mode 100644
index 0000000..7ae4e4a
--- /dev/null
+++ b/Tools/AddExtremes.py
@@ -0,0 +1,45 @@
+#
+#
+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()
diff --git a/Tools/AddFacDisto.py b/Tools/AddFacDisto.py
new file mode 100644
index 0000000..22a5572
--- /dev/null
+++ b/Tools/AddFacDisto.py
@@ -0,0 +1,193 @@
+#
+# 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()
+
--
GitLab