Merge pull request #9016 from gem/years-store

Exporting AOP, EOP curves

openquake/calculators/extract.py

@@ -752,35 +752,39 @@ def extract_agg_curves(dstore, what):
         agg_id = lst.index(','.join(tagvalues))
     else:
         agg_id = 0  # total aggregation
+    ep_fields = dstore.get_attr('aggcurves', 'ep_fields')
     if qdic['rlzs']:
         [li] = qdic['loss_type']  # loss type index
         units = dstore.get_attr('aggcurves', 'units').split()
         df = dstore.read_df('aggcurves', sel=dict(agg_id=agg_id, loss_id=li))
         rps = list(df.return_period.unique())
         P = len(rps)
         R = len(qdic['kind'])
-        arr = numpy.zeros((P, R))
+        EP = len(ep_fields)
+        arr = numpy.zeros((R, P, EP))
         for rlz in df.rlz_id.unique():
-            # NB: df may contains zeros but there are no missing periods
-            # by construction (see build_aggcurves)
-            arr[:, rlz] = df[df.rlz_id == rlz].loss
+            for ep_field_idx, ep_field in enumerate(ep_fields):
+                # NB: df may contains zeros but there are no missing periods
+                # by construction (see build_aggcurves)
+                arr[rlz, :, ep_field_idx] = df[df.rlz_id == rlz][ep_field]
     else:
         name = 'agg_curves-stats/' + lts[0]
         shape_descr = hdf5.get_shape_descr(dstore.get_attr(name, 'json'))
         rps = list(shape_descr['return_period'])
         units = dstore.get_attr(name, 'units').split()
-        arr = dstore[name][agg_id, k].T  # shape P, R
+        arr = dstore[name][agg_id, k]  # shape (P, S, EP)
     if qdic['absolute'] == [1]:
         pass
     elif qdic['absolute'] == [0]:
         evalue, = dstore['agg_values'][agg_id][lts]
         arr /= evalue
     else:
         raise ValueError('"absolute" must be 0 or 1 in %s' % what)
-    attrs = dict(shape_descr=['return_period', 'kind'] + tagnames)
-    attrs['return_period'] = rps
+    attrs = dict(shape_descr=['kind', 'return_period', 'ep_field'] + tagnames)
     attrs['kind'] = qdic['kind']
+    attrs['return_period'] = rps
     attrs['units'] = units  # used by the QGIS plugin
+    attrs['ep_field'] = ep_fields
     for tagname, tagvalue in zip(tagnames, tagvalues):
         attrs[tagname] = [tagvalue]
     if tagnames:
@@ -1186,7 +1190,7 @@ def extract_mean_rates_by_src(dstore, what):
     [imt] = qdict['imt']
     [iml] = qdict['iml']
     [site_id] = qdict.get('site_id', ['0'])
-    site_id = int(site_id)    
+    site_id = int(site_id)
     imt_id = list(oq.imtls).index(imt)
     rates = dset[site_id, imt_id]
     L1, Ns = rates.shape

openquake/calculators/post_risk.py

@@ -52,21 +52,29 @@ def save_curve_stats(dstore):
     aggcurves_df = dstore.read_df('aggcurves')
     periods = aggcurves_df.return_period.unique()
     P = len(periods)
+    ep_fields = ['loss']
+    if 'loss_aep' in aggcurves_df:
+        ep_fields.append('loss_aep')
+    if 'loss_oep' in aggcurves_df:
+        ep_fields.append('loss_oep')
+    EP = len(ep_fields)
     for lt in oq.ext_loss_types:
         loss_id = scientific.LOSSID[lt]
-        out = numpy.zeros((K + 1, S, P))
+        out = numpy.zeros((K + 1, S, P, EP))
         aggdf = aggcurves_df[aggcurves_df.loss_id == loss_id]
         for agg_id, df in aggdf.groupby("agg_id"):
             for s, stat in enumerate(stats.values()):
                 for p in range(P):
