damage_calculator_tile_version.py 10.7 KB
Newer Older
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
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
#
# 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/.

Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
19
20
21
22
23
import numpy as np
import csv
import datetime
import losslib

Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
24
25

def damageCalculator_TileVersion(
26
27
28
29
    full_ground_motion_field,
    result_filepath,
    fragility_pathname,
    exposure,
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
30
    taxonomy_conversion_path,
31
    ground_motion_filepath,
32
    geometry_source_path,
33
    exposure_type="cell",
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
34
35
36
37
    method="linear",
):

    """
38
    Returns a file "result_filepath" including damage results of a scenario earthquake using the
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
39
40
41
42
    "losslib" functions.

    Input:
    ------
43
    - full_ground_motion_field: (Numpy nd-array of shape(m,n); m records with long,
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
44
45
46
        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
47
        gmValueofTypeN.
48
49
        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
50
51
        interpolation multiple times and each time using one of the
        gmValueTypes with the format as below:
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
52
53
54
55
        [lon, lat, gmValueofType1, gmValueofType2, ..., gmValueofTypeN]
        As an example:
        [23.6875,38.3402777778,'PGA','SA(0.3)','SA(0.6)','SA(1.0)']
        Example extract:
56
        >>> full_ground_motion_field
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
57
58
59
60
61
62
63
64
65
66
        array([[
            [23.687   , 38.340,  0.179,  0.362 ,  0.252,
         0.152],
            ...,
            [23.690, 38.340,  0.177,  0.359,  0.25,
         0.151],
            ...,
            [23.693, 38.340,  0.176,  0.356,  0.248,
         0.149]
            ]])
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
67

68
    - result_filepath: (str)
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
69
70
        Result file address.
        Example extract:
71
        >>> result_filepath
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
72
        "/home/TileCalculations/DamageResult.csv"
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
73

74
    - fragility_pathname: (str)
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
75
76
        Address to the directory including all the fragility functions.
        Example extract:
77
        >>> fragility_pathname
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
78
        "/home/TileCalculations/Fragilities"
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
79

80
    - exposure: (pandas.DataFrame(series))
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
81
82
        informaion of the assets of the region of interest under hazard.
        Example extract:
83
        >>> exposure
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
84
85
86
87
88
89
90
91
92
93
94
95
96
97
                    id        lon        lat  ... admin_name    admin_ID        origin_id
        0    GDE_Ind_0  23.687500  38.340278  ...     Oropos  GR_3514913  cell_2410244527
        1    GDE_Ind_1  23.687500  38.340278  ...     Oropos  GR_3514913  cell_2410244527
        ...
        8    GDE_Ind_8  23.690278  38.340278  ...     Oropos  GR_3514913  cell_2410244528
        9    GDE_Ind_9  23.690278  38.340278  ...     Oropos  GR_3514913  cell_2410244528
        ...
        22  GDE_Ind_22  23.693056  38.340278  ...     Oropos  GR_3514913  cell_2410244529
        23  GDE_Ind_23  23.693056  38.340278  ...     Oropos  GR_3514913  cell_2410244529

    - taxonomy_conversion_path: (str)
        Address of the file including old and new taxonomy names.
        Example extract:
        >>> taxonomy_conversion_path
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
98
99
        "/home/TileCalculations/taxonomy_mapping_Europe.csv"

100
    - ground_motion_filepath: (str)
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
101
102
        Address to the ground-motion values file.
        Example extract:
103
        >>> ground_motion_filepath
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
104
        "/home/TileCalculations/shakemap1381_2.csv"
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
105

106
107
    - geometry_source_path: (str)
        Address to the file including the origin-ids and their respecive polygons
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
108
        (Could be whether each tile oan OSM building).
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
109
        Example extract:
110
        >>> geometry_source_path
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
111
112
        "/home/TileCalculations
        /M008_exposure_Attica_GDE_visual_v001_sat_27f_by_cell_reOrder.csv"
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
113

114
    - exposure_type: (string), optional
115
        {‘cell’, ‘building’}. Cell by default.
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
116
117
118

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

Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
120
121
    Output:
    -------
122
    - result_filepath: (arrays written to file)
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
123
        file containing the damage elements for each asset of the exposure.
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
124
        Example extract of the result file:
125
        geometry,,origin_id,,asset_id,,lon,,lat,,taxonomy,,gmfValue,,PoEs,,PoOs,,
126
        num_buildings,,structural_No-damage,,structural_Slight,,structural_Moderate,,
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
127
128
129
130
131
132
133
134
135
136
        structural_Extensive,,structural_Complete
        "POLYGON ((23.68611111111113 38.3388888888889, 23.6888888888889 38.3388888888889,
        23.6888888888889 38.34166666666667, 23.68611111111113 38.34166666666667,
        23.68611111111113 38.3388888888889))",,cell_2410244527,,GDE_Ind_0,,23.6875,,
        38.340277777800004,,CR/LFINF+CDM/H:2,,0.3623185006574049,,"[0.2495074860984066,
        0.008729386056721225, 0.0010441231376296645, 0.00021103727864033706]",,
        "[0.7504925139015934, 0.24077810004168537, 0.007685262919091561, 0.0008330858589893275,
        0.00021103727864033706]",,0.000629094517209,,0.0004721307257018916,,
        0.00015147218260022435,,4.834756765710135e-06,,5.24089746254536e-07,,
        1.3276239491934405e-07
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
137
    """
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
138

Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
139
140
141
142
143
    # Show the time the script begins running
    startTime = datetime.datetime.now()
    print(startTime)

    # Read inputs as Numpy arrays or Pandas data-frames
