AELO: added output ASCE-7

openquake/calculators/export/__init__.py

@@ -22,7 +22,9 @@
 
 
 DISPLAY_NAME = {
+    'asce7': 'ASCE-7 Parameters',
     'asce41': 'ASCE-41 Parameters',
+    'mag_dst_eps_sig': "(src, mag, dst, eps, sig, imt) Table",
     'asset_risk': 'Exposure + Risk',
     'gmf_data': 'Ground Motion Fields',
     'damages-rlzs': 'Asset Risk Distributions',

openquake/calculators/export/hazard.py

@@ -695,8 +695,8 @@ def export_rtgm(ekey, dstore):
     return [fname]
 
 
-@export.add(('asce41', 'csv'))
-def export_asce41(ekey, dstore):
+@export.add(('asce7', 'csv'), ('asce41', 'csv'))
+def export_asce(ekey, dstore):
     js = dstore[ekey[0]][()].decode('utf8')
     dic = json.loads(js)
     writer = writers.CsvWriter(fmt='%.5f')

openquake/calculators/postproc/compute_rtgm.py

@@ -132,6 +132,7 @@ def calc_rtgm_df(rtgm_haz, facts, oq):
     :param oq: OqParam instance
     """
     M = len(imts)
+    assert len(oq.imtls) == M
     riskCoeff, RTGM, UHGM, RTGM_max, MCE = (
         np.zeros(M), np.zeros(M), np.zeros(M), np.zeros(M), np.zeros(M))
     results = rtgmpy.BuildingCodeRTGMCalc.calc_rtgm(rtgm_haz, 'ASCE7')
@@ -214,16 +215,69 @@ def get_deterministic(prob_mce, mag_dist_eps, sigma_by_src):
     return det.to_dict(), np.array(mag_dist_eps_sig, dt)
 
 
-def get_mce(prob_mce, det_imt, DLLs):
+def get_mce_asce7(prob_mce, det_imt, DLLs, dstore):
     """
     :returns: a dictionary imt -> MCE
+    :returns: a dictionary imt -> det MCE
+    :returns: a dictionary all ASCE7 parameters
     """
+    rtgm = dstore['rtgm']
+    imts = rtgm['IMT']
+    for i, imt in enumerate(imts):
+        if imt == b'SA0P2':
+            crs = rtgm['RiskCoeff'][i]
+        elif imt == b'SA1P0':
+            cr1 = rtgm['RiskCoeff'][i]
+
     det_mce = {}
     mce = {}  # imt -> MCE
+    prob_mce_out = {}
     for i, imt in enumerate(det_imt):
         det_mce[imt] = max(det_imt[imt], DLLs[i])
-        mce[imt] = min(prob_mce[i], det_mce[imt])    
-    return mce, det_mce
+        mce[imt] = min(prob_mce[i], det_mce[imt]) 
+        prob_mce_out[imt] = prob_mce[i]
+
+    if mce['SA(0.2)'] < 0.25:
+        SS_seismicity = "Low"
+    elif mce['SA(0.2)'] <0.5:
+        SS_seismicity = "Moderate"
+    elif mce['SA(0.2)'] <1:
+        SS_seismicity = "Moderately High"
+    elif mce['SA(0.2)']  <1.5:
+        SS_seismicity = "High"
+    else:
+        SS_seismicity = "Very High"
+
+    if mce['SA(1.0)'] < 0.1:
+        S1_seismicity = "Low"
+    elif mce['SA(1.0)'] < 0.2:
+        S1_seismicity = "Moderate"
+    elif mce['SA(1.0)'] < 0.4:
+        S1_seismicity = "Moderately High"
+    elif mce['SA(1.0)']< 0.6:
+        S1_seismicity = "High"
+    else:
+        S1_seismicity = "Very High"
+
+    asce7 = {'PGA_2_50': prob_mce_out['PGA'],
+            'PGA_84th': det_mce['PGA'],
+            'PGA': mce['PGA'],
+
+            'SS_RT': prob_mce_out['SA(0.2)'],
+            'CRS': crs,
+            'SS_84th': det_mce['SA(0.2)'],
+            'SS': mce['SA(0.2)'],
+            'SS_seismicity': SS_seismicity,
+
+            'S1_RT': prob_mce_out['SA(1.0)'],
+            'CR1': cr1,
+            'S1_84th': det_mce['SA(1.0)'],
+            'S1': mce['SA(1.0)'],
+            'S1_seismicity': S1_seismicity,
+            }
+
+    return prob_mce_out, mce, det_mce, asce7
+
 
 def get_asce41(dstore, mce, facts):
     """
@@ -236,8 +290,12 @@ def get_asce41(dstore, mce, facts):
     imts = list(oq.imtls)
     sa02 = imts.index('SA(0.2)')
     sa10 = imts.index('SA(1.0)')
-    poe5_50 = poes.index(0.001025)  # NB: not existing in Japan!!
-    poe20_50 = poes.index(0.004453)  # NB: not existing in Japan!!
+    if int(oq.investigation_time) == 1:
+        poe5_50 = poes.index(0.001025)  
+        poe20_50 = poes.index(0.004453) 
+    elif int(oq.investigation_time) == 50:
+        poe5_50 = poes.index(0.05)  
+        poe20_50 = poes.index(0.2)  
 
     BSE2N_Ss = mce['SA(0.2)']
     Ss_5_50 = hmap[sa02, poe5_50] * fact['SA(0.2)']
@@ -258,14 +316,14 @@ def get_asce41(dstore, mce, facts):
             'BSE2E_Ss': BSE2E_Ss,
             'BSE1E_Ss': BSE1E_Ss,
             'Ss_20_50': Ss_20_50,
-            'BSE1E_Ss': BSE1E_Ss,
+            'BSE1N_Ss': BSE1N_Ss,
 
             'BSE2N_S1': BSE2N_S1,
             'S1_5_50': S1_5_50,
             'BSE2E_S1': BSE2E_S1,
             'BSE1E_S1': BSE1E_S1,
             'S1_20_50': S1_20_50,
-            'BSE1E_S1': BSE1E_S1}
+            'BSE1N_S1': BSE1N_S1}
 
 
 def main(dstore, csm):
