Tools to check the output and correct eventual fetch and inflows errors

Diags/HTUs_test_and_corr.pbs

+#!/bin/bash
+#
+#PBS -N HTUs_test
+#
+#PBS -j oe
+#PBS -l nodes=1:ppn=1
+#PBS -l walltime=1:00:00
+#PBS -l mem=40gb
+#PBS -l vmem=40gb
+#
+cd ${PBS_O_WORKDIR}
+#
+python3 Routing_Test_and_Correction.py
+#

Diags/Makefile

+#
+FC = gfortran
+FFLAG = -fPIC
+FDBG = -g -fbounds-check
+#
+FPY = f2py3
+PDBG =  --debug
+#
+all : hydrosuper.so
+
+clean :
+	rm -f compilation.log *.o *.so *.mod hydrosuper.so *~ 
+
+hydrosuper.so : hydrosuper.f90
+	$(FPY) -c hydrosuper.f90 $(PDBG) --fcompiler=gfortran -m hydrosuper > compilation.log

Diags/Routing_Test_and_Correction.py

+#
+#
+import sys
+sys.path.append("/home/anthony/Documents/LOCAL/RoutingPP/ROUTING_Analysis_tools/Programs")
+import numpy as np
+from netCDF4 import Dataset as NetCDFFile
+from hydrosuper import check_fetch_order, check_inflows_in_outflow, check_outflow_in_inflows, test_outflow, compute_fetch, compute_inflows
+#
+############################
+#
+# Direction of the netCDF graph file obtained from RoutingPP
+input_file = ""
+
+# Direction to save the corrected netCDF graph file
+output_file = ""
+#
+############################
+#
+def searchInlist(listname, nameFind):
+    """ Function to search a value within a list
+    listname = list
+    nameFind = value to find
+    >>> searInlist(['1', '2', '3', '5'], '5')
+    True
+    """
+    for x in listname:
+      if x == nameFind:
+        return True
+    return False
+#
+class routing_test:
+   def __init__(self, dire):
+      self.nc = NetCDFFile(dire, "r")
+      self.nbpt_glo = self.nc.variables["nbpt_glo"][:]
+
+      self.num_pt = int(np.nanmax(self.nbpt_glo))
+      self.conv_land_ji = [[np.where(self.nbpt_glo == i)[0][0], np.where(self.nbpt_glo == i)[1][0]] for i in range(1, self.num_pt+1)]
+
+      self.basin_count_land = self.var1d("routenbintobas")
+      self.max_bas_count = int(np.max(self.basin_count_land))
+
+      self.HTUoutgrid_land = self.var2d("routetogrid")
+      self.HTUoutbasin_land = self.var2d("routetobasin")
+      self.basin_area_land = self.var2d('basin_area')
+      self.fetch_land = self.var2d("fetch")
+      self.HTUinnum = self.var2d("route_innum")
+      self.max_innum = int(np.max(self.HTUinnum))
+
+      self.HTUingrid = self.var3d("route_ingrid")
+      self.HTUinbas = self.var3d("route_inbasin")
+
+   def var1d(self, varn):
+       var = self.nc.variables[varn][:]
+       output = np.array([var[j,i] for j,i in self.conv_land_ji])
+       return output
+
+   def var2d(self, varn):
+       var = self.nc.variables[varn][:]
+       output = np.array([var[:self.max_bas_count,j,i] for j,i in self.conv_land_ji])
+       return output
+
+   def var3d(self, varn):
+       var = self.nc.variables[varn][:]
+       output = np.zeros((self.num_pt,self.max_bas_count, self.max_innum))
+       for l, ji in enumerate(self.conv_land_ji):
+           j,i = ji
+           output[l,:,:] = np.transpose(var[:self.max_innum,:self.max_bas_count,j,i])
+       return output
+
+   ##################
+
+   def check_fetch_order(self):
+       check_fetch_order(nbpt = self.num_pt, nbbas = self.max_bas_count, fetch = self.fetch_land, \
+            outgrid = self.HTUoutgrid_land, outbas = self.HTUoutbasin_land, basin_count = self.basin_count_land)
+
+   def check_inflows_in_outflow(self):
+       check_inflows_in_outflow(nbpt = self.num_pt, nbbas = self.max_bas_count, outgrid = self.HTUoutgrid_land, outbas = self.HTUoutbasin_land, \
+             basin_count = self.basin_count_land, ingrid = self.HTUingrid, inbas = self.HTUinbas, innum = self.HTUinnum)
+
+   def check_outflow_in_inflows(self):
+       check_outflow_in_inflows(nbpt = self.num_pt, nbbas = self.max_bas_count, outgrid = self.HTUoutgrid_land, outbas = self.HTUoutbasin_land, \
+             basin_count = self.basin_count_land, ingrid = self.HTUingrid, inbas = self.HTUinbas, innum = self.HTUinnum)
+
+   def test_outflow(self):
+       test_outflow(nbpt = self.num_pt, nbbas = self.max_bas_count, outgrid = self.HTUoutgrid_land, outbas = self.HTUoutbasin_land, \
+             basin_count = self.basin_count_land)
+
+   def compute_fetch(self):
+       self.fetch_new = compute_fetch(nbpt = self.num_pt, nbbas = self.max_bas_count, outgrid = self.HTUoutgrid_land, outbas = self.HTUoutbasin_land, \
+             basin_count = self.basin_count_land, basin_area = self.basin_area_land)
+
+
+   def compute_inflows(self):
+       self.ingrid_new, self.inbas_new, self.innum_new = compute_inflows(nbpt = self.num_pt, nbbas = self.max_bas_count, nbin = self.max_innum,outgrid = self.HTUoutgrid_land, outbas = self.HTUoutbasin_land, \
+             basin_count = self.basin_count_land)
+
+   def full_test(self):
+     self.check_fetch_order()
+     self.check_inflows_in_outflow()
+     self.check_outflow_in_inflows()
+     self.test_outflow()
+
+   ####
+   def get_var_corr(self, X):
+       if X.ndim == 2:
+           A = np.zeros(self.nc.variables["fetch"][:].shape)
+           for k,ji in enumerate(self.conv_land_ji):
+               j,i = ji 
+               A[:self.max_bas_count,j,i] = X[k,:]
+
+       if X.ndim == 3:
+           A = np.zeros(self.nc.variables["route_ingrid"][:].shape)
+           for k,ji in enumerate(self.conv_land_ji):
+               j,i = ji     
+               A[:self.max_innum, :self.max_bas_count,j,i] = np.transpose(X[k,:,:])
+
+       return A
+
+   def map_var(self):
+       self.D_corr = {}
+       rr.compute_fetch()
+       rr.compute_inflows()
+
+       self.D_corr["fetch"] = self.get_var_corr(self.fetch_new)
+       self.D_corr["route_innum"] = self.get_var_corr(self.innum_new)
+       self.D_corr["route_ingrid"] = self.get_var_corr(self.ingrid_new)
+       self.D_corr["route_inbasin"] = self.get_var_corr(self.inbas_new)
+
+
+   def netcdf_corr(self, ofile):
+      self.map_var()
+
+      with NetCDFFile(ofile, 'w') as onew: 
+         # Dimensions
+         for dimn in self.nc.dimensions:
+             if not searchInlist(onew.dimensions,dimn):
+                 odim = self.nc.dimensions[dimn]
+                 if odim.isunlimited():
+                     onew.createDimension(dimn, None)
+                 else:
+                     onew.createDimension(dimn, len(odim))
+         # Variables
+         for varn in self.nc.variables:
+             print(varn + " ....")
+             if not searchInlist(onew.variables,varn):
+                 ovar = self.nc.variables[varn]
+                 newvar = onew.createVariable(varn, ovar.dtype, ovar.dimensions)
+                 for attrn in ovar.ncattrs():
+                     attrv = ovar.getncattr(attrn)
+                     newvar.setncattr(attrn,attrv)
+
+                 if np.prod(ovar.shape) > 300*300*100:
+                     for it in range(0,ovar.shape[0], 5):
+                          print("  ", iit, ',', eit)
+                          iit = it
+                          eit = it + 5
+                          if varn in list(self.D_corr.keys()):
+                              newvar[iit:eit,] = self.D_corr[varn][iit:eit,]
+                          else:
+                              newvar[iit:eit,] = ovar[iit:eit,]
+                     if eit != ovar.shape[0]:
+                          print("  ", iit, ',', eit)
+                          iit = eit
+                          eit = ovar.shape[0]
+                          if varn in list(self.D_corr.keys()):
+                              newvar[iit:eit,] = self.D_corr[varn][iit:eit,]
+                          else:
+                              newvar[iit:eit,] = ovar[iit:eit,]
+                 else:
+                     if varn in list(self.D_corr.keys()):
+                         newvar[:] = self.D_corr[varn][:]
+                     else:
+                         newvar[:] = ovar[:]
+
+                 onew.sync()
+          
+         # global attributes
+         for attrn in self.nc.ncattrs():
+             attrv = self.nc.getncattr(attrn)
+             onew.setncattr(attrn,attrv)
+         onew.sync()
+
+
+##################
+#
+# Class test
+if __name__ == "__main__":
+    """
+    Run the program
+    """
+    if input_file == "":
+        "Error: please inform an input file ! "
+    else:
+        if output_file == "": 
+            output_file = input_file.replace(".nc", "_corr.nc")
+        filename = input_file.split("/")[-1]
+        print("-- {0} --".format(filename))
+        print("Testing the input file")
+        rr = routing_test(input_file)
+        rr.full_test()
+    
+        print("Generation of the corrected file")
+        rr.netcdf_corr(output_file)
+        print("Testing the corrected file")
+        rr2 = routing_test(output_file)
+        rr2.full_test()
+

