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

use pipe for each reclass in step4


Signed-off-by: default avatarJulien Veyssier <eneiluj@posteo.net>
parent cf3b1677
......@@ -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')
......
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