Commit 805a7268 authored by Danijel Schorlemmer's avatar Danijel Schorlemmer Committed by Tara Evaz Zadeh
Browse files

Added the ability to deal with continuous fragility functions

parent 36396e58
Pipeline #30309 passed with stage
in 4 minutes and 4 seconds
import csv
import numpy as np
import xml.etree.ElementTree as ET
from scipy.stats import lognorm
class FragilityModel:
......@@ -9,24 +10,88 @@ class FragilityModel:
for different structural types
"""
def __init__(self, fragility_functions_filepath, intensity_measure_map_filepath):
def __init__(
self,
fragility_functions_filepath,
gmf_filepath,
intensity_measure_map_filepath=None,
taxonomy_map_filepath=None,
):
"""
Initializing the fragility model by importing the XML file with fragility functions
and creating the intensity-measure map that maps intensity types from the
ground-motion field to intensity types as given in the fragility functions.
and creating the intensity-measure map that maps intensity-measure types as given in
the fragility functions to the ones in the ground-motion field file. It then creates
a mapping dictionary between the intensity measure types and their respective
column number in the ground-motion field file for faster access to data.
This function is limited to one fragility function per building taxonomy.
This function currently only supports `discrete` fragility functions.
Another limitation is that for a building taxonomy only one fragility function
can be used.
Dictionary structure for discrete functions:
'CR/LFINF+CDL/HBET:7-20/0.0':
{'imv': array([2.225000e-03, 2.397000e-03, 2.581000e-03]),
'slight': array([0.00000e+00, 0.00000e+00, 0.00000e+00]),
'moderate': array([0.00000e+00, 0.00000e+00, 0.00000e+00]),
'extensive': array([0.00000e+00, 0.00000e+00, 0.00000e+00]),
'complete': array([0.00000e+00, 0.00000e+00, 0.00000e+00]),
'format': 'discrete',
'imt': 'SA(1)'}
Dictionary structure for a continuous function:
'CR/LFINF+CDL/HBET:7-20/0.0':
{'slight': {'mean': "1.692", 'stddev': "1.179"},
'moderate': {'mean': "1.692", 'stddev': "1.179"},
'extensive': {'mean': "1.692", 'stddev': "1.179"},
'complete': {'mean': "1.692", 'stddev': "1.179"},
'format': 'continuous',
'imt': 'SA(1)',
'maxIML': "5",
'minIML': "1e-10",
'noDamageLimit': "0.05",
'function': function object}
Mapping dictionary example:
{'PGA': 2, 'SA(0.3)': 3, 'SA(0.6)': 4, 'SA(1)': 5}
Args:
fragility_functions_filepath (str):
File path to the XML file containing fragility functions in NRML as
provided by the Global Earthquake Model.
gmf_filepath (str):
File path to the ground-motion field.
(Header line contains the intensity-measure names)
File example extract:
lon,lat,gmv_PGA,gmv_SA(0.3),gmv_SA(0.6),gmv_SA(1.0)
22.89028,36.31806,4.822530E-03,1.255729E-02,8.736364E-03,5.612176E-03
22.89028,36.32083,4.830462E-03,1.257595E-02,8.748890E-03,5.619673E-03
...
intensity_measure_map_filepath (str):
File path to a CSV file mapping the intensity types from the
ground-motion field to intensity types as given in the fragility functions.
If `None`, a default mapping will be applied.
File example extract:
in_fragility_function, in_ground_motion_field
SA(0.3),gmv_SA(0.3)
SA(0.6),gmv_SA(0.6)
SA(1),gmv_SA(1.0)
SA(1.0),gmv_SA(1.0)
PGA,gmv_PGA
taxonomy_map_filepath (str):
File path to the mapping of taxonomy strings to their relevant
fragility function names.
File example extract:
taxonomy,conversion
CR+PC/LFM+CDN/H:1,CR_LFM-CDN_H1_0
CR/LFM+CDN/H:1,CR_LFM-CDN_H1_0
...
"""
tree = ET.parse(fragility_functions_filepath)
......@@ -49,23 +114,73 @@ class FragilityModel:
self.num_limit_states = len(self.limit_states)
elif "fragilityFunction" in child.tag:
limit_states = []
# Intensity-measure values and probabilities of exceedance
imv_poe = {}
# Intensity-measure type
im_type = None
for ichild in list(child):
key, value = ichild.items()[0]
if key == "imt":
im_type = value
im_value = np.array(list(map(float, ichild.text.split())))
imv_poe["imv"] = im_value
if key == "ls":
ls = value
values = np.array(list(map(float, ichild.text.split())))
imv_poe[ls] = values
limit_states.append(value)
imv_poe["format"] = child.get("format")
imv_poe["imt"] = im_type
parameters = {}
# Import a discrete function
if child.get("format") == "discrete":
# Intensity-measure values and probabilities of exceedance
for ichild in list(child):
key, value = ichild.items()[0]
if key == "imt":
im_type = value
im_value = np.array(list(map(float, ichild.text.split())))
parameters["imv"] = im_value
if key == "ls":
ls = value
values = np.array(list(map(float, ichild.text.split())))
parameters[ls] = values
limit_states.append(value)
parameters["format"] = child.get("format")
parameters["imt"] = im_type
if im_type is None:
raise ValueError(
"No intensity measure given for fragility function %s"
% (child.get("id"))
)
if not (limit_states == self.limit_states):
raise ValueError(
"Limit states for fragility function"
+ "%s differ from global limit states %s"
% (child.get("id"), self.limit_states)
)
# Import a continuous function
elif child.get("format") == "continuous":
parameters["format"] = child.get("format")
if child.get("shape") == "logncdf":
for ichild in list(child):
if "imls" in ichild.tag:
for key, value in ichild.items():
parameters[key] = value
if key == "imt":
im_type = value
elif "params" in ichild.tag:
ls_params = {}
limit_state = None
for key, value in ichild.items():
if key == "ls":
limit_state = value
limit_states.append(value)
else:
ls_params[key] = float(value)
variance = ls_params["stddev"] ** 2.0
sigma = np.sqrt(
np.log((variance / (ls_params["mean"] ** 2.0)) + 1.0)
)
mu = (ls_params["mean"] ** 2.0) / np.sqrt(
variance + ls_params["mean"] ** 2.0
)
ls_params["function"] = lognorm(sigma, scale=mu)
parameters[limit_state] = ls_params
else:
raise ValueError(
"Unknown function form %s of fragility function %s"
% (child.get("shape"), child.get("id"))
)
# Final checks to make sure the fragility function is valid
if im_type is None:
raise ValueError(
"No intensity measure given for fragility function %s"
......@@ -77,7 +192,7 @@ class FragilityModel:
+ "%s differ from global limit states %s"
% (child.get("id"), self.limit_states)
)
self.fragility_functions[child.get("id")] = imv_poe
self.fragility_functions[child.get("id")] = parameters
# Read the intensity measure map into a dictionary
if intensity_measure_map_filepath is not None:
......@@ -87,12 +202,34 @@ class FragilityModel:
self.intensity_measure_map = {rows[0]: rows[1] for rows in reader}
else:
self.intensity_measure_map = {
"gmv_SA(0.3)": "SA(0.3)",
"gmv_SA(0.6)": "SA(0.6)",
"gmv_SA(1.0)": "SA(1)",
"gmv_PGA": "PGA",
"SA(0.3)": "gmv_SA(0.3)",
"SA(0.6)": "gmv_SA(0.6)",
"SA(1.0)": "gmv_SA(1.0)",
"SA(1)": "gmv_SA(1.0)",
"PGA": "gmv_PGA",
"MMI": "gmv_MMI",
}
# Read the taxonomy map if filename provided
if taxonomy_map_filepath is not None:
with open(taxonomy_map_filepath, mode="r") as input_file:
reader = csv.reader(input_file)
next(reader)
self.taxonomy_map = {rows[0]: rows[1] for rows in reader}
# Extracting the ground-motion types available in the ground-motion field.
with open(gmf_filepath) as gmf:
im_types = gmf.readline().strip().split(",")[2:]
im_types_dict = {}
for index, value in enumerate(im_types):
im_types_dict[value] = index + 2
# Add for each intensity measure from the fragility functions the respective
# column number to the dictionary. Hereby, using the intensity-measure map to
# translate between the name in the fragility functions and the name in
# ground-motion field.
for im_ff in self.intensity_measure_map:
self.gm_type_column_map[im_ff] = im_types_dict[self.intensity_measure_map[im_ff]]
def __repr__(self):
# Informative string representation
return "Fragility Model: %s\nLimit States: %s" % (
......@@ -104,40 +241,6 @@ class FragilityModel:
# Return the number of fragility functions in the model
return len(self.fragility_functions)
def get_imt(self, taxonomy):
return self.fragility_functions[taxonomy]
def create_gmf_column_map(self, gmf_filepath):
"""
Creates a mapping dictionary between the intensity measure types and their respective
column number in the ground-motion field file.
Mapping dictionary example:
{'PGA': 2, 'SA(0.3)': 3, 'SA(0.6)': 4, 'SA(1)': 5}
Args:
gmf_filepath (str):
File path to the ground-motion field
(Header line contains the intensity-measure names)
File example extract:
lon,lat,gmv_PGA,gmv_SA(0.3),gmv_SA(0.6),gmv_SA(1.0)
22.89028,36.31806,4.822530E-03,1.255729E-02,8.736364E-03,5.612176E-03
22.89028,36.32083,4.830462E-03,1.257595E-02,8.748890E-03,5.619673E-03
...
"""
# Extracting the ground-motion types available in the ground-motion field.
with open(gmf_filepath) as gmf:
im_types = gmf.readline().strip().split(",")[2:]
# Add for each intensity measure from the ground-motion field the respective
# column number to the dictionary. Hereby, using the intensity-measure map to
# translate between the name in the ground-motion field and the name in
# the fragility functions.
for i in range(len(im_types)):
self.gm_type_column_map[self.intensity_measure_map[im_types[i]]] = i + 2
def __call__(self, taxonomy, intensity_measure_values):
"""
Computing PoE values for all the different damage states available in the
......@@ -174,19 +277,21 @@ class FragilityModel:
Intensity-measure type used to compute the PoEs.
"""
# Apply the taxonomy mapping
tax = self.taxonomy_map.get(taxonomy, taxonomy)
# Retrieve the fragility function matching the taxonomy
ff = self.fragility_functions[taxonomy]
if ff["format"] == "discrete":
# Process a discrete fragility function
im_type = ff["imt"]
column = self.gm_type_column_map.get(im_type, None)
if column is None:
raise ValueError(
"Necessary intensity measure type %s not given in ground-motion file."
% (im_type)
)
im_value = intensity_measure_values[column]
ff = self.fragility_functions[tax]
# Retrieve the necessary intensity-measure type and value for the given taxonomy
im_type = ff["imt"]
column = self.gm_type_column_map.get(im_type, None)
if column is None:
raise ValueError(
"Necessary intensity measure type %s not given in ground-motion file."
% (im_type)
)
im_value = intensity_measure_values[column]
if ff["format"] == "discrete":
if im_value <= ff["imv"][0]:
# The intensity-measure value is smaller than or equal to the minimum
# intensity-measure value of the fragility function
......@@ -222,8 +327,28 @@ class FragilityModel:
# Interpolating PoE values for the given gm-value between x and y
# range (using interpolation)
poes.append(np.interp(im_value, x, y))
# Calculating PoOs
# pylint: disable=E1130
poos = list(-1 * np.diff(poes))
poos.append(poes[-1])
return poes[1:5], poos, im_value, im_type
# Process a continuous fragility function
else:
# If the minimum, maximum and no-damage limit intensity measures are defined
# then retrieve them and clip in input ground motion intensity measure
min_iml = float(ff.get("minIML", 0.0))
max_iml = float(ff.get("maxIML", np.inf))
# no_damage_limit = float(ff.get("noDamageLimit", 0.0))
if im_value < min_iml:
im_value = min_iml
if im_value > max_iml:
im_value = max_iml
# Prepare return value
poes = [1.0]
# Get the PoEs
for ls in self.limit_states:
# Interpolating PoE values for the given gm-value between x and y
# range (using interpolation)
poes.append(ff[ls]["function"].cdf(im_value))
# Calculating PoOs
# pylint: disable=E1130
poos = list(-1 * np.diff(poes))
poos.append(poes[-1])
return poes[1:5], poos, im_value, im_type
......@@ -29,7 +29,8 @@ import FragilityModel
def damage_calculator(
exposure_filepath,
fragility_filepath,
intensity_measure_map,
intensity_measure_map_filepath,
taxonomy_map_filepath,
gm_field_filepath,
interpolation_method="linear",
result_filepath="damage_result.csv",
......@@ -52,10 +53,12 @@ def damage_calculator(
# Read inputs.
gm_field = np.loadtxt(gm_field_filepath, delimiter=",", skiprows=1)
exposure = pd.read_csv(exposure_filepath)
esrm20_fragility_model = FragilityModel.FragilityModel(
fragility_filepath, intensity_measure_map
fragility_model = FragilityModel.FragilityModel(
fragility_filepath,
gm_field_filepath,
intensity_measure_map_filepath,
taxonomy_map_filepath,
)
esrm20_fragility_model.create_gmf_column_map(gm_field_filepath)
# Interpolate the ground-motion values for all the assets of the exposure model.
# pylint: disable=E1101
......@@ -70,7 +73,7 @@ def damage_calculator(
num_buildings = exposure["number"][asset]
# Computing probabilities of exceedance and occurrence.
poes, poos, im_value, im_type = esrm20_fragility_model(taxonomy, full_gm_field[asset])
poes, poos, im_value, im_type = fragility_model(taxonomy, full_gm_field[asset])
all_poes.append(poes)
all_poos.append(poos)
all_im_values.append(im_value)
......@@ -86,13 +89,13 @@ def damage_calculator(
exposure["PoEs"] = all_poes
exposure["PoOs"] = all_poos
exposure[esrm20_fragility_model.loss_category + "_no_damage"] = [
exposure[fragility_model.loss_category + "_no_damage"] = [
row[0] for row in damage_by_assets
]
for i in range(len(esrm20_fragility_model.limit_states)):
exposure[
esrm20_fragility_model.loss_category + "_" + esrm20_fragility_model.limit_states[i]
] = [row[(i + 1)] for row in damage_by_assets]
for i in range(len(fragility_model.limit_states)):
exposure[fragility_model.loss_category + "_" + fragility_model.limit_states[i]] = [
row[(i + 1)] for row in damage_by_assets
]
exposure.to_csv(result_filepath, index=False)
......@@ -165,6 +168,16 @@ if __name__ == "__main__":
help="Result file exists. Choose another name or set --overwrite"
+ "to overwrite the existing result file.",
)
parser.add_argument(
"-t",
"--taxonomy-map",
required=False,
type=str,
help="file path to the file that includes the mapping from complex taxonomy "
+ "strings to simple taxonomy strings as provided by the fragility "
+ "functions. If not provided, taxonomies will be passed directly to "
+ "the fragility functions.",
)
args = parser.parse_args()
# read arguments from command line.
......@@ -175,6 +188,7 @@ if __name__ == "__main__":
exposure_filepath = args.exposure
result_filepath = args.results
overwrite_result_file = args.overwrite
taxonomy_map_filepath = args.taxonomy_map
if os.path.exists(result_filepath):
if not overwrite_result_file:
......@@ -192,6 +206,7 @@ if __name__ == "__main__":
exposure_filepath,
fragility_filepath,
intensity_measure_map,
taxonomy_map_filepath,
gm_field_filepath,
interpolation_method,
result_filepath,
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment