From 19ea670125b83b1c1a5feffec55258520ce5f10d Mon Sep 17 00:00:00 2001
From: Julien Veyssier <eneiluj@posteo.net>
Date: Tue, 21 Apr 2020 05:52:04 +0200
Subject: [PATCH] use pipe for each reclass in step4
Signed-off-by: Julien Veyssier <eneiluj@posteo.net>
---
modules/hrudelin_parms_J2000.py | 111 ++++++++++++++++----------------
1 file changed, 54 insertions(+), 57 deletions(-)
diff --git a/modules/hrudelin_parms_J2000.py b/modules/hrudelin_parms_J2000.py
index 2d8696f..8c68d02 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')
--
GitLab