144
    taxonomy_to_fragility_source = csv.reader(open(taxonomy_conversion_path))
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
145
    # Skip the header
146
    next(taxonomy_to_fragility_source, None)
147
    geometry_source = csv.reader(open(geometry_source_path), delimiter=";")
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
148
    # Skip the header
149
    next(geometry_source, None)
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
150

151
152
153
154
155
156
157
    # Read each column of the exposure model
    taxonomies = exposure.taxonomy
    tot_num_buildings = exposure.number
    lons = exposure.lon
    lats = exposure.lat
    asset_ids = exposure.id
    origin_ids = exposure.origin_id
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
158
159

    # Begin Computation
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
160
161
162
    # Define a dictionary with keys as the ground-motion type and value as the column
    # number of the ground-motion type in the shakemap file.
    gmDict = {"PGA": 2, "SA(0.3)": 3, "SA(0.6)": 4, "SA(1.0)": 5, "SA(1)": 5}
163
    # Calling the function "taxonomy_to_fragility" to get a dictionary with keys as the
164
165
    # taxonomy and the values as both the fragility function name (excluding the ".csv" part)
    # and the column number of the respective ground-motion type in `ground-motion-field` file.
166
167
    taxonomy_to_fragility_map = losslib.taxonomy_to_fragility(
        gmDict, taxonomy_to_fragility_source, fragility_pathname
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
168
    )
169
170
    # Calling the function "origin_id_to_geometry" to get a dictionary with keys as the
    # origin_id and the value as the respective polygon.
171
    origin_id_to_geometry_map = losslib.origin_id_to_geometry(geometry_source, exposure_type)
172
    # Define number of columns that contain the data in the fragility function files.
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
173
174
    cls = range(1, 101)
    # Just a trick to have multiple commas between each result element, since we do
175
    # not want a single comma as the delimiter due to having geometries as a result element.
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
176
    a = [0, 2, 4, 6, 8]
177
178
179
180
181
182
    # Looping through each line of the exposure file to run the computations asset by asset.

    with open(result_filepath, "a+", newline="") as write_obj:
        csv_writer = csv.writer(write_obj)
        for asset in range(exposure.shape[0]):
            taxonomy = taxonomies.iloc[asset]
183
            fragilityfunction_filename = taxonomy_to_fragility_map[taxonomy][0] + ".csv"
184
185
186
187
188
            num_buildings = tot_num_buildings.iloc[asset]
            lon = lons.iloc[asset]
            lat = lats.iloc[asset]
            asset_id = asset_ids.iloc[asset]
            origin_id = origin_ids.iloc[asset]
189
            # In case of `exposure_type == "building"`, not only the geometry, but also the
190
191
            # respective cell ID (in which the asset is located), is extracted
            # using `origin_id_to_geometry_map`
192
            if exposure_type == "building":
193
194
195
                [geometry, respective_cell_id] = origin_id_to_geometry_map[origin_id]
            else:
                geometry = origin_id_to_geometry_map[origin_id]
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
196

197
198
            # Read fragility functions as numpy arrays.
            fragility_function = np.loadtxt(
199
200
201
                fragility_pathname + "/" + fragilityfunction_filename,
                delimiter=",",
                usecols=cls,
202
203
204
205
206
207
            )
            # Computing the ground-motion values from `full_ground_motion_field`. Please
            # note that `full_ground_motion_field` contains many identical lines (because
            # many of the assets have same locations and same location leads to the same
            # ground-motion value). This way, `full_ground_motion_field` has the same number
            # and order of assets as in the exposure model for simpler processing.
208
            gm_value = full_ground_motion_field[asset, taxonomy_to_fragility_map[taxonomy][1]]
209
210
            # Computing probabilities of exceedance and occurrence.
            [PoEs, PoOs] = losslib.get_PoEs(fragility_function, gm_value)
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
211

212
213
214
215
216
            # Compute damage by assets
            dmg_by_asset = [i * num_buildings for i in PoOs]
            for h in a:
                dmg_by_asset.insert(h, "")
            # Append results
217
            if exposure_type == "building":
218
219
220
221
222
                asset = [
                    geometry,
                    "",
                    origin_id,
                    "",
223
                    respective_cell_id,
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
                    "",
                    asset_id,
                    "",
                    lon,
                    "",
                    lat,
                    "",
                    taxonomy,
                    "",
                    gm_value,
                    "",
                    PoEs,
                    "",
                    PoOs,
                    "",
                    num_buildings,
                ]
            else:
                asset = [
                    geometry,
                    "",
                    origin_id,
                    "",
                    asset_id,
                    "",
                    lon,
                    "",
                    lat,
                    "",
                    taxonomy,
                    "",
                    gm_value,
                    "",
                    PoEs,
                    "",
                    PoOs,
                    "",
                    num_buildings,
                ]
            asset.extend(dmg_by_asset)
            csv_writer.writerow(asset)
Tara Evaz Zadeh's avatar
Tara Evaz Zadeh committed
265
    print("Execution time of the script", (datetime.datetime.now() - startTime))