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

use grass wrappers in step3, step4 and reach.py


Signed-off-by: default avatarJulien Veyssier <eneiluj@posteo.net>
parent e084502e
# -*- coding: utf-8 -*-
# coding: utf-8
import os, sys, shutil
try:
# Python 3
......@@ -7,6 +7,44 @@ except ImportError:
#DEVNULL = open('trace_step3', 'wb')
DEVNULL = open(os.devnull, 'wb')
import grass.script as grass
import subprocess
import platform
isWindows = (platform.system() == 'Windows')
################# command wrappers to automatically disable "show window" for dumb windows
if isWindows:
def getSi():
si = None
if hasattr(subprocess, 'STARTUPINFO'):
si = subprocess.STARTUPINFO()
si.dwFlags |= subprocess.STARTF_USESHOWWINDOW
return si
else:
def getSi():
return None
def grass_run_command(*args, **kwargs):
kwargs['startupinfo'] = getSi()
return grass.run_command(*args, **kwargs)
def grass_parse_command(*args, **kwargs):
kwargs['startupinfo'] = getSi()
return grass.parse_command(*args, **kwargs)
def grass_feed_command(*args, **kwargs):
kwargs['startupinfo'] = getSi()
return grass.feed_command(*args, **kwargs)
def grass_read_command(*args, **kwargs):
kwargs['startupinfo'] = getSi()
return grass.read_command(*args, **kwargs)
def grass_pipe_command(*args, **kwargs):
kwargs['startupinfo'] = getSi()
return grass.pipe_command(*args, **kwargs)
################# environment utils
def exportRasters(rMap):
for name in rMap:
......
......@@ -337,9 +337,8 @@ def main(parms_file):
if __name__ == '__main__':
from progressColors import *
from parallel import buildGrassEnv
from grassUtils import buildGrassEnv, grass_run_command, grass_parse_command, grass_feed_command
from ogr2ogr import main as ogr2ogrMain
from utils import grass_run_command, grass_parse_command, grass_feed_command
try:
from tqdm import tqdm
except Exception as e:
......@@ -362,6 +361,5 @@ if __name__ == '__main__':
pass
else:
from .progressColors import *
from .parallel import buildGrassEnv
from .ogr2ogr import main as ogr2ogrMain
from .utils import grass_run_command, grass_parse_command, grass_feed_command
\ No newline at end of file
from .grassUtils import buildGrassEnv, grass_run_command, grass_parse_command, grass_feed_command
from .ogr2ogr import main as ogr2ogrMain
\ No newline at end of file
......@@ -712,9 +712,9 @@ def main(parms_file, nbProc, generator=False):
print('---------- HRU-delin Step 2 ended ---------------------------------------------')
if __name__ == '__main__':
from parallel import buildGrassEnv, buildGrassLocation, exportRasters, importRastersInEnv
from grassUtils import buildGrassEnv, buildGrassLocation, exportRasters, importRastersInEnv,\
grass_run_command, grass_parse_command, grass_feed_command, grass_read_command
from progressColors import *
from utils import grass_run_command, grass_parse_command, grass_feed_command, grass_read_command
# check TQDM presence only if we are executed
try:
from tqdm import tqdm
......@@ -744,6 +744,6 @@ if __name__ == '__main__':
except Exception as e:
pass
else:
from .parallel import buildGrassEnv, buildGrassLocation, exportRasters, importRastersInEnv
from .progressColors import *
from .utils import grass_run_command, grass_parse_command, grass_feed_command, grass_read_command
\ No newline at end of file
from .grassUtils import buildGrassEnv, buildGrassLocation, exportRasters, importRastersInEnv,\
grass_run_command, grass_parse_command, grass_feed_command, grass_read_command
from .progressColors import *
\ No newline at end of file
This diff is collapsed.
This diff is collapsed.
......@@ -2,11 +2,17 @@
import sys, os, shutil
import time
import grass.script as grass
#import grass.script as grass
from grass.script.utils import decode, encode
import multiprocessing
from multiprocessing import Pool, cpu_count
# trick to always be able to import something next to me (QGIS context in particular)
MY_ABS_PATH=os.path.abspath(__file__)
MY_DIR=os.path.dirname(MY_ABS_PATH)
sys.path.append(MY_DIR)
from grassUtils import grass_run_command, grass_parse_command, grass_feed_command, grass_read_command
try:
# Python 3
from subprocess import DEVNULL
......@@ -25,31 +31,31 @@ def getTopo(inp_vals):
topo_dict = {}
# CHANGE NULL VALUE TO 9999
grass.run_command('r.null', map='%(streams)s' % inp_vals, null=9999, stdout=DEVNULL, stderr=DEVNULL)
grass_run_command('r.null', map='%(streams)s' % inp_vals, null=9999, stdout=DEVNULL, stderr=DEVNULL)
# GET JUNCTION OF REACH (END POINT)
grass.run_command('r.mapcalc',
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=DEVNULL, stderr=DEVNULL)
grass.run_command('r.null', map='map_a', setnull='9999,0', stdout=DEVNULL, stderr=DEVNULL)
grass.run_command('r.mapcalc',
grass_run_command('r.null', map='map_a', setnull='9999,0', stdout=DEVNULL, stderr=DEVNULL)
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=DEVNULL, stderr=DEVNULL)
grass.run_command('r.cross',
grass_run_command('r.cross',
flags='z',
input='map_a,map_b',
output='map_c',
overwrite='True', stdout=DEVNULL, stderr=DEVNULL)
# 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=DEVNULL, stderr=DEVNULL)
#grass_run_command('r.null', map='map_c', setnull='0', stdout=DEVNULL, stderr=DEVNULL)
# 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 = 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]
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 = decode(grass.read_command('r.stats', quiet=True, flags='lnN', input='map_c')).rstrip(os.linesep).split(os.linesep)
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':
......@@ -76,17 +82,17 @@ def processReach(params):
# CALC LENGTH OF REACH
# TODO check why s1 and s13 are not used
s1 = decode(grass.read_command('r.mask',
s1 = decode(grass_read_command('r.mask',
raster=streams, maskcats='%s' % id,
overwrite='True', stdout=DEVNULL, stderr=DEVNULL
)).rstrip(os.linesep).split(os.linesep)
s13 = decode(grass.read_command('r.to.vect',
s13 = decode(grass_read_command('r.to.vect',
flags='v', input='MASK', output='reach_%s' % id, type='line',
overwrite='True'
)).rstrip(os.linesep).split(os.linesep)[0]
report_length = decode(grass.read_command('v.to.db',
report_length = decode(grass_read_command('v.to.db',
flags='p', map='reach_%s' % id, option='length', columns='s'
)).rstrip(os.linesep).split(os.linesep)
try:
......@@ -96,7 +102,7 @@ def processReach(params):
length = 1.0
# SINUOUSITY
report_sin = decode(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'
)).rstrip(os.linesep).split(os.linesep)
try:
......@@ -105,7 +111,7 @@ def processReach(params):
sin = 1.0
# CALC SLOPE
range = decode(grass.read_command('r.describe',
range = decode(grass_read_command('r.describe',
flags='dr', map=dem
)).rstrip(os.linesep).split(os.linesep)[0].split('thru')
try:
......@@ -116,7 +122,7 @@ def processReach(params):
except Exception:
# CASE MIN=MAX IF REACH IS POINT
slp = 0.0
grass.run_command('g.remove', flags='f', type='raster', name='MASK', stdout=DEVNULL, stderr=DEVNULL)
grass_run_command('g.remove', flags='f', type='raster', name='MASK', stdout=DEVNULL, stderr=DEVNULL)
return (id, length, sin, slp)
def buildReachPar(streamsP, watershedP, drainageP, demP, outputP, pathP, nbProcP, generator):
......@@ -139,15 +145,15 @@ def buildReachPar(streamsP, watershedP, drainageP, demP, outputP, pathP, nbProcP
'drain': drainageP
}
grass.run_command('r.mapcalc',
grass_run_command('r.mapcalc',
expression='streams_in_ws=if((%s),%s,null())' % (watershedP, streamsP),
overwrite='True', stdout=DEVNULL, stderr=DEVNULL)
inp_vals['streams'] = 'streams_in_ws'
grass.run_command('g.region', raster='%(ws)s' % inp_vals, stdout=DEVNULL, stderr=DEVNULL)
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 = decode(grass.read_command('r.info', quiet=True, flags='g', map=demP)).rstrip(os.linesep).split(os.linesep)
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'):
......@@ -159,7 +165,7 @@ def buildReachPar(streamsP, watershedP, drainageP, demP, outputP, pathP, nbProcP
inp_vals['dem'] = demP
else:
#print('DEM as FCELL map...convert to CELL...')
grass.run_command('r.mapcalc',
grass_run_command('r.mapcalc',
expression='dem_int=int(%s)' % demP,
overwrite='True', stdout=DEVNULL, stderr=DEVNULL)
inp_vals['dem'] = 'dem_int'
......@@ -169,9 +175,9 @@ def buildReachPar(streamsP, watershedP, drainageP, demP, outputP, pathP, nbProcP
#print('=== '+str(topo_dict))
# SET NULL VALUE BACK TO DEFAULT
grass.run_command('r.null', map='%(streams)s' % inp_vals, setnull='0,9999', stdout=DEVNULL, stderr=DEVNULL)
grass_run_command('r.null', map='%(streams)s' % inp_vals, setnull='0,9999', stdout=DEVNULL, stderr=DEVNULL)
# GET REACH IDs
reachIDs = decode(grass.read_command('r.category', map='%(streams)s' % inp_vals, stdout=DEVNULL, stderr=DEVNULL)).rstrip(os.linesep).split(os.linesep)
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 = {}
......@@ -189,7 +195,7 @@ def buildReachPar(streamsP, watershedP, drainageP, demP, outputP, pathP, nbProcP
inp_vals['dem']: os.path.join(path, 'tmp', 'r.r_dem.tif'),
}
for name in rastersForWorkers:
grass.run_command('r.out.gdal',
grass_run_command('r.out.gdal',
input=name,
output=rastersForWorkers[name],
overwrite='True', stdout=DEVNULL, stderr=DEVNULL)
......@@ -213,10 +219,10 @@ def buildReachPar(streamsP, watershedP, drainageP, demP, outputP, pathP, nbProcP
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=DEVNULL, stderr=DEVNULL)
grass_run_command('g.proj', flags='c', georef=fp, stdout=DEVNULL, stderr=DEVNULL)
# import rasters in envs
for name in rastersForWorkers:
grass.run_command('r.in.gdal', flags='o', input=rastersForWorkers[name], output=name, overwrite='True', stdout=DEVNULL, stderr=DEVNULL)
grass_run_command('r.in.gdal', flags='o', input=rastersForWorkers[name], output=name, overwrite='True', stdout=DEVNULL, stderr=DEVNULL)
reach_ids = topo_dict.keys()
nbReachs = len(reach_ids)
......
# coding: utf-8
import pandas as pd
import grass.script as grass
import subprocess
import platform
isWindows = (platform.system() == 'Windows')
def myJoin(f1Path, f2Path, resultPath):
df1 = pd.read_table(f1Path, delim_whitespace=True)
df2 = pd.read_table(f2Path, delim_whitespace=True)
merge = pd.merge(df1, df2)
merge.to_csv(resultPath, header=True, index=False, sep=' ')
if isWindows:
def getSi():
si = None
if hasattr(subprocess, 'STARTUPINFO'):
si = subprocess.STARTUPINFO()
si.dwFlags |= subprocess.STARTF_USESHOWWINDOW
return si
else:
def getSi():
return None
def grass_run_command(*args, **kwargs):
kwargs['startupinfo'] = getSi()
return grass.run_command(*args, **kwargs)
def grass_parse_command(*args, **kwargs):
kwargs['startupinfo'] = getSi()
return grass.parse_command(*args, **kwargs)
def grass_feed_command(*args, **kwargs):
kwargs['startupinfo'] = getSi()
return grass.feed_command(*args, **kwargs)
def grass_read_command(*args, **kwargs):
kwargs['startupinfo'] = getSi()
return grass.read_command(*args, **kwargs)
\ No newline at end of file
merge.to_csv(resultPath, header=True, index=False, sep=' ')
\ No newline at end of file
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