-                    dfp = df[df.return_period == periods[p]]
-                    ws = weights[dfp.rlz_id.to_numpy()]
-                    ws /= ws.sum()
-                    out[agg_id, s, p] = stat(dfp.loss.to_numpy(), ws)
+                    for e, ep_field in enumerate(ep_fields):
+                        dfp = df[df.return_period == periods[p]]
+                        ws = weights[dfp.rlz_id.to_numpy()]
+                        ws /= ws.sum()
+                        out[agg_id, s, p, e] = stat(dfp[ep_field].to_numpy(),
+                                                    ws)
         stat = 'agg_curves-stats/' + lt
-        dstore.create_dset(stat, F64, (K + 1, S, P))
+        dstore.create_dset(stat, F64, (K + 1, S, P, EP))
         dstore.set_shape_descr(stat, agg_id=K+1, stat=list(stats),
-                               return_period=periods)
+                               return_period=periods, ep_fields=ep_fields)
         dstore.set_attrs(stat, units=units)
         dstore[stat][:] = out
 
@@ -198,15 +206,16 @@ def fix_dtypes(dic):
     fix_dtype(dic, F32, floatcolumns)
 
 
-def build_aggcurves(items, builder):
+def build_aggcurves(items, builder, aggregate_loss_curves_types):
     """
     :param items: a list of pairs ((agg_id, rlz_id, loss_id), losses)
     :param builder: a :class:`LossCurvesMapsBuilder` instance
     """
     dic = general.AccumDict(accum=[])
     for (agg_id, rlz_id, loss_id), data in items:
         year = data.pop('year', ())
-        curve = {kind: builder.build_curve(year, data[kind], rlz_id)
+        curve = {kind: builder.build_curve(
+                    year, kind, data[kind], aggregate_loss_curves_types, rlz_id)
                  for kind in data}
         for p, period in enumerate(builder.return_periods):
             dic['agg_id'].append(agg_id)
@@ -215,8 +224,8 @@ def build_aggcurves(items, builder):
             dic['return_period'].append(period)
             for kind in data:
                 # NB: kind be ['fatalities', 'losses'] in a scenario_damage test
-                c = curve[kind]['ep']
-                dic[kind].append(c[p])
+                for k, c in curve[kind].items():
+                    dic[k].append(c[p])
     return dic
 
 
@@ -289,13 +298,15 @@ def build_store_agg(dstore, rbe_df, num_events):
                     data['year'] = year[df.event_id.to_numpy()]
                 items.append([(agg_id, rlz_id, loss_id), data])
         dic = parallel.Starmap.apply(
-            build_aggcurves, (items, builder),
+            build_aggcurves, (items, builder, oq.aggregate_loss_curves_types),
             concurrent_tasks=oq.concurrent_tasks,
             h5=dstore.hdf5).reduce()
         fix_dtypes(dic)
+        ep_fields = ['loss' + suffix
+                     for suffix in oq.aggregate_loss_curves_types.split(',')]
         dstore.create_df('aggcurves', pandas.DataFrame(dic),
                          limit_states=' '.join(oq.limit_states),
-                         units=units)
+                         units=units, ep_fields=ep_fields)
     return aggrisk
 
 
@@ -342,13 +353,18 @@ def build_reinsurance(dstore, num_events):
                 years = year[df.index.to_numpy()]
             else:
                 years = ()
