From 49d70d0b67c9444cb5314f03427e239fbc4e5288 Mon Sep 17 00:00:00 2001
From: Anthony <anthony.schrapffer@polytechnique.fr>
Date: Mon, 4 May 2020 15:36:15 +0200
Subject: [PATCH] Tools to check the output and correct eventual fetch and
inflows errors
---
Diags/HTUs_test_and_corr.pbs | 14 ++
Diags/Makefile | 15 ++
Diags/Routing_Test_and_Correction.py | 207 +++++++++++++++++++++
Diags/hydrosuper.f90 | 259 +++++++++++++++++++++++++++
4 files changed, 495 insertions(+)
create mode 100644 Diags/HTUs_test_and_corr.pbs
create mode 100644 Diags/Makefile
create mode 100644 Diags/Routing_Test_and_Correction.py
create mode 100644 Diags/hydrosuper.f90
diff --git a/Diags/HTUs_test_and_corr.pbs b/Diags/HTUs_test_and_corr.pbs
new file mode 100644
index 0000000..81a352f
--- /dev/null
+++ b/Diags/HTUs_test_and_corr.pbs
@@ -0,0 +1,14 @@
+#!/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
+#
diff --git a/Diags/Makefile b/Diags/Makefile
new file mode 100644
index 0000000..0f90e28
--- /dev/null
+++ b/Diags/Makefile
@@ -0,0 +1,15 @@
+#
+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
diff --git a/Diags/Routing_Test_and_Correction.py b/Diags/Routing_Test_and_Correction.py
new file mode 100644
index 0000000..298d925
--- /dev/null
+++ b/Diags/Routing_Test_and_Correction.py
@@ -0,0 +1,207 @@
+#
+#
+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()
+
diff --git a/Diags/hydrosuper.f90 b/Diags/hydrosuper.f90
new file mode 100644
index 0000000..79453d7
--- /dev/null
+++ b/Diags/hydrosuper.f90
@@ -0,0 +1,259 @@
+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
+
+
+
+
--
GitLab