losslib.py 16 KB
Newer Older
1
2
#!/usr/bin/env python3

3
4
# Copyright (c) 2020-2021:
#   Helmholtz-Zentrum Potsdam Deutsches GeoForschungsZentrum GFZ
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
#
# 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/.


20
21
22
23
import numpy as np
from scipy.interpolate import griddata
import csv
from scipy import interpolate
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
24
25


26
def get_full_GMF(ground_motion_field, lons, lats, method="linear"):
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
27
28

    """
29
30
    Returns ground-motion values using 2-dimensional interpolation for the
    given longitudes and latitudes. The deafult interpolation method is linear.
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
31
32
33

    Input:
    ------
34
    - ground_motion_field: (Numpy nd-array of shape(m,n); m records with long,
35
36
37
        lat, and n-2 gm-types)
        contains locations and the ground-motion values of each location
        in the format of lon, lat, gmValueofType1, gmValueofType2, ...,
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
38
        gmValueofTypeN.
39
40
        Please note that if the ground_motion_field contains values for
        more than one ground-motion type, This function will do
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
41
42
        interpolation multiple times and each time using one of the
        gmValueTypes with the format as below:
43
44
45
46
        [lon, lat, gmValueofType1, gmValueofType2, ..., gmValueofTypeN]
        As an example:
        [23.6875,38.3402777778,'PGA','SA(0.3)','SA(0.6)','SA(1.0)']
        Example extract:
47
        >>> ground_motion_field
48
49
50
51
52
53
        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]
            ])
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
54

55
56
57
58
    - lons: (pandas.DataFrame(series) OR numpy-1darray of len a)
        longitudes of the locations we want to interpolate over.
        Example extract:
        >>> lons
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
59
60
61
        23.6875
        ...
        24.0264
62

Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
63
    - lats: (pandas.DataFrame(series) OR numpy-1darray of len a)
64
65
66
        latitudes of the locations we want to interpolate over.
        Example extract:
        >>> lats
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
67
68
        38.340278
        ...
69
70
71
72
        37.651389

    - method: (string), optional
        {‘linear’, ‘nearest’, ‘cubic’}, optional'. Linear by default.
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
73

74
75
    Output:
    -------
76
    - full_ground_motion_field: (numpy-ndarray of shape(a,n))
77
78
        an array of logitudes and latitudes along with their interpolated values
        Example extract:
79
        >>> full_ground_motion_field
80
81
82
83
84
85
86
        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]
        ])
    """

Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
87
    # points_given: the input points for interpolation
88
    points_given = np.vstack((ground_motion_field[:, 0], ground_motion_field[:, 1])).transpose()
89
    # points_todo: points to do interpolation over
90
    points_todo = np.vstack((lons, lats)).transpose()
91
    full_ground_motion_field = np.vstack((np.array(lons), np.array(lats))).transpose()
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
92
    # Griddata Interpolation
93
    # The loop changes over the columns of the ground motion field and enables
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
94
95
    # us to have interpolation for all different given ground-motion types.
    # Please note that the columns follow the order below:
96
97
    # lon,lat,gmValueofType1,gmValueofType2,etc. We want to do the interpolation
    # over all the ground motion types. Thus the range begins from third column
98
99
    # to the last column of the ground_motion_field.
    for gm_type in range(2, ground_motion_field.shape[1], 1):
100
        # gmvs_given: the input ground-motion values for interpolation.
101
        gmvs_given = ground_motion_field[:, gm_type]
102
        # gm_value_griddata : interpolated values for each ground motion type.
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
103
        gm_value_griddata = griddata(points_given, gmvs_given, points_todo, method=method)
104
105
106
107
        full_ground_motion_field = np.column_stack(
            (full_ground_motion_field, gm_value_griddata)
        )
    return full_ground_motion_field
108
109


110
def taxonomy_to_fragility(gmDict, taxonomy_to_fragility_source, fragility_pathname):
111
    """
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
112
    Creates an extended map of taxonomies to fragility function.
113
    The input map 'taxonomy_to_fragility_source' contains the mapping for each
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
114
    taxonomy to a fragility function file from which the ground-motion type is
115
    read to be written to the extended map 'taxonomy_to_fragility_map'.
116
117
118

    Input:
    ------
119
    - gmDict: (dictionary)
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
120
        key: ground-motion type; value: column number in the ground-motion
121
122
123
124
        field file
        Example extract:
        >>> gmDict
        {'PGA': 2, 'SA(0.3)': 3, 'SA(0.6)': 4, 'SA(1)': 5}
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
125

126
    - taxonomy_to_fragility_source: (csv.reader)
127
128
129
130
131
132
133
        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']
        ...
134

135
    - fragility_pathname: (string)
136
137
        directory of fragility-function files
        Example extract:
138
        >>> fragility_pathname
139
        '/home/laptop/fragilities'
140
141
142
143


    Output:
    ------
144
    - taxonomy_to_fragility_map: (Dictionary)
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
145
        containing the taxonomy to fragility function map considering the
146
        ground-motion types. It follows the format:
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
147
        {taxonomy_string: [fragility-function_filename, column of
148
149
150
151
152
        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], ... }
153
154
155
    """

    # Prepare return variable
156
    taxonomy_to_fragility_map = {}
157
    # Loop through the taxonomy-to-fragility-function map
158
    for mapping_item in taxonomy_to_fragility_source:
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
159
        # Open the fragility-function file corresponding to the taxonomy in
160
        # 'mapping_item[1]'
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
161
        fragilityFunction = list(
162
            csv.reader(open(fragility_pathname + "/" + mapping_item[1] + ".csv"))
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
163
164
        )
        # Check if already one fragility function for a given GM type has been
165
        # selected
166
        if mapping_item[0] in taxonomy_to_fragility_map:
167
168
169
            # Ignore the additional fragility function to keep everything
            # unambiguous
            pass
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
170
        # Create the entry in the extended map with adding the type of ground
171
        # motion 'fragilityFunction[0][0]'
172
173
        taxonomy_to_fragility_map[mapping_item[0]] = [
            mapping_item[1],
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
174
175
            gmDict[fragilityFunction[0][0]],
        ]
176
    return taxonomy_to_fragility_map
177
178


179
def origin_id_to_geometry(geometry_source, exposure_type):
180
    """
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
181
    Creates a dictionary of Origin-ids and their polygons.
182
    The input map 'geometry_source' contains the mapping for each
183
    Origin-id to a polygon along with its cell-id if it's a building geometry source
184
    as the 'origin_id_to_geometry_map'.
185

186
    Input:
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
187
    ------
188
    - exposure_type: (string)
189
        Either 'building' or 'cell'
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
190

191
192
    - geometry_source: (csv.reader)
        Origin-id to geometry file map with semicolon as the delimiter,
193
        following the format:
194
        ['OriginID'; 'geometry'; 'cell_ID']. 'cell_ID' would only be available if
195
        the input geometry belongs to the building data(single buildings). 'cell_ID'
196
        here points to the cell that contains the geometry.
197
        Example extract:
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
198

199
200
        if you have cell geometry source:
        >>> geometry_source
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
201
202
        cell_2410244527;POLYGON ((23.6861 38.3389,
        23.6889 38.3389, 23.6889 38.34167,
203
        23.6861 38.3417, 23.6861 38.3389))
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
204
205
        cell_2410244528;POLYGON ((23.6889 38.3389,
        23.6916 38.3389, 23.6916 38.3416,
206
        23.6889 38.3416, 23.6889 38.3389))
207
        ...
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
208

209
        if you have building geometry source:
210
        >>> geometry_source
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
211
212
213
        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
214
215
        4582013.8551,2641289.2688 4582010.1248));cell_2425278141
        OSM_452860266;POLYGON((2647123.74 4559674.98,2647130.88 4559680.21,
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
216
        2647136 4559673.23,2647138.57 4559675.12,2647143.21 4559668.8,2647133.5
217
        4559661.68,2647123.74 4559674.98));cell_2432665360
218
        ...
219

220
        please note that it doesn't matter if your input geometry_source has more
221
222
        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
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
223
224
225
        only needs to have the same format as this column has in the exposure
        file.

226

227
228
    Output:
    ------
229
    - origin_id_to_geometry_map: (Dictionary)
230
        contains the origin-id to polygons items.
231
        (If the geometry show single buildings, the dictionary also considers
232
        their cell-IDs ) following the format below:
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
233

234
235
        When geometry defines cells:
        {OriginID_string: [geometry]}
236
237

        Example extract:
238
        >>> origin_id_to_geometry_map
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
239
240
        {'cell_2410244527': 'POLYGON ((23.6861 38.33889,
        23.6889 38.3389, 23.6889 38.3416,
241
        23.6861 38.3416, 23.6861 38.3389))'
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
242
243
        ,'cell_2410244528': 'POLYGON ((23.6889 38.3389,
        23.6916 38.3389, 23.6916 38.3416,
244
245
246
247
        23.6889 38.3416, 23.6889 38.3389))',
         ... }


248
249
        In case the geometry defines buildings:
        {OriginID_string: [geometry, CellID]}
250
251

        Example extract:
252
        >>> origin_id_to_geometry_map
253
254
255
        {'OSM_517924352': ['POLYGON((2641289.2688 4582010.1248,
        2641299.6772 4582007.1292,2641304.1745 4582019.8887,
        2641294.4229 4582023.4778,2641291.5175 4582013.4736,
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
256
        2641290.5490 4582013.8551,2641289.2688 4582010.1248))',
257
258
259
260
        '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'], ...}
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
261
262


263
    Please Note that the first columns of the geometry_source does not have to
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
264
265
    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
266
267
268
    input file.
    """
    # Prepare return variable
