#!/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();