@@ -298,10 +356,12 @@ def main(dstore, csm):
         prob_mce, mag_dist_eps, sigma_by_src)
     dstore['mag_dst_eps_sig'] = mag_dst_eps_sig
     logging.info(f'{det_imt=}')
-    mce, det_mce = get_mce(prob_mce, det_imt, DLLs)
+    prob_mce_out, mce, det_mce, asce7 = get_mce_asce7(
+        prob_mce, det_imt, DLLs,dstore)
     logging.info(f'{mce=}')
     logging.info(f'{det_mce=}')
+    dstore['asce7'] = hdf5.dumps(asce7)
     asce41 = get_asce41(dstore, mce, facts)
     dstore['asce41'] = hdf5.dumps(asce41)
     logging.info(asce41)
-
+    logging.info(asce7)

openquake/engine/aelo.py

@@ -29,6 +29,11 @@
 
 CDIR = os.path.dirname(__file__)  # openquake/engine
 
+IMTLS = '''\
+{"PGA": logscale(0.005, 3.00, 25),
+ "SA(0.2)": logscale(0.005, 9.00, 25),
+ "SA(1.0)": logscale(0.005, 3.60, 25)}
+'''
 
 def get_params_from(inputs, mosaic_dir=config.directory.mosaic_dir):
     """
@@ -47,6 +52,7 @@ def get_params_from(inputs, mosaic_dir=config.directory.mosaic_dir):
         params['description'] += ' (%(lon)s, %(lat)s)' % inputs
     params['ps_grid_spacing'] = '0.'  # required for disagg_by_src
     params['pointsource_distance'] = '100.'
+    params['intensity_measure_types_and_levels'] = IMTLS
     params['disagg_by_src'] = 'true'
     params['uniform_hazard_spectra'] = 'true'
     params['use_rates'] = 'true'

openquake/engine/tests/aelo_test.py

@@ -34,8 +34,11 @@
 
 SITES = ['far -90.071 16.60'.split(), 'close -85.071 10.606'.split()]
 EXPECTED = [[0.318191, 0.661792, 0.758443], [0.763373, 1.84959, 1.28969]]
-ASCE41 = [1.5, 1.45972, 1.45972, 0.83824, 0.83824,
-          0.6, 1.00811, 0.6, 0.4, 0.57329]
+ASCE7 = ['0.78388', '0.50000', '0.50000', '1.84959', '0.95914', '1.50000',
+         '1.50000', 'Very High', '1.28969', '0.95669', '0.60000', '0.60000',
+         'Very High']
+ASCE41 = [1.5, 1.45972, 1.45972, 0.83824, 0.83824, 1., 0.6,
+          1.00811, 0.6 , 0.4, 0.57329, 0.4]
 
 
 def test_CCA():
@@ -53,6 +56,11 @@ def test_CCA():
             aac(df.RTGM, expected, atol=1E-6)
 
     if rtgmpy:
+        # check asce7 exporter
+        [fname] = export(('asce7', 'csv'), calc.datastore)
+        df = pandas.read_csv(fname, skiprows=1)
+        numpy.testing.assert_equal(df.value.to_numpy(), ASCE7)
+
         # check asce41 exporter
         [fname] = export(('asce41', 'csv'), calc.datastore)
         df = pandas.read_csv(fname, skiprows=1)