269
270
271
    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
272
    # is located in.
273
    if exposure_type == "building":
274
275
276
        # 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
277
            # selected.
278
279
            if geometry_mapping_item[0] in origin_id_to_geometry_map:
                # Ignore the additional geometry to keep everything unambiguous
280
281
                pass
            # Create the entry in the extended map with adding the cell-ID that
282
283
284
285
            # 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],
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
286
            ]
287
    else:
288
289
290
        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
291
                pass
292
293
            origin_id_to_geometry_map[geometry_mapping_item[0]] = geometry_mapping_item[1]
    return origin_id_to_geometry_map
294

295

296
297
def get_PoEs(fragility_function, gm_value):
    """
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
298
299
300
    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
301
    the fragility function).
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
302

303
304
    Input:
    ------
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
305
    - gm_value : (int OR float)
306
307
308
309
        value of the ground-motion.
        Example extract:
        >>> gm_value
            2.152E-02
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
310
311
312
313
314

    - 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.)
315
316
317
318
319
320
321
322
323
        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]
            ])
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
324

325
326
327
    Output:
    ------
    - PoEs: (list)
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
328
329
        A list with element as probabilities of exceedance.
        Example extract:
330
331
        >>> PoEs
        [0.73, 0.52, 0.32, 0.18]
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
332

333
    - PoOs: (list)
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
334
        A list with element as probabilities of occurence.
335
336
337
        Example extract:
        >>> PoOs
        [0.27, 0.21, 0.2 , 0.14, 0.18]
338
    """
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
339
340
    # Define and extract intensity measure level values, slight, moderate,
    # extensive and complete damage probabilities of exceedance from the
341
    # fragility function
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
342
343
344
345
346
    imls = fragility_function[0]
    PoEs_slight = fragility_function[1]
    PoEs_moderate = fragility_function[2]
    PoEs_extensive = fragility_function[3]
    PoEs_complete = fragility_function[4]
347
348
349
350
351
    # 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.
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
352
353
354
    if gm_value < imls[0]:
        PoEs = [1, 0, 0, 0, 0]
        PoOs = [1, 0, 0, 0, 0]
355

Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
356
357
    # in the case where our ground-motion value is bigger than the largest iml
    elif gm_value > imls[-1]:
358

Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
359
        iml_min = imls[-2]
360
        index_iml_min = -2
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
361
        iml_max = imls[-1]
362
        index_iml_max = -1
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
363
        x = [iml_min, iml_max]
364
365
        PoEs = [1]

366
        for Poes in [PoEs_slight, PoEs_moderate, PoEs_extensive, PoEs_complete]:
367
            y = [Poes[index_iml_min], Poes[index_iml_max]]
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
368
            f = interpolate.interp1d(x, y, fill_value="extrapolate")
369
            PoEs.append(f(gm_value))
370
        # Calculating PoOs
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
371
        # pylint: disable=E1130
372
        PoOs = list(-np.diff(PoEs))
373
        PoOs.append(PoEs[-1])
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
374
375

    # In the case where our ground-motion value is in the range of defined imls
376
377
    else:
        # Largest iml that is smaller than the given ground-motion value
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
378
        iml_min = imls[imls < gm_value].max()
379
        # Finding index of the iml_min
380
        index_iml_min = np.searchsorted(imls, iml_min, side="left")
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
381
382
        # Smallest iml that is larger than the given ground-motion value
        iml_max = imls[imls > gm_value].min()
383
        # Finding index of the iml_max
384
        index_iml_max = np.searchsorted(imls, iml_max, side="left")
385
386

        # Two bounds of the interpolation
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
387
        x = [iml_min, iml_max]
388
389
390
        # Prepare return value
        PoEs = [1]
        # Looping over all damage states
391
        for Poes in [PoEs_slight, PoEs_moderate, PoEs_extensive, PoEs_complete]:
392
            y = [Poes[index_iml_min], Poes[index_iml_max]]
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
393
            # interpolating PoE values for the given gm-value between x and y
394
            # range (using interpolation)
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
395
            PoEs.append(np.interp(gm_value, x, y))
396
        # Calculating PoOs
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
397
        # pylint: disable=E1130
398
        PoOs = list(-np.diff(PoEs))
399
        PoOs.append(PoEs[-1])
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
400
    return (PoEs[1:5], PoOs)