First implementation of the reading and placement of gauging stations or...

First implementation of the reading and placement of gauging stations or infrastructures. Tested on Mallorca with some imaginary structures.

GraphHydro.py

@@ -8,6 +8,13 @@ sys.path.append(localdir+'/F90subroutines')
 import numpy as np
 from netCDF4 import Dataset
 #
+import getargs
+config = getargs.SetupConfig()
+# Maximum error in the distance of the station in km
+MaxDistErr = config.getfloat("OverAll", "MaxDistErr", fallback=25.0)
+# Maximum error in the upstream area in %
+MaxUpstrErr = config.getfloat("OverAll", "MaxUpstrErr", fallback=10.0)
+#
 import routing_interface
 import NcHelp as NH
 import RPPtools as RPP
@@ -88,6 +95,25 @@ class HydroGraph :
     #
     #
     #
+    def position(self, modelgrid, locations) :
+        npt, nbas = self.routing_fetch.shape
+        if len(locations.id) > 0 :
+            for i,id in enumerate(locations.id) :
+                dist = RPP.distances(np.array(modelgrid.coordll)[:,0], np.array(modelgrid.coordll)[:,1], \
+                                    locations.lon[i], locations.lat[i])
+                # Calculation in km
+                disterr = (MaxDistErr-dist/1000.)/MaxDistErr
+                disterr[disterr < 0] = np.nan
+                # Calculation in km^2 and result in %
+                upstrerr = (MaxUpstrErr-np.abs(self.routing_fetch/1.e+6-locations.upstream[i])/locations.upstream[i]*100)/MaxUpstrErr
+                upstrerr[upstrerr<0] = np.nan
+                cost = np.tile(disterr,nbas).reshape((nbas,npt)).transpose()+upstrerr
+                if np.count_nonzero(~np.isnan(cost)) > 0 :
+                    ij, ib = np.unravel_index(np.nanargmax(cost),cost.shape)
+                    locations.addhtupos(i, ij, ib, np.nanmax(cost))
+    #
+    #
+    #
     def add_variable(self, outnf, procgrid, NCFillValue, part, coord, name, title, units, data, vtyp, orig_type = "float"):
         var = procgrid.landscatter(data.astype(vtyp), order='F')
 

Locations.py

+#
+#
+import numpy as np
+from netCDF4 import Dataset
+import getargs
+config = getargs.SetupConfig()
+#
+class Locations :
+    def __init__(self, bbox) :
+        self.id = []
+        self.lon = []
+        self.lat = []
+        self.upstream = []
+        self.type = []
+        structure_file = config.get("OverAll","StructuresFile",fallback=None)
+        if structure_file != None :
+            self.addstructures(structure_file, bbox)
+        grdc_file = config.get("OverAll","GRDCFile",fallback=None)
+        if grdc_file != None :
+            self.addgrdcstations(grdc_file, bbox)
+        #
+        self.ij = [-1]*len(self.id)
+        self.ib = [-1]*len(self.id)
+        self.cost = [np.nan]*len(self.id)
+    #
+    def addstructures(self, f, bbox) :
+        for line in open(f, 'r') :
+            if line.find("#") != 0 :
+                col = line.split()
+                lon = float(col[1]); lat = float(col[2])
+                if min(bbox[0]) < lon < max(bbox[0]) and min(bbox[1]) < lat < max(bbox[1]) :
+                    self.id.append(int(col[0]))
+                    self.lon.append(lon)
+                    self.lat.append(lat)
+                    self.upstream.append(float(col[3]))
+                    self.type.append(" ".join(col[4:]).lower())
+        return
+    #
+    def addgrdcstations(self, f, bbox) :
+        g = Dataset(f, 'r')
+        for i, id in enumerate(g.variables["number"][:]) :
+            if min(bbox[0]) < g.variables["lon"][i] < max(bbox[0]) and \
+               min(bbox[1]) < g.variables["lat"][i] < max(bbox[1]) :
+                self.id.append(id)
+                self.lon.append(g.variables["lon"][i])
+                self.lat.append(g.variables["lat"][i])
+                self.upstream.append(g.variables["area"][i])
+                self.type.append("gauging station")
+        g.close()
+        return
+    #
+    def addhtupos(self, i, ij, ib, cost) :
+        self.ij[i] = ij
+        self.ib[i] = ib
+        self.cost[i] = cost
+        return
+#

RPPtools.py

@@ -68,6 +68,17 @@ def boxit(cent, dlon, dlat, polyres) :
 def loladist(a,b) :
     return np.sqrt((a[0]-b[0])**2+(a[1]-b[1])**2)
 #
