[GRASS-SVN] r56036 - in grass-addons/grass6/raster: . r.tsunami
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Apr 29 07:03:54 PDT 2013
Author: maxi
Date: 2013-04-29 07:03:53 -0700 (Mon, 29 Apr 2013)
New Revision: 56036
Added:
grass-addons/grass6/raster/r.tsunami/
grass-addons/grass6/raster/r.tsunami/r.tsunami
grass-addons/grass6/raster/r.tsunami/r.tsunami.html
Log:
module r.tsunami
Added: grass-addons/grass6/raster/r.tsunami/r.tsunami
===================================================================
--- grass-addons/grass6/raster/r.tsunami/r.tsunami (rev 0)
+++ grass-addons/grass6/raster/r.tsunami/r.tsunami 2013-04-29 14:03:53 UTC (rev 56036)
@@ -0,0 +1,362 @@
+#!/usr/bin/python
+# -*- coding:utf-8 -*-
+############################################################################
+#
+# MODULE: r.tsunami
+#
+# AUTHOR(S): M.Molinari, M. Cannata (Earth Science Institute, http://istgeo.ist.supsi.ch)
+# Implementation of the base method developed by B. Federici
+#
+# PURPOSE: Calculation of flooded area due to tsunami event
+#
+# COPYRIGHT: (c) 2007 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: calculation of flooded area due to tsunami event
+#% keywords: tsunami, modelling
+#%End
+#%Flag
+#% key: m
+#% description: Runup values (m)
+#%End
+#%Flag
+#% key: a
+#% description: Calculation of total flooded area (mq)
+#%End
+#%Flag
+#% key: i
+#% description: Flood height values (m)
+#%End
+#%option
+#% key: dtm
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Digital Terrain Model
+#% required : yes
+#%end
+#%option
+#% key: rough
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Roughness map
+#% required : yes
+#%end
+#%option
+#% key: size
+#% type: integer
+#% options:1,3,5,7,9,11,13,15,17,19,21,23,25
+#% description: Size of the neighborhood
+#% required : yes
+#% answer : 3
+#%end
+#%option
+#% key: h
+#% type: string
+#% description: Sea depth at coastline (m)
+#% required : yes
+#%end
+#%option
+#% key: ho
+#% type: string
+#% description: Sea depth at source point (m)
+#% required : yes
+#%end
+#%option
+#% key: Ho
+#% type: string
+#% description: Wave height at source point (m)
+#% required : yes
+#%end
+#%option
+#% key: runup
+#% type: string
+#% gisprompt: new,cell,raster
+#% description:Runup map
+#% required : no
+#%end
+#%option
+#% key: flood
+#% type: string
+#% gisprompt: new,cell,raster
+#% description:Flooded area map
+#% required : yes
+#%end
+#%option
+#% key: hflood
+#% type: string
+#% gisprompt: new,cell,raster
+#% description:Flood height map
+#% required : yes
+#%end
+#%option
+#% key: x
+#% type: string
+#% description: X coordinate (sea point)
+#% required : yes
+#%end
+#%option
+#% key: y
+#% type: string
+#% description: Y coordinate (sea point)
+#% required : yes
+#%end
+
+import os,sys,subprocess,string,random,math,tempfile,time,re
+
+def main():
+
+ #### get debug level ####
+ cmdargs00=["g.gisenv","get=DEBUG"]
+ proc00 = subprocess.Popen(cmdargs00, stdout=subprocess.PIPE,stderr=subprocess.PIPE)
+ d = proc00.communicate()[0]
+
+ if d:
+ if int(d)!=0:
+ print''
+ sys.stdout.write ("Input data:\n")
+ sys.stdout.write ("Value of h: %s\n" % os.getenv("GIS_OPT_h"))
+ sys.stdout.write ("Value of ho: %s\n" % os.getenv("GIS_OPT_ho"))
+ sys.stdout.write ("Value of Ho: %s\n" % os.getenv("GIS_OPT_Ho"))
+ sys.stdout.write ("Value of size: %s\n" % os.getenv("GIS_OPT_size"))
+ sys.stdout.write ("Name of dtm: %s\n" % os.getenv("GIS_OPT_dtm"))
+ sys.stdout.write ("Name of roughness_map: %s\n" % os.getenv("GIS_OPT_rough"))
+
+ #### get inputs ####
+ size = int(os.getenv('GIS_OPT_size'))
+ h = float(os.getenv('GIS_OPT_h'))
+ ho = float(os.getenv('GIS_OPT_ho'))
+ Ho = float(os.getenv('GIS_OPT_Ho'))
+ dtm = os.getenv('GIS_OPT_dtm')
+ roughness_map = os.getenv('GIS_OPT_rough')
+ output1 = os.getenv('GIS_OPT_runup')
+ output2 = os.getenv('GIS_OPT_flood')
+ output3 = os.getenv('GIS_OPT_hflood')
+ x = float(os.getenv('GIS_OPT_x'))
+ y = float(os.getenv('GIS_OPT_y'))
+ a = os.getenv('GIS_FLAG_a')
+ m = os.getenv('GIS_FLAG_m')
+ i = os.getenv('GIS_FLAG_i')
+
+ #### setup temporary files ####
+ p=re.compile('\W')
+ slope = p.sub("",tempfile.mktemp())
+ slope_med = p.sub("",tempfile.mktemp())
+ new_value = p.sub("",tempfile.mktemp())
+ sea_mask = p.sub("",tempfile.mktemp())
+ r_coast = p.sub("",tempfile.mktemp())
+ slope_r_coast= p.sub("",tempfile.mktemp())
+ run_up= p.sub("",tempfile.mktemp())
+ run_up_ok = p.sub("",tempfile.mktemp())
+ run_up_idw = p.sub("",tempfile.mktemp())
+ run_up_idw2 = p.sub("",tempfile.mktemp())
+ b500_coast = p.sub("",tempfile.mktemp())
+ b500_runup = p.sub("",tempfile.mktemp())
+ runup = p.sub("",tempfile.mktemp())
+ inond_provv = p.sub("",tempfile.mktemp())
+ lake = p.sub("",tempfile.mktemp())
+ lake_ok = p.sub("",tempfile.mktemp())
+ lake_ok2 = p.sub("",tempfile.mktemp())
+ output_3 = p.sub("",tempfile.mktemp())
+
+ #### Estimate slope ####
+ time.sleep(10)
+ cmdargs1= ["elevation=%s"%dtm,"slope=%s"%slope, "format=percent", "--overwrite","--quiet"]
+ os.spawnvp(os.P_WAIT,"r.slope.aspect", ["r.slope.aspect"] + cmdargs1)
+
+ #### Estimate mean slope 3x3 ####
+ time.sleep(2)
+ cmdargs2=["input=%s"%slope, "output=%s"%slope_med, "method=average", "size=%s"%size, "--overwrite","--quiet"]
+ os.spawnvp(os.P_WAIT,"r.neighbors", ["r.neighbors"] + cmdargs2)
+
+ #### Extract coastline ####
+ time.sleep(2)
+ cmdargs3=["%s=if(isnull(%s),10000,%s)"%(new_value,dtm,dtm)]
+ os.spawnvp(os.P_WAIT,"r.mapcalc", ["r.mapcalc"] + cmdargs3)
+
+ time.sleep(2)
+ cmdargs4=["input=%s"%new_value, "output=%s"%sea_mask, "method=sum", "size=3", "--overwrite","--quiet"]
+ os.spawnvp(os.P_WAIT,"r.neighbors", ["r.neighbors"] + cmdargs4)
+
+ time.sleep(2)
+ cmdargs6=["%s=if(%s>10000 && %s<80000,1,null())"%(r_coast,sea_mask,sea_mask)]
+ os.spawnvp(os.P_WAIT,"r.mapcalc", ["r.mapcalc"] + cmdargs6)
+
+ #### Assign mean slope to coastline ####
+ time.sleep(2)
+ cmdargs7=["%s=if(%s>=1,%s,0)"%(slope_r_coast,r_coast,slope_med)]
+ os.spawnvp(os.P_WAIT,"r.mapcalc", ["r.mapcalc"] + cmdargs7)
+
+ #### Estimate hydraulic parameters ####
+ Co=math.sqrt(9.8 * h)
+ H= Ho * ((ho/h)**(0.25))
+ U=Co * math.sqrt((1+H/h))
+ L= 4 * 2.415 *h /(math.sqrt(3*H/h))
+ k1=U**2 * L/(19.6)
+ k2= U**2 * H/(39.2)
+
+ if d:
+ if int(d)!=0:
+ print''
+ sys.stdout.write ("Celerità dell'onda=%s\n" %Co)
+ sys.stdout.write ("Altezza d'onda a riva=%s\n" %H)
+ sys.stdout.write ("Velocità d'onda a riva=%s\n" %U)
+ sys.stdout.write ("Lunghezza di un'onda solitaria=%s\n" %L)
+ sys.stdout.write ("k1=%s\n" %k1)
+ sys.stdout.write ("k2=%s\n"%k2)
+
+ #### Estimates run-up map ####
+ cmdargs9=["%s = if(%s<=0,%s,%s+sqrt(%s*float(%s)/100+%s))" % (run_up,slope_r_coast,H, H, k1,slope_r_coast, k2)]
+ os.spawnvp(os.P_WAIT,"r.mapcalc", ["r.mapcalc"] + cmdargs9)
+
+ time.sleep(2)
+ cmdargs10=["%s = %s*10^6" %(run_up_ok,run_up)]
+ os.spawnvp(os.P_WAIT,"r.mapcalc", ["r.mapcalc"] + cmdargs10)
+
+ time.sleep(2)
+ cmdargs11=["input=%s"%run_up_ok,"output=%s"%run_up_idw, "npoints=1","--quiet"]
+ os.spawnvp(os.P_WAIT,"r.surf.idw", ["r.surf.idw"] + cmdargs11)
+
+ time.sleep(2)
+ cmdargs12=["%s=float(%s)/10^6"%(run_up_idw2,run_up_idw)]
+ os.spawnvp(os.P_WAIT,"r.mapcalc", ["r.mapcalc"] + cmdargs12)
+
+ time.sleep(2)
+ cmdargs13=["input=%s"%r_coast,"output=%s"%b500_coast, "distances=500", "units=meters","--quiet"]
+ os.spawnvp(os.P_WAIT,"r.buffer", ["r.buffer"] + cmdargs13)
+
+ time.sleep(2)
+ cmdargs14=["%s=if(%s>=1,%s,0)"%(b500_runup,b500_coast,run_up_idw2)]
+ os.spawnvp(os.P_WAIT,"r.mapcalc", ["r.mapcalc"] + cmdargs14)
+
+ opt_name='%s'%output1
+
+ if (len(opt_name)==0):
+ cmdargs15=["%s=if(%s>=0,%s,0)" % (runup,dtm,b500_runup)]
+ os.spawnvp(os.P_WAIT,"r.mapcalc", ["r.mapcalc"] + cmdargs15)
+ else:
+ cmdargs15=["%s=if(%s>=0,%s,0)" % (output1,dtm,b500_runup)]
+ os.spawnvp(os.P_WAIT,"r.mapcalc", ["r.mapcalc"] + cmdargs15)
+
+ if d:
+ if int(d)!=0:
+ print ("Run-up map...complete!\n")
+
+ #### Estimate flooded area map ####
+ time.sleep(2)
+ cmdargs16=["%s=if(float(%s)*%s-%s>0,1,0)" % (inond_provv,roughness_map,run_up_idw2,dtm)]
+ os.spawnvp(os.P_WAIT,"r.mapcalc", ["r.mapcalc"] + cmdargs16)
+
+ time.sleep(2)
+ cmdargs18=["%s=if(%s==0,10000,2)"%(lake,inond_provv)]
+ os.spawnvp(os.P_WAIT,"r.mapcalc", ["r.mapcalc"] + cmdargs18)
+
+ time.sleep(2)
+ cmdargs17=["map=%s"%lake,"null=0","--quiet"]
+ os.spawnvp(os.P_WAIT,"r.null", ["r.null"] + cmdargs17)
+
+ time.sleep(2)
+ cmdargs19=["dem=%s"%lake,"wl=100","lake=%s"%lake_ok,"xy=%s,%s"%(x,y),"--o","--quiet"]
+ os.spawnvp(os.P_WAIT,"r.lake", ["r.lake"] + cmdargs19)
+
+ time.sleep(2)
+ cmdargs22=["%s=if(isnull(%s),0,%s)"%(lake_ok2,lake_ok,lake_ok)]
+ os.spawnvp(os.P_WAIT,"r.mapcalc", ["r.mapcalc"] + cmdargs22)
+
+ time.sleep(2)
+ cmdargs23=["%s=if(%s,%s,0)"%(output2,dtm,lake_ok2)]
+ os.spawnvp(os.P_WAIT,"r.mapcalc", ["r.mapcalc"] + cmdargs23)
+
+ if d:
+ if int(d)!=0:
+ print ("Flooded area map...complete!")
+
+ #### Estimate flood height map ####
+ time.sleep(2)
+ cmdargs25=["%s=if(%s==98,%s*float(%s)-%s,0)"%(output3,output2,run_up_idw2,roughness_map,dtm)]
+ os.spawnvp(os.P_WAIT,"r.mapcalc", ["r.mapcalc"] + cmdargs25)
+ if d:
+ if int(d)!=0:
+ print ("Flood height map...complete!")
+
+ #### Estract run-up statistics according to flags ####
+ time.sleep(2)
+ if ('%s'%m=='1'):
+ cmdargs1=["r.univar","map=%s" %run_up]
+ proc1 = subprocess.Popen(cmdargs1, stdout=subprocess.PIPE,stderr=subprocess.PIPE)
+ min_run= proc1.communicate()[0].split("\n")[6]
+
+ cmdargs11=["r.univar","map=%s" %run_up]
+ proc11 = subprocess.Popen(cmdargs11, stdout=subprocess.PIPE,stderr=subprocess.PIPE)
+ max_run= proc11.communicate()[0].split("\n")[7]
+
+ cmdargs111=["r.univar","map=%s" %run_up]
+ proc111 = subprocess.Popen(cmdargs111, stdout=subprocess.PIPE,stderr=subprocess.PIPE)
+ mean_run= proc111.communicate()[0].split("\n")[9]
+
+ time.sleep(2)
+ if ('%s'%a=='1'):
+ cmdargs2=["r.stats","input=%s" %output2,"-ai","fs=space","nv=*","nsteps=255"]
+ proc2 = subprocess.Popen(cmdargs2, stdout=subprocess.PIPE,stderr=subprocess.PIPE)
+ r = proc2.communicate()[0].split("\n")[1]
+ ra = r.split(" ")[1]
+
+
+ time.sleep(2)
+ if ('%s'%i=='1'):
+ cmdargs3333=["%s=%s"%(output_3,output3)]
+ os.spawnvp(os.P_WAIT,"r.mapcalc", ["r.mapcalc"] + cmdargs3333)
+
+ cmdargs3333=["map=%s"%output_3,"setnull=0","--quiet"]
+ os.spawnvp(os.P_WAIT,"r.null", ["r.null"] + cmdargs3333)
+
+ cmdargs333=["r.univar","map=%s" %output_3]
+ proc333 = subprocess.Popen(cmdargs333, stdout=subprocess.PIPE,stderr=subprocess.PIPE)
+ mean_inond= proc333.communicate()[0].split("\n")[9]
+
+ if ('%s'%m=='1' or '%s'%a=='1' or '%s'%i=='1'):
+ print "************************************************"
+ print "************ TSUNAMI STATISTICS ****************"
+ print "************************************************"
+ if ('%s'%m=='1'):
+ print "* run-up height (m):"
+ print "* %s"%(min_run)
+ print "* %s"%(max_run)
+ print "* %s"%(mean_run)
+ print "* ------------------------------------------"
+ if ('%s'%a=='1'):
+ print "* estimated flooded surface (sqm): "
+ print "* area: %.2f"%(float(ra.strip('\'"')))
+ print "* ------------------------------------------"
+ if ('%s'%i=='1'):
+ print "* estimated flood height (m):"
+ print "* %s"%(mean_inond)
+ print "* ------------------------------------------"
+ print "************************************************"
+
+
+ #### Clean-up temporary files ####
+ time.sleep(2)
+ if (len(opt_name)==0):
+ cmdargs9=["%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s"%(slope,slope_med,r_coast,slope_r_coast,run_up,run_up_ok,run_up_idw,run_up_idw2,b500_coast,b500_runup,runup,inond_provv,lake,lake_ok,lake_ok2,new_value,sea_mask),"--quiet"]
+ os.spawnvp(os.P_WAIT,"g.remove", ["g.remove"] + cmdargs9)
+ else:
+ cmdargs99=["%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s"%(slope,slope_med,r_coast,slope_r_coast,run_up,run_up_ok,run_up_idw,run_up_idw2,b500_coast,b500_runup,inond_provv,lake,lake_ok,lake_ok2,new_value,sea_mask),"--quiet"]
+ os.spawnvp(os.P_WAIT,"g.remove", ["g.remove"] + cmdargs99)
+
+ if ('%s'%i=='1'):
+ cmdargs999=["%s"%output_3,"--quiet"]
+ os.spawnvp(os.P_WAIT,"g.remove", ["g.remove"] + cmdargs999)
+
+
+ return
+
+if __name__ == "__main__":
+ if ( len(sys.argv) <= 1 or sys.argv[1] != "@ARGS_PARSED@" ):
+ os.execvp("g.parser", [sys.argv[0]] + sys.argv)
+ else:
+ main();
Property changes on: grass-addons/grass6/raster/r.tsunami/r.tsunami
___________________________________________________________________
Added: svn:executable
+ *
Added: grass-addons/grass6/raster/r.tsunami/r.tsunami.html
===================================================================
--- grass-addons/grass6/raster/r.tsunami/r.tsunami.html (rev 0)
+++ grass-addons/grass6/raster/r.tsunami/r.tsunami.html 2013-04-29 14:03:53 UTC (rev 56036)
@@ -0,0 +1,115 @@
+<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
+<html>
+<head>
+<title>GRASS GIS: r.tsunami</title>
+<meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
+<link rel="stylesheet" href="grassdocs.css" type="text/css">
+</head>
+<body bgcolor="white">
+
+<img src="grass_logo.png" alt="GRASS logo"><hr align=center size=6 noshade>
+
+<h2>NAME</h2>
+<em><b>r.tsunami</b></em> - Calculation of flooded area due to tsunami event.
+<h2>KEYWORDS</h2>
+tsunami, modelling
+<h2>SYNOPSIS</h2>
+<b>r.tsunami</b><br>
+<b>r.tsunami help</b><br>
+<b>r.tsunami</b> <b>dtm</b>=<em>string</em> <b>rough</b>=<em>string</em> <b>size</b>=<em>value</em> <b>h</b>=<em>value</em> <b>ho</b>=<em>value</em>
+<b>Ho</b>=<em>value</em> [<b>runup</b>=<em>string</em>] <b>flood</b>=<em>string</em> <b>hflood</b>=<em>string</em> <b>x</b>=<em>value</em> <b>y</b>=<em>value</em> [--<b>overwrite</b>] [--<b>verbose</b>] [--<b>quiet</b>]
+
+<h3>Flags:</h3>
+<DL>
+<DT><b>-m</b></DT>
+<DD>Runup values [m]</DD>
+<DT><b>-a</b></DT>
+<DD>Calculation of total flooded area [mq]</DD>
+<DT><b>-i</b></DT>
+<DD>Flood height values [m]</DD>
+<DT><b>--overwrite</b></DT>
+<DD>Allow output files to overwrite existing files</DD>
+<DT><b>--verbose</b></DT>
+<DD>Verbose module output</DD>
+<DT><b>--quiet</b></DT>
+<DD>Quiet module output</DD>
+</DL>
+
+<h3>Parameters:</h3>
+<DL>
+<DT><b>dtm</b>=<em>string</em></DT>
+<DD>Digital terrain model raster map name</DD>
+
+<DT><b>rough</b>=<em>string</em></DT>
+<DD>Roughness raster map name</DD>
+
+<DT><b>size</b>=<em>value</em></DT>
+<DD>Size of the neighborhood [-]</DD>
+<DD>Options: <em> 1,3,5,7,9,11,13,15,17,19,21,23,25</em></DD>
+<DD>Default: <em>3</em></DD>
+
+<DT><b>h</b>=<em>value</em></DT>
+<DD> Sea depth at coastline [m]</DD>
+
+<DT><b>ho</b>=<em>value</em></DT>
+<DD>Sea depth at source point [m]</DD>
+
+<DT><b>Ho</b>=<em>value</em></DT>
+<DD>Wave height at source point [m]</DD>
+
+<DT><b>runup</b>=<em>string</em></DT>
+<DD>Output runup map name</DD>
+
+<DT><b>flood</b>=<em>string</em></DT>
+<DD>Output flooded area map name</DD>
+
+<DT><b>hflood</b>=<em>string</em></DT>
+<DD>Output flood height map name</DD>
+
+<DT><b>x</b>=<em>value</em></DT>
+<DD>X coordinate (sea point)</DD>
+
+<DT><b>y</b>=<em>value</em></DT>
+<DD>Y coordinate (sea point)</DD>
+</DL>
+<h2>DESCRIPTION</h2>
+
+<em>r.tsunami</em> allow to assess the maximum vertical height
+of the tsunami wave at the coast (run-up) and the subsequent
+inundation on the mainland.The procedure takes into account local
+morphological characteristics,vegetation and urbanization
+through the usage of digital terrain model and roughness raster
+maps.
+
+
+<h2>NOTES</h2>
+
+The digital terrain model raster map must contain null values at sea level.
+
+<p>
+
+
+
+<p>
+
+
+<h2>AUTHORS</h2>
+
+M.Molinari, M. Cannata (Earth Science Institute, http://istgeo.ist.supsi.ch).
+Implementation of the base method developed by B. Federici.
+
+<h2>REFERENCES</h2>
+
+<ul>
+ <li>Federici, B & Cosso, T, 2006, 'Tsunami inundation maps and damage sceneries through the GIS GRASS',<a href="http://www.foss4g2006.org">
+ proceedings of FOSS4G2006 </a>, Lausanne, Switzerland, viewed 13 August 2008. </li>
+ <li>Cannata M., Federici B., Molinari M., - "Carte di inondazione dovute a tsunami mediante il GIS GRASS: applicazione all'isola di St. Lucia, Caraibi",proceedings of <a href ="http://gislab.dirap.unipa.it/grass_meeting/index.htm">VIII Italian GRASS user meeting</a>, Palermo, Italy, viewed 13 August 2008. </li>
+</ul>
+
+
+<p><i>Last changed: $Date: 2008-05-16 12:09:06 -0700 (Fri, 16 May 2008) $</i>
+<HR>
+<P><a href="index.html">Main index</a> - <a href="raster.html">raster index</a> - <a href="full_index.html">Full index</a></P>
+<P>© 2003-2008 <a href="http://grass.osgeo.org">GRASS Development Team</a></p>
+</body>
+</html>
More information about the grass-commit
mailing list