Tutorial 2: Implementing perils in the GenMR digital template

Author: Arnaud Mignan, Mignan Risk Analytics GmbH
Version: 1.1.1
Last Updated: 2025-12-16
License: AGPL-3

Once a virtual environment has been generated (see previous tutorial), stochastic events are modelled to populate it. Table 1 lists the perils so far considered in version 1.1.1 as well as the ones to be added in the near future (v. 1.1.2).

Tab. 1. List of perils with identifiers (ID) according to Mignan (2024:tab.1.7).

ID	Peril	Peril type	Source class	Event size	Intensity class	Intensity measure	Status†
	NATURAL
AI	Asteroid impact	Primary	Point	Kinetic energy (kt)	Analytical	Overpressure (kPa)	✓
Dr	Drought						✗
EQ	Earthquake	Primary	Line	Magnitude	Analytical	Peak ground acceleration (m/s\(^2\))	✓
FF	Fluvial flood	Secondary (to RS)	Point	Peak flow	Cellular automaton	Inundation depth (m)	✓
HW	Heatwave						✗
Li	Lightning						✗
LS	Landslide	Secondary (to RS)	Diffuse	Area (km\(^2\))	Cellular automaton	Thickness (m)	✓
PI	Pest infestation						✗
RS	Rainstorm	Primary (invisible\(^*\))	Area	Rain intensity (mm/hr)	Threshold	-	✓
SS	Storm surge	Secondary (to TC)	Line	Coastal surge height (m)	Threshold	Inundation depth (m)	✓
TC	Tropical cyclone	Primary	Track	Max. wind speed (m/s)	Analytical	Max. wind speed (m/s)	✓
To	Tornado						✗
VE	Volcanic eruption	Primary	Point	Volume erupted (km\(^3\))	Analytical	Ash load (kPa)	✓
WF	Wildfire	Primary	Diffuse	Burnt area	Cellular automaton	Burnt/not burnt	✓
WS	Windstorm						✗
	TECHNOLOGICAL
BO	Blackout						✗
Ex	Explosion (industrial)	Secondary (to dg)	Point	TNT mass (kt)	Analytical	Overpressure \(P\) (kPa)	✓
	SOCIO-ECONOMIC
BI	Business interruption						✗
Sf	Public service failure						✗
SU	Social unrest						✗

\(^*\) Rainstorm is called an ‘invisible peril’ (following Mignan et al., 2014) since it does not cause any direct damage in the present framework. It is implemented as both the landslide and fluvial flood triggering mechanism.
^† Included in current version 1.1.1 (✓), planned for v.1.1.2 (✗).

This tutorial is divided into two main sections. The first introduces the definition of peril sources, an event table (containing each event’s identifier and annual rate), and a corresponding catalogue of hazard-intensity footprints. The second section covers probabilistic risk assessment, in which damage and loss footprints are generated for each event identifier (Fig. 1) and risk metrics are calculated. Before beginning these sections, the digital-template environment created in the previous tutorial is loaded.

Fig. 1. One instance of the GenMR digital template (default parameterisation) with an example of event hazard intensity footprint and loss footprint (here for a large earthquake). Simulation rendered using Rayshader (Morgan-Wall, 2022).

import numpy as np
import pandas as pd
import copy
import matplotlib.pyplot as plt

#import warnings 
#warnings.filterwarnings('ignore')   # commented, try to remove all warnings

from GenMR import perils as GenMR_perils
from GenMR import environment as GenMR_env
from GenMR import utils as GenMR_utils


# load inputs (i.e., outputs from Tutorial 1)
file_src = 'src.pkl'
file_topoLayer = 'envLayer_topo.pkl'
file_soilLayer = 'envLayer_soil.pkl'
file_urbLandLayer = 'envLayer_urbLand.pkl'
src = GenMR_utils.load_pickle2class('/io/' + file_src)
grid = copy.copy(src.grid)
topoLayer = GenMR_utils.load_pickle2class('/io/' + file_topoLayer)
soilLayer = GenMR_utils.load_pickle2class('/io/' + file_soilLayer)
urbLandLayer = GenMR_utils.load_pickle2class('/io/' + file_urbLandLayer)

GenMR_env.plot_EnvLayers([topoLayer, soilLayer, urbLandLayer], file_ext = 'jpg')

Note that the configuration of the digital template environment may be changed through Tutorial 1.

1. Stochastic event set

A list of stochastic events will be defined, consisting of both primary and secondary events. Because event interactions will be addressed in the next tutorial, only simple ad-hoc one-to-one interactions are introduced here to illustrate how secondary event footprints can be generated (e.g., a storm surge from a tropical cyclone or a landslide triggered by a rainstorm). Conditional probabilities of occurrence will be examined in Tutorial 3. In this tutorial, only primary events have assigned annual rates.

1.1. Peril source definition

Source characteristics for earthquakes, volcanic eruptions, and fluvial floods were defined in Tutorial 1, as some environmental layers depend directly on these perils. Additional peril sources are now introduced. Localised peril sources, those with no direct impact on environmental layers, are parameterised first (AI, SS, TC, and Ex). Diffuse sources (LS, WF), whose hazard-intensity footprints will be generated in Section 1.3, are characterised thereafter.

src.par['perils'] = ['AI','EQ','RS','TC','VE','WF', 'FF','LS','SS','Ex']   # IMPORTANT: primary perils come first
src.par['rdm_seed'] = 42      # None or integer (for stochastic sources: AI, TC)


# define new perils
src.par['AI'] = {'object': 'impact', 'N': 10}     # random uniform asteroid impacts (point sources)
src.par['LS'] = {'object': 'terrain'}
src.par['RS'] = {'object': 'atmosphere',          # as trigger of landslides
                 'duration': 24}                  # (hr) to calc. soil wetness that triggers LS & FF discharge
                                                  # WARNING: constant for now, to be improved

src.par['SS'] = {'object': 'coastline'}
SS_src_xi, SS_src_yi = topoLayer.coastline_coord
src.SS_char = {'x': SS_src_xi, 'y': SS_src_yi}

src.par['TC'] = {'object': 'stormtrack', 'N': 5,  # pseudo-random track sources moving from West to East
                 'npt': 10,      # number of points along track, which crosses the grid from West to East
                 'max_dev': 3,   # random uniform deviation along the y-axis at each track point
                 'bin_km': 1.}   # track spatial resolution

src.par['WF'] = {'object': 'forest'}

# WARNING: slow the first time it is called:
CI_x, CI_y = urbLandLayer.CI_refinery.centroid
src.par['Ex'] = {'object': 'refinery', 'N': 1, 'x': CI_x, 'y': CI_y}  # harbor refinery as source of explosions


GenMR_perils.plot_src(src, hillshading_z = topoLayer.z, file_ext = 'jpg')

NOTE: Hazard intensity footprints are modeled based on the sources shown in the map above. For stochastic sources (AI, TC), rerunning the previous cell will produce different configurations unless rdm_seed is fixed.

1.2. Event table generation

The stochastic event set (or event table) consists of a list of events per peril, with identifier evID, size S, and rate \(\lambda\) (lbd). Two types of perils are considered: primary perils for which the rate of occurrence must be assessed and secondary perils which will be assumed in this tutorial to occur in pair with their triggers (with rate NaN - for conditional probabilities, see the future Tutorial 3):

• Primary perils: AI, EQ, RS, TC, VE, WF.

• Secondary perils: LS and FF triggered by RS, SS triggered by TC, Ex triggered by damage.

Event size distributions are parameterised in the dictionary sizeDistr with:

• Event size incrementation according Smin, Smax and Sbin with binning in linear or log scale.

• Event size distribution modelled by one of 3 functions: power-law (powerlaw, Mignan, 2024:eq. 2.38), exponential law (exponential, Mignan, 2024:eq. 2.39), or Generalised Pareto Distribution (GPD, Mignan, 2024:eq. 2.50).

• Rates of the digital template are for most perils simply calibrated by taking the ratio of the region area to the area used to obtain event productivity parameters in the Mignan (2024) textbook (e.g., global or continental United States, CONUS). This is a simplification assuming a homogeneous spatial distribution of events but is used as a straightforward way to define rates for the present virtual region.

All calculations take place in the EventSetGenerator() class.

NOTE 1: The maximum possible event size, Smax, is defined by the user for all perils except earthquakes. For earthquakes, the maximum magnitude is constrained by the maximum fault length (src.EQ_char['srcMmax']) according to an empirical relationship (Mignan, 2024:fig. 2.6b).

NOTE 2: A floating rupture approach is used to generate earthquake ruptures on existing fault segments. See Appendix B to visualise the stochastic ruptures created by running the next cell, and refer to the Reference Manual for details on the method.

NOTE 3: The number of stochastic events is set equal to the number of sources when the sources themselves are stochastic (i.e., AI, TC). This number is kep very low for illustration purposes.

EQ_Mmax = np.max(src.EQ_char['srcMmax'])

sizeDistr = {
    'primary': ['AI','EQ','RS','TC','VE','WF'],
    'secondary': ['FF','LS','SS', 'Ex'],
    # primary perils
    'AI':{
        'Nstoch': src.par['AI']['N'], 'Sunit': 'Bolide energy $E$ (kton)',
        'Smin': 10, 'Smax': 1e3, 'Sbin': 1, 'Sscale': 'log',
        'distr': 'powerlaw', 'region': 'global', 'a0': .5677, 'b': .9       # see Mignan (2024:sect.2.3.1)
    },
    'EQ':{
        'Nstoch': 20, 'Sunit': 'Magnitude $M$',
        'Smin': 6, 'Smax': EQ_Mmax, 'Sbin': .1, 'Sscale': 'linear',
        'distr': 'exponential', 'a': 5, 'b': 1.                  # i.e., one M = a EQ per year (given b = 1)
    },
    'RS':{
        'Nstoch': 3, 'Sunit': 'Water column (mm/hr)',
        'Smin': 75, 'Smax': 125, 'Sbin': 25, 'Sscale': 'linear',
        'distr': 'GPD', 'region': 'CONUS', 'mu': 50., 'xi': -.051, 'sigma': 11.8, 'Lbdmin0': 522. 
            # see Mignan (2024:sect.2.3.2)
    },
    'TC':{
        'Nstoch': src.par['TC']['N'], 'Sunit': 'Max. wind speed (m/s)',
        'Smin': 40, 'Smax': 65, 'Sbin': 10, 'Sscale': 'linear',
        'distr': 'GPD', 'region': 'global', 'mu': 33., 'xi': -.280, 'sigma': 17.8,'Lbdmin0': 43.9  # IBTrACS fit
#        , 'Lbdmin': .1      # overwrites Lbdmin derived from Lbdmin0 to directly choose rate(>=Smin) = Lbdmin
    },
    'VE':{
        'Nstoch': 3, 'Sunit': 'Ash volume (km$^3$)',
        'Smin': 1, 'Smax': 10, 'Sbin': .5, 'Sscale': 'log',
        'distr': 'powerlaw', 'region': 'global', 'a0': -1.155, 'b': .66      # see Sect. 2.3.1
    },
    'WF':{
        'Nstoch': 10, 'Sunit': 'Burnt area (ha)',
        'Smin': 1e4, 'Smax': 1e5, 'Sbin': .5, 'Sscale': 'log',
        'distr': 'powerlaw', 'region': 'global', 'a0': 8.553, 'b': 1.23      # see Sect. 2.3.1
    },
    # secondary perils - here ad-hoc Pr(peril|trigger) = 1, see tutorial 3 for probabilistic case
    'FF':{'trigger': 'RS'},
    'LS':{'trigger': 'RS'},
    'SS':{'trigger': 'TC'},
    'Ex':{'trigger': 'dg',     # 'dg' label indicates event triggered by damage caused by any peril (tutorial 3)
          'S': 1., 'Sunit': 'Explosive yield (kton)'
    }
}

eventSet = GenMR_perils.EventSetGenerator(src, sizeDistr, GenMR_utils)
evtable, evchar = eventSet.generate()

#evtable
with pd.option_context('display.max_rows', None,
                       'display.max_columns', None):
    display(evtable)

	ID	srcID	evID	S	lbd
0	AI	impact1	AI1	10.000000	1.258426e-05
1	AI	impact2	AI2	10.000000	1.258426e-05
2	AI	impact3	AI3	10.000000	1.258426e-05
3	AI	impact4	AI4	100.000000	1.188199e-06
4	AI	impact5	AI5	100.000000	1.188199e-06
5	AI	impact6	AI6	100.000000	1.188199e-06
6	AI	impact7	AI7	100.000000	1.188199e-06
7	AI	impact8	AI8	1000.000000	1.994471e-07
8	AI	impact9	AI9	1000.000000	1.994471e-07
9	AI	impact10	AI10	1000.000000	1.994471e-07
10	EQ	fault1	EQ1	6.000000	7.692251e-03
11	EQ	fault2	EQ2	6.000000	7.692251e-03
12	EQ	fault1	EQ3	6.000000	7.692251e-03
13	EQ	fault1	EQ4	6.100000	6.110172e-03
14	EQ	fault1	EQ5	6.100000	6.110172e-03
15	EQ	fault1	EQ6	6.100000	6.110172e-03
16	EQ	fault1	EQ7	6.200000	7.280223e-03
17	EQ	fault1	EQ8	6.200000	7.280223e-03
18	EQ	fault1	EQ9	6.300000	3.855258e-03
19	EQ	fault1	EQ10	6.300000	3.855258e-03
20	EQ	fault1	EQ11	6.300000	3.855258e-03
21	EQ	fault1	EQ12	6.400000	4.593510e-03
22	EQ	fault1	EQ13	6.400000	4.593510e-03
23	EQ	fault2	EQ14	6.500000	3.648755e-03
24	EQ	fault1	EQ15	6.500000	3.648755e-03
25	EQ	fault1	EQ16	6.600000	5.796618e-03
26	EQ	fault1	EQ17	6.700000	4.604417e-03
27	EQ	fault1	EQ18	6.800000	3.657419e-03
28	EQ	fault1	EQ19	6.900000	2.905191e-03
29	EQ	fault1	EQ20	7.000000	2.307675e-03
30	RS	atmosphere	RS1	75.000000	5.249839e-03
31	RS	atmosphere	RS2	100.000000	5.006848e-04
32	RS	atmosphere	RS3	125.000000	3.425633e-05
33	TC	stormtrack1	TC1	40.000000	6.049596e-04
34	TC	stormtrack2	TC2	50.000000	3.677710e-04
35	TC	stormtrack3	TC3	60.000000	9.923821e-05
36	TC	stormtrack4	TC4	60.000000	9.923821e-05
37	TC	stormtrack5	TC5	70.000000	8.826349e-05
38	VE	volcano1	VE1	1.000000	1.793957e-06
39	VE	volcano1	VE2	3.162278	8.390967e-07
40	VE	volcano1	VE3	10.000000	3.924750e-07
41	WF	forest	WF1	10000.000000	3.625127e-02
42	WF	forest	WF2	10000.000000	3.625127e-02
43	WF	forest	WF3	10000.000000	3.625127e-02
44	WF	forest	WF4	10000.000000	3.625127e-02
45	WF	forest	WF5	10000.000000	3.625127e-02
46	WF	forest	WF6	10000.000000	3.625127e-02
47	WF	forest	WF7	31622.776602	1.759354e-02
48	WF	forest	WF8	31622.776602	1.759354e-02
49	WF	forest	WF9	31622.776602	1.759354e-02
50	WF	forest	WF10	100000.000000	1.280780e-02
51	FF	river1	FF_fromRS1	NaN	NaN
52	FF	river1	FF_fromRS2	NaN	NaN
53	FF	river1	FF_fromRS3	NaN	NaN
54	LS	terrain	LS_fromRS1	NaN	NaN
55	LS	terrain	LS_fromRS2	NaN	NaN
56	LS	terrain	LS_fromRS3	NaN	NaN
57	SS	coastline	SS_fromTC1	NaN	NaN
58	SS	coastline	SS_fromTC2	NaN	NaN
59	SS	coastline	SS_fromTC3	NaN	NaN
60	SS	coastline	SS_fromTC4	NaN	NaN
61	SS	coastline	SS_fromTC5	NaN	NaN
62	Ex	refinery	Ex_fromCIf	1.000000	NaN

1.3. Intensity footprint catalogue generation

Once the event table has been generated, the next step is to produce one mean hazard-intensity footprint per event. Aleatory uncertainty is not implemented for simplicity, and because it does not affect the objectives of this project, namely, modelling compound events in a virtual environment (see the third and final tutorial).

1.3.1. From analytical expressions & threshold models

Analytical expressions are available to calculate intensity \(I(x,y) = f_I(S,r)\) for the AI/Ex, EQ, VE and TC perils (Mignan, 2024:sec. 2.4.1). The size \(S\) is provided in evtable, while the distance-to-source \(r\) is obtained by matching each event to its corresponding source coordinates via evchar. Because generating TC footprints requires a multi-step procedure, additional details are provided in Appendix C.

The models are implemented as calc_I_*(S,r) functions in the Python file, using the equations from Mignan (2024:sec. 2.4.1):

ID	Event size \(S\)	Intensity \(I\)	Equation
AI, Ex	Explosive mass \(m_{TNT}\) (kt)	Blast overpressure \(P\) (kPa)	Eq. 2.56
EQ	Magnitude \(M\)	Peak ground acceleration PGA (m/s\(^2\))	Eq. 2.54
TC	Maximum wind speed \(v_{max}\) (m/s)	Wind speed \(v\) (m/s)	Eq. 2.24
VE	Erupted volume \(V\) (km\(^3\))	Ash load \(P\) (kPa), converted from thickness (m)	Eq. 2.61

Storm surge (SS) is the only peril in this tutorial for which hazard intensity is estimated using a threshold-based approach (Bathtub model; Mignan, 2024:sec. 2.4.2). The algorithm is intentionally simple and tailored to the virtual environment used here, where the background topography exhibits a predominantly west-east gradient (see the model_SS_Bathtub() function). As a secondary peril triggered by TC events, the SS size along the coastline is derived from a relationship of the form \(h_{max} \propto v_{max}^2\) (Mignan, 2024:fig. 2.7), where \(v_{max}\) is the wind speed along the coastline obtained from the TC footprint.

The catalogue of hazard footprints is generated via the HazardFootprintGenerator(() class and is visualised with plot_hazFootprints().

# additional atmospheric source parameters -> To move to atmospheric layer in v.1.1.2
src.par['TC']['vforward_m/s'] = 10          # forward storm speed [m/s]    
src.par['TC']['pn_mbar'] = 1005             # ambient pressure [mbar]
src.par['TC']['B_Holland'] = 2              # Holland's parameter (see Fig. 2.10 of textbook)
src.par['TC']['lat_deg'] = grid.lat_deg     # average latitude of digital template (see Tutorial 1)

# virtual region bathymetry characteristics for SS height given storm wind speed
src.par['SS']['bathy'] = 'New York harbor'  # 'generic' (Saffir-Simpson, Mignan, 2024:eq. 2.16) 
                                            # 'New York harbor' (eq. in Mignan, 2024:fig. 2.7)

hazFootprints = GenMR_perils.HazardFootprintGenerator(evtable, evchar, src, topoLayer.z)
catalog_hazFootprints = hazFootprints.generate()

generating footprints for: AI, EQ, RS, TC, VE, WF, FF, LS, SS, Ex, ... catalogue completed

# max intensity Imax in color range (see units in previous table)
# to compare similar levels of expected damage, choose a value per peril for a constant MDR in the vulnerability
# curves (see curves generated by plot_vulnFunctions())
Imax_Ex_totaldestruction = 4 * 6.89476  # kPa (here 4 psi as threshold of total destruction)

plot_Imax = {'AI': Imax_Ex_totaldestruction, 'EQ': 6, 'Ex': Imax_Ex_totaldestruction, 'LS': 2., \
             'SS': 2., 'TC': 70., 'VE': .5}

# nstoch (default = 5) defines the number of events per peril to be displayed by row
# if the number of events > nstoch, footprints are randomly sampled from all possibilities
GenMR_perils.plot_hazFootprints(catalog_hazFootprints, grid, topoLayer.z, plot_Imax)  #, nstoch = 5)

Note: Observe the TC footprint boundary effect (cyclonic eye), which is due to the track artificially starting at grid.xmin.

1.3.2. From numerical models (cellular automata)

Generation of hazard-intensity footprints via dynamical processes (here, cellular automata - CA) is handled separately because these computations are more expensive. This functionality is implemented in the DynamicHazardFootprintGenerator() class and applied to the LS, FF and WF perils. Details of the corresponding CA models can be found in the Python file within the CellularAutomaton_* classes, where * denotes the peril ID.

The following movie illustrates an example of dynamical process. It is produced by running LS_CA.write_gif() in the next cell and shows the time-evolving landslide propagation generated by the CA. The final hazard-intensity footprint included in the hazard-footprint catalogue is static and defined as \(\max_t I(x,y,t)\).

reproduce_LSmovie = False          # True to create 'movs/LS_CA_xrg38_68_yrg20_35.gif'
reproduce_FFmovie = False          # True to create 'movs/FF_CA_xrg0_100_yrg0_40.gif'

if reproduce_LSmovie:
    movie = {'create': True, 'xmin': 38, 'xmax': 68, 'ymin': 20, 'ymax': 35}   # zoom on specific hill

    evID_trigger = 'RS3'                                       # case considered for movie scenario
    S_trigger = evtable['S'][evtable['evID'] == evID_trigger].values
    hw = S_trigger * 1e-3 * src.par['RS']['duration']          # water column (m)
    wetness = hw / soilLayer.h
    wetness[wetness > 1] = 1                                   # max possible saturation
    wetness[soilLayer.h == 0] = 0                              # no soil case

    LS_CA = GenMR_perils.CellularAutomaton_LS(soilLayer, wetness, movie)
    LS_CA.run()
    LS_CA.write_gif()
    print('LS movie created (go to mov/LS_CA_*).')
else:
    print('No LS movie created.')
    
if reproduce_FFmovie:
    I_RS = evtable[evtable['evID'] == 'RS3']['S'].values[0] * 1e-3 / 3600    # (mm/hr) to (m/s)
    movie = {'create': True, 'xmin': 0, 'xmax': 100, 'ymin': 0, 'ymax': 40, 'tmin': 40} 

    FF_CA = GenMR_perils.CellularAutomaton_FF(I_RS, src, grid, topoLayer.z, movie)
    FF_CA.run()
    FF_CA.write_gif()
    print('FF movie created (go to mov/FF_CA_*).')
else:
    print('No FF movie created.')

No LS movie created.
No FF movie created.

Footprints are initialised in dynhazFootprints, an instance of the DynamicHazardFootprintGenerator() class. Running dynhazFootprints.generate() computes all dynamical intensity footprints at once. To generate only a specific subset of perils, one must use dynhazFootprints.generate(selected_perils = ['LS', 'FF', ...]). Note that if a peril has already been run, cached footprints will be loaded automatically. To force recomputation, one needs to include the option force_recompute = True in the class command (or delete the folder *_CA_frames).

# additional WF source/landuse parameters
src.par['WF'] = {'ratio_grass': .6,        # ratio of forest cells to turn to grass at initial state
             'p_lightning': .1,            # probability of lightning strike per iteration
             'rate_newtrees': 1000,        # number of new trees to add per iteration
             'nsim': 100000,               # number of iterations
             'Smin_plot': 1000,            # minimum WF footprint size to save in 'figs/WF_CA_frames'
            }

dynhazFootprints = GenMR_perils.DynamicHazardFootprintGenerator(evtable, src, soilLayer, urbLandLayer)
catalog_dynhazFootprints = dynhazFootprints.generate()
#_ = dynhazFootprints.generate(selected_perils = ['LS'], force_recompute = True)
#catalog_dynhazFootprints = dynhazFootprints.generate(selected_perils = ['FF'], force_recompute = True)
catalog_dynhazFootprints.keys()

generating footprints for:
Computing potential footprints
0 / 100000
1000 / 100000
2000 / 100000
3000 / 100000
4000 / 100000
5000 / 100000
6000 / 100000
7000 / 100000
8000 / 100000
9000 / 100000
10000 / 100000
11000 / 100000
12000 / 100000
13000 / 100000
14000 / 100000
15000 / 100000
16000 / 100000
17000 / 100000
18000 / 100000
19000 / 100000
20000 / 100000
21000 / 100000
22000 / 100000
23000 / 100000
24000 / 100000
25000 / 100000
26000 / 100000
27000 / 100000
28000 / 100000
29000 / 100000
30000 / 100000
31000 / 100000
32000 / 100000
33000 / 100000
34000 / 100000
35000 / 100000
36000 / 100000
37000 / 100000
38000 / 100000
39000 / 100000
40000 / 100000
41000 / 100000
42000 / 100000
43000 / 100000
44000 / 100000
45000 / 100000
46000 / 100000
47000 / 100000
48000 / 100000
49000 / 100000
50000 / 100000
51000 / 100000
52000 / 100000
53000 / 100000
54000 / 100000
55000 / 100000
56000 / 100000
57000 / 100000
58000 / 100000
59000 / 100000
60000 / 100000
61000 / 100000
62000 / 100000
63000 / 100000
64000 / 100000
65000 / 100000
66000 / 100000
67000 / 100000
68000 / 100000
69000 / 100000
70000 / 100000
71000 / 100000
72000 / 100000
73000 / 100000
74000 / 100000
75000 / 100000
76000 / 100000
77000 / 100000
78000 / 100000
79000 / 100000
80000 / 100000
81000 / 100000
82000 / 100000
83000 / 100000
84000 / 100000
85000 / 100000
86000 / 100000
87000 / 100000
88000 / 100000
89000 / 100000
90000 / 100000
91000 / 100000
92000 / 100000
93000 / 100000
94000 / 100000
95000 / 100000
96000 / 100000
97000 / 100000
98000 / 100000
99000 / 100000
Fetching footprints from potential footprints
WF1 (fetched from cache)
WF2 (fetched from cache)
WF3 (fetched from cache)
WF4 (fetched from cache)
WF5 (fetched from cache)
WF6 (fetched from cache)
WF7 (fetched from cache)
WF8 (fetched from cache)
WF9 (fetched from cache)
WF10 (fetched from cache)
FF_fromRS1 (loaded from cache)
FF_fromRS2 (loaded from cache)
FF_fromRS3 (loaded from cache)
LS_fromRS1 (loaded from cache)
LS_fromRS2 (loaded from cache)
LS_fromRS3 (loaded from cache)
... catalogue completed

dict_keys(['WF1', 'WF2', 'WF3', 'WF4', 'WF5', 'WF6', 'WF7', 'WF8', 'WF9', 'WF10', 'FF_fromRS1', 'FF_fromRS2', 'FF_fromRS3', 'LS_fromRS1', 'LS_fromRS2', 'LS_fromRS3'])

plot_Imax.update({'LS': 2, 'FF': 1, 'WF': 1})

catalog_dynhazFootprints = dynhazFootprints.catalog_hazFootprints
GenMR_perils.plot_hazFootprints(catalog_dynhazFootprints, grid, topoLayer.z, plot_Imax)

# merge the two catalogues
catalog_allhazFootprints = catalog_hazFootprints | catalog_dynhazFootprints

2. Probabilistic risk assessment

In this final section, loss footprints are modelled, and total losses per event are added as new column in the event table, creating an event loss table (ELT). Risk metrics are then calculated from the ELT (considering only primary events in this tutorial). Year loss tables will be addressed in the future Tutorial 3.

2.1. Vulnerability curve definition

Generic intensity-specific vulnerability functions, \(D = f_D(I)\), are implemented and apply to any building type for simplicity. The vulnerability functions (f_D()) are listed below and correspond to the equations described in Section 3.2.5 of Mignan (2024). Damage is expressed in terms of the mean damage ratio (MDR):

ID	Intensity \(I\)	Vulnerability function	Parameters
AI, Ex	Blast overpressure \(P\) [kPa]	Cum. lognormal distr. (Eq. 3.13)	\(\mu\) = log(20), \(\sigma\) = 0.4
EQ	Peak ground acceleration PGA [m/s\(^2\)]	Cum. lognormal distr. (Eq. 3.13)	\(\mu\) = log(6), \(\sigma\) = 0.6
FF, SS	Inundation depth \(h\) [m]	Square root (Eq. 3.10)	\(c\) = 0.45
LS	Landslide thickness \(h\) [m]	Weibull distr. (Eq. 3.12)	\(c_1\) = -1.671, \(c_2\) = 3.189, \(c_3\) = 1.746
VE	Ash load \(P\) [kPa] (converted from thickness [m])	Cum. lognormal distr. (Eq. 3.13)	\(\mu\) = 1.6, \(\sigma\) = 0.4
TC, To, WS	Maximum wind speed \(v_{max}\) [m/s]	Power-law (Eq. 3.14)	\(v_{thres}\) = 25.7, \(v_{half}\) = 74.7
WF	Burnt yes/no	Heaviside step function	MDR = 1 if burnt, 0 otherwise

GenMR_perils.plot_vulnFunctions()

2.2. Loss footprint catalogue generation

Damage footprints, \(D(x,y)\), are obtained by mapping hazard-intensity footprints, \(I(x,y)\), into mean damage ratio (MDR) values using the vulnerability functions defined earlier. Loss footprints, \(L(x,y)\), are calculated as the product of \(D(x,y)\) and the exposure footprint, \(\nu(x,y)\).

It should be noted that losses are computed for independent events occurring individually. For example, the total loss from a TC event could depend on losses caused by the associated SS event. Such interactions, whether at the exposure or damage level, will be addressed in Tutorial 3.

riskFootprints_all = GenMR_perils.RiskFootprintGenerator(catalog_allhazFootprints, urbLandLayer, evtable)
catalog_dmgFootprints = riskFootprints_all.catalog_dmgFootprints
catalog_lossFootprints = riskFootprints_all.catalog_lossFootprints
riskFootprints_all.ELT

Computing MDR & Loss: 100%|█████████████████████| 60/60 [00:01<00:00, 47.64it/s]

	ID	srcID	evID	S	lbd	loss
0	AI	impact1	AI1	10.0	0.000013	7.865502e+04
1	AI	impact2	AI2	10.0	0.000013	0.000000e+00
2	AI	impact3	AI3	10.0	0.000013	0.000000e+00
3	AI	impact4	AI4	100.0	0.000001	3.876696e+05
4	AI	impact5	AI5	100.0	0.000001	2.741698e+05
...	...	...	...	...	...	...
58	SS	coastline	SS_fromTC2	NaN	NaN	1.082570e+10
59	SS	coastline	SS_fromTC3	NaN	NaN	2.141890e+10
60	SS	coastline	SS_fromTC4	NaN	NaN	5.520996e+10
61	SS	coastline	SS_fromTC5	NaN	NaN	7.956129e+10
62	Ex	refinery	Ex_fromCIf	1.0	NaN	1.345011e+09

63 rows × 6 columns

nfp = len(catalog_lossFootprints)
evIDs = list(catalog_lossFootprints.keys())
_, ax = plt.subplots(nfp,3, figsize=(20,5*len(riskFootprints_all.ELT)))
for i in range(nfp):
    ax[i,0].contourf(grid.xx, grid.yy, urbLandLayer.bldg_value, cmap = 'Greens')
    ax[i,0].contourf(grid.xx, grid.yy, GenMR_env.ls.hillshade(topoLayer.z, vert_exag=.1), cmap='gray', alpha = .1)
    ax[i,0].set_xlabel('$x$ (km)')
    ax[i,0].set_ylabel('$y$ (km)')
    ax[i,0].set_title(evIDs[i])
    ax[i,0].set_aspect(1)

    ax[i,1].contourf(grid.xx, grid.yy, catalog_dmgFootprints[evIDs[i]], cmap = 'Reds', vmin = 0, vmax = 1)
    ax[i,1].contourf(grid.xx, grid.yy, GenMR_env.ls.hillshade(topoLayer.z, vert_exag=.1), cmap='gray', alpha = .1)
    ax[i,1].set_xlabel('$x$ (km)')
    ax[i,1].set_ylabel('$y$ (km)')
    ax[i,1].set_aspect(1)

    ax[i,2].contourf(grid.xx, grid.yy, catalog_lossFootprints[evIDs[i]], cmap = 'Reds')
    ax[i,2].contourf(grid.xx, grid.yy, GenMR_env.ls.hillshade(topoLayer.z, vert_exag=.1), cmap='gray', alpha = .1)
    ax[i,2].set_xlabel('$x$ (km)')
    ax[i,2].set_ylabel('$y$ (km)')
    ax[i,2].set_aspect(1)

plt.tight_layout();

# NB: adapt event size considered and vulnerability to focus on highly damaging events only (illustrative cases)

2.3. Risk metrics

Several risk metrics are derived from the ELT, including the annual average loss (AAL), exceeding probability (EP) curves, value at risk (VaR) at a given quantile \(q\), and the tail value at risk (TVaR):

q_VAR = .99
ELT = riskFootprints_all.ELT
ELT['L'] = ELT['loss']
ELT = ELT[~np.isnan(ELT['lbd'])]   # secondary events (no rate in Tutorial 2 removed)
ELT_wEP, AAL, VaRq, TVaRq, _, _ = GenMR_utils.calc_riskmetrics_fromELT(ELT, q_VAR)

print('** from ELT **')
print('AAL [$ M.]:', AAL * 1e-6)
print('VaR_q [$ B.]:', VaRq * 1e-9)
print('TVaR_q [$ B.]:', TVaRq * 1e-9)  # much higher than VaR due to very long risk tail (default parameterisation)

plt.rcParams['font.size'] = '20'
fig, ax = plt.subplots(1,1, figsize = (20,8))
plt.scatter(ELT_wEP['L'] * 1e-9, ELT_wEP['EP'], color = 'black')
plt.plot(ELT_wEP['L'] * 1e-9, ELT_wEP['EP'], color = 'darkred')
plt.axhline(y = 1 - q_VAR, linestyle = 'dotted', color = 'black')
plt.axvline(x = VaRq * 1e-9, linestyle = 'dotted', color = 'darkred')
#plt.xlim(0,20)
#plt.ylim(0,.1)
plt.xlabel('Economic losses [$ B.]')
plt.ylabel('Exceedance probability')
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False);

** from ELT **
AAL [$ M.]: 2678.5205716147916
VaR_q [$ B.]: 67.89485712334694
TVaR_q [$ B.]: 287.23171766165484

References

• Mignan A (2024), Introduction to Catastrophe Risk Modelling - A Physics-based Approach. Cambridge University Press, doi: 10.1017/9781009437370

Appendix

This appendix provides additional details that are kept outside the main notebook to maintain the clarity and flow of the main content.

A. Peril sources

Tab. A1. Full list of attributes and properties in source class src.

Attribute	ID	Description
Parameters	par	Input source parameters
Spatial grid	grid	copy of RasterGrid class
Asteroid impact source characteristics	AI_char	Stochastic point-source coordinates
Earthquake source characteristics	EQ_char	Fault source and segment coordinates, magnitudes and geometries
Fluvial flood source characteristics	FF_char	River bed coordinates
Storm surge source characteristics	SS_char	Coastline coordinates
Tropical cyclone source characteristics	TC_char	Stochastic track-source coordinates
Volcanic eruption source characteristics	VE_char	Point-source coordinates

B. Stochastic event set

B.1. Earthquake case

The coordinates of stochastic events (defined in evchar) generally match those of their sources for localised peril sources (for diffuse sources, cellular automata are applied to determine the spatial characteristics of each event). In the case of earthquakes however, a floating rupture approach is used to generate earthquake ruptures on existing fault segments. The next cell illustrates how to extract the relevant information from evchar and how to examine the resulting spatial configuration in a stochastic run of EventSetGenerator().

EQcoord = evchar[GenMR_perils.get_peril_evID(evchar['evID']) == 'EQ'].reset_index(drop = True)
EQi = evtable[evtable['ID'] == 'EQ']['evID'].values

plt.rcParams['font.size'] = '18'
fig, ax = plt.subplots(4, 5, figsize=(20,16))
ax = ax.flatten()
for i in range(sizeDistr['EQ']['Nstoch']):
    EQcoord_evID = EQcoord[EQcoord['evID'] == EQi[i]]
    ax[i].contourf(grid.xx, grid.yy, GenMR_env.ls.hillshade(topoLayer.z, vert_exag=.1), cmap='gray', alpha = .1)
    ax[i].plot(EQcoord_evID['x'], EQcoord_evID['y'], color = 'darkred', linewidth = 2)
    ax[i].set_xlim(grid.xmin, grid.xmax)
    ax[i].set_ylim(grid.ymin, grid.ymax)
    ax[i].set_xlabel('$x$ (km)')
    ax[i].set_ylabel('$y$ (km)')
    ax[i].set_title(EQi[i], pad = 10)
    ax[i].set_aspect(1)
plt.tight_layout()
plt.pause(1)   # to avoid rare case where plot not displayed on notebook
plt.show()

C. Hazard footprints

C.1. Tropical cyclone case