+# Function to compute the distance in m
+#
+def distances(lonvec, latvec, lon, lat) :
+    dLat = np.deg2rad(latvec-lat)
+    dLon = np.deg2rad(lonvec-lon)
+    lat1 = np.deg2rad(latvec)
+    lat2 = np.deg2rad(lat)
+    a = np.sin(dLat/2) * np.sin(dLat/2) + np.sin(dLon/2) * np.sin(dLon/2) * np.cos(lat1) * np.cos(lat2)
+    c = 2 * np.arctan2(np.sqrt(a), np.sqrt(1-a))
+    return c * EarthRadius
+#
 # Function to compute maximum radius of cicle around polygon
 #
 @jit(nopython = True)

RoutingPreProc.py

@@ -30,6 +30,7 @@ import HydroGrid as HG
 import Interface as IF
 import GraphHydro as GR
 from Truncate import trunc as TR
+import Locations as LO
 #
 # Logging in MPI : https://groups.google.com/forum/#!topic/mpi4py/SaNzc8bdj6U
 #
@@ -124,6 +125,20 @@ print("Rank : {0} - Basin_count After Truncate : ".format(part.rank), hsuper.bas
 
 INFO("=================== HydroGraph ====================")
 hgraph = GR.HydroGraph(nbasmax, hsuper, part, modelgrid)
+
+INFO("=================== Locate gauging stations and other structures ====================")
+loc = LO.Locations(modelgrid.box_land)
+hgraph.position(modelgrid, loc)
+
+for i, id in enumerate(loc.id) :
+    if loc.ij[i] >= 0 :
+        print(rank," --> Locations to be placed : ", id)
+        print("Locations : lon, lat, upstream : ", loc.lon[i], loc.lat[i], loc.upstream[i])
+        print("Locations : ij, ib, cost :", loc.ij[i], loc.ib[i], loc.cost[i])
+        print("Locations in HTU :", hgraph.routing_cg[loc.ij[i],loc.ib[i],:], hgraph.routing_fetch[loc.ij[i],loc.ib[i]]/1.e+6)
+        print("Locations : type : ", loc.type[i])
+
 hgraph.dumpnetcdf(OutGraphFile, gg, modelgrid, part, hydrotcst)
 
 IF.closeinterface(comm)
+

tests/Mallorca/BuildHTUs_Mallorca.pbs

@@ -4,7 +4,7 @@
 #
 #PBS -j oe
 #PBS -l nodes=1:ppn=32:Piledriver
-#PBS -l walltime=4:00:00
+#PBS -l walltime=8:00:00
 #PBS -l mem=80gb
 #PBS -l vmem=80gb
 #
@@ -17,7 +17,7 @@ source ../../Environment
 #
 # Clean-up
 #
-/bin/rm -f DocumentationInterface *.txt
+/bin/rm -f DocumentationInterface Out_*.txt Weight_*.txt
 #
 #
 # 2 Proc
@@ -40,7 +40,7 @@ fi
 #
 # 4 Proc
 #
-/bin/rm -f run.def *.txt
+/bin/rm -f run.def Out_*.txt Weight_*.txt
 cp run_E2OFD.def run.def
 mpirun -n 3 python ../../RoutingPreProc.py
 if [ $? -gt 0 ] ; then
@@ -57,7 +57,7 @@ fi
 #
 # Test MERIT on E2OFD grid
 #
-/bin/rm -f run.def *.txt
+/bin/rm -f run.def Out_*.txt Weight_*.txt
 cp run_E2OFDMERIT.def run.def
 mpirun -n 3 python ../../RoutingPreProc.py
 if [ $? -gt 0 ] ; then
@@ -74,7 +74,7 @@ fi
 #
 # More points on 32 proc
 #
-/bin/rm -f run.def *.txt
+/bin/rm -f run.def Out_*.txt Weight_*.txt
 cp run_EuroSW.def run.def
 mpirun -n 32 python ../../RoutingPreProc.py
 if [ $? -gt 0 ] ; then

tests/Mallorca/Structures.txt

+# ID   Lon(decimal)     lat      upstream area (km^2)   Nature
+1     2.7918           39.7810       763.0              Reservoir
+2     3.1861           39.5861        4.0               Run of river

tests/Mallorca/run_MEDCORDEX.def

@@ -25,3 +25,9 @@ numop = 200
 #
 GraphFile = MEDCORDEX_test_graph.nc
 #
+# File containing infrastructures to be placed.
+# Maximum error in the distance of the station in km^2
+MaxDistErr = 25.0
+# Maximum error in the upstream area in %
+MaxUpstrErr = 10.0
+StructuresFile = Structures.txt