From 1a36416f556e0b3151fc5639aea1daa0baf3ae79 Mon Sep 17 00:00:00 2001 From: Michael Rabotin <michael.rabotin@inrae.fr> Date: Mon, 28 Feb 2022 16:52:25 +0100 Subject: [PATCH] added cluster command, newest options in param, OF domain creation --- bin/create_OF_domain.py | 68 ++++++++++++++++ bin/create_config_file.py | 13 ++- bin/hru-delin_step1.sh | 87 ++++++++++---------- bin/hru-delin_step2.sh | 100 +++++++++++------------ bin/hru-delin_step3.sh | 100 +++++++++++------------ bin/hru-delin_step4.sh | 102 +++++++++++------------- modules/hrudelin_3_hrugen.py | 136 ++++++++++++++++++++++---------- modules/hrudelin_parms_J2000.py | 84 +++++++++++++++++++- 8 files changed, 438 insertions(+), 252 deletions(-) create mode 100644 bin/create_OF_domain.py diff --git a/bin/create_OF_domain.py b/bin/create_OF_domain.py new file mode 100644 index 0000000..a195dda --- /dev/null +++ b/bin/create_OF_domain.py @@ -0,0 +1,68 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +def isint(value): + try: + int(value) + return True + except ValueError: + return False + + + +print('------------- Export HRU and reach into OpenFLUID format domain.fluidx -------------') +# export in fluidx the topology of HRU and reach + +# open hru.par +#skip 5 line +#read first column (cat) and column 13 (to_poly) and column 14 (to_reach) +#test if column 12 and colum 13 not Null or not 0.0 +fo=open('domain.fluidx',"w") +fo.write("<?xml version=\"1.0\" standalone=\"yes\"?>\n") +fo.write("<openfluid>\n") +fo.write(" <domain>\n") +fo.write(" <definition>\n") +#for loop for each class unit +#read hru.par +hru_file=open("hru.par",'r') +hru_Lines=hru_file.readlines() +for i, line in enumerate( hru_Lines): + if i>4: + stripped_line = line.strip() + line_list = stripped_line.split() + if isint(line_list[0]): + + fo.write(" <unit class=\"HRU\" ID=\"%s\" pcsorder=\"1\">\n" %(line_list[0])) + # get the length of list + if len(line_list)>12: + if float(line_list[12])>0: + fo.write(" <to class=\"HRU\" ID=\"%i\" />\n" %(int(float(line_list[12])))) + if len(line_list)>13: + if float(line_list[13])>0: + fo.write(" <to class=\"Reach\" ID=\"%s\" />\n" %(int(float(line_list[13])))) + fo.write(" </unit>\n") +hru_file.close() +#for loop for each class unit +#read reach.par +reach_file=open("reach.par",'r') +reach_Lines=reach_file.readlines() +for j, line in enumerate( reach_Lines): + if j>4: + stripped_line = line.strip() + line_list = stripped_line.split() + if isint(line_list[0]): + + fo.write(" <unit class=\"Reach\" ID=\"%s\" pcsorder=\"1\">\n" %(line_list[0])) + if float(line_list[1])>0: + fo.write(" <to class=\"Reach\" ID=\"%s\" />\n" %(line_list[1])) + fo.write(" </unit>\n") +reach_file.close() + +fo.write(" </definition>\n") +fo.write(" <calendar>\n") +fo.write(" </calendar>\n") +fo.write(" </domain>\n") +fo.write("</openfluid>\n") +fo.close() + + diff --git a/bin/create_config_file.py b/bin/create_config_file.py index 34b6687..b3c99d3 100644 --- a/bin/create_config_file.py +++ b/bin/create_config_file.py @@ -47,6 +47,7 @@ file.write("relocated_gauges=\n") file.write("\n") file.write("\n") file.write("[irrigation]\n") +file.write("# yes or no\n") file.write("to_do:\n") file.write("irrigation:\n") file.write("irrig_col_name=\n") @@ -68,6 +69,7 @@ file.write("irrigation_table:\n") file.write("relocated_irrigation=\n") file.write("\n") file.write("[dams]\n") +file.write("# yes or no\n") file.write("to_do:\n") file.write("dams=\n") file.write("dams_col_name=\n") @@ -130,6 +132,7 @@ file.write("# 2nd step : hru-delin_basins\n") file.write("# ---------------------------\n") file.write("# So it's possible to specify a variable using : or = ???\n") file.write("[auto_relocation]\n") +file.write("# yes or no\n") file.write("to_do:\n") file.write("# -------- first rule\n") file.write("# surface is in percent!\n") @@ -161,6 +164,7 @@ file.write("#\n") file.write("# MNT-derived layers to be integrated in the overlay operation\n") file.write("#\n") file.write("[layer_overlay]\n") +file.write("# yes or no\n") file.write("dem:\n") file.write("slope:\n") file.write("aspect:\n") @@ -168,8 +172,14 @@ file.write("\n") file.write("# --------------------------------\n") file.write("# 4th step : hru-delin_parms_J2000\n") file.write("# --------------------------------\n") +file.write("# yes or no\n") file.write("[topology]\n") file.write("dissolve_cycle:\n") +file.write("hru_no_topology_log:\n") +file.write("OF_domain_export:\n") +file.write("[hru_param]\n") +file.write("hru_cat:\n") +file.write("hru_landuse:\n") @@ -177,5 +187,4 @@ file.write("dissolve_cycle:\n") - -file.close() \ No newline at end of file +file.close() diff --git a/bin/hru-delin_step1.sh b/bin/hru-delin_step1.sh index adb8829..1ee92c0 100755 --- a/bin/hru-delin_step1.sh +++ b/bin/hru-delin_step1.sh @@ -25,51 +25,41 @@ function printUsage(){ scriptName=`basename $MYPATH` echo "USAGE:" echo - echo "$scriptName [-d|--debug] [config-file]" + echo "$scriptName [-d|--debug] [-p|--nbprocess] [-c|--cluster] [-f|--file]" echo - echo "config-file: path to hru-delin configuration file. DEFAULT VALUE: 'hrudelin_config.cfg'" echo "-d|--debug : enable debug output" + echo "-p|--nbprocess N : use N parallel processes (default=number of available CPU cores)" + echo "-c|--cluster : mandatory for using on HIICS cluster" + echo "-f|--file : path to hru-delin configuration file. DEFAULT VALUE: 'hrudelin_config.cfg'" } -# parse arguments -ARGS=$(getopt -o dp: -l "debug,nbprocess:" -n "$0" -- "$@"); - # allow empty arg list -#if [ $? -ne 0 ] || [ $# -eq 0 ]; -#then -# printUsage -# exit -#fi +if [ $? -ne 0 ] || [ $# -eq 0 ]; +then + printUsage + exit +fi + + DEBUG="" -eval set -- "$ARGS"; - -while true; do - case "$1" in - -d|--debug) - DEBUG=1 - shift; - ;; - -p|--nbprocess) - # useless, just to avoid argument parsing errors - NBPROC=$2 - shift; - shift; - ;; - --) - break; - ;; +CLUSTER="" +NBPROCESS="" +CONFIGFILE="" + +while getopts 'dp:cf:' flag; do + case "${flag}" in + d) DEBUG='true' ;; + p) NBPROCESS=$OPTARG ;; + c) CLUSTER='true' ;; + f) CONFIGFILE="${OPTARG}" ;; + *) printUsage + exit 1;; esac -done -# getting non-option arguments : THE MAIN ONE is the config file -# only first non-option argument will be taken -shift; -while [ -n "$1" ]; do - CONFIGFILE=$1 - shift - break done + + if [ -z "$CONFIGFILE" ]; then CONFIGFILE=hrudelin_config.cfg fi @@ -90,19 +80,22 @@ CONFIGFILEPATH=`readlink -f $CONFIGFILEPATH` grassDbPath=$CONFIGFILEPATH/grass_db rc_file=$grassDbPath/grassdata/hru-delin/.grassrc -# find grass base dir -res=`find /usr/lib/grass7* -name "r.thin"` -if [ -z "$res" ]; then - echo 'Grass7x not found, did you install grass package?' - exit 1 -else - par=`dirname "$res"` - par2=`dirname "$par"` - export GISBASE="$par2" +# find grass base dir only if no calcul on cluster +if [ -z "$CLUSTER" ]; then + res=`find /usr/lib/grass7* -name "r.thin"` + if [ -z "$res" ]; then + echo 'Grass7x not found, did you install grass package?' + exit 1 + else + par=`dirname "$res"` + par2=`dirname "$par"` + export GISBASE="$par2" + fi + export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$GISBASE/lib + export PYTHONPATH=$PYTHONPATH:$GISBASE/etc/python fi + export PATH=$PATH:$GISBASE/bin:$GISBASE/scripts -export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$GISBASE/lib -export PYTHONPATH=$PYTHONPATH:$GISBASE/etc/python # DEBUG if [ ! -z "$DEBUG" ]; then echo @@ -120,7 +113,7 @@ fi # clean work environment # ---------------------- -# test id dir output exist +# test if dir output exist clean_files=`grep "files\s*:" $CONFIGFILE | cut -d ':' -f2 | sed -e 's/\s*$//' | sed -e 's/^\s*//'` diff --git a/bin/hru-delin_step2.sh b/bin/hru-delin_step2.sh index c798bdb..a3c01ab 100755 --- a/bin/hru-delin_step2.sh +++ b/bin/hru-delin_step2.sh @@ -16,6 +16,7 @@ # HRU-DELIN # 2nd step : hru-delin_basins + MYPATH=`readlink -f $0` MYDIR=`dirname $MYPATH` . $MYDIR/tools.sh @@ -24,57 +25,45 @@ function printUsage(){ scriptName=`basename $MYPATH` echo "USAGE:" echo - echo "$scriptName [-d|--debug] [-p|--nbprocess N] [config-file]" + echo "$scriptName [-d|--debug] [-p|--nbprocess] [-c|--cluster] [-f|--file]" echo - echo "config-file: path to hru-delin configuration file. DEFAULT VALUE: 'hrudelin_config.cfg'" echo "-d|--debug : enable debug output" echo "-p|--nbprocess N : use N parallel processes (default=number of available CPU cores)" + echo "-c|--cluster : mandatory for using on HIICS cluster" + echo "-f|--file : path to hru-delin configuration file. DEFAULT VALUE: 'hrudelin_config.cfg'" } -# parse arguments -ARGS=$(getopt -o dp: -l "debug,nbprocess:" -n "$0" -- "$@"); - # allow empty arg list -#if [ $? -ne 0 ] || [ $# -eq 0 ]; -#then -# printUsage -# exit -#fi +if [ $? -ne 0 ] || [ $# -eq 0 ]; +then + printUsage + exit +fi + + DEBUG="" -NBPROC="" -eval set -- "$ARGS"; - -while true; do - case "$1" in - -d|--debug) - DEBUG=1 - shift; - ;; - -p|--nbprocess) - NBPROC=$2 - shift; - shift; - ;; - --) - break; - ;; +CLUSTER="" +NBPROCESS="" +CONFIGFILE="" + +while getopts 'dp:cf:' flag; do + case "${flag}" in + d) DEBUG='true' ;; + p) NBPROCESS="${OPTARG}" ;; + c) CLUSTER='true' ;; + f) CONFIGFILE="${OPTARG}" ;; + *) printUsage + exit 1;; esac -done -# getting non-option arguments : THE MAIN ONE is the config file -# only first non-option argument will be taken -shift; -while [ -n "$1" ]; do - CONFIGFILE=$1 - shift - break done + + + if [ -z "$CONFIGFILE" ]; then CONFIGFILE=hrudelin_config.cfg fi -CONFIGFILEPATH=`dirname $CONFIGFILE` -CONFIGFILEPATH=`readlink -f $CONFIGFILEPATH` if [ ! -f "$CONFIGFILE" ]; then echo "!!! file '$CONFIGFILE' does not exist" @@ -83,30 +72,37 @@ if [ ! -f "$CONFIGFILE" ]; then exit fi -# predict env +echo "Config file: $CONFIGFILE" +CONFIGFILEPATH=`dirname $CONFIGFILE` +# get absolute path of config file directory +CONFIGFILEPATH=`readlink -f $CONFIGFILEPATH` + +## predict the environment grassDbPath=$CONFIGFILEPATH/grass_db rc_file=$grassDbPath/grassdata/hru-delin/.grassrc -# find grass base dir -res=`find /usr/lib/grass7* -name "r.thin"` -if [ -z "$res" ]; then - echo 'Grass7x not found, did you install grass package?' - exit 1 -else - par=`dirname "$res"` - par2=`dirname "$par"` - export GISBASE="$par2" +# find grass base dir only if no calcul on cluster +if [ -z "$CLUSTER" ]; then + res=`find /usr/lib/grass7* -name "r.thin"` + if [ -z "$res" ]; then + echo 'Grass7x not found, did you install grass package?' + exit 1 + else + par=`dirname "$res"` + par2=`dirname "$par"` + export GISBASE="$par2" + fi + export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$GISBASE/lib + export PYTHONPATH=$PYTHONPATH:$GISBASE/etc/python fi + export PATH=$PATH:$GISBASE/bin:$GISBASE/scripts -export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$GISBASE/lib -export PYTHONPATH=$PYTHONPATH:$GISBASE/etc/python # DEBUG if [ ! -z "$DEBUG" ]; then - echo "config file path: '$CONFIGFILE'" echo echo "execute following lines to be able to call grass tools from this terminal:" echo - echo export GISRC=$grassDbPath/grassdata/hru-delin/.grassrc + echo export GISRC=$rc_file echo export GISBASE="$par2" echo export PATH=$PATH:$GISBASE/bin:$GISBASE/scripts echo export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$GISBASE/lib @@ -162,4 +158,4 @@ fi # ----------------------------------- # exec second step of HRU-DELIN batch # ----------------------------------- -python3 $MYDIR/../modules/hrudelin_2_basins.py $CONFIGFILE $NBPROC +python3 $MYDIR/../modules/hrudelin_2_basins.py $CONFIGFILE $NBPROCESS diff --git a/bin/hru-delin_step3.sh b/bin/hru-delin_step3.sh index 06cbaba..00bb262 100755 --- a/bin/hru-delin_step3.sh +++ b/bin/hru-delin_step3.sh @@ -15,6 +15,7 @@ # HRU-DELIN # 3rd step : hru-delin_hrugen + MYPATH=`readlink -f $0` MYDIR=`dirname $MYPATH` . $MYDIR/tools.sh @@ -23,58 +24,45 @@ function printUsage(){ scriptName=`basename $MYPATH` echo "USAGE:" echo - echo "$scriptName [-d|--debug] [-p|--nbprocess N] [config-file]" + echo "$scriptName [-d|--debug] [-p|--nbprocess] [-c|--cluster] [-f|--file]" echo - echo "config-file: path to hru-delin configuration file. DEFAULT VALUE: 'hrudelin_config.cfg'" echo "-d|--debug : enable debug output" echo "-p|--nbprocess N : use N parallel processes (default=number of available CPU cores)" + echo "-c|--cluster : mandatory for using on HIICS cluster" + echo "-f|--file : path to hru-delin configuration file. DEFAULT VALUE: 'hrudelin_config.cfg'" } -# parse arguments -ARGS=$(getopt -o dp: -l "debug,nbprocess:" -n "$0" -- "$@"); - # allow empty arg list -#if [ $? -ne 0 ] || [ $# -eq 0 ]; -#then -# printUsage -# exit -#fi +if [ $? -ne 0 ] || [ $# -eq 0 ]; +then + printUsage + exit +fi + + DEBUG="" -NBPROC="" -eval set -- "$ARGS"; - -while true; do - case "$1" in - -d|--debug) - DEBUG=1 - shift; - ;; - -p|--nbprocess) - NBPROC=$2 - shift; - shift; - ;; - --) - break; - ;; +CLUSTER="" +NBPROCESS="" +CONFIGFILE="" + +while getopts 'dp:cf:' flag; do + case "${flag}" in + d) DEBUG='true' ;; + p) NBPROCESS="${OPTARG}" ;; + c) CLUSTER='true' ;; + f) CONFIGFILE="${OPTARG}" ;; + *) printUsage + exit 1;; esac -done -# getting non-option arguments : THE MAIN ONE is the config file -# only first non-option argument will be taken -shift; -while [ -n "$1" ]; do - CONFIGFILE=$1 - shift - break done + + if [ -z "$CONFIGFILE" ]; then CONFIGFILE=hrudelin_config.cfg fi -CONFIGFILEPATH=`dirname $CONFIGFILE` -CONFIGFILEPATH=`readlink -f $CONFIGFILEPATH` if [ ! -f "$CONFIGFILE" ]; then echo "!!! file '$CONFIGFILE' does not exist" @@ -83,37 +71,43 @@ if [ ! -f "$CONFIGFILE" ]; then exit fi -# predict env +echo "Config file: $CONFIGFILE" +CONFIGFILEPATH=`dirname $CONFIGFILE` +# get absolute path of config file directory +CONFIGFILEPATH=`readlink -f $CONFIGFILEPATH` + +## predict the environment grassDbPath=$CONFIGFILEPATH/grass_db rc_file=$grassDbPath/grassdata/hru-delin/.grassrc -# find grass base dir -res=`find /usr/lib/grass7* -name "r.thin"` -if [ -z "$res" ]; then - echo 'Grass7x not found, did you install grass package?' - exit 1 -else - par=`dirname "$res"` - par2=`dirname "$par"` - export GISBASE="$par2" +# find grass base dir only if no calcul on cluster +if [ -z "$CLUSTER" ]; then + res=`find /usr/lib/grass7* -name "r.thin"` + if [ -z "$res" ]; then + echo 'Grass7x not found, did you install grass package?' + exit 1 + else + par=`dirname "$res"` + par2=`dirname "$par"` + export GISBASE="$par2" + fi + export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$GISBASE/lib + export PYTHONPATH=$PYTHONPATH:$GISBASE/etc/python fi + export PATH=$PATH:$GISBASE/bin:$GISBASE/scripts -export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$GISBASE/lib -export PYTHONPATH=$PYTHONPATH:$GISBASE/etc/python # DEBUG if [ ! -z "$DEBUG" ]; then - echo "config file path: '$CONFIGFILE'" echo echo "execute following lines to be able to call grass tools from this terminal:" echo - echo export GISRC=$grassDbPath/grassdata/hru-delin/.grassrc + echo export GISRC=$rc_file echo export GISBASE="$par2" echo export PATH=$PATH:$GISBASE/bin:$GISBASE/scripts echo export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$GISBASE/lib echo export PYTHONPATH=$PYTHONPATH:$GISBASE/etc/python echo fi - # ---------------------- # clean work environment # ---------------------- @@ -161,4 +155,4 @@ fi # ----------------------------------- # exec third step of HRU-DELIN batch # ----------------------------------- -python3 $MYDIR/../modules/hrudelin_3_hrugen.py $CONFIGFILE $NBPROC +python3 $MYDIR/../modules/hrudelin_3_hrugen.py $CONFIGFILE $NBPROCESS diff --git a/bin/hru-delin_step4.sh b/bin/hru-delin_step4.sh index 1d79a43..4615fb8 100755 --- a/bin/hru-delin_step4.sh +++ b/bin/hru-delin_step4.sh @@ -16,6 +16,7 @@ # HRU-DELIN # 4th step : hru-delin_parms_J2000 + MYPATH=`readlink -f $0` MYDIR=`dirname $MYPATH` . $MYDIR/tools.sh @@ -24,58 +25,45 @@ function printUsage(){ scriptName=`basename $MYPATH` echo "USAGE:" echo - echo "$scriptName [-d|--debug] [-p|--nbprocess N] [config-file]" + echo "$scriptName [-d|--debug] [-p|--nbprocess] [-c|--cluster] [-f|--file]" echo - echo "config-file: path to hru-delin configuration file. DEFAULT VALUE: 'hrudelin_config.cfg'" echo "-d|--debug : enable debug output" echo "-p|--nbprocess N : use N parallel processes (default=number of available CPU cores)" + echo "-c|--cluster : mandatory for using on HIICS cluster" + echo "-f|--file : path to hru-delin configuration file. DEFAULT VALUE: 'hrudelin_config.cfg'" } -# parse arguments -ARGS=$(getopt -o dp: -l "debug,nbprocess:" -n "$0" -- "$@"); - # allow empty arg list -#if [ $? -ne 0 ] || [ $# -eq 0 ]; -#then -# printUsage -# exit -#fi +if [ $? -ne 0 ] || [ $# -eq 0 ]; +then + printUsage + exit +fi + + DEBUG="" -NBPROC="" -eval set -- "$ARGS"; - -while true; do - case "$1" in - -d|--debug) - DEBUG=1 - shift; - ;; - -p|--nbprocess) - NBPROC=$2 - shift; - shift; - ;; - --) - break; - ;; +CLUSTER="" +NBPROCESS="" +CONFIGFILE="" + +while getopts 'dp:cf:' flag; do + case "${flag}" in + d) DEBUG='true' ;; + p) NBPROCESS="${OPTARG}" ;; + c) CLUSTER='true' ;; + f) CONFIGFILE="${OPTARG}" ;; + *) printUsage + exit 1;; esac -done -# getting non-option arguments : THE MAIN ONE is the config file -# only first non-option argument will be taken -shift; -while [ -n "$1" ]; do - CONFIGFILE=$1 - shift - break done + + if [ -z "$CONFIGFILE" ]; then CONFIGFILE=hrudelin_config.cfg fi -CONFIGFILEPATH=`dirname $CONFIGFILE` -CONFIGFILEPATH=`readlink -f $CONFIGFILEPATH` if [ ! -f "$CONFIGFILE" ]; then echo "!!! file '$CONFIGFILE' does not exist" @@ -84,37 +72,43 @@ if [ ! -f "$CONFIGFILE" ]; then exit fi -# predict env +echo "Config file: $CONFIGFILE" +CONFIGFILEPATH=`dirname $CONFIGFILE` +# get absolute path of config file directory +CONFIGFILEPATH=`readlink -f $CONFIGFILEPATH` + +## predict the environment grassDbPath=$CONFIGFILEPATH/grass_db rc_file=$grassDbPath/grassdata/hru-delin/.grassrc -# find grass base dir -res=`find /usr/lib/grass7* -name "r.thin"` -if [ -z "$res" ]; then - echo 'Grass7x not found, did you install grass package?' - exit 1 -else - par=`dirname "$res"` - par2=`dirname "$par"` - export GISBASE="$par2" +# find grass base dir only if no calcul on cluster +if [ -z "$CLUSTER" ]; then + res=`find /usr/lib/grass7* -name "r.thin"` + if [ -z "$res" ]; then + echo 'Grass7x not found, did you install grass package?' + exit 1 + else + par=`dirname "$res"` + par2=`dirname "$par"` + export GISBASE="$par2" + fi + export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$GISBASE/lib + export PYTHONPATH=$PYTHONPATH:$GISBASE/etc/python fi + export PATH=$PATH:$GISBASE/bin:$GISBASE/scripts -export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$GISBASE/lib -export PYTHONPATH=$PYTHONPATH:$GISBASE/etc/python # DEBUG if [ ! -z "$DEBUG" ]; then - echo "config file path: '$CONFIGFILE'" echo echo "execute following lines to be able to call grass tools from this terminal:" echo - echo export GISRC=$grassDbPath/grassdata/hru-delin/.grassrc + echo export GISRC=$rc_file echo export GISBASE="$par2" - echo export PATH=$MYDIR/../modules:$PATH:$GISBASE/bin:$GISBASE/scripts + echo export PATH=$PATH:$GISBASE/bin:$GISBASE/scripts echo export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$GISBASE/lib echo export PYTHONPATH=$PYTHONPATH:$GISBASE/etc/python echo fi - # ---------------------- # clean work environment # ---------------------- @@ -150,4 +144,4 @@ rm -f topologie_* tmp/topologie_* # ----------------------------------- # exec last step of HRU-DELIN batch # ----------------------------------- -python3 $MYDIR/../modules/hrudelin_parms_J2000.py $CONFIGFILE $NBPROC +python3 $MYDIR/../modules/hrudelin_parms_J2000.py $CONFIGFILE $NBPROCESS diff --git a/modules/hrudelin_3_hrugen.py b/modules/hrudelin_3_hrugen.py index eb2a0b0..7b28530 100755 --- a/modules/hrudelin_3_hrugen.py +++ b/modules/hrudelin_3_hrugen.py @@ -44,27 +44,22 @@ import pandas as pd import multiprocessing from multiprocessing import Pool, cpu_count -from utils import isint - - - -try: - # Python 3 - from subprocess import DEVNULL -except ImportError: - DEVNULL = open(os.devnull, 'wb') - -MY_ABS_PATH=os.path.abspath(__file__) -MY_DIR=os.path.dirname(MY_ABS_PATH) +def isint(value): + try: + int(value) + return True + except ValueError: + return False def reclass(map_in, map_out, size): p = grass_pipe_command('r.stats', flags='lnNc', input=map_in) + pReclass = grass_feed_command('r.reclass', overwrite=True, input=map_in, output=map_out, rules='-') for l in p.stdout: lSpl = decode(l).rstrip(os.linesep).split() @@ -79,6 +74,7 @@ def reclass(map_in, map_out, size): def count_only(map_in, size): statLines = decode(grass_read_command('r.stats', quiet=True, flags='lnNc', input=map_in)).rstrip(os.linesep).split(os.linesep) c = 0 + for l in statLines: lSpl = l.split() if int(lSpl[1]) <= int(size): @@ -106,8 +102,20 @@ def getAreasUpTo(map_in, size): return areas +try: + # Python 3 + from subprocess import DEVNULL +except ImportError: + DEVNULL = open(os.devnull, 'wb') + +MY_ABS_PATH=os.path.abspath(__file__) +MY_DIR=os.path.dirname(MY_ABS_PATH) + + + # function launched in parallel, yeyeah def processSubbasin(params): + id, configFileDir, iRun, nbBasins, min_area, buffer_distance, tmpPath, generator, nbProc = params # define which grass env we are using current = multiprocessing.current_process() @@ -116,13 +124,16 @@ def processSubbasin(params): processN = (processN % nbProc) + 1 location = 'hru-delin_%s' % (processN) os.environ['GISRC'] = os.path.join(configFileDir, 'grass_db', 'grassdata', location, '.grassrc') - + grass_run_command('r.mask', raster='subbasins_msk', maskcats=id.rstrip('\n'), overwrite='True', stdout=DEVNULL, stderr=DEVNULL) grass_run_command('g.region', raster='clumps', stdout=DEVNULL, stderr=DEVNULL) grass_run_command('g.region', zoom='MASK', stdout=DEVNULL, stderr=DEVNULL) + # Get size of subbasin subbasin_size = int(decode(grass_read_command('r.stats', quiet=True, flags='nNc', input='MASK')).rstrip(os.linesep).split(os.linesep)[0].split()[1]) + if (subbasin_size <= min_area): + if not generator: iterable3 = tqdm( total=2, @@ -146,6 +157,7 @@ def processSubbasin(params): iterable3.update(1) iterable3.close() else: + # LOOP THAT ELIMINATES PIXELS ONLY sizes = range(2, min_area + 1) if generator: @@ -164,13 +176,18 @@ def processSubbasin(params): (processN, padLeft(str(iRun), len(str(nbBasins))), nbBasins, pad(id.strip(), 6), pad('', 1))) initmap = 'clumps' - counter_old = -1 + counter_old = 0 + while (True): + counter = count_only(initmap, 1) - if counter != counter_old: + + if counter != counter_old : reclass(initmap, 'clumps_sz1', 1) + counter_old = counter - grass_run_command('r.mapcalc', expression='base=clumps_sz1', overwrite='True', stdout=DEVNULL, stderr=DEVNULL) + grass_run_command('r.mapcalc', expression='base=clumps_sz1', overwrite='True',stdout=DEVNULL, stderr=DEVNULL) + # MERGE HORIZONTAL/VERTICAL ADJACENT PIXELS grass_run_command('r.mapcalc', expression='map_a=if((not(isnull(base))&¬(isnull(base[0,1]))),base[0,1],if((not(isnull(base))&¬(isnull(base[1,0]))),base[1,0],base))', @@ -179,39 +196,54 @@ def processSubbasin(params): input='map_a,clumps_cp', output='newmap', overwrite='True', stdout=DEVNULL, stderr=DEVNULL) grass_run_command('r.mapcalc', expression='clumps_cp=newmap', overwrite='True', stdout=DEVNULL, stderr=DEVNULL) initmap = 'newmap' + + else: + + grass_run_command('g.copy', raster='clumps,clumps_sz1') + grass_run_command('r.mapcalc', expression='newmap=clumps_cp', overwrite='True', stdout=DEVNULL, stderr=DEVNULL) break # START TO ELIMINATE REMAINING AND ISOLATED PIXELS - grass_run_command('r.mapcalc', expression='newmap_sz1=clumps_sz1', overwrite='True', stdout=DEVNULL, stderr=DEVNULL) + grass_run_command('r.mapcalc', expression='newmap_sz1=clumps_sz1', overwrite='True',stdout=DEVNULL, stderr=DEVNULL) grass_run_command('r.buffer', input='newmap_sz1', output='buffer', distances='%d'%buffer_distance, - overwrite='True', stdout=DEVNULL, stderr=DEVNULL) - grass_run_command('r.mapcalc', expression='newmap2=if((buffer==2),newmap,null())', overwrite='True', stdout=DEVNULL, stderr=DEVNULL) + overwrite='True',stdout=DEVNULL, stderr=DEVNULL) + + grass_run_command('r.mapcalc', expression='newmap2=if((buffer==2),newmap,null())', overwrite='True',stdout=DEVNULL, stderr=DEVNULL) + # there were sort/awk system calls here # now less filesystem access, cleaner code, most importantely: portable to windows # pandatana ! stats1 = decode(grass_read_command('r.stats', quiet=True, flags='cnN', input='newmap')).rstrip(os.linesep) + stats2 = decode(grass_read_command('r.stats', quiet=True, flags='nN', input='newmap2')).rstrip(os.linesep) + stio1 = StringIO(stats1) - stio2 = StringIO(stats2) - df1 = pd.read_csv(stio1, header=None, dtype={0: int, 1: int}, sep=' ') - df2 = pd.read_csv(stio2, header=None, dtype={0: int}, sep=' ') - df3 = df1.merge(df2, on=0, sort=True) + stio2 = StringIO(stats2) + + if stats2: + df2 = pd.read_csv(stio2, header=None, dtype={0: int}, sep=' ') + df3 = df1.merge(df2, on=0, sort=True) + else: + df3=df1 + + + #out3Path = os.path.join(tmpPath, 'out3_proc%s' % processN) pReclass = grass_feed_command('r.reclass', input='newmap', output='test', rules='-', overwrite='True') for index, row in df3.iterrows(): pReclass.stdin.write(encode('%s = %s %s\n' % (row[0], row[0], row[1]))) pReclass.stdin.close() pReclass.wait() - + #os.system('cat %s/out3_proc%s | r.reclass --o input=newmap output=test rules=-' % (tmpPath, processN)) grass_run_command('g.region', raster='clumps', stdout=DEVNULL, stderr=DEVNULL) grass_run_command('r.buffer', input='MASK', output='buf_mask', distances='%d'%buffer_distance, overwrite='True', stdout=DEVNULL, stderr=DEVNULL) grass_run_command('g.remove', flags='f', type='raster', name='MASK', stdout=DEVNULL, stderr=DEVNULL) grass_run_command('g.region', zoom='buf_mask', stdout=DEVNULL, stderr=DEVNULL) - + #os.system('echo "2 = 2 9999999" | r.reclass --o input=buf_mask output=buf_mask_new rules=-') pReclass = grass_feed_command('r.reclass', input='buf_mask', output='buf_mask_new', rules='-', @@ -219,20 +251,22 @@ def processSubbasin(params): pReclass.stdin.write(encode('2 = 2 9999999\n')) pReclass.stdin.close() pReclass.wait() - + grass_run_command('r.patch', input='buf_mask_new,test', output='sum_buf', overwrite='True', stdout=DEVNULL, stderr=DEVNULL) grass_run_command('r.mapcalc', expression='newmap3=if((isnull(newmap_sz1)),newmap,eval(a=@sum_buf[0,-1],b=@sum_buf[-1,0],c=@sum_buf[0,1],d=@sum_buf[1,0],m=min(a,b,c,d),if((m==int(a)),sum_buf[0,-1],if((m==int(b)),sum_buf[-1,0],if((m==int(c)),sum_buf[0,1],if((m==int(d)),sum_buf[1,0]))))))', overwrite='True', stdout=DEVNULL, stderr=DEVNULL) - + # START ELIMINATING NOW FOR AREAS GREATER THAN 1 PIXEL grass_run_command('r.mask', raster='subbasins', maskcats=id.rstrip('\n'), overwrite='True', stdout=DEVNULL, stderr=DEVNULL) for sz in iterable2: + if not generator: iterable2.set_description('[process %s] basin [%s/%s] id %s group size %s' % (processN, padLeft(str(iRun), len(str(nbBasins))), nbBasins, pad(id.strip(), 6), pad(str(sz-1), 4)) ) + if sz not in getAreasUpTo('newmap3', min_area): pass # I should create a MASK here to prevent the MASK removal to send an error later in the code... @@ -240,7 +274,9 @@ def processSubbasin(params): else: initmap = 'newmap3' counter_old = -1 + while (True): + counter = count_only(initmap, sz) if (counter != counter_old): reclass(initmap, 'clumps_sz2', sz) @@ -268,6 +304,7 @@ def processSubbasin(params): input='base_new,newmap3', output='newout', overwrite='True', stdout=DEVNULL, stderr=DEVNULL) # TODO there was a mapcalc here in v4 which was removed in v5...why? initmap = 'newmap3' + else: # ELIMINATION STOPPED break @@ -283,7 +320,7 @@ def processSubbasin(params): p.wait() pReclass.stdin.close() pReclass.wait() - + reclass('newmap3', 'newmap_sz2', sz) grass_run_command('r.buffer', input='newmap_sz2', output='buffer', distances='%d'%buffer_distance, overwrite='True', stdout=DEVNULL, stderr=DEVNULL) @@ -293,7 +330,7 @@ def processSubbasin(params): pReclass.stdin.write(encode('1 thru 2 = 1\n')) pReclass.stdin.close() pReclass.wait() - + grass_run_command('r.mapcalc', expression='b1=if((buffer_new==1&&isnull(newmap_sz2)),newmap_sz2[-1,0],null())', overwrite='True', stdout=DEVNULL, stderr=DEVNULL) grass_run_command('r.mapcalc', expression='b2=if((buffer_new==1&&isnull(newmap_sz2)),newmap_sz2[1,0],null())', overwrite='True', stdout=DEVNULL, stderr=DEVNULL) grass_run_command('r.mapcalc', expression='b3=if((buffer_new==1&&isnull(newmap_sz2)),newmap_sz2[0,-1],null())', overwrite='True', stdout=DEVNULL, stderr=DEVNULL) @@ -320,21 +357,25 @@ def processSubbasin(params): p.wait() pReclass.stdin.close() pReclass.wait() - + grass_run_command('r.patch', input='newmap_sz2_rec,newmap3', output='newout2', overwrite='True', stdout=DEVNULL, stderr=DEVNULL) + grass_run_command('r.mapcalc', expression='newmap3=newout2', overwrite='True', stdout=DEVNULL, stderr=DEVNULL) + # SOLVE STRANGE PROBLEM OF EDGED WATERSHEDS IF EXISTS grass_run_command('r.mapcalc', expression='map_a=if((not(isnull(newmap3))&&isnull(newmap3[-1,0])&&isnull(newmap3[1,0])&&'+ 'isnull(newmap3[0,-1])&&isnull(newmap3[0,1])),-1,newmap3)', overwrite='True', stdout=DEVNULL, stderr=DEVNULL) + check = int(decode(grass_read_command('r.stats', quiet=True, flags='lnN', input='map_a')).rstrip(os.linesep).split(os.linesep)[0].strip()) if (check != -1): # NOTHING TO DO pass else: + grass_run_command('r.mapcalc', expression='map_b=if((map_a==-1),if((map_a[-1,1]>-1|||map_a[1,1]>-1|||map_a[1,-1]>-1|||map_a[-1,-1]>-1),-2,null()),map_a)', overwrite='True', stdout=DEVNULL, stderr=DEVNULL) @@ -354,7 +395,7 @@ def processSubbasin(params): pReclass.stdin.write(encode('1 thru 2 = 1\n')) pReclass.stdin.close() pReclass.wait() - + grass_run_command('g.remove', flags='f', type='raster', name='MASK', stdout=DEVNULL, stderr=DEVNULL) grass_run_command('r.clump', input='map_d3', output='map_d4', @@ -362,7 +403,7 @@ def processSubbasin(params): grass_run_command('r.mapcalc', expression='map_d5=if((not(isnull(map_d4[-1,0]))&¬(isnull(map_d4[1,0]))&¬(isnull(map_d4[0,-1]))&¬(isnull(map_d4[0,1]))),map_d4,null())', overwrite='True', stdout=DEVNULL, stderr=DEVNULL) - + grass_run_command('r.statistics', base='map_d5', cover='map_c', method='max', output='map_f', overwrite='True', stdout=DEVNULL, stderr=DEVNULL) @@ -375,20 +416,23 @@ def processSubbasin(params): grass_run_command('r.mapcalc', expression='newmap3=map_h', overwrite='True', stdout=DEVNULL, stderr=DEVNULL) - + grass_run_command('g.remove', flags='f', type='raster', name='MASK', stdout=DEVNULL, stderr=DEVNULL) - + grass_run_command('g.region', raster='clumps', stdout=DEVNULL, stderr=DEVNULL) + grass_run_command('r.patch', input='result,newmap3', output='result_tmp', overwrite='True', stdout=DEVNULL, stderr=DEVNULL) + grass_run_command('r.mapcalc', expression='result=result_tmp', overwrite='True', stdout=DEVNULL, stderr=DEVNULL) - + grass_run_command('r.mapcalc', expression='newmap3=null()', overwrite='True', stdout=DEVNULL, stderr=DEVNULL) + grass_run_command('r.mapcalc', expression='clumps_cp=clumps', overwrite='True', stdout=DEVNULL, stderr=DEVNULL) @@ -428,7 +472,7 @@ def main(parms_file, nbProc, generator=False): for data in parms.items('data'): data = os.path.join(directory, data[1]) if not os.path.isfile(data): - print(data) + sys.exit('------------> ERROR : Input data not found' ) @@ -582,6 +626,7 @@ def main(parms_file, nbProc, generator=False): print('----------------------------- Extracting subbasins IDs to use as units in loop') grass_run_command('r.mask', raster='clumps', overwrite='True', stdout=DEVNULL, stderr=DEVNULL) basin_ids = decode(grass_read_command('r.stats', quiet=True, flags='nN', input='subbasins')).rstrip(os.linesep).split(os.linesep) + grass_run_command('r.mapcalc', expression='clumps_cp=clumps', overwrite='True', stdout=DEVNULL, stderr=DEVNULL) @@ -619,31 +664,42 @@ def main(parms_file, nbProc, generator=False): nbElems = len(elems) # the locks are here to prevent concurrent terminal tqdm writing # this is the interesting part, launching N processes in parallel to process basins + if generator: print('Starting subbasins loop with %s process' % nbProc) + with Pool(nbProc) as p: + params = [(id, configFileDir, i+1, nbElems, min_area, buffer_distance, tmpPath, generator, nbProc) for (i, id) in enumerate(elems)] + for i, _ in enumerate(p.imap_unordered(processSubbasin, params), 1): + yield (i/nbElems*100) + + else: + with Pool(nbProc, initializer=tqdm.set_lock, initargs=(tqdm.get_lock(),)) as p: + params = [(id, configFileDir, i+1, nbElems, min_area, buffer_distance, tmpPath, generator, nbProc) for (i, id) in enumerate(elems)] + r = list(tqdm(p.imap_unordered(processSubbasin, params), desc='[main process] Loop on subbasins [%s process] ' % nbProc, total=nbElems, unit='basin', bar_format=bar_format1 )) - print('\n' * nbProc, end='') - + + # export rasters from parallel grass environments to files for i in range(nbProc): location = 'hru-delin_%s' % (i+1) + rastersToExport = { 'result': os.path.join(tmpPath, 'step3_result_%s.tif' % (i+1)) } exportRastersFromEnv(rastersToExport, grassDbPath, location) - + # restore main grass env os.environ['GISRC'] = MAIN_GISRC @@ -660,7 +716,7 @@ def main(parms_file, nbProc, generator=False): grass_run_command('r.mapcalc', expression='result=result_tmp', overwrite='True', stdout=DEVNULL, stderr=DEVNULL) - + print('') timestr = time.strftime("%a, %d %b %Y, %H:%M:%S", time.localtime()) print('---------- HRU-delin Step 3 : Reassign category values of HRUs --- '+str(timestr)) diff --git a/modules/hrudelin_parms_J2000.py b/modules/hrudelin_parms_J2000.py index 36cfee4..fdf7231 100755 --- a/modules/hrudelin_parms_J2000.py +++ b/modules/hrudelin_parms_J2000.py @@ -43,7 +43,7 @@ from osgeo.gdalconst import * from osgeo import ogr import pandas as pd import json - +from utils import isint import multiprocessing from multiprocessing import Pool, cpu_count @@ -634,7 +634,7 @@ def main(parms_file, nbProc, generator=False): header_unit += '\t' + ranges['area'][2] sql += ',area' - if (parms.get('layer_overlay', 'dem')) == 'x': + if (parms.get('layer_overlay', 'dem')) == 'yes': dem_avg='dem_avg' dem_int='dem_int' grass_run_command('r.mapcalc', expression='dem_int=int(dem)', overwrite='True') @@ -654,7 +654,7 @@ def main(parms_file, nbProc, generator=False): else: yield 83 - if parms.get('layer_overlay', 'slope') == 'x': + if parms.get('layer_overlay', 'slope') == 'yes': slope_mult = 'slope_mult' slope_labeled='slope_labeled' slope_avg = 'slope_avg' @@ -688,7 +688,7 @@ def main(parms_file, nbProc, generator=False): else: yield 85 - if (parms.get('layer_overlay', 'aspect')) == 'x': + if (parms.get('layer_overlay', 'aspect')) == 'yes': grass_run_command('r.statistics', base='hrus', cover='aspect_rcl', method='mode', output='aspect_md', overwrite='True') grass_run_command('r.mapcalc', expression='aspect_mj=if((aspect_rcl==@aspect_md),aspect,null())', overwrite='True') grass_run_command('r.mapcalc', expression='aspect_mjj=if((aspect_mj>315.0),aspect_mj-360.0,aspect_mj)', overwrite='True') @@ -942,6 +942,14 @@ def main(parms_file, nbProc, generator=False): + if (parms.get('topology', 'hru_no_topology_log')) == 'yes': + print('------------- Export HRU with no topology in hru_with_no_toology.par -------------') + # export in csv the hru with no topology (no downstream reach and hru) + notopology_file = os.path.join(dir_results, 'hru_with_no_toology.par') + grass_run_command('v.db.select', flags='c', map='hrus_c_pnt', columns='cat', where='to_reach IS NULL AND to_poly IS NULL', file=notopology_file, overwrite=True) + + + ##optim @@ -977,6 +985,74 @@ def main(parms_file, nbProc, generator=False): p.write(footer) + if (parms.get('topology', 'OF_domain_export')) == 'yes': + print('------------- Export HRU and reach into OpenFLUID format domain.fluidx -------------') + # export in fluidx the topology of HRU and reach + + # open hru.par + #skip 5 line + #read first column (cat) and column 13 (to_poly) and column 14 (to_reach) + #test if column 12 and colum 13 not Null or not 0.0 + domain_file = os.path.join(dir_results, 'domain.fluidx') + fo=open(domain_file,"w") + fo.write("<?xml version=\"1.0\" standalone=\"yes\"?>\n") + fo.write("<openfluid>\n") + fo.write(" <domain>\n") + fo.write(" <definition>\n") + #for loop for each class unit + #read hru.par + hru_file=open(par,'r') + hru_Lines=hru_file.readlines() + for i, line in enumerate( hru_Lines): + if i>4: + stripped_line = line.strip() + line_list = stripped_line.split() + if isint(line_list[0]): + + fo.write(" <unit class=\"HRU\" ID=\"%s\" pcsorder=\"1\">\n" %(line_list[0])) + # get the length of list + if len(line_list)>12: + if float(line_list[12])>0: + fo.write(" <to class=\"HRU\" ID=\"%i\" />\n" %(int(float(line_list[12])))) + if len(line_list)>13: + if float(line_list[13])>0: + fo.write(" <to class=\"Reach\" ID=\"%s\" />\n" %(int(float(line_list[13])))) + fo.write(" </unit>\n") + hru_file.close() + #for loop for each class unit + #read reach.par + reach_file=open(reach_par,'r') + reach_Lines=reach_file.readlines() + for j, line in enumerate( reach_Lines): + if j>4: + stripped_line = line.strip() + line_list = stripped_line.split() + if isint(line_list[0]): + + fo.write(" <unit class=\"Reach\" ID=\"%s\" pcsorder=\"1\">\n" %(line_list[0])) + if float(line_list[1])>0: + fo.write(" <to class=\"Reach\" ID=\"%s\" />\n" %(line_list[1])) + fo.write(" </unit>\n") + reach_file.close() + + fo.write(" </definition>\n") + fo.write(" <calendar>\n") + fo.write(" </calendar>\n") + fo.write(" </domain>\n") + fo.write("</openfluid>\n") + fo.close() + + if (parms.get('hru_param', 'hru_cat')) == 'yes': + grass_run_command('r.out.gdal', + input=hrus, output=os.path.join(dir_results, 'hru_cat.tif'), + overwrite='True', stdout=DEVNULL, stderr=DEVNULL) + + if (parms.get('hru_param', 'hru_landuse')) == 'yes': + grass_run_command('r.mapcalc', expression='hrulanduse = if(hrus,landuse)', overwrite='True') + grass_run_command('r.out.gdal', + input='hrulanduse', output=os.path.join(dir_results, 'hru_landuse.tif'), + overwrite='True', stdout=DEVNULL, stderr=DEVNULL) + save_config = os.path.join(directory_out, 'hrudelin_config_save.cfg') #os.system('cp "%s" "%s" ' % (parms_file, save_config)) shutil.copyfile(parms_file, save_config) -- GitLab