En raison du déménagement des baies serveurs, les services gitlab.irstea.fr et mattermost.irstea.fr seront interrompus le samedi 2 octobre 2021 au matin. Ils devraient revenir à la normale dans la journée.

Unverified Commit 8523d22d authored by Julien Veyssier's avatar Julien Veyssier
Browse files

clarify and fix reclass rules generation in step1 and 2


Signed-off-by: default avatarJulien Veyssier <eneiluj@posteo.net>
parent 4ad6c679
......@@ -23,6 +23,7 @@ try:
except Exception as e:
import configparser as ConfigParser
import grass.script as grass
from grass.script.utils import decode, encode
try:
from osgeo import gdal
......@@ -201,18 +202,24 @@ def main(parms_file):
# Reclassifying rasters
# DEM
recl_dem_file = 'reclass_rules_dem'
dem_rcl = 'dem_rcl'
if (parms.get('reclass_dem', 'rules_auto_dem')) == 'yes':
pRecode = grass.feed_command('r.recode', input=dem_filled, output=dem_rcl, rules='-', quiet=True)
step_reclass = int(parms.get('reclass_dem', 'step_dem'))
alti_min = float(grass.parse_command('r.info', flags='r', map=dem_filled, stdout=DEVNULL, stderr=DEVNULL)['min'])
alti_max = float(grass.parse_command('r.info', flags='r', map=dem_filled, stdout=DEVNULL, stderr=DEVNULL)['max'])
# just a copy, still do the pipe
recl_dem_file = os.path.join(directory_out, 'reclass_rules_dem')
with open(recl_dem_file, 'w') as dem_rules:
new_max = math.floor(alti_max / step_reclass) * step_reclass + step_reclass
if alti_min < 0:
x = 1
dem_rules.write('%.1f:0:1:%d\n' % (alti_min, x))
rule = '%.1f:0:1:%d\n' % (alti_min, x)
pRecode.stdin.write(encode(rule))
dem_rules.write(rule)
new_min = 0
else:
x = 0
......@@ -220,43 +227,52 @@ def main(parms_file):
while new_min < new_max:
x += 1
# keep the decimal to make sure r.recode does not round things
dem_rules.write('%.1f:%.1f:%d:%d\n' % (new_min, new_min+step_reclass, x, x))
rule = '%.1f:%.1f:%d:%d\n' % (new_min, new_min+step_reclass, x, x)
pRecode.stdin.write(encode(rule))
dem_rules.write(rule)
new_min += step_reclass
pRecode.stdin.close()
pRecode.wait()
else:
recl_dem_file = os.path.join(configFileDir, 'reclass_rules_dem')
if not os.path.isfile(recl_dem_file):
sys.exit('------------> File reclass_rules_dem is missing')
dem_rcl = 'dem_rcl'
# r.recode does a better job with grass7, values close to interval limits are now correctly reclassified
grass.run_command('r.recode', input=dem_filled, output=dem_rcl, rules=recl_dem_file, stdout=DEVNULL, stderr=DEVNULL)
# Save the rules to directory_out
recl_dem_name = os.path.basename(recl_dem_file)
shutil.move(recl_dem_file, os.path.join(directory_out, recl_dem_name))
# r.recode does a better job with grass7, values close to interval limits are now correctly reclassified
grass.run_command('r.recode', input=dem_filled, output=dem_rcl, rules=recl_dem_file, stdout=DEVNULL, stderr=DEVNULL)
# Save the rules to directory_out
recl_dem_name = os.path.basename(recl_dem_file)
shutil.copyfile(recl_dem_file, os.path.join(directory_out, recl_dem_name))
# Slope
recl_slope_file = 'reclass_rules_slope'
slp_rcl = 'slp_rcl'
if (parms.get('reclass_slope', 'rules_auto_slope')) == 'yes':
pRecode = grass.feed_command('r.recode', input=dem_slope, output=slp_rcl, rules='-', quiet=True)
recl_slope_file = os.path.join(directory_out, 'reclass_rules_slope')
default_recl_slope_file = os.path.join(DATA_PATH, 'reclass_default_rules_slope')
slp_min = float(grass.parse_command('r.info', flags='r', map=dem_slope, stdout=DEVNULL, stderr=DEVNULL)['min'])
slp_max = float(grass.parse_command('r.info', flags='r', map=dem_slope, stdout=DEVNULL, stderr=DEVNULL)['max'])
recl_slope_file = os.path.join(DATA_PATH, 'reclass_default_rules_slope')
with open(recl_slope_file, 'r') as slp_rules:
slp_class = slp_rules.read()
with open(default_recl_slope_file, 'r') as default_slp_rules:
slp_class = default_slp_rules.read()
max_slp = str(slp_max)
slp_class_new = slp_class.replace('max_slp', max_slp)
max_class = str(int(math.ceil(slp_max)))
slp_class = slp_class_new.replace('max_class', max_class)
recl_slope_file = 'reclass_rules_slope'
# write a copy but still do the pipe
with open(recl_slope_file, 'w') as slp_rules:
slp_rules.write(slp_class)
pRecode.stdin.write(encode(slp_class))
pRecode.stdin.close()
pRecode.wait()
else:
recl_slope_file = os.path.join(configFileDir, 'reclass_rules_slope')
if not os.path.isfile(recl_slope_file):
sys.exit('------------> file reclass_rules_slope is missing')
slp_rcl='slp_rcl'
# improve some subprocess.Popen => run_command with rules=FILE
grass.run_command('r.recode', input=dem_slope, output=slp_rcl, rules=recl_slope_file, stdout=DEVNULL, stderr=DEVNULL)
recl_slope_name = os.path.basename(recl_slope_file)
shutil.move(recl_slope_file, os.path.join(directory_out, recl_slope_name))
# improve some subprocess.Popen => run_command with rules=FILE
grass.run_command('r.recode', input=dem_slope, output=slp_rcl, rules=recl_slope_file, stdout=DEVNULL, stderr=DEVNULL)
recl_slope_name = os.path.basename(recl_slope_file)
shutil.copyfile(recl_slope_file, os.path.join(directory_out, recl_slope_name))
# Aspect
recl_aspect_file = os.path.join(directory_out, 'reclass_rules_aspect')
......@@ -264,8 +280,11 @@ def main(parms_file):
recl_aspect = os.path.join(DATA_PATH, 'reclass_default_rules_aspect')
shutil.copyfile(recl_aspect, recl_aspect_file)
else:
if not os.path.isfile(recl_aspect_file):
sys.exit('------------> File %s is missing' % recl_aspect_file)
user_aspect = os.path.join(configFileDir, 'reclass_rules_aspect')
if not os.path.isfile(user_aspect):
sys.exit('------------> File %s is missing' % user_aspect)
shutil.copyfile(user_aspect, recl_aspect_file)
asp_rcl = 'asp_rcl'
grass.run_command('r.recode', input=dem_aspect, output=asp_rcl, rules=recl_aspect_file, stdout=DEVNULL, stderr=DEVNULL)
......
......@@ -530,10 +530,10 @@ def main(parms_file, nbProc, generator=False):
print('---------- HRU-delin Step 2 : Creating the mask')
grass.run_command('r.null', map=basins_rcl, setnull=0, stdout=DEVNULL, stderr=DEVNULL)
mask_tmp = 'mask_tmp'
recl_file = os.path.join(directory_out, 'reclass_masktmp_rules')
with open(recl_file, 'w') as rf:
rf.write('0 thru 10000000 = 1')
grass.run_command('r.reclass', input=basins_rcl, output=mask_tmp, overwrite='True', rules=recl_file, stdout=DEVNULL, stderr=DEVNULL)
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)
......@@ -655,17 +655,13 @@ def main(parms_file, nbProc, generator=False):
print('')
grass.run_command('g.remove', flags='f', type='raster', name='MASK', stdout=DEVNULL, stderr=DEVNULL)
reclassReachPath = os.path.join(configFileDir, 'tmp', 'reclass_reach.txt')
if os.path.isfile(reclassReachPath):
os.remove(reclassReachPath)
with open(reclassReachPath, 'w') as reclassFile:
for k in range(0, len(reach_ids_cleaned)):
reclassFile.write('%s=%s\n' % (subbasins_cleaned[k], reach_ids_cleaned[k]))
print('---------- Reclassifying (watersheds*100000+subbasins)')
grass.run_command('r.reclass',
input='cross1', output='cross1_reclassed', rules=reclassReachPath,
overwrite='True', stdout=DEVNULL, stderr=DEVNULL)
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()
print('---------- Saving \'step2_subbasins.tif\'')
grass.run_command('r.out.gdal', input='cross1_reclassed', output=os.path.join(directory_out, 'step2_subbasins_2.tif'),
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment