-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathimpact_calc_func.py
342 lines (285 loc) · 11.8 KB
/
impact_calc_func.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
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
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Useful functions for impact calculations.
@author: Pui Man (Mannie) Kam
"""
import os
import glob
import numpy as np
import pandas as pd
import json
from typing import Union, List, Tuple
from pathlib import Path
from climada.hazard import TCTracks
from climada.entity import ImpactFunc, ImpfTropCyclone, ImpactFuncSet
from climada.engine import Impact
# List of regions and the countries
iso3_to_basin = {'NA1': ['AIA', 'ATG', 'ARG', 'ABW', 'BHS', 'BRB', 'BLZ', 'BMU',
'BOL', 'CPV', 'CYM', 'CHL', 'COL', 'CRI', 'CUB', 'DMA',
'DOM', 'ECU', 'SLV', 'FLK', 'GUF', 'GRD', 'GLP', 'GTM',
'GUY', 'HTI', 'HND', 'JAM', 'MTQ', 'MEX', 'MSR', 'NIC',
'PAN', 'PRY', 'PER', 'PRI', 'SHN', 'KNA', 'LCA', 'VCT',
'SXM', 'SUR', 'TTO', 'TCA', 'URY', 'VEN', 'VGB', 'VIR'],
'NA2': ['CAN', 'USA'],
'NI': ['AFG', 'ARM', 'AZE', 'BHR', 'BGD', 'BTN', 'DJI', 'ERI',
'ETH', 'GEO', 'IND', 'IRN', 'IRQ', 'ISR', 'JOR', 'KAZ',
'KWT', 'KGZ', 'LBN', 'MDV', 'MNG', 'MMR', 'NPL', 'OMN',
'PAK', 'QAT', 'SAU', 'SOM', 'LKA', 'SYR', 'TJK', 'TKM',
'UGA', 'ARE', 'UZB', 'YEM'],
'OC1': ['ASM', 'COK', 'FJI', 'PYF', 'GUM', 'KIR', 'MHL', 'FSM',
'NRU', 'NCL', 'NIU', 'NFK', 'MNP', 'PLW', 'PNG', 'PCN',
'WSM', 'SLB', 'TLS', 'TKL', 'TON', 'TUV', 'VUT', 'WLF'],
'OC2': ['AUS', 'NZL'],
'SI': ['COM', 'COD', 'SWZ', 'MDG', 'MWI', 'MLI', 'MUS', 'MOZ',
'ZAF', 'TZA', 'ZWE'],
'WP1': ['KHM', 'IDN', 'LAO', 'MYS', 'THA', 'VNM'],
'WP2': ['PHL'],
'WP3': ['CHN'],
'WP4': ['HKG', 'JPN', 'KOR', 'MAC', 'TWN'],
'ROW': ['ALB', 'DZA', 'AND', 'AGO', 'ATA', 'AUT', 'BLR', 'BEL',
'BEN', 'BES', 'BIH', 'BWA', 'BVT', 'BRA', 'IOT', 'BRN',
'BGR', 'BFA', 'BDI', 'CMR', 'CAF', 'TCD', 'CXR', 'CCK',
'COG', 'HRV', 'CUW', 'CYP', 'CZE', 'CIV', 'DNK', 'EGY',
'GNQ', 'EST', 'FRO', 'FIN', 'FRA', 'ATF', 'GAB', 'GMB',
'DEU', 'GHA', 'GIB', 'GRC', 'GRL', 'GGY', 'GIN', 'GNB',
'HMD', 'VAT', 'HUN', 'ISL', 'IRL', 'IMN', 'ITA', 'JEY',
'KEN', 'PRK', 'XKX', 'LVA', 'LSO', 'LBR', 'LBY', 'LIE',
'LTU', 'LUX', 'MLT', 'MRT', 'MYT', 'MDA', 'MCO', 'MNE',
'MAR', 'NAM', 'NLD', 'NER', 'NGA', 'MKD', 'NOR', 'PSE',
'POL', 'PRT', 'ROU', 'RUS', 'RWA', 'REU', 'BLM', 'MAF',
'SPM', 'SMR', 'STP', 'SEN', 'SRB', 'SYC', 'SLE', 'SGP',
'SVK', 'SVN', 'SGS', 'SSD', 'ESP', 'SDN', 'SJM', 'SWE',
'CHE', 'TGO', 'TUN', 'TUR', 'UKR', 'GBR', 'UMI', 'ESH',
'ZMB', 'ALA']}
# RMSF optimised v_half for each region
v_half_per_region = {'NA1': 51.6,
'NA2': 84.1,
'NI': 41.3,
'OC1': 44.3,
'OC2': 47.4,
'SI': 40.8,
'WP1': 42.2,
'WP2': 46.7,
'WP3': 35.7,
'WP4': 93.1,
'ROW': 49.5}
def impf_set_exposed_pop(threshold: np.float64 = 32.92):
"""
Impact function set that estimate the number of people exposed
to a threshold of windspeed.
Parameters
----------
threshold : np.float64
Wind speed threshold that people are exposed to.
Default: 32.92 (Hurrane Cat. 1)
Returns
-------
impf_set : climada.entity.ImpactDuncSet
Impact function set that contains a step impact function.
"""
impf = ImpactFunc.from_step_impf((0,threshold, 100),
haz_type="TC")
impf_set = ImpactFuncSet()
impf_set.append(impf)
return(impf_set)
def impf_set_displacement(country: str):
"""
Impact function set that estimate the number of displacement. The shape of the
impact function depends on the countries and their respective region.
Details see Kam et al. (2024).
Parameters
----------
country : str
Single country in ISO3 alpha.
Returns
-------
impf_set : climada.entity.ImpactDuncSet
Impact function set that contains a displacement impact function.
"""
v_half = get_impf_v_half(country)
impf = ImpfTropCyclone.from_emanuel_usa(v_half=v_half)
impf_set = ImpactFuncSet()
impf_set.append(impf)
return(impf_set)
def get_impf_v_half(country: str):
"""
Get the impact function parameter v_half according to selected country (country_iso3).
Parameters
----------
country: str
Single country in ISO3 alpha.
Returns
-------
v_half : np.float64
The impact function parameter v_half for the respective country.
"""
# Get basin in which country_iso3 lies
basin = [key
for key, list_of_values in iso3_to_basin.items()
if country in list_of_values]
# Get impf_distr corresponding to basin
v_half = v_half_per_region[basin[0]]
return v_half
def round_to_previous_12h_utc(timestamp: pd.Timestamp):
"""
Rounding the time into 00 or 12 UTC
"""
# Ensure the timestamp is in UTC
if timestamp.tz is None:
utc_timestamp = timestamp.tz_localize('UTC')
else:
utc_timestamp = timestamp.tz_convert('UTC')
# Round down to the nearest hour
rounded = utc_timestamp.floor('H')
# Determine the previous 12-hour mark
if rounded.hour < 12:
return rounded.replace(hour=0)
else:
return rounded.replace(hour=12)
def get_forecast_times(current_timestamp: pd.Timestamp) -> Tuple[pd.Timestamp, pd.Timestamp]:
"""Get the current and previous forecast times."""
forecast_time = round_to_previous_12h_utc(current_timestamp)
previous_forecast_time = forecast_time - pd.Timedelta(hours=12)
return forecast_time, previous_forecast_time
def get_tc_wind_files(forecast_time: pd.Timestamp,
previous_forecast_time: pd.Timestamp,
tc_wind_dir: str) -> List[str]:
"""Get the list of TC wind files for the given forecast times."""
forecast_time_str = forecast_time.strftime('%Y-%m-%d_%HUTC')
file_pattern = os.path.join(tc_wind_dir, f"*{forecast_time_str}.hdf5")
tc_wind_files = glob.glob(file_pattern)
if not tc_wind_files:
print(f"No TC activities at {forecast_time_str}. Trying the previous forecast.")
forecast_time_str = previous_forecast_time.strftime('%Y-%m-%d_%HUTC')
file_pattern = os.path.join(tc_wind_dir, f"*{forecast_time_str}.hdf5")
tc_wind_files = glob.glob(file_pattern)
return forecast_time_str, tc_wind_files
def summarize_forecast(country_iso3: str,
forecast_time: str,
impact_type: str,
tc_haz: TCTracks,
tc_name: str,
impact: Impact):
"""
Summarizing forecast into a dictionary
"""
# check impact event length
## incase the TC ensemble number is less than 51
# get the array for each event
imp_at_event = _check_event_no(impact)
imp_summary_dict={
"countryISO3": country_iso3,
"hazardType": impact.haz_type,
"impactType": impact_type,
"initializationTime": forecast_time,
"eventName": tc_name,
"mean": np.mean(imp_at_event),
"median": np.median(imp_at_event),
"05perc": np.percentile(imp_at_event, 5),
"25perc": np.percentile(imp_at_event, 25),
"75perc": np.percentile(imp_at_event, 75),
"95perc": np.percentile(imp_at_event, 95),
"weatherModel": "ECMWF",
"impactUnit": "people"
}
return imp_summary_dict
def save_forecast_summary(save_dir: Union[str, Path],
forecast_summary: dict):
"""
Save the summary into a geoJSON feature collection.
"""
# Create a GeoJSON FeatureCollection structure
geojson_data = {
"type": "FeatureCollection",
"features": [
{
"type": "Feature",
"properties": forecast_summary,
"geometry": None # Set to None if there's no specific geometry
}
]
}
# Save the GeoJSON data to a file
with open(save_dir+make_save_filename(forecast_summary, save_file_type="summary"),
'w') as f:
json.dump(geojson_data, f, indent=4)
def make_save_filename(imp_summary_dict: dict,
save_file_type: str):
"""
Make a file name for saving
Parameters
----------
imp_summary_dict: dict
Impact forecast summary
save_file_type: str
indicator of which type of file is saved
Returns
-------
forecast_filename: str
File name for saving the files
"""
forecast_filename = (
f'impact-{save_file_type}_TC_ECMWF_ens_{imp_summary_dict["eventName"]}_{imp_summary_dict["initializationTime"]}'
f'_{imp_summary_dict["countryISO3"]}_{imp_summary_dict["impactType"]}.json'
)
return forecast_filename
def save_average_impact_geospatial_points(save_dir: Union[str, Path],
imp_summary_dict: dict,
impact: Impact,
include_zeros: bool = False):
"""
Save the average impact of each grid points into a geoJSON file.
Parameters
----------
save_dir: Union[str, Path],
Directory where the output is saved to.
imp_summary_dict: dict
Summary of the forecast.
impact: climada.engine.Impact
include_zeros: bool
Whether inclode grid points with impact equals to 0.
Default: False
"""
imp_gdf = impact._build_exp().gdf
if include_zeros:
imp_gdf.to_file(save_dir+make_save_filename(imp_summary_dict, save_file_type="gdf"))
else:
imp_gdf.drop(imp_gdf[imp_gdf['value'] == 0].index, inplace=True)
imp_gdf.to_file(save_dir+make_save_filename(imp_summary_dict, save_file_type="gdf"))
def save_impact_at_event(save_dir: Union[str, Path],
imp_summary_dict: dict,
impact: Impact):
"""
Save the impact of each event into a CSV file.
Parameters
----------
save_dir: Union[str, Path],
Directory where the output is saved to.
imp_summary_dict: dict
Summary of the forecast.
impact: climada.engine.Impact
"""
at_event_dict = {'ensemble_id': [], 'at_event': []}
for idx_event, event_id in enumerate(impact.event_id):
at_event_dict['ensemble_id'].append(event_id)
at_event_dict['at_event'].append(impact.at_event[idx_event])
df = pd.DataFrame(at_event_dict)
save_file_name = (
f'impact-at-event_TC_ECMWF_ens_{imp_summary_dict["eventName"]}_{imp_summary_dict["initializationTime"]}'
f'_{imp_summary_dict["countryISO3"]}_{imp_summary_dict["impactType"]}.csv'
)
df.to_csv(save_dir +save_file_name)
def _check_event_no(impact: Impact):
"""
For some TC events, there is less than 51 esemble. Hence, checke if the
impact.at_event lenth equals to 51. If not, fill the missing event as 0.
"""
imp_at_event = impact.at_event
if len(imp_at_event) == 51:
return imp_at_event
else:
no_missing_zeros = 51 - len(imp_at_event)
new_imp_at_event = np.pad(imp_at_event, pad_width=(0, no_missing_zeros),
mode='constant', constant_values=0)
return new_imp_at_event