[GRASS-SVN] r36806 - grass-addons/raster/r.inund.fluv
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Apr 20 03:15:38 EDT 2009
Author: robertomarzocchi
Date: 2009-04-20 03:15:37 -0400 (Mon, 20 Apr 2009)
New Revision: 36806
Modified:
grass-addons/raster/r.inund.fluv/find_main_channel.f90
grass-addons/raster/r.inund.fluv/r.inund.fluv
Log:
Modified: grass-addons/raster/r.inund.fluv/find_main_channel.f90
===================================================================
--- grass-addons/raster/r.inund.fluv/find_main_channel.f90 2009-04-20 03:18:22 UTC (rev 36805)
+++ grass-addons/raster/r.inund.fluv/find_main_channel.f90 2009-04-20 07:15:37 UTC (rev 36806)
@@ -41,8 +41,10 @@
integer::pixel_temp_i, pixel_temp_j, temp_x, temp_y
!real, allocatable, dimension(:,:):: quota_terreno
!character(2), allocatable, dimension(:,:):: limiti_alveo
-real(kind=8)::x, E1, N1, E1dopo, N1dopo, z, livello, E2, N2, dir_ortog, dist, E_temp, N_temp, coord_N, coord_E, quota, dist_temp, quota_max
-real(kind=8):: dist_destra, dist_sinistra, verifica_dist, pendenza, pendenza_max, quota_limite, quota_centro
+real(kind=8)::x, E1, N1, E1dopo, N1dopo, z, livello, E2, N2, dir_ortog, dist
+real(kind=8):: E_temp, N_temp, coord_N, coord_E, quota, dist_temp, quota_max
+real(kind=8):: dist_destra, dist_sinistra, verifica_dist, pendenza, pendenza_max
+real(kind=8):: quota_limite, quota_centro
character(len=40)::xx, tab, tab1
real,dimension(2)::quota_temp
real, dimension(500):: E_ds, N_ds, E_sx, N_sx
Modified: grass-addons/raster/r.inund.fluv/r.inund.fluv
===================================================================
--- grass-addons/raster/r.inund.fluv/r.inund.fluv 2009-04-20 03:18:22 UTC (rev 36805)
+++ grass-addons/raster/r.inund.fluv/r.inund.fluv 2009-04-20 07:15:37 UTC (rev 36806)
@@ -1,513 +1,515 @@
-#!/bin/sh
-############################################################################
-#
-# MODULE: r.inund.fluv v.1 for GRASS 6.3 (april 2008)
-# 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 DTM raster map
-#% 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 line 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
-#%guisection:Optional parameters & maps
-#% key: PROFILE_BOUNDARY
-#% type: string
-#% gisprompt: old_file,file,input
-#% description: Input ASCII file with water-depht for return period T > 100 years
-#% required : no
-#%end
-#%option
-#%guisection:Optional parameters & maps
-#% key: delta_z
-#% type: double
-#% gisprompt: new
-#% description: Input delta_z to find the boundaries of the main channel [default value 0.5 m]
-#% required : no
-#%end
-#%option
-#%guisection:Optional parameters & maps
-#% key: delta_x
-#% type: double
-#% gisprompt: new
-#% description: Input delta_x to find the boundaries of the main channel [default value 3.5 m]
-#% required : no
-#%end
-#%option
-#%guisection:Optional parameters & maps
-#% key: river_boundary
-#% type: string
-#% gisprompt: new,vector,vector
-#% description: Output vector boundaries of the main channel
-#% required : no
-#%end
-#%option
-#%guisection:Optional parameters & maps
-#% key: boundary_type
-#% type: string
-#% description: Format of vector boundaries of the main channel
-#% options: line,points
-#% answer: points
-#% required: no
-#%end
-#%option
-#%guisection:Optional parameters & maps
-#% key: res_B
-#% type: integer
-#% gisprompt: new,cell,raster
-#% description: Input value: resolution B [default value 10 m]
-#% required : no
-#%end
-#%option
-#%guisection:Optional parameters & maps
-#% key: res_C
-#% type: integer
-#% gisprompt: new,cell,raster
-#% description: Input value: resolution C [default value 20 m]
-#% required : no
-#%end
-#%option
-#%guisection:Optional parameters & maps
-#% key: MAP_W_S_PROFILE
-#% type: string
-#% gisprompt: new,vector,vector
-#% description: Output vector point 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
-PROG=`basename $0`
-#request control & test (only for imput map)
-if [ -z "$GIS_OPT_DTM" ] ; then
- g.message -e "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
- g.message -w "WARNING - the DTM resolution it'snt good for this application"
-else
- g.message "The DTM resolution is $res x $res meters"
-fi
-
-
-rm .resolution.txt
-
-
-if [ -z "$GIS_OPT_W_S_PROFILE" ] ; then
- g.message -e "ERROR: file with depht not specified"
- exit 1
-fi
-
-
-
-if [ -z "$GIS_OPT_RIVER" ] ; then
- g.message -e "ERROR: vector with river-axis not specified"
- exit 1
-fi
-
-
-if [ -z "$GIS_OPT_FLOODING_MAP" ] ; then
- g.message -e "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 --quiet
-
-
-#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 --quiet
-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 --quiet
-
-#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
-g.message "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
- g.message -e "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 --quiet
-
-
-#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 --quiet
-g.region nsres=$res10 ewres=$res10 save=dtm_res10 --overwrite
-g.region -p >> "$TMPDIR"/region10.txt
-
-
-g.message "STEP 1"
-###################################################################################################################
-# STEP 1 - Thiessen interpolation
-###################################################################################################################
-# don't use low resolution [default 10m]
-r.mask input=$dtm maskcats=* --quiet
-g.region region=dtm_res10
-v.surf.idw input=profilo output=quote_punti_adiacenti npoints=1 layer=1 column=dbl_3 --quiet
-g.region region=dtm_res2
-r.mapcalculator amap=quote_punti_adiacenti bmap=$dtm formula="if(A>=B)" outfile=inondazione_step1 help=- --quiet
-
-
-g.message "STEP 2"
-###################################################################################################################
-# STEP 2 - remove "lakes"
-###################################################################################################################
-r.to.vect -s input=inondazione_step1 output=inondazione_step1 feature=area --overwrite --quiet
-v.select ainput=inondazione_step1 atype=line,area alayer=1 binput=$centerline btype=line blayer=1 output=inondazione_step2 operator=overlap --overwrite --quiet
-v.to.rast input=inondazione_step2 output=inondazione_step2 use=val layer=1 value=1 rows=4096 --overwrite --quiet
-
-
-g.message "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 --quiet
-r.out.ascii -h input=quote_punti_adiacenti output="$TMPDIR"/quote_punti_adiacenti dp=2 null=0 --quiet
-
-# FORTRAN code
-g.message "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
- g.message -e "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 --quiet
-
-
-g.message "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 --quiet
-
-g.message "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
- g.message -e "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 --quiet
-r.in.ascii input="$TMPDIR"/reticolo output=reticolo 'mult=1.0 or read from header' nv=* --overwrite --quiet
-r.thin input=reticolo output=reticolo_thin iterations=200 --overwrite --quiet
-r.to.vect input=reticolo_thin output=reticolo feature=line --overwrite --quiet
-g.remove rast=reticolo_thin --quiet
-
-
-# 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=- --quiet
-r.mapcalculator amap=inondazione_step3_isnull formula='if(A,0,1)' outfile=inondazione_step3 help=- --quiet
-r.mapcalculator amap=inondazione_step2 bmap=inondazione_step3 formula='A-B' outfile=diff_inond help=- --quiet
-r.mapcalculator amap=diff_inond formula='if(A,1,0,0)' outfile=diff_inond help=- --quiet #remove the -1 value resulted from A-B subtraction
-r.null map=diff_inond setnull=0 --quiet
-r.to.vect input=diff_inond output=diff_inond1 feature=area --overwrite --quiet
-#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 --quiet
-v.in.dxf -1 input="$TMPDIR"/reticolo.dxf output=reticolo --overwrite --quiet
-
-#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 --quiet
-
-# 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 --quiet
-
-# 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 --quiet
-v.drape input=punti_reticolo_interno type=point,centroid,line,boundary,face,kernel rast=reticolo method=nearest output=punti3d_reticolo --overwrite --quiet
-#v.drape input=punti_reticolo_interno type=point,centroid,line,boundary,face,kernel rast=quote_punti_adiacenti method=nearest output=punti3d_reticolo --overwrite --quiet
-
-
-g.region region=dtm_res10
-
-#export of 3d points
-v.out.ascii input=punti3d_reticolo output="$TMPDIR"/punti3d_reticolo format=point dp=2 --quiet
-r.out.ascii -h input=diff_inond output="$TMPDIR"/diff_inond null=0 --quiet
-
-
-
-g.message "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
- g.message -e "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 --quiet
-
-# replace null cells with 0 cells
-g.region region=dtm_res2
-r.mapcalculator amap=inondazione_step5 formula='isnull(A)' outfile=inondazione_step5_isnull help=- --quiet
-r.mapcalculator amap=inondazione_step5_isnull formula='if(A,0,1)' outfile=inondazione_step5 help=- --quiet
-r.mapcalculator amap=inondazione_step5 bmap=inondazione_step3 formula='A+B' outfile=inondazione help=- --quiet
-
-
-g.message "STEP 5"
-###################################################################################################################
-# STEP 5 - EXPORT OUTPUT FILES
-###################################################################################################################
-g.rename rast=inondazione,inondazione1 --overwrite --quiet
-r.null map=inondazione1 setnull=0 --quiet
-r.to.vect -s input=inondazione1 output=inondazione1 feature=area --overwrite --quiet
-v.patch input=centerline,punti3d_reticolo output=overlap --overwrite --quiet
-v.select ainput=inondazione1 atype=area alayer=1 binput=overlap btype=line blayer=1 output=inondazione operator=overlap --overwrite --quiet
-v.to.rast input=inondazione output=$inondazione use=val layer=1 value=1 rows=4096 --overwrite --quiet
-
-#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=- --quiet
- r.mapcalculator amap=temp formula='if(A,0,1)' outfile=temp help=- --quiet
- r.mapcalculator amap=inondazione_step2 bmap=temp formula='A-B' outfile=$dubbio help=- --quiet
- r.null map=$dubbio setnull=0 --quiet
- #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 --quiet
-g.remove rast=inondazione_step1,inondazione_step2,inondazione_step3,reticolo,quote_punti_adiacenti --quiet
-g.remove rast=inondazione_step3_isnull,diff_inond,inondazione_step5,inondazione_step5_isnull,inondazione1 --quiet
-g.remove vect=inondazione_step1,inondazione_step2,reticolo --quiet
-g.remove vect=diff_inond1,reticolo_interno,inondazione1,overlap --quiet
-
-if [ "$inondazione" == "inondazione" ]
-then
- g.message " "
-else
- g.remove vect=inondazione --quiet
-fi
-
-if [ -n "$GIS_OPT_MAP_W_S_PROFILE" ]
-then
- g.rename vect=profilo,"$GIS_OPT_MAP_W_S_PROFILE" --quiet
-else
- g.remove vect=profilo --quiet
-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 --quiet
- else
- g.rename vect=limiti_alveo,"$GIS_OPT_river_boundary"
- g.remove vect=bordi_alveo --quiet
- fi
-else
- g.remove vect=bordi_alveo,limiti_alveo --quiet
-fi
-
-rm -r "$TMPDIR"
-
-
-
-exit 0
+#!/bin/sh
+############################################################################
+#
+# MODULE: r.inund.fluv v.1 for GRASS 6.3 (april 2008)
+# 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 DTM raster map
+#% 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 line 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
+#%guisection:Optional parameters & maps
+#% key: PROFILE_BOUNDARY
+#% type: string
+#% gisprompt: old_file,file,input
+#% description: Input ASCII file with water-depht for return period T > 100 years
+#% required : no
+#%end
+#%option
+#%guisection:Optional parameters & maps
+#% key: delta_z
+#% type: double
+#% gisprompt: new
+#% description: Input delta_z to find the boundaries of the main channel [default value 0.5 m]
+#% required : no
+#%end
+#%option
+#%guisection:Optional parameters & maps
+#% key: delta_x
+#% type: double
+#% gisprompt: new
+#% description: Input delta_x to find the boundaries of the main channel [default value 3.5 m]
+#% required : no
+#%end
+#%option
+#%guisection:Optional parameters & maps
+#% key: river_boundary
+#% type: string
+#% gisprompt: new,vector,vector
+#% description: Output vector boundaries of the main channel
+#% required : no
+#%end
+#%option
+#%guisection:Optional parameters & maps
+#% key: boundary_type
+#% type: string
+#% description: Format of vector boundaries of the main channel
+#% options: line,points
+#% answer: points
+#% required: no
+#%end
+#%option
+#%guisection:Optional parameters & maps
+#% key: res_B
+#% type: integer
+#% gisprompt: new,cell,raster
+#% description: Input value: resolution B [default value 10 m]
+#% required : no
+#%end
+#%option
+#%guisection:Optional parameters & maps
+#% key: res_C
+#% type: integer
+#% gisprompt: new,cell,raster
+#% description: Input value: resolution C [default value 20 m]
+#% required : no
+#%end
+#%option
+#%guisection:Optional parameters & maps
+#% key: MAP_W_S_PROFILE
+#% type: string
+#% gisprompt: new,vector,vector
+#% description: Output vector point 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
+PROG=`basename $0`
+#request control & test (only for imput map)
+if [ -z "$GIS_OPT_DTM" ] ; then
+ g.message -e "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
+ g.message -w "WARNING - the DTM resolution it'snt good for this application"
+else
+ g.message "The DTM resolution is $res x $res meters"
+fi
+
+
+rm .resolution.txt
+
+
+if [ -z "$GIS_OPT_W_S_PROFILE" ] ; then
+ g.message -e "ERROR: file with depht not specified"
+ exit 1
+fi
+
+
+
+if [ -z "$GIS_OPT_RIVER" ] ; then
+ g.message -e "ERROR: vector with river-axis not specified"
+ exit 1
+fi
+
+
+if [ -z "$GIS_OPT_FLOODING_MAP" ] ; then
+ g.message -e "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 --quiet
+
+
+#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 --quiet
+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 --quiet
+
+#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
+g.message "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
+ g.message -e "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 --quiet
+
+
+#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 --quiet
+g.region nsres=$res10 ewres=$res10 save=dtm_res10 --overwrite
+g.region -p >> "$TMPDIR"/region10.txt
+
+
+g.message "STEP 1"
+###################################################################################################################
+# STEP 1 - Thiessen interpolation
+###################################################################################################################
+# don't use low resolution [default 10m]
+r.mask input=$dtm maskcats=* --quiet
+g.region region=dtm_res10
+v.surf.idw input=profilo output=quote_punti_adiacenti npoints=1 layer=1 column=dbl_3 --quiet
+g.region region=dtm_res2
+r.mapcalculator amap=quote_punti_adiacenti bmap=$dtm formula="if(A>=B)" outfile=inondazione_step1 help=- --quiet
+
+
+g.message "STEP 2"
+###################################################################################################################
+# STEP 2 - remove "lakes"
+###################################################################################################################
+r.to.vect -s input=inondazione_step1 output=inondazione_step1 feature=area --overwrite --quiet
+v.select ainput=inondazione_step1 atype=line,area alayer=1 binput=$centerline btype=line blayer=1 output=inondazione_step2 operator=overlap --overwrite --quiet
+v.to.rast input=inondazione_step2 output=inondazione_step2 use=val layer=1 value=1 rows=4096 --overwrite --quiet
+
+
+g.message "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 --quiet
+r.out.ascii -h input=quote_punti_adiacenti output="$TMPDIR"/quote_punti_adiacenti dp=2 null=0 --quiet
+
+# FORTRAN code
+g.message "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
+ g.message -e "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 --quiet
+
+
+g.message "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 --quiet
+
+g.message "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
+ g.message -e "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 --quiet
+r.in.ascii input="$TMPDIR"/reticolo output=reticolo 'mult=1.0 or read from header' nv=* --overwrite --quiet
+r.thin input=reticolo output=reticolo_thin iterations=200 --overwrite --quiet
+r.to.vect input=reticolo_thin output=reticolo feature=line --overwrite --quiet
+g.remove rast=reticolo_thin --quiet
+
+
+# 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=- --quiet
+r.mapcalculator amap=inondazione_step3_isnull formula='if(A,0,1)' outfile=inondazione_step3 help=- --quiet
+r.mapcalculator amap=inondazione_step2 bmap=inondazione_step3 formula='A-B' outfile=diff_inond help=- --quiet
+r.mapcalculator amap=diff_inond formula='if(A,1,0,0)' outfile=diff_inond help=- --quiet #remove the -1 value resulted from A-B subtraction
+r.null map=diff_inond setnull=0 --quiet
+r.to.vect input=diff_inond output=diff_inond1 feature=area --overwrite --quiet
+#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 --quiet
+v.in.dxf -1 input="$TMPDIR"/reticolo.dxf output=reticolo --overwrite --quiet
+
+#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 --quiet
+
+# 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 --quiet
+
+# 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 --quiet
+v.drape input=punti_reticolo_interno type=point,centroid,line,boundary,face,kernel rast=reticolo method=nearest output=punti3d_reticolo --overwrite --quiet
+#v.drape input=punti_reticolo_interno type=point,centroid,line,boundary,face,kernel rast=quote_punti_adiacenti method=nearest output=punti3d_reticolo --overwrite --quiet
+
+
+g.region region=dtm_res10
+
+#export of 3d points
+v.out.ascii input=punti3d_reticolo output="$TMPDIR"/punti3d_reticolo format=point dp=2 --quiet
+r.out.ascii -h input=diff_inond output="$TMPDIR"/diff_inond null=0 --quiet
+
+
+
+g.message "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
+ g.message -e "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 --quiet
+
+# replace null cells with 0 cells
+g.region region=dtm_res2
+r.mapcalculator amap=inondazione_step5 formula='isnull(A)' outfile=inondazione_step5_isnull help=- --quiet
+r.mapcalculator amap=inondazione_step5_isnull formula='if(A,0,1)' outfile=inondazione_step5 help=- --quiet
+r.mapcalculator amap=inondazione_step5 bmap=inondazione_step3 formula='A+B' outfile=inondazione help=- --quiet
+
+
+g.message "STEP 5"
+###################################################################################################################
+# STEP 5 - EXPORT OUTPUT FILES
+###################################################################################################################
+g.rename rast=inondazione,inondazione1 --overwrite --quiet
+r.null map=inondazione1 setnull=0 --quiet
+r.to.vect -s input=inondazione1 output=inondazione1 feature=area --overwrite --quiet
+v.patch input=centerline,punti3d_reticolo output=overlap --overwrite --quiet
+v.select ainput=inondazione1 atype=area alayer=1 binput=overlap btype=line blayer=1 output=inondazione operator=overlap --overwrite --quiet
+v.to.rast input=inondazione output=$inondazione use=val layer=1 value=1 rows=4096 --overwrite --quiet
+#water level assignement
+r.mapcalculator amap=$inondazione bmap=quote_punti_adiacenti cmap=$dtm formula='A*B-C' outfile=$inondazione help=- --quiet
+
+#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=- --quiet
+ r.mapcalculator amap=temp formula='if(A,0,1)' outfile=temp help=- --quiet
+ r.mapcalculator amap=inondazione_step2 bmap=temp formula='A-B' outfile=$dubbio help=- --quiet
+ r.null map=$dubbio setnull=0 --quiet
+ #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 --quiet
+g.remove rast=inondazione_step1,inondazione_step2,inondazione_step3,reticolo,quote_punti_adiacenti --quiet
+g.remove rast=inondazione_step3_isnull,diff_inond,inondazione_step5,inondazione_step5_isnull,inondazione1 --quiet
+g.remove vect=inondazione_step1,inondazione_step2,reticolo --quiet
+g.remove vect=diff_inond1,reticolo_interno,inondazione1,overlap --quiet
+
+if [ "$inondazione" == "inondazione" ]
+then
+ g.message " "
+else
+ g.remove vect=inondazione --quiet
+fi
+
+if [ -n "$GIS_OPT_MAP_W_S_PROFILE" ]
+then
+ g.rename vect=profilo,"$GIS_OPT_MAP_W_S_PROFILE" --quiet
+else
+ g.remove vect=profilo --quiet
+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 --quiet
+ else
+ g.rename vect=limiti_alveo,"$GIS_OPT_river_boundary"
+ g.remove vect=bordi_alveo --quiet
+ fi
+else
+ g.remove vect=bordi_alveo,limiti_alveo --quiet
+fi
+
+rm -r "$TMPDIR"
+
+
+
+exit 0
More information about the grass-commit
mailing list