Commit 98fd22ba authored by Danijel Schorlemmer's avatar Danijel Schorlemmer
Browse files

Merge branch 'losslib' into 'master'

losslib

See merge request !4
parents 76b181dc 32310cd7
import pandas as pd
import numpy as np
from scipy.interpolate import griddata
import csv
from scipy import interpolate
def Get_Full_GMF(groundMotionField, 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:
------
- groundMotionField: contains locations and the ground-motion values of
each location in the format of lon,lat,gmValueofType1,gmValueofType2,etc.
Please note that if the groundMotionField contains values for more than one
ground motion type, This function will do interpolation multiple times and
each time using one of the gmValueTypes (Numpy nd-array of shape (m,n))
with the format as below.
[lon,lat,gmValueofType1,gmValueofType2...]
As an example:
[23.6875,38.3402777778,'PGA','SA(0.3)','SA(0.6)','SA(1.0)']
Example extract:
>>> groundMotionField
array([[2.289028e+01, 3.631806e+01, 4.822530e-03, 1.255729e-02,
8.736364e-03, 5.612176e-03],
[2.289028e+01, 3.632083e+01, 4.830462e-03, 1.257595e-02,
8.748890e-03, 5.619673e-03],
...,
[2.410694e+01, 3.818194e+01, 4.961839e-02, 1.078893e-01,
7.710517e-02, 4.852522e-02]])
- lons: longitudes of the locations we want to interpolate over
(pandas.DataFrame(series) OR numpy-1darray of len a)
Example extract:
>>> lons
23.687500
23.687500
...
24.026389
24.026389
- lats: latitudes of the locations we want to interpolate over
(pandas.DataFrame(series) OR numpy-1darray of len a)
Example extract:
>>> lats
38.340278
38.340278
...
37.651389
37.651389
- method: {‘linear’, ‘nearest’, ‘cubic’}, optional'. Linear by default.
Output:
------
- fullGroundMotionField: an array of logitudes and latitudes along with
their interpolated values (numpy-ndarray of shape(a,n))
Example extract:
>>> fullGroundMotionField
array([[2.36875000e+01, 3.83402778e+01, 1.79051641e-01, 3.62318501e-01,
2.52607434e-01, 1.52471069e-01],
...,
[2.40263889e+01, 3.76513889e+01, 2.03267876e-02, 4.61331339e-02,
3.06813504e-02, 1.81060137e-02]])
"""
# points_given: the input points for interpolation
points_given = np.vstack((groundMotionField[:,0], groundMotionField[:,1]))\
.transpose()
# points_todo: points to do interpolation over
points_todo = np.vstack((lons, lats)).transpose()
fullGroundMotionField = 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 groundMotionField.
for gmType in range(2, groundMotionField.shape[1], 1):
# gmvs_given: the input ground-motion values for interpolation.
gmvs_given = groundMotionField[:, gmType]
# gm_value_griddata : interpolated values for each ground motion type.
gm_value_griddata = griddata(points_given, gmvs_given, points_todo, \
method=method)
fullGroundMotionField = np.column_stack((fullGroundMotionField, \
gm_value_griddata))
return(fullGroundMotionField)
def Taxonomy_to_Fragility(gmDict, taxonomyToFragilitySource, fragilityFileDir):
"""
Creates an extended map of taxonomies to fragility function.
The input map 'taxonomyToFragilitySource' 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 'taxonomyToFragilityMap'.
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.0)': 5, 'SA(1)': 5}
- taxonomyToFragilitySource: taxonomy to fragility-function file map
(csv.reader ) 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']
...
- fragilityFileDir: directory of fragility-function files
Example extract:
>>> fragilityFileDir
'/home/laptop/fragilities'
Output:
------
- taxonomyToFragilityMap: 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
taxonomyToFragilityMap = {}
# Loop through the taxonomy-to-fragility-function map
for mappingItem in taxonomyToFragilitySource:
# Open the fragility-function file corresponding to the taxonomy in
# 'mappingitem[1]'
fragilityFunction = list(csv.reader(open(fragilityFileDir + "/" + \
mappingItem[1] + ".csv")))
# Check if already one fragility function for a given GM type has been
# selected
if mappingItem[0] in taxonomyToFragilityMap:
# 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]'
taxonomyToFragilityMap[mappingItem[0]] = [mappingItem[1], \
gmDict[fragilityFunction[0][0]]]
return(taxonomyToFragilityMap)
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),
considering the ground-motion value (in the first line of the fragility
function).
Input:
------
- gm_value : value of the ground-motion.(Either int or float)
- fragility_function: An array with first lines indicating intensity
measure levels,slight damages,moderate damages,extensive damages and
complete damages in order. (Numpy nd-array of shape (5,100) in our case.)
Example extract:
>>> fragility_function
array([[5.000000e-02, 5.217900e-02, 5.445400e-02, 5.682700e-02,
5.930400e-02, 6.188900e-02, 6.458700e-02, 6.740200e-02,
...
3.004239e+00, 3.135186e+00, 3.271841e+00, 3.414451e+00],
[0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00,
0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00,
...
7.572240e-01, 7.797630e-01, 8.010520e-01, 8.210500e-01,
8.397340e-01, 8.570950e-01, 8.731380e-01, 8.878820e-01],
[0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00,
...
3.827660e-01, 4.112880e-01, 4.402870e-01, 4.696100e-01],
[0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00,
...
1.812150e-01, 2.013660e-01, 2.228030e-01, 2.454840e-01],
[0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00,
...
9.403000e-02, 1.070540e-01, 1.213310e-01, 1.368970e-01]])
Output:
------
- PoEs: A list with element as probabilities of exceedance.
Example extract:
>>> PoEs
[1, 5.6318317919292305e-06, 0.0, 0.0, 0.0]
- PoOs: A list with element as probabilities of occurence.
Example extract:
>>> PoOs
[0.9999943681682081, 5.6318317919292305e-06, 0.0, 0.0, 0.0]
"""
# 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]
idx_iml_min = -2
iml_max=imls[-1]
idx_iml_max = -1
x=[iml_min, iml_max]
PoEs = [1]
for Poes in [PoEs_slight, PoEs_moderate, PoEs_extensive, PoEs_complete]:
y = [Poes[idx_iml_min], Poes[idx_iml_max]]
f = interpolate.interp1d(x, y, fill_value='extrapolate')
PoEs.append(f(gm_value))
# Calculating PoOs
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
idx_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
idx_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[idx_iml_min], Poes[idx_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
PoOs = list(-np.diff(PoEs))
PoOs.append(PoEs[-1])
return(PoEs[1:5],PoOs)
Markdown is supported
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