diff --git a/modules/hrudelin_parms_J2000.py b/modules/hrudelin_parms_J2000.py index 2d8696fc85481c2c0d24e5c39c7cea1b36274340..8c68d024995d6fc21d4c728a298111abbe83fdd2 100755 --- a/modules/hrudelin_parms_J2000.py +++ b/modules/hrudelin_parms_J2000.py @@ -21,6 +21,7 @@ try: except Exception as e: import configparser as ConfigParser import grass.script as grass +from grass.script.utils import decode, encode import numpy as np from osgeo import gdal @@ -110,14 +111,14 @@ def processSubbasin(params): (processN, padLeft(str(iRun), len(str(nbBasins))), nbBasins, pad(id.strip(), 6), pad('', 4))) # GET HRU AREAS #os.system("r.stats -cnN input=hrus | awk '{print $1 \"=\" $1\" \"$2'} | r.reclass --o input=hrus output=hru_sz rules=-") - # get hru areas pure python - stats = grass.read_command('r.stats', flags='cnN', input='hrus').decode('utf-8').strip().split('\n') - hruRulesPath = os.path.join(tmpPath, 'hruRules%s' % processN) - with open(hruRulesPath, 'w') as hs: - for l in stats: - spl = l.split() - hs.write('%s = %s %s\n' % (spl[0], spl[0], spl[1])) - grass.run_command('r.reclass', overwrite='True', input='hrus', output='hru_sz', rules=hruRulesPath) + p = grass.pipe_command('r.stats', flags='cnN', input='hrus') + pReclass = grass.feed_command('r.reclass', overwrite='True', input='hrus', output='hru_sz', rules='-') + for l in p.stdout: + spl = decode(l).rstrip(os.linesep).split() + pReclass.stdin.write(encode('%s = %s %s\n' % (spl[0], spl[0], spl[1]))) + p.wait() + pReclass.stdin.close() + pReclass.wait() # CREATE OVERLAY MAP FOR FURTHER PROCESSING # there is a missing category in grass7 r.cross -z result. Behaviour of -z option changed # so now we need to remove -z to get all combinations with NULL and then reclassify to put them all in cat 0 @@ -126,25 +127,23 @@ def processSubbasin(params): output='res1', overwrite='True', stdout=DEVNULL, stderr=DEVNULL) # generate reclass rules - res1Lines = grass.read_command('r.stats', flags='l', input='res1').decode('utf-8').strip().split('\n') - reclassRes1Path = os.path.join(tmpPath, 'reclassRes1_proc%s' % processN) - with open(reclassRes1Path, 'w') as res1ReclassRules: - c = 1 - for l in res1Lines: - spl = l.split() - res1Cat = spl[0] - res1Label = ' '.join(spl[1:]) - # used to be this but does not work in grass74 because result of r.cross is slightly different than grass >= 76... - #if 'NULL' in l: - if not l.startswith('*') and ('NULL' in l or 'no data' in l or l.strip() == '0'): - res1ReclassRules.write('%s = 0\n' % (res1Cat)) - else: - res1ReclassRules.write('%s = %s %s\n' % (res1Cat, c, res1Label)) - c += 1 - - grass.run_command('r.reclass', - input='res1', output='res2', rules=reclassRes1Path, - overwrite='True', stdout=DEVNULL, stderr=DEVNULL) + p = grass.pipe_command('r.stats', flags='l', input='res1') + pReclass = grass.feed_command('r.reclass', input='res1', output='res2', rules='-', overwrite='True') + c = 1 + for l in p.stdout: + spl = decode(l).rstrip(os.linesep).split() + res1Cat = spl[0] + res1Label = ' '.join(spl[1:]) + # used to be this but does not work in grass74 because result of r.cross is slightly different than grass >= 76... + #if 'NULL' in l: + if not l.startswith('*') and ('NULL' in l or 'no data' in l or l.strip() == '0'): + pReclass.stdin.write(encode('%s = 0\n' % (res1Cat))) + else: + pReclass.stdin.write(encode('%s = %s %s\n' % (res1Cat, c, res1Label))) + c += 1 + p.wait() + pReclass.stdin.close() + pReclass.wait() if not generator: iterable2.update(1) @@ -566,14 +565,15 @@ def main(parms_file, nbProc, generator=False): single_pixel_size = nres*sres #os.system("r.stats -cnN input=hrus | awk '{print $1\"=\"$2*%d}' | r.reclass --o input=hrus output='hrus_sz' rules=-" % single_pixel_size) - stats = grass.read_command('r.stats', flags='cnN', input='hrus').decode('utf-8').strip().split('\n') - hruRulesPath = os.path.join(tmpPath, 'hruRules2') - with open(hruRulesPath, 'w') as hs: - for l in stats: - spl = l.split() - n = int(spl[1]) * single_pixel_size - hs.write('%s=%s\n' % (spl[0], n)) - grass.run_command('r.reclass', overwrite='True', input='hrus', output='hrus_sz', rules=hruRulesPath) + p = grass.pipe_command('r.stats', flags='cnN', input='hrus') + pReclass = grass.feed_command('r.reclass', overwrite='True', input='hrus', output='hrus_sz', rules='-') + for l in p.stdout: + spl = decode(l).rstrip(os.linesep).split() + n = int(spl[1]) * single_pixel_size + pReclass.stdin.write(encode('%s=%s\n' % (spl[0], n))) + p.wait() + pReclass.stdin.close() + pReclass.wait() grass.run_command('db.execute', sql='ALTER TABLE hrus_c_pnt ADD area INTEGER') grass.run_command('v.what.rast', map='hrus_c_pnt', raster='hrus_sz', column='area') @@ -607,13 +607,14 @@ def main(parms_file, nbProc, generator=False): grass.run_command('r.mapcalc', expression='slope_mult=int(slope*1000.0)', overwrite='True') #os.system("r.stats -nN input=slope_mult | awk '{print $1\"=\"$1\" \"$1/1000}' | r.reclass --o input=slope_mult output=slope_labeled rules=-") - stats = grass.read_command('r.stats', flags='nN', input='slope_mult').decode('utf-8').strip().split('\n') - rulesPath = os.path.join(tmpPath, 'slopeMultRules') - with open(rulesPath, 'w') as ru: - for l in stats: - spl = l.split() - ru.write('%s = %s %s\n' % (spl[0], spl[0], int(spl[0]) / 1000)) - grass.run_command('r.reclass', overwrite='True', input='slope_mult', output='slope_labeled', rules=rulesPath) + p = grass.pipe_command('r.stats', flags='nN', input='slope_mult') + pReclass = grass.feed_command('r.reclass', overwrite='True', input='slope_mult', output='slope_labeled', rules='-') + for l in p.stdout: + spl = decode(l).rstrip(os.linesep).split() + pReclass.stdin.write(encode('%s = %s %s\n' % (spl[0], spl[0], int(spl[0]) / 1000))) + p.wait() + pReclass.stdin.close() + pReclass.wait() grass.run_command('r.statistics', flags='c',base=hrus, cover=slope_labeled, method='average', output=slope_avg, overwrite='True') grass.run_command('r.mapcalc', expression='slope_avg_cval=int(@slope_avg*1000)/1000.0', overwrite='True') @@ -700,28 +701,24 @@ def main(parms_file, nbProc, generator=False): # FROM N1 ROUTING FILE, JUST SKIP FIRST TWO LINES (HEADER) #os.system("tail -n +3 %s/topologie_n1.par | awk 'BEGIN{}{print $1\"=\"$2 >> \"%s/to_poly\";print $1\"=\"$3*-1 >> \"%s/to_reach\";}END{}'" % (tmpPath, tmpPath, tmpPath)) topoN1Path = os.path.join(tmpPath, 'topologie_n1.par') - toPolyPath = os.path.join(tmpPath, 'to_poly') - toReachPath = os.path.join(tmpPath, 'to_reach') + pToPoly = grass.feed_command('r.reclass', input='hrus', output='to_poly_map', rules='-', overwrite='True') + pToReach = grass.feed_command('r.reclass', input='hrus', output='to_reach_map', rules='-', overwrite='True') with open(topoN1Path, 'r') as topoN1: - with open(toPolyPath, 'w') as toPoly: - with open(toReachPath, 'w') as toReach: - for l in topoN1: - if not l.startswith('#'): - lSpl = l.strip().split() - toPoly.write('%s=%s\n' % (lSpl[0], lSpl[1])) - toReach.write('%s=%s\n' % (lSpl[0], -int(lSpl[2]))) + for l in topoN1: + if not l.startswith('#'): + lSpl = l.rstrip(os.linesep).split() + pToPoly.stdin.write(encode('%s=%s\n' % (lSpl[0], lSpl[1]))) + pToReach.stdin.write(encode('%s=%s\n' % (lSpl[0], -int(lSpl[2])))) + pToPoly.stdin.close() + pToPoly.wait() + pToReach.stdin.close() + pToReach.wait() if not generator: iterable3.update(1) else: yield 90 - grass.run_command('r.reclass', - input='hrus', output='to_poly_map', rules=os.path.join(tmpPath, 'to_poly'), - overwrite='True', stdout=DEVNULL, stderr=DEVNULL) - grass.run_command('r.reclass', - input='hrus', output='to_reach_map', rules=os.path.join(tmpPath, 'to_reach'), - overwrite='True', stdout=DEVNULL, stderr=DEVNULL) grass.run_command('db.execute', sql='ALTER TABLE hrus_c_pnt ADD to_poly INTEGER') grass.run_command('db.execute', sql='ALTER TABLE hrus_c_pnt ADD to_reach INTEGER') grass.run_command('v.what.rast', map='hrus_c_pnt', raster='to_poly_map', column='to_poly', overwrite='True')