diff --git a/bin/hru-delin_step2.sh b/bin/hru-delin_step2.sh index a3c01ab276fd4804da3d3a0f1b0650dc3ed2a850..1797ef99c9b8389c889d5ffedbd8c26fe85ddc86 100755 --- a/bin/hru-delin_step2.sh +++ b/bin/hru-delin_step2.sh @@ -113,49 +113,58 @@ fi # ---------------------- # clean work environment # ---------------------- -clean_files=`grep "files\s*:" $CONFIGFILE | cut -d ':' -f2 | sed -e 's/\s*$//' | sed -e 's/^\s*//'` -if [ -z "$clean_files" ]; then - clean_files=`grep "files\s*=" $CONFIGFILE | cut -d '=' -f2 | sed -e 's/\s*$//' | sed -e 's/^\s*//'` -fi +#clean_files=`grep "files\s*:" $CONFIGFILE | cut -d ':' -f2 | sed -e 's/\s*$//' | sed -e 's/^\s*//'` +#if [ -z "$clean_files" ]; then +# clean_files=`grep "files\s*=" $CONFIGFILE | cut -d '=' -f2 | sed -e 's/\s*$//' | sed -e 's/^\s*//'` +#fi -if [ -z "$clean_files" ]; then - echo "------------> ERROR : Output FILE Directory not provided !" - exit 1 -fi +#if [ -z "$clean_files" ]; then +# echo "------------> ERROR : Output FILE Directory not provided !" +# exit 1 +#fi # is the path absolute? -if [[ "$clean_files" = /* ]]; then - rm -f $clean_files/step2* $clean_files/step3* -else - rm -f $CONFIGFILEPATH/$clean_files/step2* $CONFIGFILEPATH/$clean_files/step3* -fi - -clean_results=`grep "results\s*:" $CONFIGFILE | cut -d ':' -f2 | sed -e 's/\s*$//' | sed -e 's/^\s*//'` -if [ -z "$clean_results" ]; then - clean_results=`grep "results\s*=" $CONFIGFILE | cut -d '=' -f2 | sed -e 's/\s*$//' | sed -e 's/^\s*//'` -fi - -if [ -z "$clean_results" ]; then - echo "------------> ERROR : Output RESULTS Directory not provided !" - exit 1 -fi +#if [[ "$clean_files" = /* ]]; then +# rm -f $clean_files/step2* $clean_files/step3* +#else +# rm -f $CONFIGFILEPATH/$clean_files/step2* $CONFIGFILEPATH/$clean_files/step3* +#fi + +#clean_results=`grep "results\s*:" $CONFIGFILE | cut -d ':' -f2 | sed -e 's/\s*$//' | sed -e 's/^\s*//'` +#if [ -z "$clean_results" ]; then +# clean_results=`grep "results\s*=" $CONFIGFILE | cut -d '=' -f2 | sed -e 's/\s*$//' | sed -e 's/^\s*//'` +#fi + +#if [ -z "$clean_results" ]; then +# echo "------------> ERROR : Output RESULTS Directory not provided !" +# exit 1 +#fi # is the path absolute? -if [[ "$clean_results" = /* ]]; then - rm -rf $clean_results - mkdir $clean_results -else - rm -rf $CONFIGFILEPATH/$clean_results - mkdir $CONFIGFILEPATH/$clean_results -fi +#if [[ "$clean_results" = /* ]]; then +# rm -rf $clean_results +# mkdir $clean_results +#else +# rm -rf $CONFIGFILEPATH/$clean_results +# mkdir $CONFIGFILEPATH/$clean_results +#fi # test if mkdir works -if [ $? -ne 0 ] ; then - echo "------------> ERROR : Impossible to create Output directory !" - exit 1 -fi +#if [ $? -ne 0 ] ; then +# echo "------------> ERROR : Impossible to create Output directory !" +# exit 1 +#fi # ----------------------------------- # exec second step of HRU-DELIN batch # ----------------------------------- -python3 $MYDIR/../modules/hrudelin_2_basins.py $CONFIGFILE $NBPROCESS +#python3 $MYDIR/../modules/hrudelin_2_basins.py $CONFIGFILE $NBPROCESS + +python3 $MYDIR/../modules/hrudelin_2_1_env_relocate.py $CONFIGFILE $NBPROCESS +python3 $MYDIR/../modules/hrudelin_2_2_derivation_watershed.py $CONFIGFILE $NBPROCESS +python3 $MYDIR/../modules/hrudelin_2_3_watershed_reclassification.py $CONFIGFILE $NBPROCESS +python3 $MYDIR/../modules/hrudelin_2_4_mask.py $CONFIGFILE $NBPROCESS +python3 $MYDIR/../modules/hrudelin_2_5_cut_stream_subbassins.py $CONFIGFILE $NBPROCESS +python3 $MYDIR/../modules/hrudelin_2_6_parrallele.py $CONFIGFILE $NBPROCESS +python3 $MYDIR/../modules/hrudelin_2_7_isolate_pixel.py $CONFIGFILE $NBPROCESS + diff --git a/modules/hrudelin_1_init.py b/modules/hrudelin_1_init.py index 27b52d997e45281243e49b5d9b8b0ca53e29d95b..538da85b02242a18dd401abbdb2ed754dc80f0bc 100755 --- a/modules/hrudelin_1_init.py +++ b/modules/hrudelin_1_init.py @@ -8,6 +8,8 @@ # AUTHOR(S): adapted from GRASS-HRU (ILMS) - JENA University # by IRSTEA - Christine Barachet, # Julien Veyssier +# Michael Rabotin +# Florent Veillon # PURPOSE: Prepare files for the next steps of HRU Delineation # shapefile : selected gauges # rasters : bounded DEM and others (geology, soils, landuse) @@ -631,6 +633,35 @@ if __name__ == '__main__': print('or if you don\'t have superuser access:\n') print('%spip3 install tqdm%s\n' % (COLOR_GREEN, COLOR_RESET)) sys.exit(1) + try: + from rastertodataframe import raster_to_dataframe + except Exception as e: + print('!!! %s rastertodataframe python3 module not found !!%s\n' % (COLOR_RED, COLOR_RESET)) + sys.exit(1) + try: + import geopandas as gpd + except Exception as e: + print('!!! %s geopandas python3 module not found !!%s\n' % (COLOR_RED, COLOR_RESET)) + print('On Debian/Ubuntu/Linux Mint you can install it with:\n') + print('%ssudo apt install python3-geopandas%s\n' % (COLOR_GREEN, COLOR_RESET)) + print('or if you don\'t have superuser access:\n') + print('%spip3 install geopandas%s\n' % (COLOR_GREEN, COLOR_RESET)) + sys.exit(1) + try: + import rtree + except Exception as e: + print('!!! %s rtree python3 module not found !!%s\n' % (COLOR_RED, COLOR_RESET)) + print('On Debian/Ubuntu/Linux Mint you can install it with:\n') + print('%ssudo apt install python3-rtree%s\n' % (COLOR_GREEN, COLOR_RESET)) + print('or if you don\'t have superuser access:\n') + print('%spip3 install rtree%s\n' % (COLOR_GREEN, COLOR_RESET)) + sys.exit(1) + try: + import pygeos + except Exception as e: + print('!!! %s pygeos python3 module not found !!%s\n' % (COLOR_RED, COLOR_RESET)) + print('%spip install pygeos%s\n' % (COLOR_GREEN, COLOR_RESET)) + sys.exit(1) parms_file = 'hrudelin_config.cfg' if len(sys.argv) > 1: @@ -639,7 +670,7 @@ if __name__ == '__main__': main(parms_file) try: - os.system('notify-send "hru-delin step 1 complete"') + os.system('notify-send "hru-delin-6-2 step 1 complete"') except Exception as e: pass else: diff --git a/modules/hrudelin_2_1_env_relocate.py b/modules/hrudelin_2_1_env_relocate.py new file mode 100755 index 0000000000000000000000000000000000000000..c5a9bd4710b19ff4565b1feb03fa125c59e27cfd --- /dev/null +++ b/modules/hrudelin_2_1_env_relocate.py @@ -0,0 +1,356 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + + +############################################################################ +# +# MODULE: hru-delin_basins.py +# AUTHOR(S): adapted from GRASS-HRU (ILMS) - JENA University +# by IRSTEA - Christine Barachet, +# Julien Veyssier +# Michael Rabotin +# Florent Veillon +# PURPOSE: 1. Relocates the gauges on the reaches +# 2. Calculates watersheds at the gauges +# +# +# COPYRIGHT: (C) 2020 UR RIVERLY - INRAE +# +# This program is free software under the GNU General Public +# License (>=v2). Read the file LICENSE that comes with +# HRU-DELIN for details. +# +############################################################################# + + + + +# to keep python2 compatibility +from __future__ import print_function +import string, os, sys, glob, types, time, platform +import numpy as np +try: + import ConfigParser +except Exception as e: + import configparser as ConfigParser +#import grass.script as grass +from grass.script.utils import decode, encode +import struct, math, csv, shutil + +from osgeo import gdal +from osgeo.gdalnumeric import * +from osgeo.gdalconst import * +from osgeo import ogr + +import multiprocessing +from multiprocessing import Pool, cpu_count + +from utils import isint, write_log +from reach import snapping_points_to_reaches, cut_streams_at_points +from reach import updateAttributeTable, processReachStats + +MY_ABS_PATH=os.path.abspath(__file__) +MY_DIR=os.path.dirname(MY_ABS_PATH) + +try: + # Python 3 + from subprocess import DEVNULL +except ImportError: + DEVNULL = open(os.devnull, 'wb') + +import pandas as pd + +''' + + MAIN + +''' +def main(parms_file, nbProc, generator=False): + print("-----------------------------------------------------------------------------------") + print('---------- HRU-delin Step 2-1 started ---------------------------------------------') + print("-----------------------------------------------------------------------------------") + + """OUTPUT files + - gauges_reloc.csv + - dams_reloc.csv + - step2_dams_for_watersheds.shp + - step2_gauges_for_watersheds.shp + """ + + #################### + #GRASS ENVIRONNEMENT + #################### + + configFileDir = os.path.dirname(parms_file) + # create main env + buildGrassEnv(os.path.join(configFileDir, 'grass_db'), 'hru-delin') + os.environ['GISRC'] = os.path.join(configFileDir, 'grass_db', 'grassdata', 'hru-delin', '.grassrc') + # Get parameters from configuration file + parms = ConfigParser.ConfigParser(allow_no_value=True) + tmpPath = os.path.join(configFileDir, 'tmp') + if not os.path.isdir(tmpPath): + os.mkdir(tmpPath) + parms.read(parms_file) + directory_out = parms.get('dir_out', 'files') + # manage absolute and relative paths + if not os.path.isabs(directory_out): + directory_out = os.path.join(configFileDir, directory_out) + + + + # test parameters from configuration file + # if auto_relocation == yes, test int value for surface_tolerance_1 and distance_tolerance_1 + if (parms.get('auto_relocation', 'to_do')) == 'yes': + + if not isint(parms.get('auto_relocation', 'surface_tolerance_1')): + sys.exit('------------> ERROR : Surface_tolerance_1 value not provided or is not integer' ) + if not isint(parms.get('auto_relocation', 'distance_tolerance_1')): + sys.exit('------------> ERROR : Distance_tolerance_1 value not provided or is not integer' ) + + + ## test if basin min size is valid + if not isint(parms.get('basin_min_size', 'size')): + sys.exit('------------> ERROR : Basin min size value not provided or is not integer' ) + + + ####### + #GAUGES + ####### + + # Get the shape of gauges + gauges_file = parms.get('gauges', 'relocated_gauges') + if gauges_file == '': + gauges_file = os.path.join(directory_out, 'gauges_selected.shp') + else: + if ogr.Open(gauges_file) is None: + sys.exit('------------> ERROR : Relocated Gauges file not found') + + + gauges_in = ogr.Open(gauges_file) + gauges_lyr = gauges_in.GetLayer() + + # Set the new shape + gauges_reloc_name = 'step2_gauges_for_watersheds' + gauges_reloc_file = os.path.join(directory_out, gauges_reloc_name + '.shp') + driver = ogr.GetDriverByName('ESRI Shapefile') + if os.path.exists(gauges_reloc_file): + driver.DeleteDataSource(gauges_reloc_file) + gauges_reloc_shp = driver.CreateDataSource(gauges_reloc_file) + gauges_reloc_lyr = gauges_reloc_shp.CopyLayer(gauges_lyr, gauges_reloc_name) + + # Relocation of the gauges + if (parms.get('auto_relocation', 'to_do')) == 'yes': + print('---------- HRU-delin Step 2-1 : Relocation of the gauges') + + gauges_area_col_name = parms.get('gauges', 'gauges_area_col_name') + gauges_col_name = parms.get('gauges', 'gauges_col_name') + snapping_points_to_reaches(parms, directory_out,gauges_col_name,gauges_area_col_name,gauges_reloc_lyr,'gauges_reloc.csv','gauges') + + gauges_reloc_shp.ExecuteSQL('REPACK ' + gauges_reloc_lyr.GetName()) + gauges_reloc_shp.Destroy() + + ##### + #DAMS + ##### + + # relocation of dams if provided + dams_reloc_file=0 + if str(parms.get('dams', 'to_do')) == 'yes': + # Get the shape of dams + dams_file = parms.get('dams', 'relocated_dams') + if dams_file == '': + dams_file = os.path.join(directory_out, 'dams_selected.shp') + else: + if ogr.Open(dams_file) is None: + sys.exit('------------> ERROR : Relocated dams file not found') + + + dams_in = ogr.Open(dams_file) + dams_lyr = dams_in.GetLayer() + + # Set the new shape + dams_reloc_name = 'step2_dams_for_watersheds' + dams_reloc_file = os.path.join(directory_out, dams_reloc_name + '.shp') + driver = ogr.GetDriverByName('ESRI Shapefile') + if os.path.exists(dams_reloc_file): + driver.DeleteDataSource(dams_reloc_file) + dams_reloc_shp = driver.CreateDataSource(dams_reloc_file) + dams_reloc_lyr = dams_reloc_shp.CopyLayer(dams_lyr, dams_reloc_name) + + # Relocation of the dams + if (parms.get('auto_relocation', 'to_do')) == 'yes': + print('---------- HRU-delin Step 2-1 : Relocation of the dams') + + dams_area_col_name = parms.get('dams', 'dams_area_col_name') + dams_col_name = parms.get('dams', 'dams_col_name') + snapping_points_to_reaches(parms, directory_out,dams_col_name,dams_area_col_name,dams_reloc_lyr,'dams_reloc.csv','dams') + + dams_reloc_shp.ExecuteSQL('REPACK ' + dams_reloc_lyr.GetName()) + dams_reloc_shp.Destroy() + + + + + # Import drain raster + drain_layer = os.path.join(directory_out, 'step1_drain.tif') + drain_wk = 'drain_wk' + grass_run_command('r.in.gdal', flags='o', input=drain_layer, output=drain_wk, overwrite='True', stdout=DEVNULL, stderr=DEVNULL) + grass_run_command('g.proj', flags='p', georef=drain_layer, stdout=DEVNULL, stderr=DEVNULL) + grass_run_command('g.region', flags='sp', raster=drain_wk, stdout=DEVNULL, stderr=DEVNULL) + + # Watersheds derivation + basins = 'basins' + rasters_list = ['basins'] + + # r.cross doesn't accept more than 30 layers + max_rasters = 29 + i = 0 + grass_run_command('r.mapcalc', expression='basins=null()', overwrite=True, stdout=DEVNULL, stderr=DEVNULL) + #gauges + gaugesDs = ogr.Open(gauges_reloc_file) + gauges_reloc_lyr = gaugesDs.GetLayer() + nb_gauges = gauges_reloc_lyr.GetFeatureCount() + #dams + nb_dams=0 + + + list_point=[] + if str(parms.get('dams', 'to_do')) == 'yes': + dams_col_name = parms.get('dams', 'dams_col_name') + damsDs = ogr.Open(dams_reloc_file) + dams_reloc_lyr = damsDs.GetLayer() + nb_dams = dams_reloc_lyr.GetFeatureCount() + + for dam in dams_reloc_lyr: + geom = dam.GetGeometryRef() + dam_x, dam_y = geom.GetX(), geom.GetY() + dam_ID=int(dam.GetField(dams_col_name)) + tuple_dam=(dam_x,dam_y,dam_ID) + list_point.append(tuple_dam) + + nb_points=nb_gauges+nb_dams + + #test if nb_points still has feature + if nb_points == 0: + print("Error on gauges (and dams if provided) relocated layer : 0 features found ") + print("Maybe check basin min size parameter") + sys.exit() + + gauges_col_name = parms.get('gauges', 'gauges_col_name') + for gauge in gauges_reloc_lyr: + geom = gauge.GetGeometryRef() + gauge_x, gauge_y = geom.GetX(), geom.GetY() + gauge_ID=int(gauge.GetField(gauges_col_name)) + tuple_gauge=(gauge_x,gauge_y,gauge_ID) + list_point.append(tuple_gauge) + + #export of list_point + df_list_point = pd.DataFrame(list_point) + df_list_point.to_csv(os.path.join(directory_out,"list_point.csv"),index=False, header=False) + + #TEST EXIST FILES AND FILL FILES + print('---------- HRU-delin Step 2-1 : Test of existing and completed files') + #GAUGES + + #gauges_reloc + gauges_reloc_path = os.path.join(directory_out, 'gauges_reloc.csv') + if (parms.get('auto_relocation', 'to_do')) == 'yes': + if os.stat(gauges_reloc_path).st_size == 0: + print('--------------- gauges_reloc.csv is empty or nonexistent') + else: + pd_gauges_reloc = pd.read_csv(gauges_reloc_path) + len_gauges_reloc = len(pd_gauges_reloc) + if len_gauges_reloc > 0 : + print('--------------- gauges_reloc.csv is created and it has ', len_gauges_reloc, " lines") + else : + print("--------------- gauges_reloc.csv is created but it empty") + + + #gauges_selected + gauges_selected_path = os.path.join(directory_out, 'gauges_selected.shp') + + if os.stat(gauges_selected_path).st_size == 0: + print('--------------- gauges_selected.shp is empty or nonexistent') + else: + datasource_gauges = ogr.Open(gauges_selected_path) + layer_gauges = datasource_gauges.GetLayer() + featureCount_gauges = layer_gauges.GetFeatureCount() + if featureCount_gauges > 0 : + print('--------------- gauges_selected.shp is created and it has ', featureCount_gauges, " features") + else : + print("--------------- gauges_selected.shp is created but it empty") + + + #DAMS + if str(parms.get('dams', 'to_do')) == 'yes': + #dams_reloc + dams_reloc_path = os.path.join(directory_out, 'dams_reloc.csv') + + if os.stat(dams_reloc_path).st_size == 0: + print('--------------- dams_reloc.csv is empty or nonexistent') + else: + pd_dams_reloc = pd.read_csv(dams_reloc_path) + len_dams_reloc = len(pd_dams_reloc) + if len_dams_reloc > 0 : + print('--------------- dams_reloc.csv is created and it has ', len_dams_reloc, " lines") + else : + print("--------------- dams_reloc.csv is created but it empty") + + #dams_selected + dams_selected_path = os.path.join(directory_out, 'dams_selected.shp') + + if os.stat(dams_selected_path).st_size == 0: + print('--------------- dams_selected.shp is empty or nonexistent') + else: + datasource_dams = ogr.Open(dams_selected_path) + layer_dams = datasource_dams.GetLayer() + featureCount_dams = layer_dams.GetFeatureCount() + if featureCount_dams > 0 : + print('--------------- dams_selected.shp is created and it has ', featureCount_dams, " features") + else : + print("--------------- dams_selected.shp is created but it empty") + + + print('---------- HRU-delin Step 2-1 ended ---------------------------------------------') + + if generator: + yield 10 + + +#MAIN _PARAM +if __name__ == '__main__': + from grassUtils import buildGrassEnv, buildGrassLocation, exportRasters, importRastersInEnv,\ + grass_run_command, grass_parse_command, grass_feed_command, grass_read_command, grass_pipe_command + from progressColors import * + # check TQDM presence only if we are executed + try: + from tqdm import tqdm + except Exception as e: + print('!! %stqdm module not found%s\n' % (COLOR_RED, COLOR_RESET)) + sys.exit(1) + + parms_file = 'hrudelin_config.cfg' + nbProcArg = '' + if len(sys.argv) > 1: + parms_file = sys.argv[1] + if len(sys.argv) > 2: + nbProcArg = sys.argv[2] + + # determine how many processes we can launch + if str(nbProcArg).isnumeric() and int(nbProcArg) > 0: + nbProc = int(nbProcArg) + else: + nbProc = cpu_count() + + # main is a generator but we don't use it here + for pc in main(parms_file, nbProc, False): + pass + + try: + os.system('notify-send "hru-delin-6-2 step 2-1 complete"') + except Exception as e: + pass +else: + from .grassUtils import buildGrassEnv, buildGrassLocation, exportRasters, importRastersInEnv,\ + grass_run_command, grass_parse_command, grass_feed_command, grass_read_command, grass_pipe_command + from .progressColors import * diff --git a/modules/hrudelin_2_2_derivation_watershed.py b/modules/hrudelin_2_2_derivation_watershed.py new file mode 100755 index 0000000000000000000000000000000000000000..a99ac2bcd3a075b1103f65e762cee3121ce208b2 --- /dev/null +++ b/modules/hrudelin_2_2_derivation_watershed.py @@ -0,0 +1,203 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + + +############################################################################ +# +# MODULE: hru-delin_basins.py +# AUTHOR(S): adapted from GRASS-HRU (ILMS) - JENA University +# by IRSTEA - Christine Barachet, +# Julien Veyssier +# Michael Rabotin +# Florent Veillon +# PURPOSE: 1. Relocates the gauges on the reaches +# 2. Calculates watersheds at the gauges +# +# +# COPYRIGHT: (C) 2020 UR RIVERLY - INRAE +# +# This program is free software under the GNU General Public +# License (>=v2). Read the file LICENSE that comes with +# HRU-DELIN for details. +# +############################################################################# + + + + +# to keep python2 compatibility +from __future__ import print_function +import string, os, sys, glob, types, time, platform +import numpy as np +try: + import ConfigParser +except Exception as e: + import configparser as ConfigParser +#import grass.script as grass +from grass.script.utils import decode, encode +import struct, math, csv, shutil + +from osgeo import gdal +from osgeo.gdalnumeric import * +from osgeo.gdalconst import * +from osgeo import ogr + +import multiprocessing +from multiprocessing import Pool, cpu_count + +from utils import isint, write_log +from reach import snapping_points_to_reaches, cut_streams_at_points +from reach import updateAttributeTable, processReachStats + +MY_ABS_PATH=os.path.abspath(__file__) +MY_DIR=os.path.dirname(MY_ABS_PATH) + +try: + # Python 3 + from subprocess import DEVNULL +except ImportError: + DEVNULL = open(os.devnull, 'wb') + +import pandas as pd +from rastertodataframe import raster_to_dataframe + +''' + + MAIN + +''' +def main(parms_file, nbProc, generator=False): + + """OUTPUT files + - basins.tif + """ + print(" ") + print('---------- HRU-delin Step 2-2 started ---------------------------------------------') + print("-----------------------------------------------------------------------------------") + + configFileDir = os.path.dirname(parms_file) + parms = ConfigParser.ConfigParser(allow_no_value=True) + tmpPath = os.path.join(configFileDir, 'tmp') + if not os.path.isdir(tmpPath): + os.mkdir(tmpPath) + parms.read(parms_file) + directory_out = parms.get('dir_out', 'files') + # manage absolute and relative paths + if not os.path.isabs(directory_out): + directory_out = os.path.join(configFileDir, directory_out) + + #open list_point.csv + file = pd.read_csv(os.path.join(directory_out,"list_point.csv"),header=None) + #convert to tuples + list_point = list(file.itertuples(index=False, name=None)) + nb_points = len(list_point) + #Set Grass environnement + os.environ['GISRC'] = os.path.join(configFileDir, 'grass_db', 'grassdata', 'hru-delin', '.grassrc') + # Import drain raster + drain_layer = os.path.join(directory_out, 'step1_drain.tif') + drain_wk = 'drain_wk' + grass_run_command('r.in.gdal', flags='o', input=drain_layer, output=drain_wk, overwrite='True', stdout=DEVNULL, stderr=DEVNULL) + grass_run_command('g.proj', flags='p', georef=drain_layer, stdout=DEVNULL, stderr=DEVNULL) + grass_run_command('g.region', flags='sp', raster=drain_wk, stdout=DEVNULL, stderr=DEVNULL) + + i = 0 + basins = 'basins' + rasters_list = ['basins'] + max_rasters = 29 + + + + print('---------- HRU-delin Step 2-2 : r.water.outlet') + + for points in list_point: + + grass_run_command('r.water.outlet', + input=drain_wk, output='basin_tmp%d' % i, + coordinates='%s,%s' % (points[0], points[1]), + overwrite='True', + stdout=DEVNULL, stderr=DEVNULL + ) + + + rasters_list.append('basin_tmp%d' % i) + i += 1 + + if i % max_rasters == 0 or i == nb_points: + grass_run_command('r.cross', input=','.join(rasters_list), output=basins, overwrite='True', stdout=DEVNULL, stderr=DEVNULL) + grass_run_command('g.remove', flags='f', type='raster', pattern='basin_tmp*', stdout=DEVNULL, stderr=DEVNULL) + rasters_list = ['basins'] + + if generator: + yield 5 + + # with grass7 watershed ids begin at 0 + # empty areas value is already NULL, no need to set it + grass_run_command('r.mapcalc', expression='basins=basins+1', overwrite='True', stdout=DEVNULL, stderr=DEVNULL) + #TODO remove this line + grass_run_command('r.out.gdal', input='basins', output=os.path.join(directory_out, 'basins.tif'), + overwrite='True', stdout=DEVNULL, stderr=DEVNULL) + + + + #TEST EXIST FILES AND FILL FILES + print('---------- HRU-delin Step 2-2 : Test of existing and completed files') + + #basins.tif + tif_path = os.path.join(directory_out, 'basins.tif') + + if os.stat(tif_path).st_size == 0: + print('--------------- basins.tif is empty or nonexistent') + else: + tif_gdal = gdal.Open(tif_path) + tif_band = tif_gdal.GetRasterBand(1) + (min_tif,max_tif) = tif_band.ComputeRasterMinMax(True) + + if min_tif > 0 and max_tif <= nb_points: + print('--------------- basins.tif is created and it has',nb_points, "watersheds") + else : + print("--------------- basins.tif is created but it empty") + + print('---------- HRU-delin Step 2-2 ended ---------------------------------------------') + + + + + + + +if __name__ == '__main__': + from grassUtils import buildGrassEnv, buildGrassLocation, exportRasters, importRastersInEnv,\ + grass_run_command, grass_parse_command, grass_feed_command, grass_read_command, grass_pipe_command + from progressColors import * + # check TQDM presence only if we are executed + try: + from tqdm import tqdm + except Exception as e: + print('!! %stqdm module not found%s\n' % (COLOR_RED, COLOR_RESET)) + sys.exit(1) + + parms_file = 'hrudelin_config.cfg' + nbProcArg = '' + if len(sys.argv) > 1: + parms_file = sys.argv[1] + if len(sys.argv) > 2: + nbProcArg = sys.argv[2] + + # determine how many processes we can launch + if str(nbProcArg).isnumeric() and int(nbProcArg) > 0: + nbProc = int(nbProcArg) + else: + nbProc = cpu_count() + + # main is a generator but we don't use it here + for pc in main(parms_file, nbProc, False): + pass + + try: + os.system('notify-send "hru-delin-6-2 step 2-2 complete"') + except Exception as e: + pass +else: + from .grassUtils import buildGrassEnv, buildGrassLocation, exportRasters, importRastersInEnv,\ + grass_run_command, grass_parse_command, grass_feed_command, grass_read_command, grass_pipe_command + from .progressColors import * diff --git a/modules/hrudelin_2_3_watershed_reclassification.py b/modules/hrudelin_2_3_watershed_reclassification.py new file mode 100755 index 0000000000000000000000000000000000000000..3b005f1c1c48b17778067419e2bda6bd31d1c782 --- /dev/null +++ b/modules/hrudelin_2_3_watershed_reclassification.py @@ -0,0 +1,190 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + + +############################################################################ +# +# MODULE: hru-delin_basins.py +# AUTHOR(S): adapted from GRASS-HRU (ILMS) - JENA University +# by IRSTEA - Christine Barachet, +# Julien Veyssier +# Michael Rabotin +# Florent Veillon +# PURPOSE: 1. Relocates the gauges on the reaches +# 2. Calculates watersheds at the gauges +# +# +# COPYRIGHT: (C) 2020 UR RIVERLY - INRAE +# +# This program is free software under the GNU General Public +# License (>=v2). Read the file LICENSE that comes with +# HRU-DELIN for details. +# +############################################################################# + + + + +# to keep python2 compatibility +from __future__ import print_function +import string, os, sys, glob, types, time, platform +import numpy as np +try: + import ConfigParser +except Exception as e: + import configparser as ConfigParser +#import grass.script as grass +from grass.script.utils import decode, encode +import struct, math, csv, shutil + +from osgeo import gdal +from osgeo.gdalnumeric import * +from osgeo.gdalconst import * +from osgeo import ogr + +import multiprocessing +from multiprocessing import Pool, cpu_count + +from utils import isint, write_log +from reach import snapping_points_to_reaches, cut_streams_at_points +from reach import updateAttributeTable, processReachStats + +MY_ABS_PATH=os.path.abspath(__file__) +MY_DIR=os.path.dirname(MY_ABS_PATH) + +try: + # Python 3 + from subprocess import DEVNULL +except ImportError: + DEVNULL = open(os.devnull, 'wb') + +import pandas as pd +from rastertodataframe import raster_to_dataframe + +''' + + MAIN + +''' +def main(parms_file, nbProc, generator=False): + + """OUTPUT files + - step2_watersheds.tif + """ + + print(" ") + print('---------- HRU-delin Step 2-3 started ---------------------------------------------') + print("-----------------------------------------------------------------------------------") + + configFileDir = os.path.dirname(parms_file) + parms = ConfigParser.ConfigParser(allow_no_value=True) + tmpPath = os.path.join(configFileDir, 'tmp') + if not os.path.isdir(tmpPath): + os.mkdir(tmpPath) + parms.read(parms_file) + directory_out = parms.get('dir_out', 'files') + # manage absolute and relative paths + if not os.path.isabs(directory_out): + directory_out = os.path.join(configFileDir, directory_out) + + #open list_point.csv + file = pd.read_csv(os.path.join(directory_out,"list_point.csv"),header=None) + #convert to tuples + list_point = list(file.itertuples(index=False, name=None)) + #Set Grass environnement + os.environ['GISRC'] = os.path.join(configFileDir, 'grass_db', 'grassdata', 'hru-delin', '.grassrc') + # Import drain raster + drain_layer = os.path.join(directory_out, 'step1_drain.tif') + drain_wk = 'drain_wk' + grass_run_command('r.in.gdal', flags='o', input=drain_layer, output=drain_wk, overwrite='True', stdout=DEVNULL, stderr=DEVNULL) + grass_run_command('g.proj', flags='p', georef=drain_layer, stdout=DEVNULL, stderr=DEVNULL) + grass_run_command('g.region', flags='sp', raster=drain_wk, stdout=DEVNULL, stderr=DEVNULL) + # Import basins raster + basins_tif = os.path.join(directory_out, 'basins.tif') + basins = "basins" + + + basins_rcl = 'basins_rcl' + + #create a GRASS basins raster + grass_run_command('r.in.gdal', flags='o', input=basins_tif, output=basins, overwrite='True', stdout=DEVNULL, stderr=DEVNULL) + + pRecode = grass_feed_command('r.recode', input=basins, output=basins_rcl, overwrite='True', rules='-', quiet=True) + + + for points in list_point: + out = decode(grass_read_command('r.what', map=basins, coordinates='%s,%s' % (points[0], points[1]))) + cat_grass = out.rstrip(os.linesep).split('|')[3].strip() + pRecode.stdin.write(encode('%s:%s:%s\n' % (cat_grass, cat_grass, points[2]))) + pRecode.stdin.close() + pRecode.wait() + + grass_run_command('r.out.gdal', + input=basins_rcl, + type='UInt16', + output=os.path.join(directory_out, 'step2_watersheds.tif'), + overwrite='True', stdout=DEVNULL, stderr=DEVNULL) + + + #TEST EXIST FILES AND FILL FILES + print('---------- HRU-delin Step 2-3 : Test of existing and completed files') + + #basins.tif + tif_path = os.path.join(directory_out, 'step2_watersheds.tif') + + if os.stat(tif_path).st_size == 0: + print('--------------- step2_watersheds.tif is empty or nonexistent') + else: + tif_gdal = gdal.Open(tif_path) + tif_band = tif_gdal.GetRasterBand(1) + (min_tif,max_tif) = tif_band.ComputeRasterMinMax(True) + + if min_tif > 0 : + print('--------------- step2_watersheds.tif is created') + else : + print("--------------- step2_watersheds.tif is created but it empty") + + + if generator: + yield 10 + + + print('---------- HRU-delin Step 2-3 ended ---------------------------------------------') + + +if __name__ == '__main__': + from grassUtils import buildGrassEnv, buildGrassLocation, exportRasters, importRastersInEnv,\ + grass_run_command, grass_parse_command, grass_feed_command, grass_read_command, grass_pipe_command + from progressColors import * + # check TQDM presence only if we are executed + try: + from tqdm import tqdm + except Exception as e: + print('!! %stqdm module not found%s\n' % (COLOR_RED, COLOR_RESET)) + sys.exit(1) + + parms_file = 'hrudelin_config.cfg' + nbProcArg = '' + if len(sys.argv) > 1: + parms_file = sys.argv[1] + if len(sys.argv) > 2: + nbProcArg = sys.argv[2] + + # determine how many processes we can launch + if str(nbProcArg).isnumeric() and int(nbProcArg) > 0: + nbProc = int(nbProcArg) + else: + nbProc = cpu_count() + + # main is a generator but we don't use it here + for pc in main(parms_file, nbProc, False): + pass + + try: + os.system('notify-send "hru-delin-6-2 step 2-3 complete"') + except Exception as e: + pass +else: + from .grassUtils import buildGrassEnv, buildGrassLocation, exportRasters, importRastersInEnv,\ + grass_run_command, grass_parse_command, grass_feed_command, grass_read_command, grass_pipe_command + from .progressColors import * diff --git a/modules/hrudelin_2_4_mask.py b/modules/hrudelin_2_4_mask.py new file mode 100755 index 0000000000000000000000000000000000000000..889b0f9d71b1820a9d3adbbe69ba9a3d62b98c58 --- /dev/null +++ b/modules/hrudelin_2_4_mask.py @@ -0,0 +1,201 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + + +############################################################################ +# +# MODULE: hru-delin_basins.py +# AUTHOR(S): adapted from GRASS-HRU (ILMS) - JENA University +# by IRSTEA - Christine Barachet, +# Julien Veyssier +# Michael Rabotin +# Florent Veillon +# PURPOSE: 1. Relocates the gauges on the reaches +# 2. Calculates watersheds at the gauges +# +# +# COPYRIGHT: (C) 2020 UR RIVERLY - INRAE +# +# This program is free software under the GNU General Public +# License (>=v2). Read the file LICENSE that comes with +# HRU-DELIN for details. +# +############################################################################# + + + + +# to keep python2 compatibility +from __future__ import print_function +import string, os, sys, glob, types, time, platform +import numpy as np +try: + import ConfigParser +except Exception as e: + import configparser as ConfigParser +#import grass.script as grass +from grass.script.utils import decode, encode +import struct, math, csv, shutil + +from osgeo import gdal +from osgeo.gdalnumeric import * +from osgeo.gdalconst import * +from osgeo import ogr + +import multiprocessing +from multiprocessing import Pool, cpu_count + +from utils import isint, write_log +from reach import snapping_points_to_reaches, cut_streams_at_points +from reach import updateAttributeTable, processReachStats + +MY_ABS_PATH=os.path.abspath(__file__) +MY_DIR=os.path.dirname(MY_ABS_PATH) + +try: + # Python 3 + from subprocess import DEVNULL +except ImportError: + DEVNULL = open(os.devnull, 'wb') + +import pandas as pd +from rastertodataframe import raster_to_dataframe + +''' + + MAIN + +''' +def main(parms_file, nbProc, generator=False): + + """OUTPUT files + - step2_mask.tif + """ + print(" ") + print('---------- HRU-delin Step 2-4 started ---------------------------------------------') + print("-----------------------------------------------------------------------------------") + + configFileDir = os.path.dirname(parms_file) + parms = ConfigParser.ConfigParser(allow_no_value=True) + tmpPath = os.path.join(configFileDir, 'tmp') + if not os.path.isdir(tmpPath): + os.mkdir(tmpPath) + parms.read(parms_file) + directory_out = parms.get('dir_out', 'files') + # manage absolute and relative paths + if not os.path.isabs(directory_out): + directory_out = os.path.join(configFileDir, directory_out) + + #Set Grass environnement + os.environ['GISRC'] = os.path.join(configFileDir, 'grass_db', 'grassdata', 'hru-delin', '.grassrc') + # Import drain raster + drain_layer = os.path.join(directory_out, 'step1_drain.tif') + drain_wk = 'drain_wk' + grass_run_command('r.in.gdal', flags='o', input=drain_layer, output=drain_wk, overwrite='True', stdout=DEVNULL, stderr=DEVNULL) + grass_run_command('g.proj', flags='p', georef=drain_layer, stdout=DEVNULL, stderr=DEVNULL) + grass_run_command('g.region', flags='sp', raster=drain_wk, stdout=DEVNULL, stderr=DEVNULL) + + #open list_point.csv + file = pd.read_csv(os.path.join(directory_out,"list_point.csv"),header=None) + #convert to tuples + list_point = list(file.itertuples(index=False, name=None)) + # Get parameters from configuration file + parms = ConfigParser.ConfigParser(allow_no_value=True) + tmpPath = os.path.join(configFileDir, 'tmp') + if not os.path.isdir(tmpPath): + os.mkdir(tmpPath) + parms.read(parms_file) + + basins_rcl_tif = os.path.join(directory_out, 'step2_watersheds.tif') + basins_rcl = "basins_rcl" + + grass_run_command('r.in.gdal', flags='o', input=basins_rcl_tif, output=basins_rcl, overwrite='True', stdout=DEVNULL, stderr=DEVNULL) + + + grass_run_command('r.null', map=basins_rcl, setnull=0, stdout=DEVNULL, stderr=DEVNULL) + mask_tmp = 'mask_tmp' + pReclass = grass_feed_command('r.reclass', input=basins_rcl, output=mask_tmp, overwrite='True', rules='-') + pReclass.stdin.write(encode('0 thru 10000000 = 1\n')) + pReclass.stdin.close() + pReclass.wait() + + # TODO check that (was removed before) + #grass_run_command('r.null', map=mask_tmp, null=0, setnull=1, stdout=DEVNULL, stderr=DEVNULL) + grass_run_command('r.out.gdal', input=mask_tmp, type='UInt16', output=os.path.join(directory_out, 'step2_mask.tif'), + overwrite='True', stdout=DEVNULL, stderr=DEVNULL) + + # exit if mask is too small and reloc is disabled + inGd = gdal.Open(os.path.join(directory_out, 'step2_mask.tif')) + inband1 = inGd.GetRasterBand(1) + indata1 = BandReadAsArray(inband1) + inNodata = inband1.GetNoDataValue() + + nbGoodValues = np.count_nonzero(indata1 != inNodata) + xl = inGd.RasterXSize + yl = inGd.RasterYSize + # mask represents less than 0.1 % of area + maskProportion = nbGoodValues / (xl * yl) * 100 + maskIsTooSmall = maskProportion < 0.1 + if parms.get('auto_relocation', 'to_do') != 'yes' and maskIsTooSmall: + sys.exit('!!! Relocation is disabled and you are getting a very small mask. Please consider enabling gauge relocation.') + + # Cutting the streams at gauges + listReachID=cut_streams_at_points(directory_out,list_point,'step1_streams.tif','step2_streams_new.tif') + + #update attribute table with ReachID + gauges_col_name = parms.get('gauges', 'gauges_col_name') + updateAttributeTable(directory_out,listReachID,'step2_gauges_for_watersheds.shp',gauges_col_name) + + if str(parms.get('dams', 'to_do')) == 'yes': + dams_col_name = parms.get('dams', 'dams_col_name') + updateAttributeTable(directory_out,listReachID,'step2_dams_for_watersheds.shp',dams_col_name) + + if generator: + yield 20 + + print('---------- HRU-delin Step 2-4 ended ---------------------------------------------') + + if generator: + yield 15 + + + + + + +if __name__ == '__main__': + from grassUtils import buildGrassEnv, buildGrassLocation, exportRasters, importRastersInEnv,\ + grass_run_command, grass_parse_command, grass_feed_command, grass_read_command, grass_pipe_command + from progressColors import * + # check TQDM presence only if we are executed + try: + from tqdm import tqdm + except Exception as e: + print('!! %stqdm module not found%s\n' % (COLOR_RED, COLOR_RESET)) + sys.exit(1) + + parms_file = 'hrudelin_config.cfg' + nbProcArg = '' + if len(sys.argv) > 1: + parms_file = sys.argv[1] + if len(sys.argv) > 2: + nbProcArg = sys.argv[2] + + # determine how many processes we can launch + if str(nbProcArg).isnumeric() and int(nbProcArg) > 0: + nbProc = int(nbProcArg) + else: + nbProc = cpu_count() + + # main is a generator but we don't use it here + for pc in main(parms_file, nbProc, False): + pass + + try: + os.system('notify-send "hru-delin-6-2 step 2-4 complete"') + except Exception as e: + pass +else: + from .grassUtils import buildGrassEnv, buildGrassLocation, exportRasters, importRastersInEnv,\ + grass_run_command, grass_parse_command, grass_feed_command, grass_read_command, grass_pipe_command + from .progressColors import * diff --git a/modules/hrudelin_2_5_cut_stream_subbassins.py b/modules/hrudelin_2_5_cut_stream_subbassins.py new file mode 100755 index 0000000000000000000000000000000000000000..c175b2644b0a1af12bfe2376c7abdf66770e2893 --- /dev/null +++ b/modules/hrudelin_2_5_cut_stream_subbassins.py @@ -0,0 +1,162 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + + +############################################################################ +# +# MODULE: hru-delin_basins.py +# AUTHOR(S): adapted from GRASS-HRU (ILMS) - JENA University +# by IRSTEA - Christine Barachet, +# Julien Veyssier +# Michael Rabotin +# Florent Veillon +# PURPOSE: 1. Relocates the gauges on the reaches +# 2. Calculates watersheds at the gauges +# +# +# COPYRIGHT: (C) 2020 UR RIVERLY - INRAE +# +# This program is free software under the GNU General Public +# License (>=v2). Read the file LICENSE that comes with +# HRU-DELIN for details. +# +############################################################################# + + + + +# to keep python2 compatibility +from __future__ import print_function +import string, os, sys, glob, types, time, platform +import numpy as np +try: + import ConfigParser +except Exception as e: + import configparser as ConfigParser +#import grass.script as grass +from grass.script.utils import decode, encode +import struct, math, csv, shutil + +from osgeo import gdal +from osgeo.gdalnumeric import * +from osgeo.gdalconst import * +from osgeo import ogr + +import multiprocessing +from multiprocessing import Pool, cpu_count + +from utils import isint, write_log +from reach import snapping_points_to_reaches, cut_streams_at_points +from reach import updateAttributeTable, processReachStats + +MY_ABS_PATH=os.path.abspath(__file__) +MY_DIR=os.path.dirname(MY_ABS_PATH) + +try: + # Python 3 + from subprocess import DEVNULL +except ImportError: + DEVNULL = open(os.devnull, 'wb') + +import pandas as pd +from rastertodataframe import raster_to_dataframe + +''' + + MAIN + +''' +def main(parms_file, nbProc, generator=False): + + """OUTPUT files + + """ + print(" ") + print('---------- HRU-delin Step 2-5 started ---------------------------------------------') + print("-----------------------------------------------------------------------------------") + + configFileDir = os.path.dirname(parms_file) + parms = ConfigParser.ConfigParser(allow_no_value=True) + tmpPath = os.path.join(configFileDir, 'tmp') + if not os.path.isdir(tmpPath): + os.mkdir(tmpPath) + parms.read(parms_file) + directory_out = parms.get('dir_out', 'files') + # manage absolute and relative paths + if not os.path.isabs(directory_out): + directory_out = os.path.join(configFileDir, directory_out) + #Set Grass environnement + os.environ['GISRC'] = os.path.join(configFileDir, 'grass_db', 'grassdata', 'hru-delin', '.grassrc') + + # Import drain raster + drain_layer = os.path.join(directory_out, 'step1_drain.tif') + drain_wk = 'drain_wk' + grass_run_command('r.in.gdal', flags='o', input=drain_layer, output=drain_wk, overwrite='True', stdout=DEVNULL, stderr=DEVNULL) + grass_run_command('g.proj', flags='p', georef=drain_layer, stdout=DEVNULL, stderr=DEVNULL) + grass_run_command('g.region', flags='sp', raster=drain_wk, stdout=DEVNULL, stderr=DEVNULL) + + #open list_point.csv + file = pd.read_csv(os.path.join(directory_out,"list_point.csv"),header=None) + #convert to tuples + list_point = list(file.itertuples(index=False, name=None)) + # Get parameters from configuration file + parms = ConfigParser.ConfigParser(allow_no_value=True) + tmpPath = os.path.join(configFileDir, 'tmp') + if not os.path.isdir(tmpPath): + os.mkdir(tmpPath) + parms.read(parms_file) + + # Cutting the streams at gauges + listReachID=cut_streams_at_points(directory_out,list_point,'step1_streams.tif','step2_streams_new.tif') + + #update attribute table with ReachID + gauges_col_name = parms.get('gauges', 'gauges_col_name') + updateAttributeTable(directory_out,listReachID,'step2_gauges_for_watersheds.shp',gauges_col_name) + + if str(parms.get('dams', 'to_do')) == 'yes': + dams_col_name = parms.get('dams', 'dams_col_name') + updateAttributeTable(directory_out,listReachID,'step2_dams_for_watersheds.shp',dams_col_name) + + if generator: + yield 20 + + print('---------- HRU-delin Step 2-5 ended ---------------------------------------------') + + + +if __name__ == '__main__': + from grassUtils import buildGrassEnv, buildGrassLocation, exportRasters, importRastersInEnv,\ + grass_run_command, grass_parse_command, grass_feed_command, grass_read_command, grass_pipe_command + from progressColors import * + # check TQDM presence only if we are executed + try: + from tqdm import tqdm + except Exception as e: + print('!! %stqdm module not found%s\n' % (COLOR_RED, COLOR_RESET)) + sys.exit(1) + + parms_file = 'hrudelin_config.cfg' + nbProcArg = '' + if len(sys.argv) > 1: + parms_file = sys.argv[1] + if len(sys.argv) > 2: + nbProcArg = sys.argv[2] + + # determine how many processes we can launch + if str(nbProcArg).isnumeric() and int(nbProcArg) > 0: + nbProc = int(nbProcArg) + else: + nbProc = cpu_count() + + # main is a generator but we don't use it here + for pc in main(parms_file, nbProc, False): + pass + + try: + os.system('notify-send "hru-delin-6-2 step 2-5 complete"') + except Exception as e: + pass +else: + from .grassUtils import buildGrassEnv, buildGrassLocation, exportRasters, importRastersInEnv,\ + grass_run_command, grass_parse_command, grass_feed_command, grass_read_command, grass_pipe_command + from .progressColors import * diff --git a/modules/hrudelin_2_6_parrallele.py b/modules/hrudelin_2_6_parrallele.py new file mode 100755 index 0000000000000000000000000000000000000000..79dc97e0866fdd1cd669094fa53738148b246bd1 --- /dev/null +++ b/modules/hrudelin_2_6_parrallele.py @@ -0,0 +1,374 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + + +############################################################################ +# +# MODULE: hru-delin_basins.py +# AUTHOR(S): adapted from GRASS-HRU (ILMS) - JENA University +# by IRSTEA - Christine Barachet, +# Julien Veyssier +# Michael Rabotin +# Florent Veillon +# PURPOSE: 1. Relocates the gauges on the reaches +# 2. Calculates watersheds at the gauges +# +# +# COPYRIGHT: (C) 2020 UR RIVERLY - INRAE +# +# This program is free software under the GNU General Public +# License (>=v2). Read the file LICENSE that comes with +# HRU-DELIN for details. +# +############################################################################# + + + + +# to keep python2 compatibility +from __future__ import print_function +import string, os, sys, glob, types, time, platform +import numpy as np +try: + import ConfigParser +except Exception as e: + import configparser as ConfigParser +#import grass.script as grass +from grass.script.utils import decode, encode +import struct, math, csv, shutil + +from osgeo import gdal +from osgeo.gdalnumeric import * +from osgeo.gdalconst import * +from osgeo import ogr + +import multiprocessing +from multiprocessing import Pool, cpu_count + +from utils import isint, write_log +from reach import snapping_points_to_reaches, cut_streams_at_points +from reach import updateAttributeTable, processReachStats + +MY_ABS_PATH=os.path.abspath(__file__) +MY_DIR=os.path.dirname(MY_ABS_PATH) + +try: + # Python 3 + from subprocess import DEVNULL +except ImportError: + DEVNULL = open(os.devnull, 'wb') + +import pandas as pd +from rastertodataframe import raster_to_dataframe + +''' + + MAIN + +''' +def main(parms_file, nbProc, generator=False): + + """OUTPUT files + - step2_mask.tif + """ + print(" ") + print('---------- HRU-delin Step 2-6 started ---------------------------------------------') + print("-----------------------------------------------------------------------------------") + + configFileDir = os.path.dirname(parms_file) + parms = ConfigParser.ConfigParser(allow_no_value=True) + tmpPath = os.path.join(configFileDir, 'tmp') + if not os.path.isdir(tmpPath): + os.mkdir(tmpPath) + parms.read(parms_file) + directory_out = parms.get('dir_out', 'files') + # manage absolute and relative paths + if not os.path.isabs(directory_out): + directory_out = os.path.join(configFileDir, directory_out) + + #Set Grass environnement + os.environ['GISRC'] = os.path.join(configFileDir, 'grass_db', 'grassdata', 'hru-delin', '.grassrc') + + #open list_point.csv + #file = pd.read_csv(os.path.join(directory_out,"list_point.csv")) + #convert to tuples + #list_point = list(file.itertuples(index=False, name=None)) + # Get parameters from configuration file + #parms = ConfigParser.ConfigParser(allow_no_value=True) + tmpPath = os.path.join(configFileDir, 'tmp') +# if not os.path.isdir(tmpPath): +# os.mkdir(tmpPath) +# parms.read(parms_file) + + # Import drain raster + print('---------- Importing raster \'step1_drain.tif\'') + drain_layer = os.path.join(directory_out, 'step1_drain.tif') + drain_wk = 'drain_wk' + grass_run_command('r.in.gdal', flags='o', input=drain_layer, output=drain_wk, overwrite='True', stdout=DEVNULL, stderr=DEVNULL) + + + print('---------- Importing raster \'step1_dem_reclass.tif\'') + dem_recl = os.path.join(directory_out, 'step1_dem_reclass.tif') + grass_run_command('g.proj', flags='c', georef=dem_recl, stdout=DEVNULL, stderr=DEVNULL) + + print('---------- Importing raster \'step1_subbasins.tif\'') + subbasins = os.path.join(directory_out, 'step1_subbasins.tif') + grass_run_command('r.in.gdal', flags='o', input=subbasins, output='subbasins', overwrite='True', stdout=DEVNULL, stderr=DEVNULL) + + print('---------- Importing raster \'step2_watersheds.tif\'') + basins = os.path.join(directory_out, 'step2_watersheds.tif') + grass_run_command('r.in.gdal', flags='o', input=basins, output='watersheds', overwrite='True', stdout=DEVNULL, stderr=DEVNULL) + + print('---------- Importing raster \'step2_streams_new.tif\'') + reachraster = os.path.join(directory_out, 'step2_streams_new.tif') + grass_run_command('r.in.gdal', flags='o', input=reachraster, output='reachraster', overwrite='True', stdout=DEVNULL, stderr=DEVNULL) + + print('---------- Compute cross product \'watersheds*100000+subbasins\'') + grass_run_command('r.mapcalc', expression='cross1=watersheds*100000+subbasins', overwrite='True', stdout=DEVNULL, stderr=DEVNULL) + + print('---------- Compute cross product cross1*reachraster') + grass_run_command('r.mapcalc', expression='cross2=if(reachraster!=0, cross1, 0)', overwrite='True', stdout=DEVNULL, stderr=DEVNULL) +# + if generator: + yield 20 + + print('---------- Setting nulls in \'cross1*reachraster\' and \'reachraster\'') + grass_run_command('r.null', map='cross2', setnull=0, stdout=DEVNULL, stderr=DEVNULL) + grass_run_command('r.null', map='reachraster', setnull=0, stdout=DEVNULL, stderr=DEVNULL) + + print('---------- Saving \'step2_subbasins.tif\' (watersheds*100000+subbasins)') + grass_run_command('r.out.gdal', input='cross1', output=os.path.join(directory_out, 'step2_subbasins.tif'), + overwrite='True', stdout=DEVNULL, stderr=DEVNULL) + + print('---------- Computing links between (watersheds*100000+subbasins) and reach ids') + reach_ids = decode(grass_read_command('r.stats', quiet=True, flags='nN', input='reachraster')).rstrip(os.linesep).split(os.linesep) + reach_ids_cleaned = [] + subbasins_cleaned = [] + n_reach = len(reach_ids) + print('') + # main loop + + # export rasters that are necessary for parallel environments + rastersForWorkers = { + 'reachraster': os.path.join(tmpPath, 'step2_reachraster.tif'), + 'cross1': os.path.join(tmpPath, 'step2_cross1.tif'), + } + exportRasters(rastersForWorkers) + + # save main grass env which is being overriden later + MAIN_GISRC = os.environ['GISRC'] + + # build the environments and load exported rasters in each of them + grassDbPath = os.path.join(configFileDir, 'grass_db') + for i in range(nbProc): + location = 'hru-delin_%s' % (i+1) + buildGrassLocation(grassDbPath, location) + # set projection + # TODO test with a raster we want to pass to // + os.environ['GISRC'] = os.path.join(grassDbPath, 'grassdata', location, '.grassrc') + dem_recl = os.path.join(directory_out, 'step1_dem_reclass.tif') + grass_run_command('g.proj', flags='c', georef=dem_recl, stdout=DEVNULL, stderr=DEVNULL) + + importRastersInEnv(rastersForWorkers, grassDbPath, location) + + nbReachs = len(reach_ids) + if generator: + print('Starting reach loop with %s process' % nbProc) + with Pool(nbProc) as p: + params = [(id, configFileDir, nbProc) for (i, id) in enumerate(reach_ids)] + results = [] + for i, _ in enumerate(p.imap_unordered(processReachStats, params), 1): + results.append(_) + loopProgress = i/nbReachs*100 + globalProgress = 25 + (loopProgress/100*60) + yield globalProgress + else: + # this is the interesting part, launching N processes in parallel to process basins + # the locks are here to prevent concurrent terminal tqdm writing + with Pool(nbProc, initializer=tqdm.set_lock, initargs=(tqdm.get_lock(),)) as p: + params = [(id, configFileDir, nbProc) for (i, id) in enumerate(reach_ids)] + results = list(tqdm(p.imap_unordered(processReachStats, params), + desc='[main process] get reach id => subbasins id [%s process] ' % nbProc, + total=nbReachs, + unit='reach', + bar_format=bar_format1 + )) + + # merge results + for r in results: + reach_ids_cleaned.append(r[0]) + subbasins_cleaned.append(r[1]) + + # restore main grass env + os.environ['GISRC'] = MAIN_GISRC + + print('') + grass_run_command('g.remove', flags='f', type='raster', name='MASK', stdout=DEVNULL, stderr=DEVNULL) + + print('---------- Reclassifying (watersheds*100000+subbasins)') + pReclass = grass_feed_command('r.reclass', input='cross1', output='cross1_reclassed', rules='-', overwrite='True') + for k in range(0, len(reach_ids_cleaned)): + pReclass.stdin.write(encode('%s=%s\n' % (subbasins_cleaned[k], reach_ids_cleaned[k]))) + pReclass.stdin.close() + pReclass.wait() + + + if generator: + yield 90 + + print('---------- Saving \'step2_subbasins_2.tif\'') + try: + grass_run_command('r.out.gdal', input='cross1_reclassed', output=os.path.join(directory_out, 'step2_subbasins_2.tif'), + overwrite='True', stdout=DEVNULL, stderr=DEVNULL) + except: + sys.exit('------------> ERROR : Too many reclass category; check the value of ID of gauges and/or dams') + + print('---------- Creating vector layers from raster layers ... ') + grass_run_command('r.to.vect', flags='v', quiet=True, input='subbasins', output='subbasins_vector', type='area', overwrite='True') + grass_run_command('r.to.vect', flags='v', quiet=True, input='watersheds', output='watersheds_vector', type='area', overwrite='True') + ## we could clean the vector data in case isolated pixels in edges produce false small areas + ## (with no category because it comes from an isolated nodata pixel in watershed.tif) + #rasterWsheds = gdal.Open(basins) + #ulx, xres, xskew, uly, yskew, yres = rasterWsheds.GetGeoTransform() + #pixelArea = abs(xres) * abs(yres) + #grass_run_command('v.clean', + # #flags='v', + # quiet=True, + # input='watersheds_vector', + # type='area', + # tool='rmarea', + # threshold=pixelArea, + # output='watersheds_vector_clean', + # overwrite='True') + grass_run_command('r.to.vect', flags='v', quiet=True, input='cross1', output='cross1_vector', type='area', overwrite='True') + grass_run_command('r.to.vect', flags='v', quiet=True, input='cross1_reclassed', output='cross1_reclassed_vector', type='area', overwrite='True') + + if generator: + yield 95 + + # to make sure there is a projection in exported files (grass78 complains) + grass_run_command('g.proj', flags='c', georef=drain_layer, stdout=DEVNULL, stderr=DEVNULL) + grass_run_command('v.out.ogr', + #flags='c', + quiet=True, + overwrite='True', + input='subbasins_vector', type='area', format='ESRI_Shapefile', output=os.path.join(directory_out, 'step2_step1_subbasins.shp')) + # we avoid c flag to skip areas with no category (result of vectorisation error) + grass_run_command('v.out.ogr', + #flags='c', + quiet=True, + overwrite='True', + input='watersheds_vector', type='area', format='ESRI_Shapefile', output=os.path.join(directory_out, 'step2_step2_watersheds.shp')) + grass_run_command('v.out.ogr', + #flags='c', + quiet=True, + overwrite='True', + input='cross1_vector', type='area', format='ESRI_Shapefile', output=os.path.join(directory_out, 'step2_step2_subbasins.shp')) + grass_run_command('v.out.ogr', + #flags='c', + quiet=True, + overwrite='True', + input='cross1_reclassed_vector', type='area', format='ESRI_Shapefile', output=os.path.join(directory_out, 'step2_step2_subbasins_2.shp')) + + + #TEST EXIST FILES AND FILL FILES + print('---------- HRU-delin Step 2-2 : Test of existing and completed files') + + #step2_step1_subbasins.shp + shp_subb_path = os.path.join(directory_out, 'step2_step1_subbasins.shp') + + if os.stat(shp_subb_path).st_size == 0: + print('--------------- step2_step1_subbasins.shp is empty or nonexistent') + else: + datasource_subb1 = ogr.Open(shp_subb_path) + layer_datasource_subb1 = datasource_subb1.GetLayer() + featureCount_datasource_subb1 = layer_datasource_subb1.GetFeatureCount() + if featureCount_datasource_subb1 > 0 : + print('--------------- step2_step1_subbasins.shp is created and it has ', featureCount_datasource_subb1, " features") + else : + print("--------------- step2_step1_subbasins.shp is created but it empty") + + #step2_step2_watersheds.shp + shp_watershed_path = os.path.join(directory_out, 'step2_step2_watersheds.shp') + + if os.stat(shp_watershed_path).st_size == 0: + print('--------------- step2_step2_watersheds.shp is empty or nonexistent') + else: + datasource_watershed2 = ogr.Open(shp_watershed_path) + layer_datasource_watershed2 = datasource_watershed2.GetLayer() + featureCount_datasource_watershed2 = layer_datasource_watershed2.GetFeatureCount() + if featureCount_datasource_watershed2 > 0 : + print('--------------- step2_step2_watersheds.shp is created and it has ', featureCount_datasource_watershed2, " features") + else : + print("--------------- step2_step2_watersheds.shp is created but it empty") + + #step2_step2_subbasins.shp + shp_subb2_path = os.path.join(directory_out, 'step2_step2_subbasins.shp') + + if os.stat(shp_subb2_path).st_size == 0: + print('--------------- step2_step2_subbasins.shp is empty or nonexistent') + else: + datasource_subb2 = ogr.Open(shp_subb2_path) + layer_datasource_subb2 = datasource_subb2.GetLayer() + featureCount_datasource_subb2 = layer_datasource_subb2.GetFeatureCount() + if featureCount_datasource_subb2 > 0 : + print('--------------- step2_step2_subbasins.shp is created and it has ', featureCount_datasource_subb2, " features") + else : + print("--------------- step2_step2_subbasins.shp is created but it empty") + + #step2_step2_subbasins_2.shp + shp_subb2_2_path = os.path.join(directory_out, 'step2_step2_subbasins_2.shp') + + if os.stat(shp_subb2_2_path).st_size == 0: + print('--------------- step2_step2_subbasins.shp is empty or nonexistent') + else: + datasource_subb2_2 = ogr.Open(shp_subb2_2_path) + layer_datasource_subb2_2 = datasource_subb2_2.GetLayer() + featureCount_datasource_subb2_2 = layer_datasource_subb2_2.GetFeatureCount() + if featureCount_datasource_subb2_2 > 0 : + print('--------------- step2_step2_subbasins_2.shp is created and it has ', featureCount_datasource_subb2_2, " features") + else : + print("--------------- step2_step2_subbasins_2.shp is created but it empty") + + + + print('---------- HRU-delin Step 2 ended ---------------------------------------------') + + + +if __name__ == '__main__': + from grassUtils import buildGrassEnv, buildGrassLocation, exportRasters, importRastersInEnv,\ + grass_run_command, grass_parse_command, grass_feed_command, grass_read_command, grass_pipe_command + from progressColors import * + # check TQDM presence only if we are executed + try: + from tqdm import tqdm + except Exception as e: + print('!! %stqdm module not found%s\n' % (COLOR_RED, COLOR_RESET)) + sys.exit(1) + + parms_file = 'hrudelin_config.cfg' + nbProcArg = '' + if len(sys.argv) > 1: + parms_file = sys.argv[1] + if len(sys.argv) > 2: + nbProcArg = sys.argv[2] + + # determine how many processes we can launch + if str(nbProcArg).isnumeric() and int(nbProcArg) > 0: + nbProc = int(nbProcArg) + else: + nbProc = cpu_count() + + # main is a generator but we don't use it here + for pc in main(parms_file, nbProc, False): + pass + + try: + os.system('notify-send "hru-delin-6-2 step 2-6 complete"') + except Exception as e: + pass +else: + from .grassUtils import buildGrassEnv, buildGrassLocation, exportRasters, importRastersInEnv,\ + grass_run_command, grass_parse_command, grass_feed_command, grass_read_command, grass_pipe_command + from .progressColors import * diff --git a/modules/hrudelin_2_7_isolate_pixel.py b/modules/hrudelin_2_7_isolate_pixel.py new file mode 100755 index 0000000000000000000000000000000000000000..f3837933a204f7eb0f66b9afd13885d3bb7baae9 --- /dev/null +++ b/modules/hrudelin_2_7_isolate_pixel.py @@ -0,0 +1,236 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + + +############################################################################ +# +# MODULE: hru-delin_basins.py +# AUTHOR(S): adapted from GRASS-HRU (ILMS) - JENA University +# by IRSTEA - Christine Barachet, +# Julien Veyssier +# Michael Rabotin +# Florent Veillon +# PURPOSE: 1. Relocates the gauges on the reaches +# 2. Calculates watersheds at the gauges +# +# +# COPYRIGHT: (C) 2020 UR RIVERLY - INRAE +# +# This program is free software under the GNU General Public +# License (>=v2). Read the file LICENSE that comes with +# HRU-DELIN for details. +# +############################################################################# + + + + +# to keep python2 compatibility +from __future__ import print_function +import string, os, sys, glob, types, time, platform +import numpy as np +try: + import ConfigParser +except Exception as e: + import configparser as ConfigParser +#import grass.script as grass +from grass.script.utils import decode, encode +import struct, math, csv, shutil + +from osgeo import gdal +from osgeo.gdalnumeric import * +from osgeo.gdalconst import * +from osgeo import ogr + +import multiprocessing +from multiprocessing import Pool, cpu_count + +from utils import isint, write_log +from reach import snapping_points_to_reaches, cut_streams_at_points +from reach import updateAttributeTable, processReachStats + +MY_ABS_PATH=os.path.abspath(__file__) +MY_DIR=os.path.dirname(MY_ABS_PATH) + +try: + # Python 3 + from subprocess import DEVNULL +except ImportError: + DEVNULL = open(os.devnull, 'wb') + +import pandas as pd +#pd.options.mode.chained_assignment = None +import geopandas as gpd +from rastertodataframe import raster_to_dataframe +import rtree +import pygeos + +import warnings +warnings.simplefilter(action='ignore', category=FutureWarning) + + +''' + + MAIN + +''' +def main(parms_file, nbProc, generator=False): + + """OUTPUT files + + """ + print(" ") + print('---------- HRU-delin Step 2-7 started ---------------------------------------------') + print("-----------------------------------------------------------------------------------") + + configFileDir = os.path.dirname(parms_file) + parms = ConfigParser.ConfigParser(allow_no_value=True) + tmpPath = os.path.join(configFileDir, 'tmp') + if not os.path.isdir(tmpPath): + os.mkdir(tmpPath) + parms.read(parms_file) + directory_out = parms.get('dir_out', 'files') + # manage absolute and relative paths + if not os.path.isabs(directory_out): + directory_out = os.path.join(configFileDir, directory_out) + #Set Grass environnement + os.environ['GISRC'] = os.path.join(configFileDir, 'grass_db', 'grassdata', 'hru-delin', '.grassrc') + + # Import step2_subbasins_2.tif + step2_subbasins_layer = os.path.join(directory_out, 'step2_subbasins_2.tif') + step2_subbasins_wk = 'step2_subbasins' + grass_run_command('r.in.gdal', flags='o', input=step2_subbasins_layer, output=step2_subbasins_wk, overwrite='True', stdout=DEVNULL, stderr=DEVNULL) + grass_run_command('g.proj', flags='p', georef=step2_subbasins_layer, stdout=DEVNULL, stderr=DEVNULL) + grass_run_command('g.region', flags='sp', raster=step2_subbasins_wk, stdout=DEVNULL, stderr=DEVNULL) + + #STEP 1 : Raster to vector + print('---------- Creating vector layers from raster layers ... ') + grass_run_command('r.to.vect', flags='v', quiet=True, input='step2_subbasins', output='step2_subbasins_vector', type='area', overwrite='True') + grass_run_command('v.out.ogr',quiet=True,input='step2_subbasins_vector', type='area', format='ESRI_Shapefile', output=os.path.join(directory_out, 'step2_subbasins_2_vector.shp'), overwrite='True') + + #STEP 2 : Search shape with same ID + print('---------- Search shape with same ID ... ') + subbasins_vector = gpd.read_file(os.path.join(directory_out, 'step2_subbasins_2_vector.shp')) + subbasins_vector_sameID = subbasins_vector[subbasins_vector.groupby(['cat'])['geometry'].transform('nunique') > 1] + print(subbasins_vector_sameID) + #subbasins_vector_sameID['index1'] = subbasins_vector_sameID.index + subbasins_vector_sameID = subbasins_vector_sameID.rename_axis('index1').reset_index() + + + print(subbasins_vector_sameID) + + #STEP 3 : Area calculation + print('---------- Area calculation ... ') + #create new column with area of each shape + subbasins_vector_sameID = subbasins_vector_sameID.assign(area=subbasins_vector_sameID.area) + print(subbasins_vector_sameID) + + #STEP 4 : Identification of smaller layer for each ID + print('---------- Identification of smaller layer for each ID ... ') + #Extraction of ID and row line of layer with smaller area + pixel_sameID = subbasins_vector_sameID.groupby('cat', as_index=False)['area'].idxmin() + #Set new column with correspondance or not (True/False) + subbasins_vector_sameID["single"] = subbasins_vector_sameID.index1.isin(pixel_sameID.area) + #Subset df with isolate pixel + single_pixels = subbasins_vector_sameID.loc[subbasins_vector_sameID['single'] == True,:] + #print(single_pixels) + #Export single pixels to .shp + single_pixels.to_file(os.path.join(directory_out,"step2_single_pixels.shp")) + + #STEP 5 : Make buffer around single pixels + print('---------- Make buffer ... ') + single_pixels_shp = gpd.read_file(os.path.join(directory_out, 'step2_single_pixels.shp')) + buffer_pixels = single_pixels_shp.buffer(10) + #Export buffer pixels to .shp + buffer_pixels.to_file(os.path.join(directory_out,"step2_buffer_pixels.shp")) + + #STEP 6 : Intersection between buffer pixel and subbasins_vector + print('---------- Make intersection ... ') + #buffer pixel + buffer_pixels_shp = gpd.read_file(os.path.join(directory_out, 'step2_buffer_pixels.shp')) + #subbasin vector + subbasins_vector_shp = subbasins_vector + intersection = gpd.overlay(subbasins_vector_shp, buffer_pixels_shp, how='intersection') + #print(intersection) + + #STEP 7 : identification of single error pixel + print('---------- Make identification ... ') + #group by FID (1 FID for part of buffer) + test = intersection.groupby('FID') + + #set empty geopandasdataframe + df_error_pixel = gpd.GeoDataFrame() + #group by iteration + for name, group in test: + lenght_same_ID = len(group) + lenght_unique = len(group['cat'].unique()) +# print(lenght_same_ID ) +# print(lenght_unique ) + #if the buffer intersect only one shape with the same id + if lenght_same_ID ==2 and lenght_unique == 1: + group = group.assign(area=group.area) + #print(group) + group = group[group.area == group.area.max()] + df_error_pixel = df_error_pixel.append(group) + + #df_error_pixel["ID_gauges"] = np.nan + #print(df_error_pixel) + + print('---------- Isolate pixels ? ... ') + if len(df_error_pixel) >=1: + print('---------------- Yes : ') + for index, row in df_error_pixel.iterrows(): + print("---------------- ID of subbassins with isolate pixel : ", row["cat"]) + df_error_pixel.to_file(os.path.join(directory_out,"step2_error_pixels_subbasins.shp")) + else : + print('---------------- No') + +# #STEP 8 : Research for the nearest gauges +# print('---------- Research for the nearest gauges ... ') +# gauges_selected_newID = gpd.read_file(os.path.join(directory_out, 'gauges_selected.shp')) +# #Recherche du plus proche voisin +# #ATTENTION : formatage du fichier d'origine gauge : mettre ID +# test = gpd.sjoin_nearest(df_error_pixel, gauges_selected_newID) +# print(test.gid) + + + print('---------- HRU-delin Step 2-7 ended ---------------------------------------------') + + + +if __name__ == '__main__': + from grassUtils import buildGrassEnv, buildGrassLocation, exportRasters, importRastersInEnv,\ + grass_run_command, grass_parse_command, grass_feed_command, grass_read_command, grass_pipe_command + from progressColors import * + # check TQDM presence only if we are executed + try: + from tqdm import tqdm + except Exception as e: + print('!! %stqdm module not found%s\n' % (COLOR_RED, COLOR_RESET)) + sys.exit(1) + + parms_file = 'hrudelin_config.cfg' + nbProcArg = '' + if len(sys.argv) > 1: + parms_file = sys.argv[1] + if len(sys.argv) > 2: + nbProcArg = sys.argv[2] + + # determine how many processes we can launch + if str(nbProcArg).isnumeric() and int(nbProcArg) > 0: + nbProc = int(nbProcArg) + else: + nbProc = cpu_count() + + # main is a generator but we don't use it here + # for pc in main(parms_file, nbProc, False): + # pass + + try: + os.system('notify-send "hru-delin-6-2 step 2-7 complete"') + except Exception as e: + pass +else: + from .grassUtils import buildGrassEnv, buildGrassLocation, exportRasters, importRastersInEnv,\ + grass_run_command, grass_parse_command, grass_feed_command, grass_read_command, grass_pipe_command + from .progressColors import * diff --git a/modules/hrudelin_3_hrugen.py b/modules/hrudelin_3_hrugen.py index 5354dae645cc998c452a106c332c2d6ab35fd05c..e5cb17c5b332e936de86eced26e3e78df5132d7a 100755 --- a/modules/hrudelin_3_hrugen.py +++ b/modules/hrudelin_3_hrugen.py @@ -7,6 +7,8 @@ # AUTHOR(S): adapted from GRASS-HRU (ILMS) - JENA University # by IRSTEA - Christine Barachet, # Julien Veyssier +# Michael Rabotin +# Florent Veillon # # PURPOSE: overlay of all selected layers # construction of HRUs @@ -160,7 +162,7 @@ def processSubbasin(params): # LOOP THAT ELIMINATES PIXELS ONLY sizes = range(2, min_area + 1) - print("ligne 163") + if generator: iterable2 = sizes else: @@ -175,10 +177,10 @@ def processSubbasin(params): ) iterable2.set_description('[process %s] basin [%s/%s] id %s single pixels %s' % (processN, padLeft(str(iRun), len(str(nbBasins))), nbBasins, pad(id.strip(), 6), pad('', 1))) - print("ligne 178") + initmap = 'clumps' counter_old = 0 - print("ligne 181") + while (True): counter = count_only(initmap, 1) @@ -209,7 +211,7 @@ def processSubbasin(params): grass_run_command('r.buffer', input='newmap_sz1', output='buffer', distances='%d'%buffer_distance, overwrite='True',stdout=DEVNULL, stderr=DEVNULL) - print("ligne 212") + 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 @@ -229,7 +231,7 @@ def processSubbasin(params): else: df3=df1 - print("ligne 232") + #out3Path = os.path.join(tmpPath, 'out3_proc%s' % processN) pReclass = grass_feed_command('r.reclass', input='newmap', output='test', rules='-', overwrite='True') @@ -252,7 +254,7 @@ def processSubbasin(params): pReclass.stdin.write(encode('2 = 2 9999999\n')) pReclass.stdin.close() pReclass.wait() - print("ligne 255") + grass_run_command('r.patch', input='buf_mask_new,test', output='sum_buf', overwrite='True', stdout=DEVNULL, stderr=DEVNULL) grass_run_command('r.mapcalc', @@ -263,12 +265,12 @@ def processSubbasin(params): grass_run_command('r.mask', raster='subbasins', maskcats=id.rstrip('\n'), overwrite='True', stdout=DEVNULL, stderr=DEVNULL) for sz in iterable2: - print("ligne 266") + 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)) ) - print("ligne 271") + 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... @@ -289,19 +291,23 @@ def processSubbasin(params): expression='map_d=if((not(isnull(base))),eval(a=base[0,-1],b=base[-1,0],c=base[0,1],d=base[1,0],i=base,max(a,b,c,d,i)))', overwrite='True', stdout=DEVNULL, stderr=DEVNULL) grass_run_command('r.statistics', base='base', cover='map_d', method='mode', output='map1', overwrite='True', stdout=DEVNULL, stderr=DEVNULL) - print("ligne 292") + #os.system("r.stats -lnN input=map1 | tail -n +2 | awk '{print $1 \" = \" $2}' | r.reclass --o input=base output=base_new rules=-") p = grass_pipe_command('r.stats', quiet=True, flags='lnN', input='map1') pReclass = grass_feed_command('r.reclass', overwrite=True, input='base', output='base_new', rules='-') # skip first line "0 0" - l = p.stdout.readline() + + # l = p.stdout.readline() + + + for l in p.stdout: lSpl = decode(l).rstrip(os.linesep).split() pReclass.stdin.write(encode('%s = %s\n' % (lSpl[0], lSpl[1]))) p.wait() pReclass.stdin.close() pReclass.wait() - print("ligne 304") + grass_run_command('r.patch', 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? @@ -492,27 +498,28 @@ def main(parms_file, nbProc, generator=True): data_list = [] mask_list = [] - if parms.get('layer_overlay', 'dem') == 'x': + if parms.get('layer_overlay', 'dem') == 'yes': print('----------------------------- Importing raster \'step1_dem_reclass.tif\'') dem = 'dem' grass_run_command('r.in.gdal', flags='o', input=dem_recl, output=dem, overwrite='True', stdout=DEVNULL, stderr=DEVNULL) data_list.append('dem') mask_list.append('dem_msk') - if parms.get('layer_overlay', 'slope') == 'x': + if parms.get('layer_overlay', 'slope') == 'yes': print('----------------------------- Importing raster \'step1_slope_reclass.tif\'') slp_recl = os.path.join(directory_out, 'step1_slope_reclass.tif') grass_run_command('r.in.gdal', flags='o', input=slp_recl, output='slp', overwrite='True', stdout=DEVNULL, stderr=DEVNULL) data_list.append('slp') mask_list.append('slp_msk') - if parms.get('layer_overlay', 'aspect') == 'x': + if parms.get('layer_overlay', 'aspect') == 'yes': print('----------------------------- Importing raster \'step1_aspect_reclass.tif\'') asp_recl = os.path.join(directory_out, 'step1_aspect_reclass.tif') grass_run_command('r.in.gdal', flags='o', input=asp_recl, output='asp', overwrite='True', stdout=DEVNULL, stderr=DEVNULL) data_list.append('asp') mask_list.append('asp_msk') + # read sub-basin raster created in step 1 print('----------------------------- Importing raster \'step2_subbasins_2.tif\'') # new raster I. H. @@ -784,7 +791,7 @@ if __name__ == '__main__': res = list(main(parms_file, nbProc, False)) try: - os.system('notify-send "hru-delin step 3 complete"') + os.system('notify-send "hru-delin-6-2 step 3 complete"') except Exception as e: pass else: diff --git a/modules/hrudelin_parms_J2000.py b/modules/hrudelin_parms_J2000.py old mode 100644 new mode 100755 index bda4c0ffc7824ac9c5c6696212215633fa9b8875..13a804dd3022c2f1262465394c80982ef79612ab --- a/modules/hrudelin_parms_J2000.py +++ b/modules/hrudelin_parms_J2000.py @@ -8,6 +8,8 @@ # AUTHOR(S): adapted from GRASS-HRU (ILMS) - JENA University # by IRSTEA - Christine Barachet, # Julien Veyssier +# Michael Rabotin +# Florent Veillon # # PURPOSE: calculates the topology # generates parameters files for J2000 @@ -30,7 +32,7 @@ from __future__ import print_function import string, os, sys, time, types, shutil import zipfile -import pyGate +import pythonGate as pyGate try: import ConfigParser except Exception as e: @@ -493,7 +495,7 @@ def main(parms_file, nbProc, generator=False): # Find maximum destination rate (need to be updated because of additiona 2 line header # output adapted to requirements for Bjoern tool (see next) - if (parms.get('topology', 'dissolve_cycle')) == 'y': + if (parms.get('topology', 'dissolve_cycle')) == 'yes': print('------------- Dissolve cycles -------------', time.strftime('%a, %d %b %Y, %H:%M:%S', time.localtime())) # adaptation N1 topoToHruPath = os.path.join(tmpPath, 'topologie_to_hru.par') @@ -1059,15 +1061,15 @@ def main(parms_file, nbProc, generator=False): ## running irrigation programs in R using pyGate by Theo L. - execFile = pyGate.ExecFile("/home/michael.rabotin/1_HYBV/HRU_DELIN/hru-delin-dev/rScript/", "MDR_areaselect_irrigated_HRUs.r", "Rscript ") - sndStub = pyGate.SenderStub(execFile) - dictionaire = pyGate.Dictionary("irrigation_files") - dictionaire.addParameter("output_dir", dir_results) - dictionaire.addParameter("hruFile", parms.get("irrigation_analysis", "HRU_file")) #passe le chemin depuis le fichier de configuration du fichier HRU.dbf - dictionaire.addParameter("cantonFile", parms.get("irrigation_analysis", "cantons_file")) #passe le chemin depuis le fichier de configuration du fichier cantons.dbf + # execFile = pyGate.ExecFile("/home/michael.rabotin/1_HYBV/HRU_DELIN/hru-delin-dev/rScript/", "MDR_areaselect_irrigated_HRUs.r", "Rscript ") + #sndStub = pyGate.SenderStub(execFile) + # dictionaire = pyGate.Dictionary("irrigation_files") + # dictionaire.addParameter("output_dir", dir_results) + # dictionaire.addParameter("hruFile", parms.get("irrigation_analysis", "HRU_file")) #passe le chemin depuis le fichier de configuration du fichier HRU.dbf + # dictionaire.addParameter("cantonFile", parms.get("irrigation_analysis", "cantons_file")) #passe le chemin depuis le fichier de configuration du fichier cantons.dbf - sndStub.dictionaries.append(dictionaire) - sndStub.run("/home/michael.rabotin/1_HYBV/HRU_DELIN/hru-delin-dev/modules/pyGate/", "gate.exe") + # sndStub.dictionaries.append(dictionaire) + # sndStub.run("/home/michael.rabotin/1_HYBV/HRU_DELIN/hru-delin-dev/gateway-main/gateway/build/", "gate.exe") if not generator: @@ -1110,7 +1112,7 @@ if __name__ == '__main__': res = list(main(parms_file, nbProc, False)) try: - os.system('notify-send "hru-delin step 4 complete"') + os.system('notify-send "hru-delin-6-2 step 4 complete"') except Exception as e: pass else: