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.

Commit ef4467e9 authored by Rabotin Michael's avatar Rabotin Michael
Browse files

Added irrigation and dams management

parent 232ac43f
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
############################################################################
#
# MODULE: create_config_file.py
# AUTHOR(S): Michael Rabotin
# PURPOSE: Create empty config file
#
#
# COPYRIGHT: (C) 2021 UR RIVERLY - INRAE
#
# This program is free software under the GNU General Public
# License (>=v2). Read the file LICENSE that comes with
# HRU-DELIN for details.
#
#############################################################################
file= open("hrudelin_config.cfg","a")
file.write("# -----------\n")
file.write("# environment\n")
file.write("# -----------\n")
file.write("\n")
file.write("[dir_in]\n")
file.write("dir:\n")
file.write("\n")
file.write("[dem]\n")
file.write("dem:\n")
file.write("\n")
file.write("\n")
file.write("[data]\n")
file.write("hgeo:\n")
file.write("landuse:\n")
file.write("soil:\n")
file.write("\n")
file.write("\n")
file.write("[gauges]\n")
file.write("gauges:\n")
file.write("#for watershed ID, used for identification of watersheds\n")
file.write("gauges_col_name=\n")
file.write("# drained surface\n")
file.write("gauges_area_col_name=\n")
file.write("relocated_gauges=\n")
file.write("\n")
file.write("\n")
file.write("[irrigation]\n")
file.write("to_do:\n")
file.write("irrigation:\n")
file.write("irrig_col_name=\n")
file.write("#for irrig_col_type, 2 for groundwater and 3 for surface")
file.write("irrig_col_type=\n")
file.write("\n")
file.write("# you can indicate a minimum surface value for an HRU to be a GU:\n")
file.write("#minimum surface can be null, global (irrig_surf_min_GU) or spatialized (irrig_col_min_GU)\n")
file.write("irrig_surf_min_GU=\n")
file.write("irrig_col_min_GU=\n")
file.write("#you can specify a maximum distance search for GU (default is 5000 m\n")
file.write("irrig_distance_GU=5000\n")
file.write("\n")
file.write("irrigation_sector:\n")
file.write("irrig_sector_col_name=\n")
file.write("irrig_col_sau_irr=\n")
file.write("irrig_col_dom_sau_irr=\n")
file.write("irrigation_table:\n")
file.write("relocated_irrigation=\n")
file.write("\n")
file.write("[dams]\n")
file.write("to_do:\n")
file.write("dams=\n")
file.write("dams_col_name=\n")
file.write("dams_smax=\n")
file.write("dams_s0=\n")
file.write("#drained surface\n")
file.write("dams_area_col_name=\n")
file.write("relocated_dams=\n")
file.write("\n")
file.write("[dir_out]\n")
file.write("files:\n")
file.write("results:\n")
file.write("# -------------------------\n")
file.write("# 1st step : hru-delin_init\n")
file.write("# -------------------------\n")
file.write("\n")
file.write("[surface]\n")
file.write("#selection: total -> full dem\n")
file.write("# polygon -> polygon: name of the shapefile\n")
file.write("# coords -> give the coords upper left (west and north) and lower right (east and south)\n")
file.write("selection:\n")
file.write("polygon:\n")
file.write("west:\n")
file.write("north:\n")
file.write("east:\n")
file.write("south:\n")
file.write("\n")
file.write("\n")
file.write("[demfill]\n")
file.write("#\n")
file.write("# if demfill = yes : depressionless DEM will be generated\n")
file.write("# no : no action on input DEM \n")
file.write("#\n")
file.write("demfill:\n")
file.write("\n")
file.write("#\n")
file.write("# if rules_auto_* = yes : rules will be calculated by the module\n")
file.write("# if no : fill the corresponding file (reclass_default_rules_*)\n")
file.write("#\n")
file.write("[reclass_dem]\n")
file.write("rules_auto_dem:\n")
file.write("step_dem:\n")
file.write("\n")
file.write("[reclass_slope]\n")
file.write("rules_auto_slope:\n")
file.write("\n")
file.write("[reclass_aspect]\n")
file.write("rules_auto_aspect:\n")
file.write("\n")
file.write("[basin_min_size]\n")
file.write("# minimum size of calculated watersheds (r.watershed)\n")
file.write("# number of pixels \n")
file.write("# size = N = SURFACE_km2 / ( RES_km2^2 )\n")
file.write("# ex: S = 10km2, RES = 200m = 0.2 km ==> N = 250 pixels\n")
file.write("# S = 20km2, RES = 90m = 0.09 km ==> N = 2469 pixels\n")
file.write("size=\n")
file.write("\n")
file.write("# ---------------------------\n")
file.write("# 2nd step : hru-delin_basins\n")
file.write("# ---------------------------\n")
file.write("# So it's possible to specify a variable using : or = ???\n")
file.write("[auto_relocation]\n")
file.write("to_do:\n")
file.write("# -------- first rule\n")
file.write("# surface is in percent!\n")
file.write("# distance is in pixels: N = D / RES\n")
file.write("# example: for 3km distance with a 50m DEM, the number of pixel is: 3000/50 = 33 pixels\n")
file.write("surface_tolerance_1=\n")
file.write("distance_tolerance_1=\n")
file.write("# -------- second rule \n")
file.write("# second rule with a distance tolerance of 6km and a surface tolerance of 30%\n")
file.write("surface_tolerance_2=\n")
file.write("distance_tolerance_2=120\n")
file.write("\n")
file.write("# unit = 1 : m , = 2 : km\n")
file.write("area_unit=\n")
file.write("\n")
file.write("\n")
file.write("\n")
file.write("# ---------------------------\n")
file.write("# 3rd step : hru-delin_hrugen\n")
file.write("# ---------------------------\n")
file.write("\n")
file.write("[hrus_min_surface]\n")
file.write("# there, this is in pixel so pay attention to the DEM resolution \n")
file.write("# same as for bassin_min_size: N = SURFACE_km2 / ( RES_km2^2 ) # see 'size' parameter in step 1 for examples\n")
file.write("# 2 km2 = 247 pixels\n")
file.write("surface=\n")
file.write("\n")
file.write("#\n")
file.write("# MNT-derived layers to be integrated in the overlay operation\n")
file.write("#\n")
file.write("[layer_overlay]\n")
file.write("dem:\n")
file.write("slope:\n")
file.write("aspect:\n")
file.write("\n")
file.write("# --------------------------------\n")
file.write("# 4th step : hru-delin_parms_J2000\n")
file.write("# --------------------------------\n")
file.write("[topology]\n")
file.write("dissolve_cycle:\n")
file.close()
\ No newline at end of file
......@@ -28,6 +28,7 @@
from __future__ import print_function
import string, os, shutil, sys, glob, types
import math
import json
try:
import ConfigParser
except Exception as e:
......@@ -101,6 +102,69 @@ def get_coords(parms, xmin_slect, ymax_slect, xmax_slect, ymin_slect):
return xmin_slect, ymax_slect, xmax_slect, ymin_slect
def is_valid_shp(shplayer):
if not os.path.isfile(shplayer):
msg="------------> ERROR : Input "+ shplayer+" not found"
sys.exit(msg)
if ogr.Open(shplayer) is None:
msg="------------> ERROR : Input "+ shplayer+" is not a valid shapefile"
sys.exit(msg)
def is_valid_geometry(Name,GeometryType):
NameDs=ogr.Open(Name)
NameLayer=NameDs.GetLayer()
NameLayer_defn=NameLayer.GetLayerDefn()
if ogr.GeometryTypeToName(NameLayer_defn.GetGeomType()) != GeometryType:
msg="------------> ERROR : Input "+ Name+" is not "+ GeometryType
sys.exit(msg)
def is_column_exist(Name,colname):
NameDs=ogr.Open(Name)
NameLayer=NameDs.GetLayer()
NameLayer_defn=NameLayer.GetLayerDefn()
FieldFound=False
for n in range(NameLayer_defn.GetFieldCount()):
fdefn=NameLayer_defn.GetFieldDefn(n)
if fdefn.name == colname:
FieldFound=True
if not FieldFound:
msg="------------> ERROR : column "+ colname+ " not found in Input "+Name
sys.exit(msg)
def is_column_valid(Name,colname):
NameDs=ogr.Open(Name)
NameLayer=NameDs.GetLayer()
listvalue= [ feature.GetField(colname) for feature in NameLayer ]
if None in listvalue or (min(listvalue) == 0):
msg="ERROR : in column "+ colname+ " in Input "+Name +" null or zero value detected"
sys.exit(msg)
def is_column_value_unique(Name,colname):
NameDs=ogr.Open(Name)
NameLayer=NameDs.GetLayer()
listvalue= [ feature.GetField(colname) for feature in NameLayer ]
if(len(set(listvalue)) != len(listvalue)):
msg="ERROR : in column "+ colname+ " in Input "+Name +" non unique value detected"
sys.exit(msg)
def are_columns_value_unique(Name1,colname1,Name2,colname2):
NameDs1=ogr.Open(Name1)
NameLayer1=NameDs1.GetLayer()
listvalue1= [ feature.GetField(colname1) for feature in NameLayer1 ]
NameDs2=ogr.Open(Name2)
NameLayer2=NameDs2.GetLayer()
listvalue2= [ feature.GetField(colname2) for feature in NameLayer2 ]
listvalue1.extend(listvalue2)
if(len(set(listvalue1)) != len(listvalue1)):
msg="ERROR : in column "+ colname1+ " in Input "+Name1 +" and in column "+colname2+" in Input "+Name2+" duplicate value !"
sys.exit(msg)
......@@ -121,10 +185,11 @@ except ImportError:
def main(parms_file):
print('---------- HRU-delin Step 1 started ---------------------------------------------')
#test gdal version; must be at least >=3
if int(gdal.VersionInfo('VERSION_NUM'))< 3000000:
sys.exit('------------> ERROR : GDAL version mut be at least >= 3.0')
configFileDir = os.path.dirname(parms_file)
# create main env
buildGrassEnv(os.path.join(configFileDir, 'grass_db'), 'hru-delin')
......@@ -150,57 +215,197 @@ def main(parms_file):
directory_out = os.path.join(configFileDir, directory_out)
## Test if dem exist
dem = os.path.join(directory, str(parms.get('files_in', 'dem')))
dem = os.path.join(directory, str(parms.get('dem', 'dem')))
if not os.path.isfile(dem):
sys.exit('------------> ERROR : Input Dem not found')
## Test if gauges exist
gauges = os.path.join(directory, str(parms.get('files_in', 'gauges')))
if not os.path.isfile(gauges):
sys.exit('------------> ERROR : Input Gauges not found')
msg="------------> ERROR : Input Dem "+ dem+" not found"
sys.exit(msg)
## if hgeon, landuse and soil provided, test if exist
for data in parms.items('data'):
data = os.path.join(directory, data[1])
if not os.path.isfile(data):
print(data)
sys.exit('------------> ERROR : Input data not found' )
## Test if gauges exist
gauges = os.path.join(directory, str(parms.get('gauges', 'gauges')))
is_valid_shp(gauges)
## Test if gauges is point or multipoint
GaugesDs=ogr.Open(gauges)
GaugesLayer=GaugesDs.GetLayer()
GaugesLayer_defn=GaugesLayer.GetLayerDefn()
if ogr.GeometryTypeToName(GaugesLayer_defn.GetGeomType()) != 'Point':
sys.exit('------------> ERROR : Input Gauges is not Point Geometry (multi Point not allowed)')
is_valid_geometry(gauges,'Point')
## Test gauges_col_name
gauges_col_name = parms.get('gauges', 'gauges_col_name')
if gauges_col_name == '':
msg="------------> ERROR : CHECK gauges_col_name PARAMETER"
sys.exit(msg)
## In gauges, test if gauges_col_name exist and if no null value
is_column_exist(gauges,gauges_col_name)
is_column_valid(gauges,gauges_col_name)
is_column_value_unique(gauges,gauges_col_name)
## Test gauges_area_col_name
gauges_area_col_name = parms.get('gauges', 'gauges_area_col_name')
if gauges_area_col_name == '':
msg="------------> ERROR : CHECK gauges_area_col_name PARAMETER"
sys.exit(msg)
## In gauges, test if gauges_area_col_name exist and if no null value
is_column_exist(gauges,gauges_area_col_name)
is_column_valid(gauges,gauges_area_col_name)
gauges_file = parms.get('gauges', 'relocated_gauges')
if gauges_file != '':
is_valid_shp(gauges_file)
## Test col_name
col_name = parms.get('for_watershed_id', 'col_name')
if col_name == '':
sys.exit('------------> ERROR : CHECK col_name PARAMETER')
## if irrigation is yes, test input data
if str(parms.get('irrigation', 'to_do')) == 'yes':
irrig = os.path.join(directory, str(parms.get('irrigation', 'irrigation')))
is_valid_shp(irrig)
## Test if irrigation is point or multipoint
is_valid_geometry(irrig,'Point')
## Test irrig_col_name
irrig_col_name = parms.get('irrigation', 'irrig_col_name')
if irrig_col_name == '':
sys.exit('------------> ERROR : CHECK irrig_col_name PARAMETER in Irrigation')
## In irrigation, test if irrig_col_name exists and if no null value
is_column_exist(irrig,irrig_col_name)
is_column_valid(irrig,irrig_col_name)
is_column_value_unique(irrig,irrig_col_name)
## Test irrig_col_type
irrig_col_type = parms.get('irrigation', 'irrig_col_type')
if irrig_col_type == '':
sys.exit('------------> ERROR : CHECK irrig_col_type PARAMETER in Irrigation')
## In irrigation, test if irrig_col_type exists
is_column_exist(irrig,irrig_col_type)
#if min_surf_GU is provided, test value:
irrig_surf_min_GU = parms.get('irrigation', 'irrig_surf_min_GU')
if irrig_surf_min_GU != '':
try:
isinstance(float(irrig_surf_min_GU),float)
except:
sys.exit('------------> ERROR : irrig_surf_min_GU PARAMETER is not integer or double')
irrig_col_min_GU = parms.get('irrigation', 'irrig_col_min_GU')
if irrig_col_min_GU != '':
is_column_exist(irrig,irrig_col_min_GU)
irrig_distance_GU = parms.get('irrigation', 'irrig_distance_GU')
try:
isinstance(float(irrig_distance_GU),float)
except:
sys.exit('------------> ERROR : irrig_distance_GU PARAMETER is not integer or double')
if irrig_distance_GU<0:
sys.exit('------------> ERROR : irrig_distance_GU PARAMETER cannot be negative')
irrig_sector = os.path.join(directory, str(parms.get('irrigation', 'irrigation_sector')))
is_valid_shp(irrig_sector)
## Test if irrig_sector is polyon or multipolygon
is_valid_geometry(irrig_sector,'Polygon')
## Test irrig_sector_col_name
irrig_sector_col_name = parms.get('irrigation', 'irrig_sector_col_name')
if irrig_sector_col_name == '':
sys.exit('------------> ERROR : CHECK irrig_sector_col_name PARAMETER in Irrigation')
## In irrig_sector, test if irrig_sector_col_name exists and if no null value
is_column_exist(irrig_sector,irrig_sector_col_name)
is_column_valid(irrig_sector,irrig_sector_col_name)
## Test irrig_sector_col_SAU_IRR
irrig_sector_col_SAU = parms.get('irrigation', 'irrig_sector_col_sau_irr')
if irrig_sector_col_SAU == '':
sys.exit('------------> ERROR : CHECK irrig_sector_col_sau_irr PARAMETER in Irrigation')
## In irrig_sector, test if irrig_sector_col_SAU_IRR exists
is_column_exist(irrig_sector,irrig_sector_col_SAU)
## Test irrig_sector_col_DOM_SAU_IRR
irrig_sector_col_DOM_SAU = parms.get('irrigation', 'irrig_sector_col_dom_sau_irr')
if irrig_sector_col_DOM_SAU == '':
sys.exit('------------> ERROR : CHECK irrig_sector_col_dom_sau_irr PARAMETER in Irrigation')
## In irrig_sector, test if irrig_sector_col_DOM_SAU_IRR exists
is_column_exist(irrig_sector,irrig_sector_col_DOM_SAU)
## In gauges, test if col_name exist and if no null value
NoColName=False
for n in range(GaugesLayer_defn.GetFieldCount()):
fdefn=GaugesLayer_defn.GetFieldDefn(n)
if fdefn.name == col_name:
NoColName=True
if NoColName == False:
sys.exit('------------> ERROR : col_name PARAMETER not found in gauges input file')
irrigation_file = parms.get('irrigation', 'relocated_irrigation')
if irrigation_file != '':
is_valid_shp(irrigation_file)
listvalue= [ feature.GetField(col_name) for feature in GaugesLayer ]
if None in listvalue:
sys.exit('------------> ERROR : in col_name in gauges file, null value detected')
if min(listvalue) == 0:
sys.exit('------------> ERROR : in col_name in gauges file, zero value detected')
## if dams is yes, test input data
if str(parms.get('dams', 'to_do')) == 'yes':
dams = os.path.join(directory, str(parms.get('dams', 'dams')))
is_valid_shp(dams)
## if hgeon, landuse and soil provided, test if exist
for data in parms.items('data'):
data = os.path.join(directory, data[1])
if not os.path.isfile(data):
print(data)
sys.exit('------------> ERROR : Input data not found' )
## Test if dams is point or multipoint
is_valid_geometry(dams,'Point')
## Test dams_col_name
dams_col_name = parms.get('dams', 'dams_col_name')
if dams_col_name == '':
sys.exit('------------> ERROR : CHECK dams_col_name PARAMETER in Irrigation')
## In dams, test if dams_col_name exists and if no null value
is_column_exist(dams,dams_col_name)
is_column_valid(dams,dams_col_name)
is_column_value_unique(dams,dams_col_name)
are_columns_value_unique(gauges,gauges_col_name,dams,dams_col_name)
## Test dams_smax
dams_smax = parms.get('dams', 'dams_smax')
if dams_smax == '':
sys.exit('------------> ERROR : CHECK dams_smax PARAMETER in dams')
## In dams, test if dams_smax exists
is_column_exist(dams,dams_smax)
## Test dams_s0
dams_s0 = parms.get('dams', 'dams_s0')
if dams_s0 == '':
sys.exit('------------> ERROR : CHECK dams_s0 PARAMETER in dams')
## In dams, test if dams_s0 exists
is_column_exist(dams,dams_s0)
# Test dams_area_col_name
dams_area_col_name = parms.get('dams', 'dams_area_col_name')
if dams_area_col_name == '':
sys.exit('------------> ERROR : CHECK dams_area_col_name PARAMETER in dams')
## In dams, test if dams_area_col_name exists
is_column_exist(dams,dams_area_col_name)
dams_file = parms.get('dams', 'relocated_dams')
if dams_file != '':
is_valid_shp(dams_file)
#test if dams_col_name and gauges_col_name have no common value
## if irrig_rast provided, test if is valid
if str(parms.get('irrigation', 'irrig_rast')) != '':
irr_rast_test = os.path.join(directory, str(parms.get('irrigation', 'irrig_rast')))
if not os.path.isfile(irr_rast_test):
sys.exit('------------> ERROR : Irrigation raster not found' )
## if rules_auto_dem is yes, test if step_dem valid
if (parms.get('reclass_dem', 'rules_auto_dem')) == 'yes':
......@@ -217,7 +422,7 @@ def main(parms_file):
dem_name = os.path.basename(parms.get('files_in', 'dem')).split('.')[0]
dem_name = os.path.basename(parms.get('dem', 'dem')).split('.')[0]
selection = (parms.get('surface', 'selection'))
xmin_slect, ymax_slect, xmax_slect, ymin_slect = get_raster_bounds(dem, xmin_slect, ymax_slect, xmax_slect, ymin_slect)
......@@ -244,18 +449,30 @@ def main(parms_file):
dem_cut = os.path.join(directory_out, 'step1_dem_cut.tif')
cutting_raster(dem, dem_cut, xmin_slect, ymax_slect, xmax_slect, ymin_slect)
irr_rast = str(parms.get('irrigation', 'irrig_rast'))
if irr_rast != '':
irr_rast = os.path.join(directory, irr_rast)
layer_cut = os.path.join(directory_out, 'step1_irrigation.tif')
cutting_raster(irr_rast, layer_cut, xmin_slect, ymax_slect, xmax_slect, ymin_slect)
# Cutting up of the gauges shapefile
gauges_shp = os.path.join(directory, parms.get('files_in', 'gauges'))
gauges_shp = os.path.join(directory, parms.get('gauges', 'gauges'))
gauges_selected = os.path.join(directory_out, 'gauges_selected.shp')
## this ensures intersection with DEM
ogr2ogrMain(['', '-spat', str(xmin_slect), str(ymin_slect), str(xmax_slect), str(ymax_slect), gauges_selected, gauges_shp])
#if irrigation, cutting irrigation and cantons vector
if str(parms.get('irrigation', 'to_do')) == 'yes':
irrigation_shp = os.path.join(directory, parms.get('irrigation', 'irrigation'))
irrigation_selected = os.path.join(directory_out, 'irrigation_selected.shp')
## this ensures intersection with DEM
ogr2ogrMain(['', '-spat', str(xmin_slect), str(ymin_slect), str(xmax_slect), str(ymax_slect), irrigation_selected, irrigation_shp])
cantons_shp = os.path.join(directory, parms.get('irrigation', 'irrigation_sector'))
cantons_selected = os.path.join(directory_out, 'irrig_sector_selected.shp')
## this ensures intersection with DEM
ogr2ogrMain(['', '-spat', str(xmin_slect), str(ymin_slect), str(xmax_slect), str(ymax_slect), cantons_selected, cantons_shp])
if str(parms.get('dams', 'to_do')) == 'yes':
dams_shp = os.path.join(directory, parms.get('dams', 'dams'))
dams_selected = os.path.join(directory_out, 'dams_selected.shp')
## this ensures intersection with DEM
ogr2ogrMain(['', '-spat', str(xmin_slect), str(ymin_slect), str(xmax_slect), str(ymax_slect), dams_selected, dams_shp])
# Cutting up and save other rasters (landuse, soils, geology, ...)
for data in parms.items('data'):
xmin_slect, ymax_slect, xmax_slect, ymin_slect = get_raster_bounds(os.path.join(directory, data[1]), xmin_slect, ymax_slect, xmax_slect, ymin_slect)
......
This diff is collapsed.
......@@ -42,6 +42,7 @@ from osgeo import gdal
from osgeo.gdalconst import *
from osgeo import ogr
import pandas as pd
import json
import multiprocessing
from multiprocessing import Pool, cpu_count
......@@ -225,8 +226,7 @@ def main(parms_file, nbProc, generator=False):
parms.read(parms_file)
directory_out = parms.get('dir_out', 'files')
if not os.path.isdir(directory_out):
sys.exit('------------> ERROR : Out Directoy is not valid')
## if hgeon, landuse and soil provided, test if exist
......@@ -238,10 +238,10 @@ def main(parms_file, nbProc, generator=False):
sys.exit('------------> ERROR : Input data not found' )
## if irrig_rast provided, test if is valid
if str(parms.get('irrigation', 'irrig_rast')) != '':
irr_rast_test = os.path.join(directory, str(parms.get('irrigation', 'irrig_rast')))
if not os.path.isfile(irr_rast_test):
sys.exit('------------> ERROR : Irrigation raster not found' )
#if str(parms.get('irrigation', 'irrig_rast')) != '':
# irr_rast_test = os.path.join(directory, str(parms.get('irrigation', 'irrig_rast')))
# if not os.path.isfile(irr_rast_test):
# sys.exit('------------> ERROR : Irrigation raster not found' )
# manage absolute and relative paths
......@@ -252,6 +252,8 @@ def main(parms_file, nbProc, generator=False):
if not os.path.isabs(dir_results):
dir_results = os.path.join(configFileDir, dir_results)