Skip navigation.
Istituto Scienze della Terra - SUPSI

Tsunami Risk Lab

This Lab proposes a simplified approach for tsunami risk mapping within OSGeo (Open Source Geospatial) software. In particular this lab uses the Open Source GRASS GIS for risk analysis. 

Working Region

In order to proceed with analysis and calculation, GRASS requires to set the computational region. For this tutorial the region settings are loaded with the following command:
g.region n=4530358 s=4529364 e=542604 w=541054 res=2 

Hazard Assessment

Definition of event statistics

The first task to be addressed is the definition of the return periods and severity of tsunami events in your region. The required information is therefore:
  1. the definition of return periods T associated with high, medium and low probability
  2. the definition of parameter values associated with high, medium and low event severity
The tool proposed in this lab to estimate the tsunami hazard is the r.tsunami GRASS script that requires as inputs the height of the generated wave (Ho) and the sea depth (ho) in the tsunami source point.
This information (expectation that every T years a wave high Ho is generated in a point with sea depth of ho) is generally retrieved by means of literature review, tectonic informations and historical data analysis (a useful source is Tsunami Data at NGDC).
For the aim of this lab we selected the following values after Camilleri, 2006 [1]:
Return period
ho [m]
Ho [m] h [m]
100 700 0.1 2
250 700 0.25 2
500 700 0.5 2

Preparation of input data

Before the simulation model can be run, all the required inputs have to be prepared. In particular the r.tsunami model needs:
  • digital terrain model raster map (dtm)
  • roughness raster map (rough)
  • height of the generated wave in the source point (ho)
  • depth of the sea in the source point (Ho)
  • sea depth near to the coast (h)
  • output flooded area map name (flood)
  • output flood height map name (hflood)
  • coordinates of a point on the sea (x) and (y)
Digital surface model: already available in the demo dataset.
Roughness map: this map can be derived from the landuse map by means of reclassification. Each category can be associated with a roughness values (a sort of water friction parameter).
Corine values
Roughness value
521, 121 1
122, 111, 141
142, 112, 133, 123 0.8
322, 243, 323 0.5
In order to build the new roughness map we have to use the following map algebra command:
r.mapcalc 'roughness = \
if(corine==521||corine==121,1, \
if(corine==122||corine==111||corine==141,0.9, \
if(corine==142||corine==112||corine==133 || corine==123,0.8, \
if( corine==322 || corine==243 || corine==323, 0.5, null() ))))'
Tsunami parameters: for the scope of this lab we have already selected the values ho, Ho and h that characterize events.
Buildings: already available in the demo data, we just need to copy the layer for our own use:
g.copy vect=buildings@PERMANENT,build
Sea point coordinates: just display the dtm, and query the map over the see. For example we could use the following values:
542218 4529858

Simulation of events

At this point all the pieces of information required to simulate tsunami events for the given return periods are known. The r.tsunami command can be executed:
r.tsunami dtm=dtm rough=roughness h=2 ho=700 Ho=0.1 runup=runup100 flood=f100 hflood=hf100 \
x=542218 y=4529858 -mai --o
r.tsunami dtm=dtm rough=roughness h=2 ho=700 Ho=0.25 runup=runup250 flood=f250 hflood=hf250 \
x=542218 y=4529858 -mai --o
r.tsunami dtm=dtm rough=roughness h=2 ho=700 Ho=0.5 runup=runup500 flood=f500 hflood=hf500 \
x=542218 y=4529858 -mai --o 

We obtain for each period the following raster maps:

runup = map of the wave runup height
flood = 2 values map denoting flooded and not flooded areas
hflood = flooding height map

Hazard classification map generation

To derive classified hazard map, the intensity class boundaries values have to be defined. The selected values (1, 3) has been derived from the Indian Government [2].
Those values have to be used as input of the r.hazard.class command to extract the hazard maps for low, medium, and high return periods, as weel as the total hazard maps:
r.hazard.class -c l_hflood=hf100 m_hflood=hf250 h_hflood=hf500 low_ub=1 med_ub=3 \
l_class=low m_class=med h_class=high t_class=tot

Exposition Assessment

The exposed elements are estimated by means of layer intersection between the object layer and the hazard layer. In this case the event intensity layers (hf100, hf250, hf500) are raster maps while the object map (buildings) is a vector.
For this propose a new command that samples the statistical parameter of the cells intersected by each single feature has been developed: v.sample.class (see man page for more details).
Given a vector map, a raster map, a prefix and one or more modes (minimum, maximum, average, sum) this command creates new columns named prefix_mode containing the required statistical values.
The next step is therefore to sample the three intensity parameter maps using the maximum value (for the worst case):
v.sample.class input=build raster=hf100 pre=h100 mode=max 
v.sample.class input=build raster=hf250 pre=h250 mode=max
v.sample.class input=build raster=hf500 pre=h500 mode=max
Now the vector map has three new columns named: h100_max, h250_max, h500_max.
Exposed elements for each return period are identified by buildings with  f100_max, f250_max, and f500_max respectively different from 0.

Vulnerability assessment

Vulnerability can be estimated by means of past events data analysis studies. For the case of inundation the U.S. Army Engineer Corps [3] individuates a serie of functions that link the flood height with the damage expressed as percentage of building value based on building typology :
1) Residential
Element Vulnerability function (damage as % of value) Condition
flood > 0
flood > 0
2) Agricultural
Element Vulnerability function (damage as % of value) Condition
0 < flood <= 4
flood > 4
0 < flood <= 4
content 0.5 flood > 4

3) Industrial 

