Commit 00e3be66 authored by Tara Evaz Zadeh's avatar Tara Evaz Zadeh
Browse files

added Get_Full_GMF.py,taxonomyToFragility.py,OriginId_to_Polygon.py and...

added Get_Full_GMF.py,taxonomyToFragility.py,OriginId_to_Polygon.py and get_PoEs.py as a module called losslib.py
parent 76b181dc
def Get_Full_GMF(groundMotionField,lons,lats,method):
"""
Input:
------
- groundMotionField: contains locations and their ground-motion values in
the format of lon,lat,gmValueofType1,gmValueofType2.... Please note that
if there are more than 3 columns, 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:
[lon,lat,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
>>> lons
0 23.687500
1 23.687500
2 23.687500
3 23.687500
4 23.687500
...
125014 24.026389
125015 24.026389
125016 24.026389
125017 24.026389
125018 24.026389
Name: lon, Length: 125019, dtype: float64
- 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
>>> lats
0 38.340278
1 38.340278
2 38.340278
3 38.340278
4 38.340278
...
125014 37.651389
125015 37.651389
125016 37.651389
125017 37.651389
125018 37.651389
Name: lat, Length: 125019, dtype: float64
- method: {‘linear’, ‘nearest’, ‘cubic’}, optional'
Output:
------
- fullGroundMotionField: an array of logitudes and latitudes along with
their interpolated values (numpy-ndarray of shape(a,n))
creates ground-motion values using 2-dimensional interpolation for the
given lons and lats.
Example extract:
>>> fullGroundMotionField
array([[2.36875000e+01, 3.83402778e+01, 1.79051641e-01, 3.62318501e-01,
2.52607434e-01, 1.52471069e-01],
[2.36875000e+01, 3.83402778e+01, 1.79051641e-01, 3.62318501e-01,
2.52607434e-01, 1.52471069e-01],
[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],
[2.40263889e+01, 3.76513889e+01, 2.03267876e-02, 4.61331339e-02,
3.06813504e-02, 1.81060137e-02],
[2.40263889e+01, 3.76513889e+01, 2.03267876e-02, 4.61331339e-02,
3.06813504e-02, 1.81060137e-02]])
"""
import pandas as pd
import numpy as np
from scipy.interpolate import griddata
# 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.
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)
import csv
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)
import csv
def OriginId_to_Polygon(PolygonSource,exposureType):
"""
Creates an extended map of Origin-id to polygons.
The input map 'PolygonSource' contains the mapping for each
Origin-id to a polygon along with its Cell-id if it's an OBM polygon source
as the 'OriginIdToPolygonMap'.
Input:
------
- exposureType: string. either 'OBM' or 'cell'.
- PolygonSource: origin-id to polygons file map with semicolon as the
delimiter (csv.reader) following the format:
['OriginID'; 'Polygon'; 'cell_ID']. The 'cell_ID' would only be available
if your input polygon belongs to the OBM data(single buildings).
Example file extract:
if you have cell polygon source:
>>> PolygonSource
cell_2410244527;POLYGON ((23.68611111111113 38.3388888888889,
23.6888888888889 38.3388888888889, 23.6888888888889 38.34166666666667,
23.68611111111113 38.34166666666667, 23.68611111111113 38.3388888888889))
cell_2410244528;POLYGON ((23.6888888888889 38.3388888888889,
23.69166666666666 38.3388888888889, 23.69166666666666 38.34166666666667,
23.6888888888889 38.34166666666667, 23.6888888888889 38.3388888888889))
...
if you have OBM polygon source:
>>> PolygonSource
OSM_517924352;POLYGON((2641289.26886243 4582010.12482994,2641299.67723482
4582007.12923406,2641304.17454225 4582019.8887829,2641294.42295486
4582023.47785041,2641291.51751615 4582013.47368108,2641290.54903658
4582013.85519583,2641289.26886243 4582010.12482994));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 polygonSource 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:
------
- OriginIdToPolygonMap: Dictionary containing the origin-id to polygons
(If the polygons show single buildings, the dictionary also considers their
cell-IDs ) following the format below:
When polygons define cells:
{OriginID_string: [Polygon]}
Example extract:
>>> OriginIdToPolygonMap
{'cell_2410244527': 'POLYGON ((23.68611111111113 38.3388888888889,
23.6888888888889 38.3388888888889, 23.6888888888889 38.34166666666667,
23.68611111111113 38.34166666666667, 23.68611111111113 38.3388888888889))'
,'cell_2410244528': 'POLYGON ((23.6888888888889 38.3388888888889,
23.69166666666666 38.3388888888889, 23.69166666666666 38.34166666666667,
23.6888888888889 38.34166666666667, 23.6888888888889 38.3388888888889))',
... }
In case the polygons define buildings:
{OriginID_string: [Polygon, CellID]}
Example extract:
>>> OriginIdToPolygonMap
{'OSM_517924352': ['POLYGON((2641289.26886243 4582010.12482994,
2641299.67723482 4582007.12923406,2641304.17454225 4582019.8887829,
2641294.42295486 4582023.47785041,2641291.51751615 4582013.47368108,
2641290.54903658 4582013.85519583,2641289.26886243 4582010.12482994))',
'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'], ... }
"""
# Prepare return variable
OriginIdToPolygonMap = {}
# to know if the polygons refers to single buildings instead of cells.
# In this case we need to also know the cell_ID that this building polygon
# is located in.
if exposureType.startswith('OBM'):
# Loop through the originId_to_polygon map
for polygonMappingItem in PolygonSource:
# Check if already one polygon for a given OriginId has been
# selected.
if polygonMappingItem[0] in OriginIdToPolygonMap:
# Ignore the additional polygon to keep everything unambiguous
pass
# Create the entry in the extended map with adding the cell-ID that
# the building belongs to ('polygonMappingItem[2]')
OriginIdToPolygonMap[polygonMappingItem[0]] = [polygonMappingItem[1],polygonMappingItem[2]]
else:
for polygonMappingItem in PolygonSource:
if polygonMappingItem[0] in OriginIdToPolygonMap:
# implement your duplicate row(polygonMappingItem) handling here
pass
OriginIdToPolygonMap[polygonMappingItem[0]] = polygonMappingItem[1]
return(OriginIdToPolygonMap)
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 excedance.
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]
"""
import numpy as np
from scipy import interpolate
#can be done using spline
#https://stackoverflow.com/questions/2745329/how-to-make-scipy-interpolate-give-an-extrapolated-result-beyond-the-input-range
# 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=[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))
PoOs = list()
for i in [1,2,3,4]:
t = np.subtract(PoEs[i-1], PoEs[i])
PoOs.append(t)
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 [PoE_slight,PoE_moderate,PoE_extensive,PoE_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))
# Prepare return value
PoOs = list()
# Looping over all 4 PoE damage state values
for i in [1,2,3,4]:
# Calculating PoOs
t = np.subtract(PoEs[i-1], PoEs[i])
PoOs.append(t)
PoOs.append(PoEs[-1])
return(PoEs,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