-            curve = {col: builder.build_curve(years, df[col].to_numpy(), rlzid)
+            curve = {col: builder.build_curve(
+                        years, col, df[col].to_numpy(),
+                        oq.aggregate_loss_curves_types,
+                        rlzid)
                      for col in columns}
             for p, period in enumerate(builder.return_periods):
                 dic['rlz_id'].append(rlzid)
                 dic['return_period'].append(period)
                 for col in curve:
-                    dic[col].append(curve[col]['ep'][p])
+                    for k, c in curve[col].items():
+                        dic[k].append(c[p])
+
     dstore.create_df('reinsurance-avg_portfolio', pandas.DataFrame(avg),
                      units=dstore['cost_calculator'].get_units(
                          oq.loss_types))

openquake/commonlib/oqvalidation.py

@@ -61,6 +61,14 @@
   Example: *aggregate_by = region, taxonomy*.
   Default: empty list
 
+aggregate_loss_curves_types:
+  Used for event-based risk and damage calculations, to estimate the aggregated
+  loss Exceedance Probability (EP) only or to also calculate (if possible) the
+  Occurrence Exceedance Probability (OEP) and/or the Aggregate Exceedance
+  Probability (AEP).
+  Example: *aggregate_loss_curves_types = ,_oep,_aep*.
+  Default: ,_oep,_aep
+
 reaggregate_by:
   Used to perform additional aggregations in risk calculations. Takes in
   input a proper subset of the tags in the aggregate_by option.
@@ -903,6 +911,13 @@ class OqParam(valid.ParamSet):
     hazard_imtls = {}
     override_vs30 = valid.Param(valid.positivefloat, None)
     aggregate_by = valid.Param(valid.namelists, [])
+    aggregate_loss_curves_types = valid.Param(
+        valid.Choice('',
+                     ',_oep',
+                     ',_aep',
+                     ',_oep,_aep',
+                     ',_aep,_oep'),
+        ',_oep,_aep')
     reaggregate_by = valid.Param(valid.namelist, [])
     amplification_method = valid.Param(
         valid.Choice('convolution', 'kernel'), 'convolution')

openquake/qa_tests_data/event_based_damage/case_11/expected/aggcurves-county-taxonomy.csv

@@ -1,17 +1,17 @@
-#,,,,,"generated_by='OpenQuake engine 3.17.0-git237fc3d19f', start_date='2023-05-25T11:44:59', checksum=572001388, risk_investigation_time=1.0, num_events=4, effective_time=200.0, limit_states='ds1 ds2 ds3 ds4'"
-county,taxonomy,return_period,loss_type,losses_value,losses_ratio
-Contra Costa,tax1,50,structural,9.41111E+02,6.27408E-03
-Contra Costa,tax1,100,structural,9.41111E+02,6.27408E-03
-Contra Costa,tax1,200,structural,9.41111E+02,6.27408E-03
-Contra Costa,tax3,50,structural,6.07736E+02,3.79835E-03
-Contra Costa,tax3,100,structural,6.07736E+02,3.79835E-03
-Contra Costa,tax3,200,structural,6.07736E+02,3.79835E-03
-Marin,tax1,50,structural,9.41111E+02,6.27408E-03
-Marin,tax1,100,structural,9.41111E+02,6.27408E-03
-Marin,tax1,200,structural,9.41111E+02,6.27408E-03
-Solano,tax1,50,structural,8.78371E+02,6.27408E-03
-Solano,tax1,100,structural,8.78371E+02,6.27408E-03
-Solano,tax1,200,structural,4.77793E+03,3.41281E-02
-Solano,tax2,50,structural,1.77647E+03,7.40194E-03
-Solano,tax2,100,structural,1.77647E+03,7.40194E-03
-Solano,tax2,200,structural,1.86057E+03,7.75239E-03
+#,,,,,,,,,"generated_by='OpenQuake engine 3.18.0-gitc771776529', start_date='2023-09-18T08:35:29', checksum=210043009, risk_investigation_time=1.0, num_events=4, effective_time=200.0, limit_states='ds1 ds2 ds3 ds4'"
+county,taxonomy,return_period,loss_type,loss_value,loss_ratio,loss_aep_value,loss_aep_ratio,loss_oep_value,loss_oep_ratio
+Contra Costa,tax1,50,structural,9.41111E+02,6.27408E-03,9.41111E+02,6.27408E-03,9.41111E+02,6.27408E-03
+Contra Costa,tax1,100,structural,9.41111E+02,6.27408E-03,9.41111E+02,6.27408E-03,9.41111E+02,6.27408E-03
+Contra Costa,tax1,200,structural,9.41111E+02,6.27408E-03,9.41111E+02,6.27408E-03,9.41111E+02,6.27408E-03
+Contra Costa,tax3,50,structural,6.07736E+02,3.79835E-03,6.07736E+02,3.79835E-03,6.07736E+02,3.79835E-03
+Contra Costa,tax3,100,structural,6.07736E+02,3.79835E-03,6.07736E+02,3.79835E-03,6.07736E+02,3.79835E-03
+Contra Costa,tax3,200,structural,6.07736E+02,3.79835E-03,6.07736E+02,3.79835E-03,6.07736E+02,3.79835E-03
+Marin,tax1,50,structural,9.41111E+02,6.27408E-03,9.41111E+02,6.27408E-03,9.41111E+02,6.27408E-03
+Marin,tax1,100,structural,9.41111E+02,6.27408E-03,9.41111E+02,6.27408E-03,9.41111E+02,6.27408E-03
+Marin,tax1,200,structural,9.41111E+02,6.27408E-03,9.41111E+02,6.27408E-03,9.41111E+02,6.27408E-03
+Solano,tax1,50,structural,8.78371E+02,6.27408E-03,8.78371E+02,6.27408E-03,8.78371E+02,6.27408E-03
+Solano,tax1,100,structural,8.78371E+02,6.27408E-03,8.78371E+02,6.27408E-03,8.78371E+02,6.27408E-03
+Solano,tax1,200,structural,4.77793E+03,3.41281E-02,4.77793E+03,3.41281E-02,4.77793E+03,3.41281E-02
+Solano,tax2,50,structural,1.77647E+03,7.40194E-03,1.77647E+03,7.40194E-03,1.77647E+03,7.40194E-03
+Solano,tax2,100,structural,1.77647E+03,7.40194E-03,1.77647E+03,7.40194E-03,1.77647E+03,7.40194E-03
+Solano,tax2,200,structural,1.86057E+03,7.75239E-03,1.86057E+03,7.75239E-03,1.86057E+03,7.75239E-03

openquake/qa_tests_data/event_based_damage/case_13/expected/aggcurves-.csv

@@ -1,5 +1,5 @@
-#,,,"generated_by='OpenQuake engine 3.13.0-git4df12b7ee8', start_date='2021-11-25T05:30:27', checksum=2445880374, risk_investigation_time=50.0, num_events=42, effective_time=1000.0, limit_states='LS1 LS2'"
-return_period,loss_type,losses_value,losses_ratio
+#,,,"generated_by='OpenQuake engine 3.18.0-gitc771776529', start_date='2023-09-18T08:35:31', checksum=361246452, risk_investigation_time=50.0, num_events=42, effective_time=1000.0, limit_states='LS1 LS2'"
+return_period,loss_type,loss_value,loss_ratio
 50,structural,2.35340E+02,4.60547E-04
 100,structural,2.43929E+04,4.77356E-02
 200,structural,6.38974E+04,1.25044E-01

openquake/qa_tests_data/event_based_damage/case_13/expected/b_damages.txt

@@ -1,4 +1,4 @@
-                rlz_id  return_period   losses
+                rlz_id  return_period     loss
 agg_id loss_id                                
 0      3             0             50      0.0
        3             0            100      8.0

openquake/qa_tests_data/event_based_damage/case_13/expected/c_damages.txt

@@ -1,4 +1,4 @@
-                rlz_id  return_period   losses
+                rlz_id  return_period     loss
 agg_id loss_id                                
 0      3             0             50      0.0
        3             0            100      8.0

openquake/qa_tests_data/event_based_damage/case_13/expected/d_damages.txt

@@ -1,4 +1,4 @@
-                rlz_id  return_period   losses
+                rlz_id  return_period     loss
 agg_id loss_id                                
 0      3             0             50      0.0
        3             0            100    125.0

openquake/qa_tests_data/event_based_damage/case_14/expected/aggcurves-NAME.csv

@@ -1,5 +1,5 @@
-#,,,,"generated_by='OpenQuake engine 3.17.0-git237fc3d19f', start_date='2023-05-25T11:45:13', checksum=2973522147, risk_investigation_time=1.0, num_events=481, effective_time=500.0, limit_states='slight moderate extreme complete'"
-NAME_1,return_period,loss_type,losses_value,losses_ratio
+#,,,,"generated_by='OpenQuake engine 3.18.0-gitc771776529', start_date='2023-09-18T08:35:57', checksum=1185162647, risk_investigation_time=1.0, num_events=481, effective_time=500.0, limit_states='slight moderate extreme complete'"
+NAME_1,return_period,loss_type,loss_value,loss_ratio
 Central,2,structural,0.00000E+00,0.00000E+00
 Central,5,structural,0.00000E+00,0.00000E+00
 Central,10,structural,1.35847E+08,6.38471E-03

openquake/qa_tests_data/event_based_damage/case_15/expected/aggcurves-.csv

@@ -1,7 +1,7 @@
-#,,,,,,"generated_by='OpenQuake engine 3.18.0-git3f93e21a24', start_date='2023-07-14T04:55:44', checksum=154880459, risk_investigation_time=1.0, num_events=16828, effective_time=10000.0, limit_states='moderate complete'"
-return_period,loss_type,rlz_id,fatalities_value,fatalities_ratio,losses_value,losses_ratio
-500,structural,0,2.78624E+02,5.45253E-02,1.26430E+06,1.06486E-01
-500,structural,1,4.24553E+02,8.30827E-02,1.81472E+06,1.52845E-01
-500,structural,2,4.95023E+02,9.68733E-02,1.85553E+06,1.56282E-01
-500,structural,3,6.23238E+02,1.21964E-01,2.61005E+06,2.19831E-01
-500,structural,4,3.60764E+02,7.05995E-02,1.91788E+06,1.61533E-01
+#,,,,,,,,,,"generated_by='OpenQuake engine 3.18.0-gitc771776529', start_date='2023-09-18T08:49:35', checksum=2353089205, risk_investigation_time=1.0, num_events=16828, effective_time=10000.0, limit_states='moderate complete'"
+return_period,loss_type,rlz_id,fatalities_value,fatalities_ratio,loss_value,loss_ratio,loss_aep_value,loss_aep_ratio,loss_oep_value,loss_oep_ratio
+500,structural,0,2.78624E+02,5.45253E-02,1.26430E+06,1.06486E-01,1.26430E+06,1.06486E-01,1.26430E+06,1.06486E-01
+500,structural,1,4.24553E+02,8.30827E-02,1.81472E+06,1.52845E-01,1.81472E+06,1.52845E-01,1.81472E+06,1.52845E-01
+500,structural,2,4.95023E+02,9.68733E-02,1.85553E+06,1.56282E-01,1.85553E+06,1.56282E-01,1.85553E+06,1.56282E-01
+500,structural,3,6.23238E+02,1.21964E-01,2.61005E+06,2.19831E-01,2.61005E+06,2.19831E-01,2.61005E+06,2.19831E-01
+500,structural,4,3.60764E+02,7.05995E-02,1.91788E+06,1.61533E-01,1.91788E+06,1.61533E-01,1.91788E+06,1.61533E-01

openquake/qa_tests_data/event_based_risk/case_1/expected/agg_curves3.csv

@@ -1,17 +1,45 @@
-+---------------+---------+--------+----------+---------+
-| return_period | kind    | policy | taxonomy | value   |
-+---------------+---------+--------+----------+---------+
-| 120           | rlz-000 | A      | RC       | 0.03578 |
-+---------------+---------+--------+----------+---------+
-| 240           | rlz-000 | A      | RC       | 0.10181 |
-+---------------+---------+--------+----------+---------+
-| 240           | rlz-001 | A      | RC       | 0.03919 |
-+---------------+---------+--------+----------+---------+
-| 480           | rlz-000 | A      | RC       | 0.23541 |
-+---------------+---------+--------+----------+---------+
-| 480           | rlz-001 | A      | RC       | 0.11574 |
-+---------------+---------+--------+----------+---------+
-| 960           | rlz-000 | A      | RC       | 0.54103 |
-+---------------+---------+--------+----------+---------+
-| 960           | rlz-001 | A      | RC       | 0.13667 |
-+---------------+---------+--------+----------+---------+
++---------+---------------+----------+--------+----------+---------+
+| kind    | return_period | ep_field | policy | taxonomy | value   |
++---------+---------------+----------+--------+----------+---------+
+| rlz-000 | 120           | loss     | A      | RC       | 0.03578 |
++---------+---------------+----------+--------+----------+---------+
+| rlz-000 | 120           | loss_oep | A      | RC       | 0.03578 |
++---------+---------------+----------+--------+----------+---------+
+| rlz-000 | 120           | loss_aep | A      | RC       | 0.03578 |
++---------+---------------+----------+--------+----------+---------+
+| rlz-000 | 240           | loss     | A      | RC       | 0.10181 |
++---------+---------------+----------+--------+----------+---------+
+| rlz-000 | 240           | loss_oep | A      | RC       | 0.10181 |
++---------+---------------+----------+--------+----------+---------+
+| rlz-000 | 240           | loss_aep | A      | RC       | 0.10181 |
++---------+---------------+----------+--------+----------+---------+
+| rlz-000 | 480           | loss     | A      | RC       | 0.23541 |
++---------+---------------+----------+--------+----------+---------+
+| rlz-000 | 480           | loss_oep | A      | RC       | 0.23541 |
++---------+---------------+----------+--------+----------+---------+
+| rlz-000 | 480           | loss_aep | A      | RC       | 0.23541 |
++---------+---------------+----------+--------+----------+---------+
+| rlz-000 | 960           | loss     | A      | RC       | 0.54103 |
++---------+---------------+----------+--------+----------+---------+
+| rlz-000 | 960           | loss_oep | A      | RC       | 0.54103 |
++---------+---------------+----------+--------+----------+---------+
+| rlz-000 | 960           | loss_aep | A      | RC       | 0.54103 |
++---------+---------------+----------+--------+----------+---------+
+| rlz-001 | 240           | loss     | A      | RC       | 0.03919 |
++---------+---------------+----------+--------+----------+---------+
+| rlz-001 | 240           | loss_oep | A      | RC       | 0.03919 |
++---------+---------------+----------+--------+----------+---------+
+| rlz-001 | 240           | loss_aep | A      | RC       | 0.03919 |
++---------+---------------+----------+--------+----------+---------+
+| rlz-001 | 480           | loss     | A      | RC       | 0.11574 |
++---------+---------------+----------+--------+----------+---------+
+| rlz-001 | 480           | loss_oep | A      | RC       | 0.11574 |
++---------+---------------+----------+--------+----------+---------+
+| rlz-001 | 480           | loss_aep | A      | RC       | 0.11574 |
++---------+---------------+----------+--------+----------+---------+
+| rlz-001 | 960           | loss     | A      | RC       | 0.13667 |
++---------+---------------+----------+--------+----------+---------+
+| rlz-001 | 960           | loss_oep | A      | RC       | 0.13667 |
++---------+---------------+----------+--------+----------+---------+
+| rlz-001 | 960           | loss_aep | A      | RC       | 0.13667 |
++---------+---------------+----------+--------+----------+---------+