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

use grass wrappers in step2


Signed-off-by: default avatarJulien Veyssier <eneiluj@posteo.net>
parent 7240850c
......@@ -22,7 +22,8 @@ try:
import ConfigParser
except Exception as e:
import configparser as ConfigParser
import grass.script as grass
# use wrappers now !
#import grass.script as grass
from grass.script.utils import decode, encode
try:
......
......@@ -19,7 +19,7 @@ try:
import ConfigParser
except Exception as e:
import configparser as ConfigParser
import grass.script as grass
#import grass.script as grass
from grass.script.utils import decode, encode
import struct, math, csv, shutil
......@@ -409,9 +409,9 @@ def processReach(params):
cleanId = id.strip()
grass.run_command('r.mask', raster='reachraster', maskcats='%s' % cleanId, overwrite='True', stdout=DEVNULL, stderr=DEVNULL)
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'))
grass_run_command('r.mask', raster='reachraster', maskcats='%s' % cleanId, overwrite='True', stdout=DEVNULL, stderr=DEVNULL)
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(os.linesep))
'''
......@@ -468,9 +468,9 @@ def main(parms_file, nbProc, generator=False):
# Import drain raster
drain_layer = os.path.join(directory_out, 'step1_drain.tif')
drain_wk = 'drain_wk'
grass.run_command('r.in.gdal', flags='o', input=drain_layer, output=drain_wk, overwrite='True', stdout=DEVNULL, stderr=DEVNULL)
grass.run_command('g.proj', flags='p', georef=drain_layer, stdout=DEVNULL, stderr=DEVNULL)
grass.run_command('g.region', flags='sp', raster=drain_wk, stdout=DEVNULL, stderr=DEVNULL)
grass_run_command('r.in.gdal', flags='o', input=drain_layer, output=drain_wk, overwrite='True', stdout=DEVNULL, stderr=DEVNULL)
grass_run_command('g.proj', flags='p', georef=drain_layer, stdout=DEVNULL, stderr=DEVNULL)
grass_run_command('g.region', flags='sp', raster=drain_wk, stdout=DEVNULL, stderr=DEVNULL)
# Watersheds derivation
basins = 'basins'
......@@ -478,7 +478,7 @@ def main(parms_file, nbProc, generator=False):
# r.cross doesn't accept more than 30 layers
max_rasters = 29
i = 0
grass.run_command('r.mapcalc', expression='basins=null()', overwrite=True, stdout=DEVNULL, stderr=DEVNULL)
grass_run_command('r.mapcalc', expression='basins=null()', overwrite=True, stdout=DEVNULL, stderr=DEVNULL)
ds = ogr.Open(gauges_reloc_file)
reloc_lyr = ds.GetLayer()
nb_gauges = reloc_lyr.GetFeatureCount()
......@@ -488,7 +488,7 @@ def main(parms_file, nbProc, generator=False):
geom = gauge.GetGeometryRef()
gauge_x, gauge_y = geom.GetX(), geom.GetY()
basin_wk = 'basin_wk'
grass.run_command('r.water.outlet',
grass_run_command('r.water.outlet',
input=drain_wk, output='basin_tmp%d' % i,
coordinates='%s,%s' % (gauge_x, gauge_y),
overwrite='True',
......@@ -497,31 +497,31 @@ def main(parms_file, nbProc, generator=False):
rasters_list.append('basin_tmp%d' % i)
i += 1
if i % max_rasters == 0 or i == nb_gauges:
grass.run_command('r.cross', input=','.join(rasters_list), output=basins, overwrite='True', stdout=DEVNULL, stderr=DEVNULL)
grass.run_command('g.remove', flags='f', type='raster', pattern='basin_tmp*', stdout=DEVNULL, stderr=DEVNULL)
grass_run_command('r.cross', input=','.join(rasters_list), output=basins, overwrite='True', stdout=DEVNULL, stderr=DEVNULL)
grass_run_command('g.remove', flags='f', type='raster', pattern='basin_tmp*', stdout=DEVNULL, stderr=DEVNULL)
rasters_list = ['basins']
# with grass7 watershed ids begin at 0
# empty areas value is already NULL, no need to set it
grass.run_command('r.mapcalc', expression='basins=basins+1', overwrite='True', stdout=DEVNULL, stderr=DEVNULL)
grass_run_command('r.mapcalc', expression='basins=basins+1', overwrite='True', stdout=DEVNULL, stderr=DEVNULL)
print('---------- HRU-delin Step 2 : Reclassifying the watersheds')
################
basins_rcl = 'basins_rcl'
pRecode = grass.feed_command('r.recode', input=basins, output=basins_rcl, overwrite='True', rules='-', quiet=True)
pRecode = grass_feed_command('r.recode', input=basins, output=basins_rcl, overwrite='True', rules='-', quiet=True)
ds = ogr.Open(gauges_reloc_file)
reloc_lyr = ds.GetLayer()
for gauge in reloc_lyr:
geom = gauge.GetGeometryRef()
gauge_x, gauge_y = geom.GetX(), geom.GetY()
out = decode(grass.read_command('r.what', map=basins, coordinates='%s,%s' % (gauge_x, gauge_y)))
out = decode(grass_read_command('r.what', map=basins, coordinates='%s,%s' % (gauge_x, gauge_y)))
cat_grass = out.rstrip(os.linesep).split('|')[3].strip()
pRecode.stdin.write(encode('%s:%s:%s\n' % (cat_grass, cat_grass, int(gauge.GetField(col_name)))))
pRecode.stdin.close()
pRecode.wait()
grass.run_command('r.out.gdal',
grass_run_command('r.out.gdal',
input=basins_rcl,
type='UInt16',
output=os.path.join(directory_out, 'step2_watersheds.tif'),
......@@ -529,16 +529,16 @@ def main(parms_file, nbProc, generator=False):
# Create mask = surface of all watersheds
print('---------- HRU-delin Step 2 : Creating the mask')
grass.run_command('r.null', map=basins_rcl, setnull=0, stdout=DEVNULL, stderr=DEVNULL)
grass_run_command('r.null', map=basins_rcl, setnull=0, stdout=DEVNULL, stderr=DEVNULL)
mask_tmp = 'mask_tmp'
pReclass = grass.feed_command('r.reclass', input=basins_rcl, output=mask_tmp, overwrite='True', rules='-')
pReclass = grass_feed_command('r.reclass', input=basins_rcl, output=mask_tmp, overwrite='True', rules='-')
pReclass.stdin.write(encode('0 thru 10000000 = 1\n'))
pReclass.stdin.close()
pReclass.wait()
# TODO check that (was removed before)
#grass.run_command('r.null', map=mask_tmp, null=0, setnull=1, stdout=DEVNULL, stderr=DEVNULL)
grass.run_command('r.out.gdal', input=mask_tmp, type='UInt16', output=os.path.join(directory_out, 'step2_mask.tif'),
#grass_run_command('r.null', map=mask_tmp, null=0, setnull=1, stdout=DEVNULL, stderr=DEVNULL)
grass_run_command('r.out.gdal', input=mask_tmp, type='UInt16', output=os.path.join(directory_out, 'step2_mask.tif'),
overwrite='True', stdout=DEVNULL, stderr=DEVNULL)
# exit if mask is too small and reloc is disabled
......@@ -567,36 +567,36 @@ def main(parms_file, nbProc, generator=False):
print('---------- HRU-delin Step 2 : Computing subbasinXwatershed raster layer')
dem_recl = os.path.join(directory_out, 'step1_dem_reclass.tif')
grass.run_command('g.proj', flags='c', georef=dem_recl, stdout=DEVNULL, stderr=DEVNULL)
grass_run_command('g.proj', flags='c', georef=dem_recl, stdout=DEVNULL, stderr=DEVNULL)
print('---------- Importing raster \'step1_subbasins.tif\'')
subbasins = os.path.join(directory_out, 'step1_subbasins.tif')
grass.run_command('r.in.gdal', flags='o', input=subbasins, output='subbasins', overwrite='True', stdout=DEVNULL, stderr=DEVNULL)
grass_run_command('r.in.gdal', flags='o', input=subbasins, output='subbasins', overwrite='True', stdout=DEVNULL, stderr=DEVNULL)
print('---------- Importing raster \'step2_watersheds.tif\'')
basins = os.path.join(directory_out, 'step2_watersheds.tif')
grass.run_command('r.in.gdal', flags='o', input=basins, output='watersheds', overwrite='True', stdout=DEVNULL, stderr=DEVNULL)
grass_run_command('r.in.gdal', flags='o', input=basins, output='watersheds', overwrite='True', stdout=DEVNULL, stderr=DEVNULL)
print('---------- Importing raster \'step2_streams_new.tif\'')
reachraster = os.path.join(directory_out, 'step2_streams_new.tif')
grass.run_command('r.in.gdal', flags='o', input=reachraster, output='reachraster', overwrite='True', stdout=DEVNULL, stderr=DEVNULL)
grass_run_command('r.in.gdal', flags='o', input=reachraster, output='reachraster', overwrite='True', stdout=DEVNULL, stderr=DEVNULL)
print('---------- Compute cross product \'watersheds*100000+subbasins\'')
grass.run_command('r.mapcalc', expression='cross1=watersheds*100000+subbasins', overwrite='True', stdout=DEVNULL, stderr=DEVNULL)
grass_run_command('r.mapcalc', expression='cross1=watersheds*100000+subbasins', overwrite='True', stdout=DEVNULL, stderr=DEVNULL)
print('---------- Compute cross product cross1*reachraster')
grass.run_command('r.mapcalc', expression='cross2=if(reachraster!=0, cross1, 0)', overwrite='True', stdout=DEVNULL, stderr=DEVNULL)
grass_run_command('r.mapcalc', expression='cross2=if(reachraster!=0, cross1, 0)', overwrite='True', stdout=DEVNULL, stderr=DEVNULL)
print('---------- Setting nulls in \'cross1*reachraster\' and \'reachraster\'')
grass.run_command('r.null', map='cross2', setnull=0, stdout=DEVNULL, stderr=DEVNULL)
grass.run_command('r.null', map='reachraster', setnull=0, stdout=DEVNULL, stderr=DEVNULL)
grass_run_command('r.null', map='cross2', setnull=0, stdout=DEVNULL, stderr=DEVNULL)
grass_run_command('r.null', map='reachraster', setnull=0, stdout=DEVNULL, stderr=DEVNULL)
print('---------- Saving \'step2_subbasins.tif\' (watersheds*100000+subbasins)')
grass.run_command('r.out.gdal', input='cross1', output=os.path.join(directory_out, 'step2_subbasins.tif'),
grass_run_command('r.out.gdal', input='cross1', output=os.path.join(directory_out, 'step2_subbasins.tif'),
overwrite='True', stdout=DEVNULL, stderr=DEVNULL)
print('---------- Computing links between (watersheds*100000+subbasins) and reach ids')
reach_ids = decode(grass.read_command('r.stats', quiet=True, flags='nN', input='reachraster')).rstrip(os.linesep).split(os.linesep)
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)
......@@ -622,7 +622,7 @@ def main(parms_file, nbProc, generator=False):
# TODO test with a raster we want to pass to //
os.environ['GISRC'] = os.path.join(grassDbPath, 'grassdata', location, '.grassrc')
dem_recl = os.path.join(directory_out, 'step1_dem_reclass.tif')
grass.run_command('g.proj', flags='c', georef=dem_recl, stdout=DEVNULL, stderr=DEVNULL)
grass_run_command('g.proj', flags='c', georef=dem_recl, stdout=DEVNULL, stderr=DEVNULL)
importRastersInEnv(rastersForWorkers, grassDbPath, location)
......@@ -656,28 +656,28 @@ def main(parms_file, nbProc, generator=False):
os.environ['GISRC'] = MAIN_GISRC
print('')
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)
print('---------- Reclassifying (watersheds*100000+subbasins)')
pReclass = grass.feed_command('r.reclass', input='cross1', output='cross1_reclassed', rules='-', overwrite='True')
pReclass = grass_feed_command('r.reclass', input='cross1', output='cross1_reclassed', rules='-', overwrite='True')
for k in range(0, len(reach_ids_cleaned)):
pReclass.stdin.write(encode('%s=%s\n' % (subbasins_cleaned[k], reach_ids_cleaned[k])))
pReclass.stdin.close()
pReclass.wait()
print('---------- Saving \'step2_subbasins.tif\'')
grass.run_command('r.out.gdal', input='cross1_reclassed', output=os.path.join(directory_out, 'step2_subbasins_2.tif'),
grass_run_command('r.out.gdal', input='cross1_reclassed', output=os.path.join(directory_out, 'step2_subbasins_2.tif'),
overwrite='True', stdout=DEVNULL, stderr=DEVNULL)
print('---------- Creating vector layers from raster layers ... ')
grass.run_command('r.to.vect', flags='v', quiet=True, input='subbasins', output='subbasins_vector', type='area', overwrite='True')
grass.run_command('r.to.vect', flags='v', quiet=True, input='watersheds', output='watersheds_vector', type='area', overwrite='True')
grass_run_command('r.to.vect', flags='v', quiet=True, input='subbasins', output='subbasins_vector', type='area', overwrite='True')
grass_run_command('r.to.vect', flags='v', quiet=True, input='watersheds', output='watersheds_vector', type='area', overwrite='True')
## we could clean the vector data in case isolated pixels in edges produce false small areas
## (with no category because it comes from an isolated nodata pixel in watershed.tif)
#rasterWsheds = gdal.Open(basins)
#ulx, xres, xskew, uly, yskew, yres = rasterWsheds.GetGeoTransform()
#pixelArea = abs(xres) * abs(yres)
#grass.run_command('v.clean',
#grass_run_command('v.clean',
# #flags='v',
# quiet=True,
# input='watersheds_vector',
......@@ -686,25 +686,25 @@ def main(parms_file, nbProc, generator=False):
# threshold=pixelArea,
# output='watersheds_vector_clean',
# overwrite='True')
grass.run_command('r.to.vect', flags='v', quiet=True, input='cross1', output='cross1_vector', type='area', overwrite='True')
grass.run_command('r.to.vect', flags='v', quiet=True, input='cross1_reclassed', output='cross1_reclassed_vector', type='area', overwrite='True')
grass_run_command('r.to.vect', flags='v', quiet=True, input='cross1', output='cross1_vector', type='area', overwrite='True')
grass_run_command('r.to.vect', flags='v', quiet=True, input='cross1_reclassed', output='cross1_reclassed_vector', type='area', overwrite='True')
# to make sure there is a projection in exported files (grass78 complains)
grass.run_command('g.proj', flags='c', georef=drain_layer, stdout=DEVNULL, stderr=DEVNULL)
grass.run_command('v.out.ogr',
grass_run_command('g.proj', flags='c', georef=drain_layer, stdout=DEVNULL, stderr=DEVNULL)
grass_run_command('v.out.ogr',
#flags='c',
quiet=True,
input='subbasins_vector', type='area', format='ESRI_Shapefile', output=os.path.join(directory_out, 'step2_step1_subbasins.shp'))
# we avoid c flag to skip areas with no category (result of vectorisation error)
grass.run_command('v.out.ogr',
grass_run_command('v.out.ogr',
#flags='c',
quiet=True,
input='watersheds_vector', type='area', format='ESRI_Shapefile', output=os.path.join(directory_out, 'step2_step2_watersheds.shp'))
grass.run_command('v.out.ogr',
grass_run_command('v.out.ogr',
#flags='c',
quiet=True,
input='cross1_vector', type='area', format='ESRI_Shapefile', output=os.path.join(directory_out, 'step2_step2_subbasins.shp'))
grass.run_command('v.out.ogr',
grass_run_command('v.out.ogr',
#flags='c',
quiet=True,
input='cross1_reclassed_vector', type='area', format='ESRI_Shapefile', output=os.path.join(directory_out, 'step2_step2_subbasins_2.shp'))
......@@ -714,6 +714,7 @@ def main(parms_file, nbProc, generator=False):
if __name__ == '__main__':
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
# check TQDM presence only if we are executed
try:
from tqdm import tqdm
......@@ -744,4 +745,5 @@ if __name__ == '__main__':
pass
else:
from .parallel import buildGrassEnv, buildGrassLocation, exportRasters, importRastersInEnv
from .progressColors import *
\ No newline at end of file
from .progressColors import *
from .utils import grass_run_command, grass_parse_command, grass_feed_command, grass_read_command
\ No newline at end of file
......@@ -660,6 +660,7 @@ def main(parms_file, nbProc, generator=False):
if __name__ == '__main__':
from parallel import buildGrassEnv, buildGrassLocation, exportRasters, importRastersInEnv, exportRastersFromEnv
from progressColors import *
from utils import grass_run_command, grass_parse_command, grass_feed_command
try:
from tqdm import tqdm
......@@ -689,3 +690,4 @@ if __name__ == '__main__':
else:
from .parallel import buildGrassEnv, buildGrassLocation, exportRasters, importRastersInEnv, exportRastersFromEnv
from .progressColors import *
from .utils import grass_run_command, grass_parse_command, grass_feed_command
\ No newline at end of file
......@@ -32,4 +32,8 @@ def grass_parse_command(*args, **kwargs):
def grass_feed_command(*args, **kwargs):
kwargs['startupinfo'] = getSi()
return grass.feed_command(*args, **kwargs)
\ No newline at end of file
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
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