Diags/hydrosuper.f90

+SUBROUTINE check_fetch_order(nbpt, nbbas, fetch, outgrid, outbas, basin_count)
+
+  IMPLICIT NONE
+
+  !! INPUT VARIABLES
+  INTEGER, INTENT(in) :: nbpt, nbbas
+  REAL,DIMENSION(nbpt, nbbas), INTENT(in) :: fetch
+  INTEGER, DIMENSION(nbpt, nbbas), INTENT(in) :: outgrid, outbas
+  INTEGER, DIMENSION(nbpt), INTENT(in) :: basin_count
+
+  !! LOCAL VARIABLES
+  INTEGER :: ig, ib, ibmx
+  REAL :: f1, f2
+  
+
+  WRITE(*,*) " *** Checking fetch order errors *** "
+  DO ig=1,nbpt
+    ibmx = basin_count(ig)
+    DO ib=1,ibmx
+      IF (outbas(ig,ib) .LE. 35) THEN
+        f1 = fetch(ig,ib)
+        f2 = fetch(outgrid(ig,ib), outbas(ig,ib))
+        IF (f1 .GT. f2) THEN 
+          WRITE(*,*) "*ERROR*", ig, ib, f1, f2
+          WRITE(*,*) "outflows : ", outgrid(ig,ib), outbas(ig,ib)
+        END IF
+      END IF
+    END DO
+  END DO 
+
+END SUBROUTINE check_fetch_order
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+SUBROUTINE check_inflows_in_outflow(nbpt, nbbas, nbin, outgrid, outbas, basin_count, ingrid, inbas, innum)
+
+  IMPLICIT NONE
+
+  !! INPUT VARIABLES
+  INTEGER, INTENT(in) :: nbpt, nbbas, nbin
+  INTEGER, DIMENSION(nbpt, nbbas), INTENT(in) :: outgrid, outbas, innum
+  INTEGER, DIMENSION(nbpt, nbbas, nbin), INTENT(in) :: ingrid, inbas
+  INTEGER, DIMENSION(nbpt), INTENT(in) :: basin_count
+
+  !! LOCAL VARIABLES
+  INTEGER :: ig, ib, ibmx, og, ob, test
+  INTEGER :: iin, iinmx, ii, limiin
+
+
+  WRITE(*,*) " *** Checking upstream HTUs is in the inflows *** "
+
+  DO ig=1,nbpt
+    ibmx = basin_count(ig)
+    DO ib=1,ibmx
+      og = outgrid(ig,ib)
+      ob = outbas(ig,ib)
+      IF ( ob .LE. 35) THEN
+        iinmx = innum(og,ob)
+        test = 0
+        ! Look through the inflows
+        DO iin=1,iinmx
+          IF (test .EQ. 0) THEN
+              IF ((ingrid(og,ob,iin) .EQ. ig) .AND. (inbas(og,ob,iin) .EQ. ib)) test = 1
+              IF (ingrid(og,ob,iin) .EQ. 0) THEN
+                test = -1
+                limiin = iin-1
+              END IF
+          END IF
+        ENDDO
+        ! Error : not present
+        IF (test .EQ. 0) THEN 
+          WRITE(*,*) "Error at", ig, ib
+          WRITE(*,*) "Outflow:", og, ob, "with", iinmx, "inflows"
+          DO ii=1,iinmx
+            WRITE(*,*) ii, ingrid(og,ob,ii), inbas(og,ob,ii)
+          END DO
+        END IF
+        ! Error, less inflows than innum (and not present)
+        IF (test .EQ. -1) WRITE(*,*) "innum does not correspond (iin, iinmax)", limiin, iinmx
+      END IF 
+
+    END DO
+  END DO 
+
+END SUBROUTINE check_inflows_in_outflow
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+
+SUBROUTINE check_outflow_in_inflows(nbpt, nbbas, nbin, outgrid, outbas, basin_count, ingrid, inbas, innum)
+
+  IMPLICIT NONE
+
+  !! INPUT VARIABLES
+  INTEGER, INTENT(in) :: nbpt, nbbas, nbin
+  INTEGER, DIMENSION(nbpt, nbbas), INTENT(in) :: outgrid, outbas, innum
+  INTEGER, DIMENSION(nbpt, nbbas, nbin), INTENT(in) :: ingrid, inbas
+  INTEGER, DIMENSION(nbpt), INTENT(in) :: basin_count
+
+  !! LOCAL VARIABLES
+  INTEGER :: ig, ib, ibmx, test, ing, inb
+  INTEGER :: iinmx, ii
+
+
+  WRITE(*,*) " *** Checking if the inflows corresponds to the upstream outflows  *** "
+
+  DO ig=1,nbpt
+    ibmx = basin_count(ig)
+    DO ib=1,ibmx
+      iinmx = innum(ig,ib)
+      IF ( iinmx .GT. 0) THEN
+        test = 0
+        DO ii=1,iinmx
+          ing = ingrid(ig,ib,ii)
+          inb = inbas(ig,ib,ii)
+          IF ((outgrid(ing,inb) .NE. ig) .OR. (outbas(ing,inb) .NE. ib)) THEN
+            WRITE(*,*) "Error at", ig, ib
+            WRITE(*,*) "outflow->", outgrid(ing,inb), outbas(ing,inb)
+          END IF
+        ENDDO
+      END IF 
+    END DO
+  END DO 
+
+END SUBROUTINE check_outflow_in_inflows
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+
+SUBROUTINE test_outflow(nbpt, nbbas, outgrid, outbas, basin_count)
+
+  IMPLICIT NONE
+  ! Changes to do : 
+  ! Optimize by keeping in memory the points already checked 
+
+  !! INPUT VARIABLES
+  INTEGER, INTENT(in) :: nbpt, nbbas
+  INTEGER, DIMENSION(nbpt, nbbas), INTENT(in) :: outgrid, outbas
+  INTEGER, DIMENSION(nbpt), INTENT(in) :: basin_count
+
+  !! LOCAL VARIABLES
+  INTEGER :: ig, ib, ibmx, og, ob, ind
+  INTEGER :: it
+
+  WRITE(*,*) " *** Checking all HTUs goes to ocean *** "
+
+  DO ig=1,nbpt
+    ibmx = basin_count(ig)
+    DO ib=1,ibmx
+      ind = 0
+      og = outgrid(ig,ib)
+      ob = outbas(ig,ib)
+      DO WHILE (ob .LE. 35)
+        it = outgrid(og,ob)
+        ob = outbas(og,ob)
+        og = it
+        ind = ind+1
+        IF (ind .GE. 5000) THEN
+          og = 0
+          WRITE(*,*) "ERROR", ig, ib
+        END IF
+      END DO
+    END DO
+  END DO
+
+
+END SUBROUTINE test_outflow
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+
+
+SUBROUTINE compute_fetch(nbpt, nbbas, outgrid, outbas, basin_count, basin_area, fetch_new)
+
+  IMPLICIT NONE
+
+  !! INPUT VARIABLES
+  INTEGER, INTENT(in) :: nbpt, nbbas
+  INTEGER, DIMENSION(nbpt, nbbas), INTENT(in) :: outgrid, outbas
+  REAL, DIMENSION(nbpt, nbbas), INTENT(in) :: basin_area
+  INTEGER, DIMENSION(nbpt), INTENT(in) :: basin_count
+
+  !! OUTPUT VARIABLES
+  REAL, DIMENSION(nbpt,nbbas), INTENT(out) :: fetch_new
+
+  !! LOCAL VARIABLES
+  INTEGER :: ig, ib, ibmx, og, ob, ind
+  INTEGER :: it
+  REAL :: area
+
+  WRITE(*,*) " *** Computing the fetch *** "
+
+  fetch_new(:,:) = 0
+
+  DO ig=1,nbpt
+    ibmx = basin_count(ig)
+    DO ib=1,ibmx
+      ind = 0
+      area = basin_area(ig,ib)
+      fetch_new(ig,ib) = fetch_new(ig,ib) + area
+
+      og = outgrid(ig,ib)
+      ob = outbas(ig,ib)
+      DO WHILE (ob .LE. 35)
+        fetch_new(og,ob) = fetch_new(og,ob) + area
+        it = outgrid(og,ob)
+        ob = outbas(og,ob)
+        og = it
+      END DO
+    END DO
+  END DO
+
+
+END SUBROUTINE compute_fetch
+
+
+
+SUBROUTINE compute_inflows(nbpt, nbbas, nbin, outgrid, outbas, basin_count, ingrid_new, inbas_new, innum_new)
+
+  IMPLICIT NONE
+
+  !! INPUT VARIABLES
+  INTEGER, INTENT(in) :: nbpt, nbbas, nbin
+  INTEGER, DIMENSION(nbpt, nbbas), INTENT(in) :: outgrid, outbas
+  INTEGER, DIMENSION(nbpt), INTENT(in) :: basin_count
+
+  INTEGER, DIMENSION(nbpt, nbbas), INTENT(OUT) :: innum_new
+  INTEGER, DIMENSION(nbpt, nbbas, nbin), INTENT(OUT) :: ingrid_new, inbas_new
+
+  !! LOCAL VARIABLES
+  INTEGER :: ig, ib, ibmx, og, ob
+  INTEGER :: iin
+
+  innum_new(:,:)  = 0
+  ingrid_new(:,:,:) = 0
+  inbas_new(:,:,:)  = 0
+
+  WRITE(*,*) " *** Computing the inflows  *** "
+
+  DO ig=1,nbpt
+    ibmx = basin_count(ig)
+    DO ib=1,ibmx
+      og = outgrid(ig,ib)
+      ob = outbas(ig,ib)
+      IF ( ob .LE. 35) THEN
+        innum_new(og,ob) = innum_new(og,ob) +1
+        iin = innum_new(og,ob)
+        ingrid_new(og,ob,iin) = ig
+        inbas_new(og,ob,iin) = ib
+      END IF
+    END DO
+  END DO 
+
+END SUBROUTINE compute_inflows
+
+
+
+