|
13 | 13 | Define functions to handle impact_yearsets |
14 | 14 | """ |
15 | 15 |
|
16 | | -import copy |
17 | 16 | import logging |
18 | 17 |
|
19 | 18 | import numpy as np |
20 | 19 | from numpy.random import default_rng |
21 | 20 |
|
22 | 21 | import climada.util.dates_times as u_dt |
| 22 | +from climada.engine import Impact |
23 | 23 |
|
24 | 24 | LOGGER = logging.getLogger(__name__) |
25 | 25 |
|
@@ -127,19 +127,23 @@ def impact_yearset_from_sampling_vect( |
127 | 127 | # compute impact per sampled_year |
128 | 128 | imp_per_year = compute_imp_per_year(imp, sampling_vect) |
129 | 129 |
|
130 | | - # copy imp object as basis for the yimp object |
131 | | - yimp = copy.deepcopy(imp) |
132 | | - |
133 | 130 | if correction_fac: # adjust for sampling error |
134 | | - imp_per_year = imp_per_year.astype(float) / calculate_correction_fac( |
135 | | - imp_per_year, imp |
136 | | - ) |
| 131 | + correction_factor = calculate_correction_fac(imp_per_year, imp) |
| 132 | + imp_per_year = imp_per_year.astype(float) * correction_factor |
| 133 | + LOGGER.info("The correction factor is %s.", correction_factor) |
| 134 | + |
| 135 | + # create yearly impact object |
| 136 | + yimp = Impact( |
| 137 | + at_event=imp_per_year, |
| 138 | + date=u_dt.str_to_date([str(date) + "-01-01" for date in sampled_years]), |
| 139 | + event_name=np.arange(1, len(sampled_years) + 1), |
| 140 | + event_id=np.arange(1, len(sampled_years) + 1), |
| 141 | + unit=imp.unit, |
| 142 | + haz_type=imp.haz_type, |
| 143 | + frequency=np.ones(len(sampled_years)) / len(sampled_years), |
| 144 | + ) |
137 | 145 |
|
138 | | - # save calculations in yimp |
139 | | - yimp.at_event = imp_per_year |
140 | | - yimp.event_id = np.arange(1, len(sampled_years) + 1) |
141 | | - yimp.date = u_dt.str_to_date([str(date) + "-01-01" for date in sampled_years]) |
142 | | - yimp.frequency = np.full_like(yimp.at_event, 1.0 / len(sampled_years), dtype=float) |
| 146 | + yimp.aai_agg = np.sum(yimp.frequency * yimp.at_event) |
143 | 147 |
|
144 | 148 | return yimp |
145 | 149 |
|
@@ -301,10 +305,9 @@ def calculate_correction_fac(imp_per_year, imp): |
301 | 305 | The correction factor is calculated as imp_eai/yimp_eai |
302 | 306 | """ |
303 | 307 |
|
304 | | - yimp_eai = np.sum(imp_per_year) / len(imp_per_year) |
| 308 | + yimp_eai = np.mean(imp_per_year) |
305 | 309 | imp_eai = np.sum(imp.frequency * imp.at_event) |
306 | 310 | correction_factor = imp_eai / yimp_eai |
307 | | - LOGGER.info("The correction factor amounts to %s", (correction_factor - 1) * 100) |
308 | 311 |
|
309 | 312 | # if correction_factor > 0.1: |
310 | 313 | # tex = raw_input("Do you want to exclude small events?") |
|
0 commit comments