#!/bin/sh ############################################################################ # # MODULE: r.inund.fluv v.1 for GRASS 6.2 # AUTHOR(S): Roberto Marzocchi, Bianca Federici, Domenico Sguerso # PURPOSE: Creates a fluvial inundation map given an high-resolution dtm and a water surface profile # COPYRIGHT: (C) 2008 by the GRASS Development Team # # This program is free software under the GNU General Public # License (>=v2). Read the file COPYING that comes with GRASS # for details. # ############################################################################# #%Module #% description: Creates a fluvial inundation map given an high-resolution dtm and a water surface profile #% keywords: Automatic procedure to compute a fluvial inundation map #%End #%option #% key: DTM #% type: string #% gisprompt: old,cell,raster #% description: Input raster map (DTM or DEM) #% required : yes #%end #%option #% key: W_S_PROFILE #% type: string #% gisprompt: old_file,file,input #% description: Input ASCII file of the water surface profile #% required : yes #%end #%option #% key: RIVER #% type: string #% gisprompt: old,vector,vector #% description: Input vector map of river-axis #% required : yes #%end #%option #% key: FLOODING_MAP #% type: string #% gisprompt: new,cell,raster #% description: Output: name of flooding map #% required : yes #%end #%option #% key: DOUBT_MAP #% type: string #% gisprompt: new,cell,raster #% description: Output: name of doubful surface areas. #% required : no #%end #%option #% key: PROFILE_T100 #% type: string #% gisprompt: old_file,file,input #% description: Input ASCII file with water-depht for return period T > 100 years #% required : no #%end #%option #% key: delta_y #% type: double #% gisprompt: new #% description: Input delta_y to find the boundaries of the main channel [default value 0.5 m] #% required : no #%end #%option #% key: delta_x #% type: double #% gisprompt: new #% description: Input delta_x to find the boundaries of the main channel [default value DTM_resolution * 1.5 m ] #% required : no #%end #%option #% key: river_boundary #% type: string #% gisprompt: new,vector,vector #% description: Output vector boundaries of the main channel #% required : no #%end #%option #% key: boundary_type #% type: string #% description: Format of vector boundaries of the main channel #% options: line,points #% answer: points #% required: no #%end #%option #% key: res_B #% type: integer #% gisprompt: new,cell,raster #% description: Input value: resolution B [default value 10 m] #% required : no #%end #%option #% key: res_C #% type: integer #% gisprompt: new,cell,raster #% description: Input value: resolution C [default value 20 m] #% required : no #%end #%option #% key: MAP_W_S_PROFILE #% type: string #% gisprompt: new,vector,vector #% description: Output vector map of the water surface profile #% required : no #%end if [ -z "$GISBASE" ] then echo "" echo "You must be in GRASS GIS to run this program" echo "" exit 1 fi if [ "$1" != "@ARGS_PARSED@" ] ; then exec g.parser "$0" "$@" fi $GRASS_VERBOSE #request control & test (only for imput map) if [ -z "$GIS_OPT_DTM" ] ; then echo "ERROR: DTM not specified" exit 1 fi #g.findfile element=cell file="$GIS_OPT_DTM" > /dev/null # test if a map exists #if [ $? -eq 0 ]; then # echo "Imput map not found" 1>&2 # exit 1 #fi g.region save=dtm_res2 --overwrite g.region rast="$GIS_OPT_DTM" save=verifica g.region -p > .resolution.txt g.region region=dtm_res2 exec 6<&0 exec < .resolution.txt read a1 b1 read a2 b2 read a3 b3 read a4 b4 read a5 b5 read a6 b6 read a7 b7 read a8 b8 read a9 b9 exec 0<&6 6<&- res=$b9 if [ "$b9" -gt 5 ] then echo "WARNING - the DTM resolution it'snt good for this application" else echo "The DTM resolution is $res x $res meters" fi rm .resolution.txt if [ -z "$GIS_OPT_W_S_PROFILE" ] ; then echo "ERROR: file with depht not specified" exit 1 fi if [ -z "$GIS_OPT_RIVER" ] ; then echo "ERROR: vector with river-axis not specified" exit 1 fi if [ -z "$GIS_OPT_FLOODING_MAP" ] ; then echo "ERROR: name of flooding map not specified" exit 1 fi #optional control: assignement default value if [ -n "$GIS_OPT_res_B" ] then res10="$GIS_OPT_res_B" else res10=10 fi if [ -n "$GIS_OPT_res_C" ] then res20="$GIS_OPT_res_C" else res20=20 fi if [ -n "$GIS_PROFILE_T100" ] then quote100="$GIS_OPT_PROFILE_T100" else quote100="$GIS_OPT_W_S_PROFILE" fi if [ -n "$GIS_OPT_delta_x" ] then deltax="$GIS_OPT_delta_x" else deltax=$(echo "$res * 1.5" | bc) fi if [ -n "$GIS_OPT_delta_y" ] then deltay="$GIS_OPT_delta_y" else deltay="0.5" fi g.remove rast=MASK #make a temporary directory TMPDIR="`g.tempfile pid=$$`" if [ $? -ne 0 ] || [ -z "$TMPDIR" ] ; then g.message "Unable to create temporary files" exit 1 fi rm -r -f "$TMPDIR" mkdir "$TMPDIR" cp "$GISBASE/etc/fortran_code"/find_main_channel "$TMPDIR" cp "$GISBASE/etc/fortran_code"/clean_inundation "$TMPDIR" cp "$GISBASE/etc/fortran_code"/2d_path "$TMPDIR" cp "$GISBASE/etc/fortran_code"/correction_from_path "$TMPDIR" #variables dtm="$GIS_OPT_DTM" cp "$GIS_OPT_W_S_PROFILE" "$TMPDIR"/profilo v.in.ascii -z input="$GIS_OPT_W_S_PROFILE" output=profilo format=point fs=" " skip=0 x=1 y=2 z=3 cat=0 --overwrite centerline="$GIS_OPT_RIVER" inondazione="$GIS_OPT_FLOODING_MAP" #calculation of the "boundaries" of the main channel g.region nsres=$res ewres=$res save=dtm_res2 --overwrite g.region -p >> "$TMPDIR"/"region2.txt" r.out.ascii -h input=$dtm output="$TMPDIR"/"dtm2x2" null=0 #fortran code run cd cp "$quote100" "$TMPDIR"/quote_100 cd "$TMPDIR" echo $res $deltax $deltay > parameter.txt echo "'"region2.txt"'" "'"dtm2x2"'" "'"limiti_alveo_vect"'" "'"quote_100"'" > nomefile.txt echo "Fortran code that automatically find bed river" ./find_main_channel rm nomefile.txt rm parameter.txt #TEST FORTRAN CODE RUN if [ -z "limiti_alveo_vect" ] then echo "ERROR IN FORTRAN CODE - see fortran error message and try understand the problem" exit 1 fi cd v.in.ascii --o input="$TMPDIR"/limiti_alveo_vect output=limiti_alveo format=point fs=' ' skip=0 x=1 y=2 cat=0 #export of 20x20 dtm, region20.txt and region10.txt g.region nsres=$res20 ewres=$res20 save=dtm_res20 --overwrite g.region -p >> "$TMPDIR"/region20.txt r.out.ascii -h input=$dtm output="$TMPDIR"/dtm20x20 null=0 g.region nsres=$res10 ewres=$res10 save=dtm_res10 --overwrite g.region -p >> "$TMPDIR"/region10.txt echo "STEP 1" ################################################################################################################### # STEP 1 - Thiessen interpolation ################################################################################################################### # don't use low resolution [default 10m] r.mask input=$dtm maskcats=* g.region region=dtm_res10 v.surf.idw input=profilo output=quote_punti_adiacenti npoints=1 layer=1 column=dbl_3 g.region region=dtm_res2 r.mapcalculator amap=quote_punti_adiacenti bmap=$dtm formula="if(A>=B)" outfile=inondazione_step1 help=- echo "STEP 2" ################################################################################################################### # STEP 2 - remove "lakes" ################################################################################################################### r.to.vect -s input=inondazione_step1 output=inondazione_step1 feature=area --overwrite v.select ainput=inondazione_step1 atype=line,area alayer=1 binput=$centerline btype=line blayer=1 output=inondazione_step2 operator=overlap --overwrite v.to.rast input=inondazione_step2 output=inondazione_step2 use=val layer=1 value=1 rows=4096 --overwrite echo "STEP 3" ################################################################################################################### # STEP 3 - looking at the minimum path from wet pixel to the channel axis # if the dtm elevation is higher than the water elevation in the pixel ==> pixel becomes dry # if the dtm elevation is always lower than the water elevation in the pixel ==> pixel remains wet ################################################################################################################### # export data #20x20 g.region region=dtm_res20 r.out.ascii -h input=inondazione_step2 output="$TMPDIR"/inondazione_step2 null=0 r.out.ascii -h input=quote_punti_adiacenti output="$TMPDIR"/quote_punti_adiacenti dp=2 null=0 # FORTRAN code echo "CODE FORTRAN OF STEP 3" cd cd "$TMPDIR" echo "'"inondazione_step2"'" "'"profilo"'" "'"inondazione_nuova.txt"'" "'"quote_punti_adiacenti"'" "'"region20.txt"'" "'"region2.txt"'" "'"dtm2x2"'" > nomefile.txt ./clean_inundation rm nomefile.txt #TEST FORTRAN CODE RUN if [ -z "inondazione_nuova.txt" ] then echo "ERROR IN FORTRAN CODE - see fortran error message and try understand the problem" exit 1 fi cd # output file r.in.ascii input="$TMPDIR"/inondazione_nuova.txt output=inondazione_step3 'mult=1.0 or read from header' nv=0 --overwrite echo "STEP 4" ################################################################################################################### # STEP 4 - calcolation of the 2d_path outside the main channel ################################################################################################################### r.out.ascii -h input=inondazione_step3 output="$TMPDIR"/inondazione_step3 null=0 echo "2D MODEL OF DIFFUSION OF FLOODING INONDATION OUTSIDE THE MAIN CHANEL" cd cd "$TMPDIR" echo "'"profilo"'" "'"inondazione_step3"'" "'"reticolo"'" "'"region20.txt"'" "'"dtm20x20"'" "'"limiti_alveo_vect"'" "'"region2.txt"'" "'"dtm2x2"'" "'"bordi_alveo"'" > nomefile.txt ./2d_path rm nomefile.txt #TEST FORTRAN CODE RUN if [ -z "reticolo" ] then echo "ERROR IN FORTRAN CODE - see fortran error message and try understand the problem" exit 1 fi cd #r.in.ascii input="$TMPDIR"/bordi_alveo output=bordi_alveo 'mult=1.0 or read from header' nv=0 --overwrite #r.thin input=bordi_alveo output=bordi_alveo_thin iterations=200 --overwrite #r.to.vect input=bordi_alveo_thin output=bordi_alveo feature=line --overwrite v.in.ascii -n input="$TMPDIR"/bordi_alveo output=bordi_alveo format=standard fs=' ' skip=0 x=1 y=2 z=0 cat=0 r.in.ascii input="$TMPDIR"/reticolo output=reticolo 'mult=1.0 or read from header' nv=* --overwrite r.thin input=reticolo output=reticolo_thin iterations=200 --overwrite r.to.vect input=reticolo_thin output=reticolo feature=line --overwrite g.remove rast=reticolo_thin # 10x10 resolution ---> 1 pò più lento ma più preciso! g.region region=dtm_res10 # difference between maps of step1 and step2 r.mapcalculator amap=inondazione_step3 formula='isnull(A)' outfile=inondazione_step3_isnull help=- r.mapcalculator amap=inondazione_step3_isnull formula='if(A,0,1)' outfile=inondazione_step3 help=- r.mapcalculator amap=inondazione_step2 bmap=inondazione_step3 formula='A-B' outfile=diff_inond help=- r.mapcalculator amap=diff_inond formula='if(A,1,0,0)' outfile=diff_inond help=- #remove the -1 value resulted from A-B subtraction r.null map=diff_inond setnull=0 r.to.vect input=diff_inond output=diff_inond1 feature=area --overwrite #to use the v.select command the line must have a layer associated ==> export of vector map in dxf format and then import the dxf file v.out.dxf input=reticolo output="$TMPDIR"/reticolo.dxf v.in.dxf -1 input="$TMPDIR"/reticolo.dxf output=reticolo --overwrite #only 2d_path overlapped the surfaces to look v.select ainput=reticolo atype=line alayer=1 binput=diff_inond1 btype=area blayer=1 output=reticolo_interno operator=overlap --overwrite # creation of points by 2d path outside the main channel v.to.points input=reticolo_interno type=point,line,boundary,centroid output=punti_reticolo_interno llayer=1 dmax=30 --overwrite # convert 2d points to 3d points with the elevation of the water surface profile at the starting point of the path in the channel axis (or with the elevation of the nearest water surface profile) g.region region=dtm_res20 v.drape input=punti_reticolo_interno type=point,centroid,line,boundary,face,kernel rast=reticolo method=nearest output=punti3d_reticolo --overwrite #v.drape input=punti_reticolo_interno type=point,centroid,line,boundary,face,kernel rast=quote_punti_adiacenti method=nearest output=punti3d_reticolo --overwrite g.region region=dtm_res10 #export of 3d points v.out.ascii input=punti3d_reticolo output="$TMPDIR"/punti3d_reticolo format=point dp=2 r.out.ascii -h input=diff_inond output="$TMPDIR"/diff_inond null=0 echo "FORTRAN CODE OF STEP 4" cd cd "$TMPDIR" echo "'"diff_inond"'" "'"punti3d_reticolo"'" "'"correzione"'" "'"region10.txt"'" "'"dtm2x2"'" "'"region2.txt"'" > nomefile.txt ./correction_from_path rm nomefile.txt #TEST FORTRAN CODE RUN if [ -z "correzione" ] then echo "ERROR IN FORTRAN CODE - see fortran error message and try understand the problem" exit 1 fi cd r.in.ascii input="$TMPDIR"/correzione output=inondazione_step5 'mult=1.0 or read from header' nv=0 --overwrite # replace null cells with 0 cells g.region region=dtm_res2 r.mapcalculator amap=inondazione_step5 formula='isnull(A)' outfile=inondazione_step5_isnull help=- r.mapcalculator amap=inondazione_step5_isnull formula='if(A,0,1)' outfile=inondazione_step5 help=- r.mapcalculator amap=inondazione_step5 bmap=inondazione_step3 formula='A+B' outfile=inondazione help=- echo "STEP 5" ################################################################################################################### # STEP 5 - EXPORT OUTPUT FILES ################################################################################################################### g.rename rast=inondazione,inondazione1 --overwrite r.null map=inondazione1 setnull=0 --overwrite r.to.vect -s input=inondazione1 output=inondazione1 feature=area --overwrite v.patch input=centerline,punti3d_reticolo output=overlap --overwrite v.select ainput=inondazione1 atype=area alayer=1 binput=overlap btype=line blayer=1 output=inondazione operator=overlap --overwrite v.to.rast input=inondazione output=$inondazione use=val layer=1 value=1 rows=4096 --overwrite #COLOURS OF INUNDATION RASTER MAP rm -f ".colore_inondazione.file" echo '1 blue' >> ".colore_inondazione.file" echo '0 grey' >> ".colore_inondazione.file" echo 'nv white' >> ".colore_inondazione.file" echo 'end' >> ".colore_inondazione.file" cat ".colore_inondazione.file" | r.colors map=$inondazione color=rules rm ".colore_inondazione.file" if [ -n "$GIS_OPT_DOUBT_MAP" ] then # It's possible to create a map with doubful surface areas dubbio="$GIS_OPT_DOUBT_MAP" r.mapcalculator amap=$inondazione formula='isnull(A)' outfile=temp help=- r.mapcalculator amap=temp formula='if(A,0,1)' outfile=temp help=- r.mapcalculator amap=inondazione_step2 bmap=temp formula='A-B' outfile=$dubbio help=- r.null map=$dubbio setnull=0 #COLOURS OF DOUBFUL RASTER MAP rm -f .colore_dubbio.file echo '1 red' >> .colore_dubbio.file echo '0 grey' >> .colore_dubbio.file echo 'nv white' >> .colore_dubbio.file echo 'end' >> .colore_dubbio.file cat .colore_dubbio.file | r.colors map=$dubbio color=rules rm .colore_dubbio.file g.remove rast=temp fi g.remove rast=MASK g.remove rast=inondazione_step1,inondazione_step2,inondazione_step3,reticolo,quote_punti_adiacenti g.remove rast=inondazione_step3_isnull,diff_inond,inondazione_step5,inondazione_step5_isnull,inondazione1 g.remove vect=inondazione_step1,inondazione_step2,reticolo g.remove vect=diff_inond1,reticolo_interno,inondazione1,overlap if [ "$inondazione" = "inondazione" ] then echo " " else g.remove vect=inondazione fi if [ -n "$GIS_OPT_MAP_W_S_PROFILE" ] then g.rename vect=profilo,"$GIS_OPT_MAP_W_S_PROFILE" else g.remove vect=profilo fi if [ -n "$GIS_OPT_river_boundary" ] then if [ "$GIS_OPT_boundary_type" = "line" ] then g.rename vect=bordi_alveo,"$GIS_OPT_river_boundary" g.remove vect=limiti_alveo else g.rename vect=limiti_alveo,"$GIS_OPT_river_boundary" g.remove vect=bordi_alveo fi else g.remove vect=bordi_alveo,limiti_alveo fi rm -r "$TMPDIR" exit 0