Element Vulnerability function (damage as % of value) Condition
0 < flood <= 5
structure 0.5 flood > 5
0 < flood <= 1.2
content 1 flood > 1.2

Risk Assessment

The first step is to estimate values of objects and of their content. This can be done by associating a value per square meter as a function of the building type. Then the risk associated at each return period can be estimated as expected annual damage by means of vulnerability function application.
To calculate these values first we add eight numerical columns to the object layer (2 for the global content and structure values and 6 for the T100, T250 and T500 content and structure damages):
v.db.addcol map=build columns="sval double, cval double"
v.db.addcol map=build columns="sar100 double, car100 double"
v.db.addcol map=build columns="sar250 double, car250 double"
v.db.addcol map=build columns="sar500 double, car500 double" 


Value of exposed elements

The value of houses can be associated to building type following market prices, while the content value can be estimated as a percentage of the structure value; the market prices used in the study are derived from the database of property selling rates available on-line at site of Agenzia del Territorio (The Italian Cadaster Agency) while the percentage values are provided by De Lotto and Testa, 1999 [4].
Type Structure value [€/m2]
Content value [% of structure value]
Residential 1700 50
Industrial 500 350
Agricultural 800 50
We use the v.db.calc command to estimate:
1) Structure value:
v.db.calc input=build field=sval exp='1700*[AREA]' \

v.db.calc input=build field=sval exp='500*[AREA]' \

v.db.calc input=build field=sval exp='800*[AREA]' \
2) Content value:
v.db.calc input=build field=cval exp='0.5*[sval]' \
where="type='residential' or type='agricultural'"

v.db.calc input=build field=cval exp='3.5*[sval]' \


Risk(T500) assessment

The expected annual damage for each return period can be assessed  applying the vulnerability functions to the exposed elements.The following commands show how to estimate it for a return period T=500.

1) Residential buildings


v.db.calc input=build field=sar500 exp='0.72*(1-math.exp(-0.1332*[h500_max])*math.exp(-0.0983))*[sval]/500' where="type='residential' and h500_max>0"

v.db.calc input=build field=car500 exp='0.7102*(1-math.exp(-0.3736*[h500_max])*math.exp(-0.0595))*[cval]/500' where="type='residential' and h500_max>0"

2) Agricultural buildings

#d%(agr,str) h<4
v.db.calc input=build field=sar500 exp='((4.35+18.20*[h500_max]+0.09*math.pow([h500_max],2)-0.94734*math.pow([h500_max],3)+0.12038*math.pow([h500_max],4))/100)*[sval]/500' where="type='agricultural' and h500_max<=4 and h500_max>0"

#d%(agr,con) h<4
v.db.calc input=build field=car500 exp='((5.532+33.164*[h500_max]+26.347*math.pow([h500_max],2)-27.765*math.pow([h500_max],3)+6.117*math.pow([h500_max],4))/100)*[cval]/500' where="type='agricultural' and h500_max<=4 and h500_max>0"

#d%(agr,str) h>4
v.db.calc input=build field=sar500 exp='0.5*[sval]/500' where="type='agricultural' and h500_max>4"

#d%(agr,con) h>4
v.db.calc input=build field=car500 exp='0.5*[cval]/500' where="type='agricultural' and h500_max>4"

3) Industrial buildings

#d%(ind,str) h<5
v.db.calc input=build field=sar500 exp='((2.157+20.354*[h500_max]-11.745*math.pow([h500_max],2)+3.507*math.pow([h500_max],3)-0.320*math.pow([h500_max],4))/100)*[sval]/500' where="type='industrial' and h500_max<=5 and h500_max>0"

#d%(ind,str) h>5
v.db.calc input=build field=sar500 exp='0.5*[sval]/500' where="type='industrial' and h500_max>5"

#d%(ind,con) h<1.2
v.db.calc input=build field=car500 exp='((-0.114+228.784*[h500_max]-235.076*math.pow([h500_max],2)+94.057*math.pow([h500_max],3))/100)*[cval]/500' where="type='industrial' and h500_max<=1.2 and h500_max>0"

#d%(ind,con) h>1.2
v.db.calc input=build field=car500 exp='[cval]*1/500' where="type='industrial' and h500_max>1.2"
The users are encouraged to finish the tutorial by applying the correct functions to the correct fields and estimating the Risk(T250) and Risk(T100).

Total Risk assessment (EAD, Expected Annual Damage)

The total risk is then the sum of the three risks, as a combination of probability.
Total Risk = Risk(T100) + Risk(T250) + Risk(T500)

 Therefore we apply in sequence the 2 commands:

v.db.addcol map=build columns="tot_risk double"
v.db.update map=build layer=1 column=tot_risk value=car100+sar100+car250+sar250+car500+sar500

[1]  Camilleri, HD 2006, 'Tsunami construction risks in the Mediterranean - outlining Malta’s scenario', Disaster Prevention  and Management, vol.15, no. 1., pp 146-162.
[2]  National Disaster management Division, Government of India 2005, 'Prevention/Protection and Mitigation from risk of Tsunami Disaster', strategy paper, viewed 14 August 2008, <>.
[3]  U.S. Army Corps of Engineers, The Hydrologic Engineering Center 1988, 'Incorporating the Effects of Mudflows into Flood Insurance Studies on Alluvial Fans', Special Projects Report 86-4, August 1988.
[4]  De Lotto P & Testa G 1999, 'Dambreak risk management and social economic impacts: a simplified method of flood damage estimation', proceedings of the 3rd CADAM workshop, Milan, May 1999.