While modelling the AI, EQ, and VE intensity footprints from analytical expressions (Sect. 1.3.1) is straightforward, this is not the case for the TC peril, which involves a time-dependent intensity field. First, because the footprint depends on a moving track, a hazard-intensity field must be computed at each time increment, and the final static footprint is then defined as \(\max_t I(x,y,t)\). Second, the windfield is composed of several components, a tangential wind vector and the storm’s forward-motion vector, that must be combined to obtain the total wind amplitude. For illustration, intensity footprints at a given instant \(t\) are displayed below, together with their associated windfields, shown both without and with the effect of storm forward motion.

evID = 'TC5'     # choose a TC event
loc_x = 25.      # location of time snapshot (km)


track_evID = evchar[evchar['evID'] == evID]
t_snapshot = np.where(track_evID['x'] > loc_x)[0][0]

I_sym_t, I_asym_t, vtot_x, vtot_y, vtan_x, vtan_y = hazFootprints.get_TC_timeshot(evID, t_snapshot)


plt.rcParams['font.size'] = '18'
fig, ax = plt.subplots(1, 2, figsize=(20,10))

Imax = plot_Imax['TC']
I_plt = np.copy(I_sym_t)
I_plt[I_plt >= Imax] = Imax
ax[0].contourf(grid.xx, grid.yy, I_plt, cmap = 'Reds', levels = np.linspace(0, Imax, 100), vmin=33, vmax=Imax)
ax[0].plot(track_evID['x'], track_evID['y'], color = GenMR_utils.col_peril('TC'))
ax[0].plot(src.SS_char['x'], src.SS_char['y'], color = 'black', linestyle = 'dashed')
ax[0].scatter(track_evID['x'].values[t_snapshot], track_evID['y'].values[t_snapshot], marker = '+', \
              color = 'black', s = 200)
delta = 50
for i in range(int(grid.xx.shape[0]/delta)):
    for j in range(int(grid.xx.shape[1]/delta)):
        ax[0].arrow(grid.xx[delta*i,delta*j], grid.yy[delta*i,delta*j], 
            .05 * vtan_x[delta*i,delta*j], .05* vtan_y[delta*i,delta*j], head_width = 1, color = 'white')
ax[0].set_xlabel('$x$ (km)')
ax[0].set_ylabel('$y$ (km)')
ax[0].set_title('Static storm', pad = 20)
ax[0].set_xlim(grid.xmin, grid.xmax)
ax[0].set_ylim(grid.ymin, grid.ymax)
ax[0].set_aspect(1)

I_plt = np.copy(I_asym_t)
I_plt[I_plt >= Imax] = Imax
ax[1].contourf(grid.xx, grid.yy, I_plt, cmap = 'Reds', levels = np.linspace(0, Imax, 100), vmin=33, vmax=Imax)
ax[1].plot(track_evID['x'], track_evID['y'], color = GenMR_utils.col_peril('TC'))
ax[1].plot(src.SS_char['x'], src.SS_char['y'], color = 'black', linestyle = 'dashed')
ax[1].scatter(track_evID['x'].values[t_snapshot], track_evID['y'].values[t_snapshot], marker = '+', \
              color = 'black', s = 200)
delta = 50
for i in range(int(grid.xx.shape[0]/delta)):
    for j in range(int(grid.xx.shape[1]/delta)):
        ax[1].arrow(grid.xx[delta*i,delta*j], grid.yy[delta*i,delta*j], 
            .05 * vtot_x[delta*i,delta*j], .05* vtot_y[delta*i,delta*j], head_width = 1, color = 'white')
ax[1].set_xlabel('$x$ (km)')
ax[1].set_ylabel('$y$ (km)')
ax[1].set_title('Storm in motion', pad = 20)
ax[1].set_xlim(grid.xmin, grid.xmax)
ax[1].set_ylim(grid.ymin, grid.ymax)
ax[1].set_aspect(1)
plt.tight_layout();