Skip to content

GitLab

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 29a156ac authored by Julien Veyssier's avatar Julien Veyssier
Browse files

fix all decode, remove r.reach.par 


Signed-off-by: default avatarJulien Veyssier <eneiluj@posteo.net>
parent 3da564c7
Hide whitespace changes
Inline Side-by-side
Showing with 36 additions and 382 deletions
modules/hrudelin_2_basins.py
View file @ 29a156ac
......@@ -409,7 +409,7 @@ def processReach(params):
cleanId = id.strip()
grass.run_command('r.mask', raster='reachraster', maskcats='%s' % cleanId, overwrite='True', stdout=DEVNULL, stderr=DEVNULL)
tmp_val = grass.read_command('r.stats', quiet=True, flags='nN', input='cross1').decode('utf-8').strip().split('\n')
tmp_val = decode(grass.read_command('r.stats', quiet=True, flags='nN', input='cross1')).rstrip(os.linesep).split(os.linesep)
return (cleanId, tmp_val[0].rstrip('\n'))
'''
......@@ -595,7 +595,7 @@ def main(parms_file, nbProc, generator=False):
overwrite='True', stdout=DEVNULL, stderr=DEVNULL)
print('---------- Computing links between (watersheds*100000+subbasins) and reach ids')
reach_ids = grass.read_command('r.stats', quiet=True, flags='nN', input='reachraster').decode('utf-8').strip().split('\n')
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)
......
modules/hrudelin_3_hrugen.py
View file @ 29a156ac
......@@ -44,7 +44,7 @@ def reclass(map_in, map_out, size):
pReclass.wait()
def count_only(map_in, size):
statLines = grass.read_command('r.stats', quiet=True, flags='lnNc', input=map_in).decode('utf-8').strip().split('\n')
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()
......@@ -56,7 +56,7 @@ def count_only(map_in, size):
def getAreasUpTo(map_in, size):
areas = []
statLists = []
statLines = grass.read_command('r.stats', quiet=True, flags='lnNc', input=map_in).decode('utf-8').strip().split('\n')
statLines = decode(grass.read_command('r.stats', quiet=True, flags='lnNc', input=map_in)).rstrip(os.linesep).split(os.linesep)
for l in statLines:
statLists.append(l.strip().split())
# sort on second column
......@@ -96,7 +96,7 @@ def processSubbasin(params):
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(grass.read_command('r.stats', quiet=True, flags='nNc', input='MASK').decode('utf-8').strip().split('\n')[0].split()[1])
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(
......@@ -113,7 +113,7 @@ def processSubbasin(params):
iterable3.update(1)
grass.run_command('r.statistics', base='MASK', cover='clumps', method='max', output='hruid',
overwrite='True', stdout=DEVNULL, stderr=DEVNULL)
hruid = int(grass.read_command('r.stats', quiet=True, flags='lnNc', input='hruid').decode('utf-8').strip().split('\n')[0].split()[1])
hruid = int(decode(grass.read_command('r.stats', quiet=True, flags='lnNc', input='hruid')).rstrip(os.linesep).split(os.linesep)[0].split()[1])
grass.run_command('r.mapcalc',
expression='newmap3=if(MASK,%d,null())'%hruid,
overwrite='True', stdout=DEVNULL, stderr=DEVNULL)
......@@ -165,8 +165,8 @@ def processSubbasin(params):
# there were sort/awk system calls here
# now less filesystem access, cleaner code, most importantely: portable to windows
# pandatana !
stats1 = grass.read_command('r.stats', quiet=True, flags='cnN', input='newmap').decode('utf-8').strip()
stats2 = grass.read_command('r.stats', quiet=True, flags='nN', input='newmap2').decode('utf-8').strip()
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)
......@@ -305,7 +305,7 @@ def processSubbasin(params):
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(grass.read_command('r.stats', quiet=True, flags='lnN', input='map_a').decode('utf-8').strip().split('\n')[0].strip())
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
......@@ -522,7 +522,7 @@ def main(parms_file, nbProc, generator=False):
#buffer_distance = int(math.ceil(float(os.popen("g.region -m").readlines()[4].split("=")[1])))
#buffer_distance_ew = int(math.ceil(float(os.popen("g.region -m").readlines()[5].split("=")[1])))
# g.region result parsing is now safer, it does not depend on line ordering
regionOutput = grass.read_command('g.region', flags='m').decode('utf-8').strip().split('\n')
regionOutput = decode(grass.read_command('g.region', flags='m')).rstrip(os.linesep).split(os.linesep)
buffer_distance = None
buffer_distance_ew = None
for ro in regionOutput:
......@@ -540,7 +540,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 = grass.read_command('r.stats', quiet=True, flags='nN', input='subbasins').decode('utf-8').strip().split('\n')
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)
......
modules/hrudelin_parms_J2000.py
View file @ 29a156ac
......@@ -183,7 +183,7 @@ def processSubbasin(params):
#os.system("r.stats -lnN input=res2,acc_sum,acc_max | awk -v BASIN=%s -v TTOHRU=%s_proc%s -v TBFL=%s_proc%s -f %s" %
# (new_id, topoToHruTmpPath, processN, topoBflTmpPath, processN, AWK_cmd_1))
# nextgen mzfk
stats = grass.read_command('r.stats', flags='lnN', input='res2,acc_sum,acc_max').decode('utf-8').strip().split('\n')
stats = decode(grass.read_command('r.stats', flags='lnN', input='res2,acc_sum,acc_max')).rstrip(os.linesep).split(os.linesep)
topoToHruTmpPath = os.path.join(tmpPath, 'topologie_to_hru.par_tmp_proc%s' % processN)
topoBflTmpPath = os.path.join(tmpPath, 'topologie_bfl.par_tmp_proc%s' % processN)
formatNm(new_id, topoToHruTmpPath, topoBflTmpPath, stats)
......@@ -342,14 +342,14 @@ def main(parms_file, nbProc, generator=False):
# Subbasins loop
print('------------- Retrieving subbasin ids -----------')
subbasins_msk = 'subbasins_msk'
print('r.patch ')
#print('r.patch ')
grass.run_command('r.patch',
input='mask_wk,subbasins', output=subbasins_msk,
overwrite='True', stdout=DEVNULL, stderr=DEVNULL)
print('r.null ')
#print('r.null ')
grass.run_command('r.null', map=subbasins_msk, setnull=0, stdout=DEVNULL, stderr=DEVNULL)
print('r.stats ')
basin_ids = grass.read_command('r.stats', quiet=True, flags='nN', input='subbasins_msk').decode('utf-8').strip().split('\n')
#print('r.stats ')
basin_ids = decode(grass.read_command('r.stats', quiet=True, flags='nN', input='subbasins_msk')).rstrip(os.linesep).split(os.linesep)
print('%s subbasins were found:' % len(basin_ids))
print('%s, %s, %s ... %s, %s, %s ' % (
......@@ -521,7 +521,7 @@ def main(parms_file, nbProc, generator=False):
grass.run_command('g.region', zoom=hrus)
grass.run_command('r.to.vect', input=hrus, output='hrus_v', type='area', overwrite='True')
rawTables = grass.read_command('db.tables').decode('utf-8').strip().split('\n')
rawTables = decode(grass.read_command('db.tables')).rstrip(os.linesep).split(os.linesep)
tables = [t.strip() for t in rawTables]
if 'hrus_v_cat2' in tables:
grass.run_command('db.droptable', flags='f', table='hrus_v_cat2')
......@@ -551,7 +551,7 @@ def main(parms_file, nbProc, generator=False):
sql = 'cat'
# versatile way to get resolution (does not depend on output lines order)
regionOutput = grass.read_command('g.region', flags='m').decode('utf-8').strip().split('\n')
regionOutput = decode(grass.read_command('g.region', flags='m')).rstrip(os.linesep).split(os.linesep)
nres = None
sres = None
for ro in regionOutput:
......
modules/r.reach.par deleted 100755 → 0
View file @ 3da564c7
#!/usr/bin/python3
#%Module
#% description: Creates a stream network parameter file.
#% keywords: raster, model, streams, reaches, parameter
#%End
#%option
#% key: streams
#% type: string
#% description: raster map of stream network
#% required: yes
#%end
#%option
#% key: watershed
#% type: string
#% description: raster map of watershed to limit calculation to
#% required: yes
#%end
#%option
#% key: drainage
#% type: string
#% description: raster map of drainage direction
#% required: yes
#%end
#%option
#% key: dem
#% type: string
#% description: elevation data as integer raster map
#% required: yes
#%end
#%option
#% key: path
#% type: string
#% description: path to create grass parallel env
#% required: yes
#%end
#%option
#% key: nbcores
#% type: string
#% description: number of parallel processes to launch
#% required: no
#%end
#%option
#% key: output
#% type: string
#% description: output parameter file
#% required: yes
#%end
import sys, os, shutil
import subprocess
import time
import grass.script as grass
from tqdm import tqdm
import multiprocessing
from multiprocessing import Pool, cpu_count
# KNOWN BUGS:
# - IF ONLY ONE REACH EXISTS, REACH.PAR IS EMPTY
def getOutletReachID(topo):
for i in topo.values():
if i not in topo.keys():
return i
def getTopo():
Nul = open(os.devnull, 'w')
topo_dict = {}
# CHANGE NULL VALUE TO 9999
grass.run_command('r.null', map='%(streams)s' % inp_vals, null=9999, stdout=Nul, stderr=Nul)
# GET JUNCTION OF REACH (END POINT)
grass.run_command('r.mapcalc',
expression='map_a=if(!(isnull(%(ws)s)),eval(m=min(%(streams)s[-1,1],%(streams)s[1,1],%(streams)s[1,-1],%(streams)s[-1,-1],%(streams)s[0,1],%(streams)s[1,0],%(streams)s[-1,0],%(streams)s[0,-1]),if((m!=%(streams)s),%(streams)s,null())))' % inp_vals,
overwrite='True', stdout=Nul, stderr=Nul)
grass.run_command('r.null', map='map_a', setnull='9999,0', stdout=Nul, stderr=Nul)
grass.run_command('r.mapcalc',
expression='map_b=if((map_a),if((%(drain)s==1&&%(streams)s[-1,1]!=map_a),%(streams)s[-1,1],if((%(drain)s==2&&%(streams)s[-1,0]!=map_a),%(streams)s[-1,0],if((%(drain)s==3&&%(streams)s[-1,-1]!=map_a),%(streams)s[-1,-1],if((%(drain)s==4&&%(streams)s[0,-1]!=map_a),%(streams)s[0,-1],if((%(drain)s==5&&%(streams)s[1,-1]!=map_a),%(streams)s[1,-1],if((%(drain)s==6&&%(streams)s[1,0]!=map_a),%(streams)s[1,0],if((%(drain)s==7&&%(streams)s[1,1]!=map_a),%(streams)s[1,1],if((%(drain)s==8&&%(streams)s[0,1]!=map_a),%(streams)s[0,1],null())))))))))' % inp_vals,
overwrite='True', stdout=Nul, stderr=Nul)
grass.run_command('r.cross',
flags='z',
input='map_a,map_b',
output='map_c',
overwrite='True', stdout=Nul, stderr=Nul)
# don't null 0 because now it's a valid category (new behaviour of r.cross -z)
# and the old cat zero is not there anymore
#grass.run_command('r.null', map='map_c', setnull='0', stdout=Nul, stderr=Nul)
# this is a fix for grass74 which fails to produce categories for the first line of r.cross result
# so we take it manually
a = grass.read_command('r.stats', quiet=True, flags='lnN', input='map_a').decode('utf-8').strip().split('\n')[0]
b = grass.read_command('r.stats', quiet=True, flags='lnN', input='map_b').decode('utf-8').strip().split('\n')[0]
topo = grass.read_command('r.stats', quiet=True, flags='lnN', input='map_c').decode('utf-8').strip().split('\n')
for relation in topo:
# this only happens with grass74
if relation.strip() == '0':
from_node = a.strip()
to_node = b.strip()
topo_dict[from_node] = to_node
else:
from_node = int(relation.split(' ')[2].rstrip(';'))
to_node = int(relation.split(' ')[4])
topo_dict[from_node] = to_node
outletReachID = getOutletReachID(topo_dict)
topo_dict[outletReachID] = 0
return topo_dict
def processReach(params):
id, path, streams, dem, nbProc = params
# define which grass env we are using
current = multiprocessing.current_process()
processN = int(current._identity[0])
# because multiprocess does not reset workers id numbers and we still use [1, N]
processN = (processN % nbProc) + 1
location = 'hru-delin_%s' % (processN)
os.environ['GISRC'] = os.path.join(path, 'grass_db', 'grassdata', location, '.grassrc')
# CALC LENGTH OF REACH
# TODO check why s1 and s13 are not used
s1 = grass.read_command('r.mask',
raster=streams, maskcats='%s' % id,
overwrite='True', stdout=Nul, stderr=Nul
).decode('utf-8').split('\n')
s13 = grass.read_command('r.to.vect',
flags='v', input='MASK', output='reach_%s' % id, type='line',
overwrite='True'
).decode('utf-8').split('\n')[0]
report_length = grass.read_command('v.to.db',
flags='p', map='reach_%s' % id, option='length', columns='s'
).decode('utf-8').split('\n')
try:
length = round(float(report_length[1].split('|')[1]))
except Exception:
# REACH IS A POINT STUPIDLY
length = 1.0
# SINUOUSITY
report_sin = grass.read_command('v.to.db',
flags='p', map='reach_%s' % id, option='sinuous', columns='s'
).decode('utf-8').split('\n')
try:
sin = round(float(report_sin[1].split('|')[1]), 2)
except Exception:
sin = 1.0
# CALC SLOPE
range = grass.read_command('r.describe',
flags='dr', map=dem
).decode('utf-8').split('\n')[0].split('thru')
try:
min = int(range[0])
max = int(range[1])
#slp_dict[id] = round((max-min)/len_dict[id],4)*100.0
slp = round((max-min) / length * 100.0, 4)
except Exception:
# CASE MIN=MAX IF REACH IS POINT
slp = 0.0
grass.run_command('g.remove', flags='f', type='raster', name='MASK', stdout=Nul, stderr=Nul)
return (id, length, sin, slp)
def buildGrassLocation(grassDbPath, location):
gisdbase = os.path.join(grassDbPath, 'grassdata')
mapset = 'PERMANENT'
if os.path.exists(os.path.join(grassDbPath, 'grassdata', location)):
shutil.rmtree(os.path.join(grassDbPath, 'grassdata', location))
os.mkdir(os.path.join(gisdbase, location))
os.mkdir(os.path.join(gisdbase, location, mapset))
content = '''proj: 1
zone: 0
north: 1
south: 0
east: 1
west: 0
cols: 1
rows: 1
e-w resol: 1
n-s resol: 1
top: 1
bottom: 0
cols3: 1
rows3: 1
depths: 1
e-w resol3: 1
n-s resol3: 1
t-b resol: 1
'''
with open(os.path.join(gisdbase, location, mapset, 'WIND'), 'w') as w:
w.write(content)
with open(os.path.join(gisdbase, location, mapset, 'DEFAULT_WIND'), 'w') as dw:
dw.write(content)
rc_file = os.path.join(gisdbase, location, '.grassrc')
with open(rc_file, 'w') as rc:
rc.write('GISDBASE: %s\n' % (os.path.abspath(gisdbase)))
rc.write('LOCATION_NAME: %s\n' % location)
rc.write('MAPSET: %s\n' % mapset)
def main():
#DEFINE RANGES FOR HEADER
ranges = {
'ID': [0, 999999, 'n/a'],
'to-reach': [0, 999999, 'n/a'],
'length': [0, 99999, 'm'],
'slope': [0, 90, '%'],
'rough': [0, 9999, 'n/a'],
'width': [0, 9999, 'm'],
'sinuous': [1, 999999, 'n/a']
}
Nul = open(os.devnull, 'w')
global inp_vals
inp_vals = {
'streams': os.getenv('GIS_OPT_STREAMS'),
'ws': os.getenv('GIS_OPT_WATERSHED'),
'drain': os.getenv('GIS_OPT_DRAINAGE')
}
grass.run_command('r.mapcalc',
expression='streams_in_ws=if((%s),%s,null())' % (os.getenv('GIS_OPT_WATERSHED'), os.getenv('GIS_OPT_STREAMS')),
overwrite='True', stdout=Nul, stderr=Nul)
inp_vals['streams'] = 'streams_in_ws'
grass.run_command('g.region', raster='%(ws)s' % inp_vals, stdout=Nul, stderr=Nul)
# -t option of r.info does not exist anymore and it was sketchy anyway
# now we get maptype in a safer way by checking every output lines of r.info
infos = grass.read_command('r.info', quiet=True, flags='g', map=os.getenv('GIS_OPT_DEM')).decode('utf-8').strip().split('\n')
mapType = None
for info in infos:
if info.startswith('datatype'):
mapType = info.split('=')[1].strip()
if mapType == None:
sys.exit('impossible to get map type')
if (mapType == 'CELL'):
#print('DEM already as CELL map...ok!')
inp_vals['dem'] = os.getenv('GIS_OPT_DEM')
else:
#print('DEM as FCELL map...convert to CELL...')
grass.run_command('r.mapcalc',
expression='dem_int=int(%s)' % os.getenv('GIS_OPT_DEM'),
overwrite='True', stdout=Nul, stderr=Nul)
inp_vals['dem'] = 'dem_int'
# GET TOPOLOGY INFORMATION
topo_dict = getTopo()
#print('=== '+str(topo_dict))
# SET NULL VALUE BACK TO DEFAULT
grass.run_command('r.null', map='%(streams)s' % inp_vals, setnull='0,9999', stdout=Nul, stderr=Nul)
# GET REACH IDs
reachIDs = grass.read_command('r.category', map='%(streams)s' % inp_vals, stdout=Nul, stderr=Nul).decode('utf-8').split('\n')
#print('reach ids: %s' % reachIDs)
len_dict = {}
slp_dict = {}
sin_dict = {}
rough_const = '30'
width_const = '1'
########### main loop
#for reach in topo_dict.keys():
path = os.getenv('GIS_OPT_PATH')
# export rasters that are necessary for parallel environments
rastersForWorkers = {
inp_vals['streams']: os.path.join(path, 'tmp', 'r.r_streams.tif'),
inp_vals['dem']: os.path.join(path, 'tmp', 'r.r_dem.tif'),
}
for name in rastersForWorkers:
grass.run_command('r.out.gdal',
input=name,
output=rastersForWorkers[name],
overwrite='True', stdout=Nul, stderr=Nul)
# save main grass env which is being overriden later
MAIN_GISRC = os.environ['GISRC']
# determine number of processes
nbProcArg = os.getenv('GIS_OPT_NBCORES')
if str(nbProcArg).isnumeric() and int(nbProcArg) > 0:
nbProc = int(nbProcArg)
else:
nbProc = cpu_count()
# build the environments and load exported rasters in each of them
grassDbPath = os.path.join(path, 'grass_db')
for i in range(nbProc):
location = 'hru-delin_%s' % (i+1)
buildGrassLocation(grassDbPath, location)
os.environ['GISRC'] = os.path.join(grassDbPath, 'grassdata', location, '.grassrc')
# set projection
fp = os.path.join(path, 'tmp', 'r.r_dem.tif')
grass.run_command('g.proj', flags='c', georef=fp, stdout=Nul, stderr=Nul)
# import rasters in envs
for name in rastersForWorkers:
grass.run_command('r.in.gdal', flags='o', input=rastersForWorkers[name], output=name, overwrite='True', stdout=Nul, stderr=Nul)
reach_ids = topo_dict.keys()
nbReachs = len(reach_ids)
print()
COLOR_GREEN='\x1b[32m'
STYLE_BRIGHT='\x1b[1m'
STYLE_RESET='\x1b[0m'
l_bar1 = '{desc}%s%s{percentage:3.0f}%%%s|' % (STYLE_BRIGHT, COLOR_GREEN, STYLE_RESET)
r_bar1 = '| %s%s{n_fmt}/{total_fmt}%s [{elapsed}<{remaining}, {rate_fmt}{postfix}]' % (STYLE_BRIGHT, COLOR_GREEN, STYLE_RESET)
bar_format1='%s%s{bar}%s%s' % (l_bar1, COLOR_GREEN, STYLE_RESET, r_bar1)
# the locks are here to prevent concurrent terminal tqdm writing
# this is the interesting part, launching N processes in parallel to process basins
with Pool(nbProc, initializer=tqdm.set_lock, initargs=(tqdm.get_lock(),)) as p:
params = [(id, path, inp_vals['streams'], inp_vals['dem'], nbProc) for (i, id) in enumerate(reach_ids)]
results = list(tqdm(p.imap_unordered(processReach, params),
desc='[main process] Loop on reaches [%s process] ' % nbProc,
total=nbReachs,
unit='reach',
bar_format=bar_format1
))
print()
# merge results
for r in results:
id = r[0]
len_dict[id] = r[1]
sin_dict[id] = r[2]
slp_dict[id] = r[3]
# restore main grass env
os.environ['GISRC'] = MAIN_GISRC
# CREATE .PAR FILE
with open(os.getenv('GIS_OPT_OUTPUT'), 'w') as outfile:
timestr = time.strftime('%a, %d %b %Y, %H:%M:%S', time.localtime())
outfile.write('# reach.par created at %s, r.reach.par called by HRU-delin\n' % timestr)
outfile.write('ID\tto-reach\tlength\tslope\tsinuosity\trough\twidth\n')
outfile.write(str(ranges['ID'][0])+'\t'+str(ranges['to-reach'][0])+'\t'+str(ranges['length'][0])+'\t'+str(ranges['slope'][0])+'\t'+str(ranges['sinuous'][0])+'\t'+str(ranges['rough'][0])+'\t'+str(ranges['width'][0])+'\n')
outfile.write(str(ranges['ID'][1])+'\t'+str(ranges['to-reach'][1])+'\t'+str(ranges['length'][1])+'\t'+str(ranges['slope'][1])+str(ranges['sinuous'][1])+'\t'+str(ranges['rough'][1])+'\t'+str(ranges['width'][1])+'\n')
outfile.write(ranges['ID'][2]+'\t'+ranges['to-reach'][2]+'\t'+ranges['length'][2]+'\t'+ranges['slope'][2]+'\t'+ranges['sinuous'][2]+'\t'+ranges['rough'][2]+'\t'+ranges['width'][2]+'\n')
for i in topo_dict.keys():
outfile.write(str(i)+'\t'+str(topo_dict[i])+'\t'+str(len_dict[i])+'\t'+str(slp_dict[i])+'\t'+str(sin_dict[i])+'\t'+rough_const+'\t'+width_const+'\n')
outfile.write('# end of reach.par')
if __name__ == '__main__':
Nul = open(os.devnull, 'w')
if len(sys.argv) <= 1 or sys.argv[1] != '@ARGS_PARSED@':
os.execvp('g.parser', [sys.argv[0]] + sys.argv)
else:
main()
modules/reach.py
View file @ 29a156ac
......@@ -3,6 +3,7 @@
import sys, os, shutil
import time
import grass.script as grass
from grass.script.utils import decode, encode
import multiprocessing
from multiprocessing import Pool, cpu_count
......@@ -45,10 +46,10 @@ def getTopo(inp_vals):
# this is a fix for grass74 which fails to produce categories for the first line of r.cross result
# so we take it manually
a = grass.read_command('r.stats', quiet=True, flags='lnN', input='map_a').decode('utf-8').strip().split('\n')[0]
b = grass.read_command('r.stats', quiet=True, flags='lnN', input='map_b').decode('utf-8').strip().split('\n')[0]
a = decode(grass.read_command('r.stats', quiet=True, flags='lnN', input='map_a')).rstrip(os.linesep).split(os.linesep)[0]
b = decode(grass.read_command('r.stats', quiet=True, flags='lnN', input='map_b')).rstrip(os.linesep).split(os.linesep)[0]
topo = grass.read_command('r.stats', quiet=True, flags='lnN', input='map_c').decode('utf-8').strip().split('\n')
topo = decode(grass.read_command('r.stats', quiet=True, flags='lnN', input='map_c')).rstrip(os.linesep).split(os.linesep)
for relation in topo:
# this only happens with grass74
if relation.strip() == '0':
......@@ -75,17 +76,19 @@ def processReach(params):
# CALC LENGTH OF REACH
# TODO check why s1 and s13 are not used
s1 = grass.read_command('r.mask',
s1 = decode(grass.read_command('r.mask',
raster=streams, maskcats='%s' % id,
overwrite='True', stdout=DEVNULL, stderr=DEVNULL
).decode('utf-8').split('\n')
s13 = grass.read_command('r.to.vect',
)).rstrip(os.linesep).split(os.linesep)
s13 = decode(grass.read_command('r.to.vect',
flags='v', input='MASK', output='reach_%s' % id, type='line',
overwrite='True'
).decode('utf-8').split('\n')[0]
report_length = grass.read_command('v.to.db',
)).rstrip(os.linesep).split(os.linesep)[0]
report_length = decode(grass.read_command('v.to.db',
flags='p', map='reach_%s' % id, option='length', columns='s'
).decode('utf-8').split('\n')
)).rstrip(os.linesep).split(os.linesep)
try:
length = round(float(report_length[1].split('|')[1]))
except Exception:
......@@ -93,18 +96,18 @@ def processReach(params):
length = 1.0
# SINUOUSITY
report_sin = grass.read_command('v.to.db',
report_sin = decode(grass.read_command('v.to.db',
flags='p', map='reach_%s' % id, option='sinuous', columns='s'
).decode('utf-8').split('\n')
)).rstrip(os.linesep).split(os.linesep)
try:
sin = round(float(report_sin[1].split('|')[1]), 2)
except Exception:
sin = 1.0
# CALC SLOPE
range = grass.read_command('r.describe',
range = decode(grass.read_command('r.describe',
flags='dr', map=dem
).decode('utf-8').split('\n')[0].split('thru')
)).rstrip(os.linesep).split(os.linesep)[0].split('thru')
try:
min = int(range[0])
max = int(range[1])
......@@ -144,7 +147,7 @@ def buildReachPar(streamsP, watershedP, drainageP, demP, outputP, pathP, nbProcP
grass.run_command('g.region', raster='%(ws)s' % inp_vals, stdout=DEVNULL, stderr=DEVNULL)
# -t option of r.info does not exist anymore and it was sketchy anyway
# now we get maptype in a safer way by checking every output lines of r.info
infos = grass.read_command('r.info', quiet=True, flags='g', map=demP).decode('utf-8').strip().split('\n')
infos = decode(grass.read_command('r.info', quiet=True, flags='g', map=demP)).rstrip(os.linesep).split(os.linesep)
mapType = None
for info in infos:
if info.startswith('datatype'):
......@@ -168,7 +171,7 @@ def buildReachPar(streamsP, watershedP, drainageP, demP, outputP, pathP, nbProcP
# SET NULL VALUE BACK TO DEFAULT
grass.run_command('r.null', map='%(streams)s' % inp_vals, setnull='0,9999', stdout=DEVNULL, stderr=DEVNULL)
# GET REACH IDs
reachIDs = grass.read_command('r.category', map='%(streams)s' % inp_vals, stdout=DEVNULL, stderr=DEVNULL).decode('utf-8').split('\n')
reachIDs = decode(grass.read_command('r.category', map='%(streams)s' % inp_vals, stdout=DEVNULL, stderr=DEVNULL)).rstrip(os.linesep).split(os.linesep)
#print('reach ids: %s' % reachIDs)
len_dict = {}
slp_dict = {}
......