From a68474ae146bb07c7330fd276ff3209bb3e9c3fb Mon Sep 17 00:00:00 2001 From: "manuela.villani@globalquakemodel.org" Date: Wed, 27 Sep 2023 16:33:43 +0100 Subject: [PATCH 1/7] adding asce7 --- .../calculators/postproc/compute_rtgm.py | 68 +++++++++++++++++-- 1 file changed, 62 insertions(+), 6 deletions(-) diff --git a/openquake/calculators/postproc/compute_rtgm.py b/openquake/calculators/postproc/compute_rtgm.py index 643c0e52d70b..2aa10c26382f 100644 --- a/openquake/calculators/postproc/compute_rtgm.py +++ b/openquake/calculators/postproc/compute_rtgm.py @@ -209,16 +209,66 @@ def get_deterministic(prob_mce, mag_dist_eps, sigma_by_src): return det.to_dict() -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 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]) + + 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['PGA'], + 'PGA_84th': det_mce['PGA'], + 'PGA': mce['PGA'], + + 'SS_RT': prob_mce['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['SA(1.0)'], + 'CR1': cr1, + 'S1_84th': det_mce['SA(1.0)'], + 'S1': mce['SA(1.0)'], + 'S1_seismicity': S1_seismicity, + } + + return mce, det_mce, asce7 def get_asce41(dstore, mce, facts): """ @@ -231,8 +281,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)'] @@ -291,9 +345,11 @@ def main(dstore, csm): dstore, csm, imts, imls_disagg) det_imt = get_deterministic(prob_mce, mag_dist_eps, sigma_by_src) logging.info(f'{det_imt=}') - mce, det_mce = get_mce(prob_mce, det_imt, DLLs) + mce, det_mce, asce7 = get_mce(prob_mce, det_imt, DLLs) logging.info(f'{mce=}') logging.info(f'{det_mce=}') asce41 = get_asce41(dstore, mce, facts) + dstore.create_df('asce41', asce41) + dstore.create_df('asce7', asce7) logging.info(asce41) From 6b9857d53ec0a425253d26e5238420393f68f826 Mon Sep 17 00:00:00 2001 From: "manuela.villani@globalquakemodel.org" Date: Wed, 27 Sep 2023 16:48:24 +0100 Subject: [PATCH 2/7] adding asce7 --- openquake/calculators/postproc/compute_rtgm.py | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/openquake/calculators/postproc/compute_rtgm.py b/openquake/calculators/postproc/compute_rtgm.py index 2aa10c26382f..caeef1701353 100644 --- a/openquake/calculators/postproc/compute_rtgm.py +++ b/openquake/calculators/postproc/compute_rtgm.py @@ -215,6 +215,8 @@ def get_mce_asce7(prob_mce, det_imt, DLLs, dstore): :returns: a dictionary imt -> det MCE :returns: a dictionary all ASCE7 parameters """ + + rtgm = dstore['rtgm'] imts = rtgm['IMT'] for i, imt in enumerate(imts): @@ -251,17 +253,17 @@ def get_mce_asce7(prob_mce, det_imt, DLLs, dstore): else: S1_seismicity = "Very High" - asce7 = {'PGA_2_50': prob_mce['PGA'], + asce7 = {'PGA_2_50': prob_mce[0], 'PGA_84th': det_mce['PGA'], 'PGA': mce['PGA'], - 'SS_RT': prob_mce['SA(0.2)'], + 'SS_RT': prob_mce[1], 'CRS': crs, 'SS_84th': det_mce['SA(0.2)'], 'SS': mce['SA(0.2)'], 'SS_seismicity': SS_seismicity, - 'S1_RT': prob_mce['SA(1.0)'], + 'S1_RT': prob_mce[2], 'CR1': cr1, 'S1_84th': det_mce['SA(1.0)'], 'S1': mce['SA(1.0)'], @@ -345,11 +347,11 @@ def main(dstore, csm): dstore, csm, imts, imls_disagg) det_imt = get_deterministic(prob_mce, mag_dist_eps, sigma_by_src) logging.info(f'{det_imt=}') - mce, det_mce, asce7 = get_mce(prob_mce, det_imt, DLLs) + mce, det_mce, asce7 = get_mce_asce7(prob_mce, det_imt, DLLs,dstore) logging.info(f'{mce=}') logging.info(f'{det_mce=}') asce41 = get_asce41(dstore, mce, facts) - dstore.create_df('asce41', asce41) - dstore.create_df('asce7', asce7) + logging.info(asce41) + logging.info(asce7) From 143f067c7e8122cc61c4b3b9ac26610843066d56 Mon Sep 17 00:00:00 2001 From: "manuela.villani@globalquakemodel.org" Date: Wed, 27 Sep 2023 17:30:08 +0100 Subject: [PATCH 3/7] adding asce7 --- .../calculators/postproc/compute_rtgm.py | 20 ++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/openquake/calculators/postproc/compute_rtgm.py b/openquake/calculators/postproc/compute_rtgm.py index caeef1701353..fe276486ff9b 100644 --- a/openquake/calculators/postproc/compute_rtgm.py +++ b/openquake/calculators/postproc/compute_rtgm.py @@ -210,13 +210,13 @@ def get_deterministic(prob_mce, mag_dist_eps, sigma_by_src): 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): @@ -227,11 +227,13 @@ def get_mce_asce7(prob_mce, det_imt, DLLs, dstore): 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]) + mce[imt] = min(prob_mce[i], det_mce[imt]) + prob_mce_out[imt] = prob_mce[i] - if mce['SA(0.2)']<0.25: + if mce['SA(0.2)'] < 0.25: SS_seismicity = "Low" elif mce['SA(0.2)'] <0.5: SS_seismicity = "Moderate" @@ -253,24 +255,24 @@ def get_mce_asce7(prob_mce, det_imt, DLLs, dstore): else: S1_seismicity = "Very High" - asce7 = {'PGA_2_50': prob_mce[0], + asce7 = {'PGA_2_50': prob_mce_out['PGA'], 'PGA_84th': det_mce['PGA'], 'PGA': mce['PGA'], - 'SS_RT': prob_mce[1], + '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[2], + '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 mce, det_mce, asce7 + return prob_mce_out, mce, det_mce, asce7 def get_asce41(dstore, mce, facts): """ @@ -347,11 +349,11 @@ def main(dstore, csm): dstore, csm, imts, imls_disagg) det_imt = get_deterministic(prob_mce, mag_dist_eps, sigma_by_src) logging.info(f'{det_imt=}') - mce, det_mce, asce7 = get_mce_asce7(prob_mce, det_imt, DLLs,dstore) + 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=}') asce41 = get_asce41(dstore, mce, facts) - + logging.info(asce41) logging.info(asce7) From 68d77ed3a4ddfa7e2272d4a93f1e2c38f49e090e Mon Sep 17 00:00:00 2001 From: "manuela.villani@globalquakemodel.org" Date: Thu, 28 Sep 2023 09:31:58 +0100 Subject: [PATCH 4/7] correcting typo in asce 41 --- openquake/calculators/postproc/compute_rtgm.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/openquake/calculators/postproc/compute_rtgm.py b/openquake/calculators/postproc/compute_rtgm.py index fe276486ff9b..fec62792e35a 100644 --- a/openquake/calculators/postproc/compute_rtgm.py +++ b/openquake/calculators/postproc/compute_rtgm.py @@ -231,7 +231,7 @@ def get_mce_asce7(prob_mce, det_imt, DLLs, dstore): 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]) - prob_mce_out[imt] = prob_mce[i] + prob_mce_out[imt] = prob_mce[i] # I think we should have the same number of IMTs everywhere as in the job.ini. For this first year we strat with 3 and then they'll become many later. if mce['SA(0.2)'] < 0.25: SS_seismicity = "Low" @@ -311,14 +311,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): From f3211f7956e38f9ce0c5499300d56100911ad8c0 Mon Sep 17 00:00:00 2001 From: Michele Simionato Date: Thu, 28 Sep 2023 11:45:48 +0200 Subject: [PATCH 5/7] Overridden oq.imtls with the AELO values --- openquake/engine/aelo.py | 6 ++++++ openquake/engine/tests/aelo_test.py | 4 ++-- 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/openquake/engine/aelo.py b/openquake/engine/aelo.py index 28b64828cd85..1c3579e71dbf 100644 --- a/openquake/engine/aelo.py +++ b/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' diff --git a/openquake/engine/tests/aelo_test.py b/openquake/engine/tests/aelo_test.py index 1178e9bd9644..f54a88a2d623 100644 --- a/openquake/engine/tests/aelo_test.py +++ b/openquake/engine/tests/aelo_test.py @@ -34,8 +34,8 @@ 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] +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(): From efd7e72647fb15995e13fe59d1def62cd5685a66 Mon Sep 17 00:00:00 2001 From: Michele Simionato Date: Thu, 28 Sep 2023 11:58:55 +0200 Subject: [PATCH 6/7] Fixed test --- openquake/calculators/export/__init__.py | 2 ++ openquake/calculators/export/hazard.py | 4 ++-- openquake/calculators/postproc/compute_rtgm.py | 5 +++-- openquake/engine/tests/aelo_test.py | 10 +++++++++- 4 files changed, 16 insertions(+), 5 deletions(-) diff --git a/openquake/calculators/export/__init__.py b/openquake/calculators/export/__init__.py index 7d885189563d..2802365bed4d 100644 --- a/openquake/calculators/export/__init__.py +++ b/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', diff --git a/openquake/calculators/export/hazard.py b/openquake/calculators/export/hazard.py index 09789c9d2fbb..265be72a239f 100644 --- a/openquake/calculators/export/hazard.py +++ b/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') diff --git a/openquake/calculators/postproc/compute_rtgm.py b/openquake/calculators/postproc/compute_rtgm.py index fe898824b7a3..11ebb282b080 100644 --- a/openquake/calculators/postproc/compute_rtgm.py +++ b/openquake/calculators/postproc/compute_rtgm.py @@ -356,11 +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=}') - prob_mce_out, mce, det_mce, asce7 = get_mce_asce7(prob_mce, det_imt, DLLs,dstore) + 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) - diff --git a/openquake/engine/tests/aelo_test.py b/openquake/engine/tests/aelo_test.py index f54a88a2d623..5800bc7d8f5b 100644 --- a/openquake/engine/tests/aelo_test.py +++ b/openquake/engine/tests/aelo_test.py @@ -34,7 +34,10 @@ 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, 1., 0.6, +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] @@ -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) From 16e6c1e0a64f0b2987eaaa1a664a4461422588a9 Mon Sep 17 00:00:00 2001 From: Michele Simionato Date: Thu, 28 Sep 2023 12:37:22 +0200 Subject: [PATCH 7/7] Cleanup --- openquake/calculators/postproc/compute_rtgm.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/openquake/calculators/postproc/compute_rtgm.py b/openquake/calculators/postproc/compute_rtgm.py index 11ebb282b080..96693bb878fa 100644 --- a/openquake/calculators/postproc/compute_rtgm.py +++ b/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') @@ -215,13 +216,11 @@ def get_deterministic(prob_mce, mag_dist_eps, sigma_by_src): 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): @@ -236,7 +235,7 @@ def get_mce_asce7(prob_mce, det_imt, DLLs, dstore): 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]) - prob_mce_out[imt] = prob_mce[i] # I think we should have the same number of IMTs everywhere as in the job.ini. For this first year we strat with 3 and then they'll become many later. + prob_mce_out[imt] = prob_mce[i] if mce['SA(0.2)'] < 0.25: SS_seismicity = "Low" @@ -279,6 +278,7 @@ def get_mce_asce7(prob_mce, det_imt, DLLs, dstore): return prob_mce_out, mce, det_mce, asce7 + def get_asce41(dstore, mce, facts): """ :returns: a dictionary with the ASCE-41 parameters