#!/usr/bin/env python3
# Copyright (c) 2020-2021:
# Helmholtz-Zentrum Potsdam Deutsches GeoForschungsZentrum GFZ
#
# This program is free software: you can redistribute it and/or modify it
# under the terms of the GNU Affero General Public License as published by
# the Free Software Foundation, either version 3 of the License, or (at
# your option) any later version.
#
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Affero
# General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with this program. If not, see http://www.gnu.org/licenses/.
import numpy as np
from scipy.interpolate import griddata
import csv
from scipy import interpolate
def get_full_GMF(ground_motion_field, lons, lats, method="linear"):
"""
Returns ground-motion values using 2-dimensional interpolation for the
given longitudes and latitudes. The deafult interpolation method is linear.
Input:
------
- ground_motion_field: (Numpy nd-array of shape(m,n); m records with long,
lat, and n-2 gm-types)
contains locations and the ground-motion values of each location
in the format of lon, lat, gmValueofType1, gmValueofType2, ...,
gmValueofTypeN.
Please note that if the ground_motion_field contains values for
more than one ground-motion type, This function will do
interpolation multiple times and each time using one of the
gmValueTypes with the format as below:
[lon, lat, gmValueofType1, gmValueofType2, ..., gmValueofTypeN]
As an example:
[23.6875,38.3402777778,'PGA','SA(0.3)','SA(0.6)','SA(1.0)']
Example extract:
>>> ground_motion_field
array([
[2.289e+01, 3.631e+01, 4.822e-03, 1.255e-02, 8.736e-03, 5.612e-03],
[2.289e+01, 3.632e+01, 4.830e-03, 1.257e-02, 8.748e-03, 5.619e-03],
...,
[2.410e+01, 3.818e+01, 4.961e-02, 1.078e-01, 7.710e-02, 4.852e-02]
])
- lons: (pandas.DataFrame(series) OR numpy-1darray of len a)
longitudes of the locations we want to interpolate over.
Example extract:
>>> lons
23.6875
...
24.0264
- lats: (pandas.DataFrame(series) OR numpy-1darray of len a)
latitudes of the locations we want to interpolate over.
Example extract:
>>> lats
38.340278
...
37.651389
- method: (string), optional
{‘linear’, ‘nearest’, ‘cubic’}, optional'. Linear by default.
Output:
-------
- full_ground_motion_field: (numpy-ndarray of shape(a,n))
an array of logitudes and latitudes along with their interpolated values
Example extract:
>>> full_ground_motion_field
array([
[2.36e+01, 3.83e+01, 1.79e-01, 3.62e-01, 2.52e-01, 1.52e-01],
...,
[2.40e+01, 3.76e+01, 2.03e-02, 4.61e-02, 3.07e-02, 1.81e-02]
])
"""
# points_given: the input points for interpolation
points_given = np.vstack((ground_motion_field[:, 0], ground_motion_field[:, 1])).transpose()
# points_todo: points to do interpolation over
points_todo = np.vstack((lons, lats)).transpose()
full_ground_motion_field = np.vstack((np.array(lons), np.array(lats))).transpose()
# Griddata Interpolation
# The loop changes over the columns of the ground motion field and enables
# us to have interpolation for all different given ground-motion types.
# Please note that the columns follow the order below:
# lon,lat,gmValueofType1,gmValueofType2,etc. We want to do the interpolation
# over all the ground motion types. Thus the range begins from third column
# to the last column of the ground_motion_field.
for gm_type in range(2, ground_motion_field.shape[1], 1):
# gmvs_given: the input ground-motion values for interpolation.
gmvs_given = ground_motion_field[:, gm_type]
# gm_value_griddata : interpolated values for each ground motion type.
gm_value_griddata = griddata(points_given, gmvs_given, points_todo, method=method)
full_ground_motion_field = np.column_stack(
(full_ground_motion_field, gm_value_griddata)
)
return full_ground_motion_field
def taxonomy_to_fragility(gmDict, taxonomy_to_fragility_source, fragility_pathname):
"""
Creates an extended map of taxonomies to fragility function.
The input map 'taxonomy_to_fragility_source' contains the mapping for each
taxonomy to a fragility function file from which the ground-motion type is
read to be written to the extended map 'taxonomy_to_fragility_map'.
Input:
------
- gmDict: (dictionary)
key: ground-motion type; value: column number in the ground-motion
field file
Example extract:
>>> gmDict
{'PGA': 2, 'SA(0.3)': 3, 'SA(0.6)': 4, 'SA(1)': 5}
- taxonomy_to_fragility_source: (csv.reader)
taxonomy to fragility-function file map following the format:
[taxonomy_string, fragility-function_filename, 'weight'], [...]
Example file extract:
['CR/LDUAL+CDM/HBET:6-/11.0', 'CR_LDUAL-DUM_H6', '1']
['CR/LDUAL+CDM/HBET:6-/5.0', 'CR_LDUAL-DUL_H6', '1']
['CR/LDUAL+CDM/HBET:6-/SOS/11.0', 'CR_LDUAL-DUL_H6', '1']
...
- fragility_pathname: (string)
directory of fragility-function files
Example extract:
>>> fragility_pathname
'/home/laptop/fragilities'
Output:
------
- taxonomy_to_fragility_map: (Dictionary)
containing the taxonomy to fragility function map considering the
ground-motion types. It follows the format:
{taxonomy_string: [fragility-function_filename, column of
ground-motion_type in ground-motion_field file]}
Example extract:
{'CR/LDUAL+CDM/HBET:6-/11.0': ['CR_LDUAL-DUM_H6', 4],
'CR/LDUAL+CDM/HBET:6-/5.0': ['CR_LDUAL-DUL_H6', 4],
'CR/LDUAL+CDM/HBET:6-/SOS/11.0': ['CR_LDUAL-DUL_H6', 4], ... }
"""
# Prepare return variable
taxonomy_to_fragility_map = {}
# Loop through the taxonomy-to-fragility-function map
for mapping_item in taxonomy_to_fragility_source:
# Open the fragility-function file corresponding to the taxonomy in
# 'mapping_item[1]'
fragilityFunction = list(
csv.reader(open(fragility_pathname + "/" + mapping_item[1] + ".csv"))
)
# Check if already one fragility function for a given GM type has been
# selected
if mapping_item[0] in taxonomy_to_fragility_map:
# Ignore the additional fragility function to keep everything
# unambiguous
pass
# Create the entry in the extended map with adding the type of ground
# motion 'fragilityFunction[0][0]'
taxonomy_to_fragility_map[mapping_item[0]] = [
mapping_item[1],
gmDict[fragilityFunction[0][0]],
]
return taxonomy_to_fragility_map
def origin_id_to_geometry(geometry_source, exposure_type):
"""
Creates a dictionary of Origin-ids and their polygons.
The input map 'geometry_source' contains the mapping for each
Origin-id to a polygon along with its cell-id if it's a building geometry source
as the 'origin_id_to_geometry_map'.
Input:
------
- exposure_type: (string)
Either 'building' or 'cell'
- geometry_source: (csv.reader)
Origin-id to geometry file map with semicolon as the delimiter,
following the format:
['OriginID'; 'geometry'; 'cell_ID']. 'cell_ID' would only be available if
the input geometry belongs to the building data(single buildings). 'cell_ID'
here points to the cell that contains the geometry.
Example extract:
if you have cell geometry source:
>>> geometry_source
cell_2410244527;POLYGON ((23.6861 38.3389,
23.6889 38.3389, 23.6889 38.34167,
23.6861 38.3417, 23.6861 38.3389))
cell_2410244528;POLYGON ((23.6889 38.3389,
23.6916 38.3389, 23.6916 38.3416,
23.6889 38.3416, 23.6889 38.3389))
...
if you have building geometry source:
>>> geometry_source
OSM_517924352;POLYGON((2641289.2688 4582010.1248,2641299.6772
4582007.1292,2641304.1745 4582019.8887,2641294.4229
4582023.4778,2641291.5175 4582013.4736,2641290.5490
4582013.8551,2641289.2688 4582010.1248));cell_2425278141
OSM_452860266;POLYGON((2647123.74 4559674.98,2647130.88 4559680.21,
2647136 4559673.23,2647138.57 4559675.12,2647143.21 4559668.8,2647133.5
4559661.68,2647123.74 4559674.98));cell_2432665360
...
please note that it doesn't matter if your input geometry_source has more
columns than needed. But the first columns should be as sited above.
The first column does not necessarily have to start with OSM or cell. It
only needs to have the same format as this column has in the exposure
file.
Output:
------
- origin_id_to_geometry_map: (Dictionary)
contains the origin-id to polygons items.
(If the geometry show single buildings, the dictionary also considers
their cell-IDs ) following the format below:
When geometry defines cells:
{OriginID_string: [geometry]}
Example extract:
>>> origin_id_to_geometry_map
{'cell_2410244527': 'POLYGON ((23.6861 38.33889,
23.6889 38.3389, 23.6889 38.3416,
23.6861 38.3416, 23.6861 38.3389))'
,'cell_2410244528': 'POLYGON ((23.6889 38.3389,
23.6916 38.3389, 23.6916 38.3416,
23.6889 38.3416, 23.6889 38.3389))',
... }
In case the geometry defines buildings:
{OriginID_string: [geometry, CellID]}
Example extract:
>>> origin_id_to_geometry_map
{'OSM_517924352': ['POLYGON((2641289.2688 4582010.1248,
2641299.6772 4582007.1292,2641304.1745 4582019.8887,
2641294.4229 4582023.4778,2641291.5175 4582013.4736,
2641290.5490 4582013.8551,2641289.2688 4582010.1248))',
'cell_2425278141'],
'OSM_452860266': ['POLYGON((2647123.74 4559674.98,2647130.88 4559680.21,
2647136 4559673.23,2647138.57 4559675.12,2647143.21 4559668.8,
2647133.5 4559661.68 ,2647123.74 4559674.98))', 'cell_2432665360'], ...}
Please Note that the first columns of the geometry_source does not have to
have the same format as either "OSM_517924352" or "cell_2410244527", but
the id have to have the same format as "origin_id" column of the exposure
input file.
"""
# Prepare return variable
origin_id_to_geometry_map = {}
# to know if the geometry refers to single buildings instead of cells.
# In this case we need to also know the cell_ID that this building geometry
# is located in.
if exposure_type == "building":
# Loop through the origin_id_to_geometry map
for geometry_mapping_item in geometry_source:
# Check if already one geometry for a given OriginId has been
# selected.
if geometry_mapping_item[0] in origin_id_to_geometry_map:
# Ignore the additional geometry to keep everything unambiguous
pass
# Create the entry in the extended map with adding the cell-ID that
# The building belongs to ('geometry_mapping_item[2]')
origin_id_to_geometry_map[geometry_mapping_item[0]] = [
geometry_mapping_item[1],
geometry_mapping_item[2],
]
else:
for geometry_mapping_item in geometry_source:
if geometry_mapping_item[0] in origin_id_to_geometry_map:
# Implement your duplicate row(geometry_mapping_item) handling here
pass
origin_id_to_geometry_map[geometry_mapping_item[0]] = geometry_mapping_item[1]
return origin_id_to_geometry_map
def get_PoEs(fragility_function, gm_value):
"""
This function interpolates PoE values (for slight,moderate,extensive and
complete damages (Lines 2, 3, 4 and 5 of the fragility function
respectively)), considering the ground-motion value (in the first line of
the fragility function).
Input:
------
- gm_value : (int OR float)
value of the ground-motion.
Example extract:
>>> gm_value
2.152E-02
- fragility_function: (Numpy nd-array of shape (5,100) in our case;
An array of 5 records with values of Intensity measure levels, values
of slight, moderate, extensive and complete damage probabilities of
exceedance, as the first to 5th records respectively.)
Example extract:
>>> fragility_function
array([
[5.0e-02, 5.2-02, 5.4-02, 5.6e-02, ..., 3.414451e+00],
[0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00, ..., 8.878820e-01],
[0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00, ..., 4.696100e-01],
[0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00, ..., 2.454840e-01],
[0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00, ..., 1.368970e-01]
])
Output:
------
- PoEs: (list)
A list with element as probabilities of exceedance.
Example extract:
>>> PoEs
[0.73, 0.52, 0.32, 0.18]
- PoOs: (list)
A list with element as probabilities of occurence.
Example extract:
>>> PoOs
[0.27, 0.21, 0.2 , 0.14, 0.18]
"""
# Define and extract intensity measure level values, slight, moderate,
# extensive and complete damage probabilities of exceedance from the
# fragility function
imls = fragility_function[0]
PoEs_slight = fragility_function[1]
PoEs_moderate = fragility_function[2]
PoEs_extensive = fragility_function[3]
PoEs_complete = fragility_function[4]
# If the ground-motion value is smaller than the smallest intesity measure
# level defined in the fragility-funtion, we have no damages, therefore the
# Probability of exceeding any of the damage states is zero and the
# probaility of occurence for no-damage, slight-damage, moderate-damage,
# extensive-damage, complete-damage is equal to 1,0,0,0,0 respectively.
if gm_value < imls[0]:
PoEs = [1, 0, 0, 0, 0]
PoOs = [1, 0, 0, 0, 0]
# in the case where our ground-motion value is bigger than the largest iml
elif gm_value > imls[-1]:
iml_min = imls[-2]
index_iml_min = -2
iml_max = imls[-1]
index_iml_max = -1
x = [iml_min, iml_max]
PoEs = [1]
for Poes in [PoEs_slight, PoEs_moderate, PoEs_extensive, PoEs_complete]:
y = [Poes[index_iml_min], Poes[index_iml_max]]
f = interpolate.interp1d(x, y, fill_value="extrapolate")
PoEs.append(f(gm_value))
# Calculating PoOs
# pylint: disable=E1130
PoOs = list(-np.diff(PoEs))
PoOs.append(PoEs[-1])
# In the case where our ground-motion value is in the range of defined imls
else:
# Largest iml that is smaller than the given ground-motion value
iml_min = imls[imls < gm_value].max()
# Finding index of the iml_min
index_iml_min = np.searchsorted(imls, iml_min, side="left")
# Smallest iml that is larger than the given ground-motion value
iml_max = imls[imls > gm_value].min()
# Finding index of the iml_max
index_iml_max = np.searchsorted(imls, iml_max, side="left")
# Two bounds of the interpolation
x = [iml_min, iml_max]
# Prepare return value
PoEs = [1]
# Looping over all damage states
for Poes in [PoEs_slight, PoEs_moderate, PoEs_extensive, PoEs_complete]:
y = [Poes[index_iml_min], Poes[index_iml_max]]
# interpolating PoE values for the given gm-value between x and y
# range (using interpolation)
PoEs.append(np.interp(gm_value, x, y))
# Calculating PoOs
# pylint: disable=E1130
PoOs = list(-np.diff(PoEs))
PoOs.append(PoEs[-1])
return (PoEs[1:5], PoOs)