The structure of the data in the routing-graph file has changed for the locations to be monitored.

There is now a list of stations with all the indices of where they are to be found n the routing graph and the monitoring variables.

GraphHydro.py

@@ -164,7 +164,7 @@ class HydroGraph :
         cornerind=[0,2,4,6]
         nbcorners = len(cornerind)
         #
-        nlmax = locations.maxlocations(self.nbasmax, part)
+        nlmax, nblocated = locations.maxlocations(self.nbasmax, part)
         #
         if part.rank == 0 :
             outnf=Dataset(filename, 'w', format='NETCDF4_CLASSIC')
@@ -175,7 +175,8 @@ class HydroGraph :
             outnf.createDimension('z', self.nbasmax)
             outnf.createDimension('bnd', nbcorners)
             outnf.createDimension('inflow', self.max_inflow)
-            outnf.createDimension('stations', nlmax)
+            outnf.createDimension('stnperhtu', nlmax)
+            outnf.createDimension('locations', nblocated)
         else :
             outnf = None
         #
@@ -258,7 +259,7 @@ class HydroGraph :
         #
         # Build and save a map of all the structures located on the graph
         #
-        locations.addtonetcdf(self.nbpt, self.nbasmax, outnf, procgrid, NCFillValue, part, ('stations','z','y','x'), vtyp)
+        locations.addtonetcdf(self.nbpt, self.nbasmax, outnf, procgrid, part, ('y','x'), ('locations',))
         #
         # Close the file
         #

Locations.py

@@ -122,33 +122,44 @@ class Locations :
             u,n = np.unique(ghl[~np.isnan(ghl)], return_counts=True)
             if len(u) > 0 :
                 nlmax = max(n)
+                # Maximum number of HTU to monitor per grid
+                iu = [int(ii) for ii in u]
+                uhtu,nhtu = np.unique(iu, return_counts=True)
+                htumonmax=max(nhtu)
             else :
                 nlmax = 1
+                htumonmax = 1
         else :
             nlmax = None
+            htumonmax = None
         nlmax = part.bcast(nlmax)
-        return nblocated, nlmax
+        htumonmax = part.bcast(htumonmax)
+        return nblocated, nlmax, htumonmax 
     #
     # Function to get the maximum number of locations per HTU
     #
     def maxlocations(self, nbasmax, part) :
         self.cleanup()
-        self.nblocated, self.nlmax = self.unique(nbasmax, part)
+        self.nblocated, self.nlmax, self.htumonmax = self.unique(nbasmax, part)
         INFO("Out of a total of "+str(self.gnbloc)+" stations "+str(self.nblocated)+" could be located on the HTU graph.")
         INFO("Maximum number of stations per HTU : "+str(self.nlmax))
-        return self.nlmax
+        INFO("Maximum number of HTU to monitor per grid : "+str(self.htumonmax))
+        return self.nlmax, self.nblocated
     #
     # Function to place all IDs of structures onto a map and remove duplicated locations
     #
-    def addtonetcdf(self, nbpt, nbasmax, outnf, procgrid, NCFillValue, part, coord, vtyp) :
+    def addtonetcdf(self, nbpt, nbasmax, outnf, procgrid, part, spacecoord, loccoord) :
         #
-        name = "locations"
-        title = "Locations of stations or structures"
-        units = "-"
+        locname = "locations"
+        loctitle = "Basic information of located stations or structures"
+        locunits = "-"
+        monname = "HTUmonitor"
+        montitle = "Index of HTU to be monitored in the model"
+        monunits = "-"
         #
-        locmap = np.zeros((nbpt,nbasmax,self.nlmax), dtype=vtyp, order='F')
+        locmap = np.zeros((nbpt,nbasmax,self.nlmax), dtype=np.float64, order='F')
         locmap[:,:,:] = np.nan
-        costmap = np.zeros((nbpt,nbasmax,self.nlmax), dtype=vtyp, order='F')
+        costmap = np.zeros((nbpt,nbasmax,self.nlmax), dtype=np.float64, order='F')
         costmap[:,:,:] = np.nan
         #
         for i,lid in enumerate(self.lid) :
@@ -161,6 +172,11 @@ class Locations :
         gl = part.gather(procgrid.landscatter(locmap, order='F'))
         gc = part.gather(procgrid.landscatter(costmap, order='F'))
         #
+        htumon = np.zeros((self.htumonmax,part.njg,part.nig), dtype=int)
+        htumon[:,:,:] = -1
+        monnb = np.zeros((part.njg,part.nig), dtype=int)
+        locationlist = np.zeros((5,self.nblocated), dtype=int)
+        #
         # Eliminate locations with lower cost
         #
         if part.rank == 0:
@@ -170,16 +186,46 @@ class Locations :
                 maxcost = np.max(gc[gl == lid])
                 gl[(gl == lid) & (gc < maxcost)] = np.nan
                 gc[(gl == lid) & (gc < maxcost)] = np.nan
+            #
+            # Build the list of all IDs which are now unique
+            #
+            u = np.unique(gl[~np.isnan(gl)])
+            for i,uid in enumerate(u) :
+                f = np.where(gl == uid)
+                if np.any(htumon[:,f[2][0],f[3][0]] == f[1][0]+1) :
+                    # HTU already monitored
+                    imom = np.where(htumon[:,f[2][0],f[3][0]] == f[1][0]+1)[0][0]
+                else :
+                    # Add HTU to those monitored
+                    imon = monnb[f[2][0],f[3][0]]
+                    htumon[imon,f[2][0],f[3][0]] = f[1][0]+1
+                    monnb[f[2][0],f[3][0]] += 1
+                # Information saved.
+                locationlist[0,i] = uid
+                locationlist[1,i] = f[3][0]+1
+                locationlist[2,i] = f[2][0]+1
+                locationlist[3,i] = f[1][0]+1
+                locationlist[4,i] = imon+1
         #
         # Add to NetCDF file
         #
         if part.rank == 0:
-            ncvar = outnf.createVariable(name, vtyp, coord, fill_value=NCFillValue)
-            ncvar.title = title
-            ncvar.units = units
-            gl[np.isnan(gl)] = NCFillValue
-        else :
-            ncvar = np.zeros((1,)*len(coord))
-        ncvar[:,:,:,:] = gl
+            # Create new dimension
+            monmax = np.max(monnb)
+            moncoord = ('htumon',)
+            outnf.createDimension(moncoord[0], monmax)
+            infocoord = ('locinfo',)
+            outnf.createDimension(infocoord[0], 5)
+            # Create variable for the list of stations
+            locvar = outnf.createVariable(locname, 'i4',  infocoord+loccoord)
+            locvar.title = loctitle
+            locvar.units = locunits
+            locvar.description = "ID, i-index, j-index, HTU-index, monitoring-index"
+            locvar[:,:] = locationlist
+            # Create variable for the indices of the HTU to monitor
+            monvar = outnf.createVariable(monname, 'i4', moncoord+spacecoord)
+            monvar.title = montitle
+            monvar.units= monunits
+            monvar[:,:,:] = htumon[0:monmax,:,:]
         return
 #