official-apps/SeismicHazardForecasting
10
0
Fork 1
Code Issues Pull Requests Activity

change search interval for PGV to be different than that for PGA/SA

Browse Source
This commit is contained in:
ftong
2025-09-09 10:56:35 +02:00
parent a8233950e1
commit 65759b86f1
1 changed files with 14 additions and 10 deletions
Download Patch File Download Diff File Expand all files Collapse all files

24
src/seismic_hazard_forecasting.py
View File

@@ -88,9 +88,7 @@ def main(catalog_file, mc_file, pdf_file, m_file, m_select, mag_label, mc, m_max
else: else:
logger.setLevel(logging.INFO) logger.setLevel(logging.INFO)
# temporary hard-coded configuration exclude_low_fxy = True # skip low probability areas of the map
# exclude_low_fxy = False
exclude_low_fxy = True
thresh_fxy = 1e-3 # minimum fxy value (location PDF) needed to do PGA estimation (to skip low probability areas); also should scale according to number of grid points thresh_fxy = 1e-3 # minimum fxy value (location PDF) needed to do PGA estimation (to skip low probability areas); also should scale according to number of grid points
# log user selections # log user selections
@@ -125,10 +123,6 @@ verbose: {verbose}")
logger.info("No magnitude label of catalog specified, therefore try Mw by default") logger.info("No magnitude label of catalog specified, therefore try Mw by default")
mag_label = 'Mw' mag_label = 'Mw'
# if cat_label == None:
# print("No magnitude label of catalog specified, therefore try 'Catalog' by default")
# cat_label='Catalog'
time, mag, lat, lon, depth = read_mat_cat(catalog_file, mag_label=mag_label, catalog_label='Catalog') time, mag, lat, lon, depth = read_mat_cat(catalog_file, mag_label=mag_label, catalog_label='Catalog')
# check for null magnitude values # check for null magnitude values
@@ -459,15 +453,19 @@ verbose: {verbose}")
if use_pp: # use dask parallel computing if use_pp: # use dask parallel computing
pbar = ProgressBar() pbar = ProgressBar()
pbar.register() pbar.register()
# iter = range(0,len(distances))
iter = indices iter = indices
iml_grid_raw = [] # raw ground motion grids iml_grid_raw = [] # raw ground motion grids
for imt in products: for imt in products:
logger.info(f"Estimating {imt}") logger.info(f"Estimating {imt}")
if imt == "PGV":
IMT_max = 200 # search interval max for velocity (cm/s)
else:
IMT_max = 2.0 # search interval max for acceleration (g)
imls = [dask.delayed(compute_IMT_exceedance)(rx_lat[i], rx_lon[i], distances[i].flatten(), fr, p, lambdas, imls = [dask.delayed(compute_IMT_exceedance)(rx_lat[i], rx_lon[i], distances[i].flatten(), fr, p, lambdas,
forecast_len, lambdas_perc, m_range, m_pdf, m_cdf, model, forecast_len, lambdas_perc, m_range, m_pdf, m_cdf, model,
log_level=logging.DEBUG, imt=imt, IMT_min=0.0, IMT_max=2.0, rx_label=i, log_level=logging.DEBUG, imt=imt, IMT_min=0.0, IMT_max=IMT_max, rx_label=i,
rtol=0.1, use_cython=True) for i in iter] rtol=0.1, use_cython=True) for i in iter]
iml = dask.compute(*imls) iml = dask.compute(*imls)
iml_grid_raw.append(list(iml)) iml_grid_raw.append(list(iml))
@@ -476,11 +474,17 @@ verbose: {verbose}")
iml_grid_raw = [] iml_grid_raw = []
iter = indices iter = indices
for imt in products: for imt in products:
if imt == "PGV":
IMT_max = 200 # search interval max for velocity (cm/s)
else:
IMT_max = 2.0 # search interval max for acceleration (g)
iml = [] iml = []
for i in iter: for i in iter:
iml_i = compute_IMT_exceedance(rx_lat[i], rx_lon[i], distances[i].flatten(), fr, p, lambdas, forecast_len, iml_i = compute_IMT_exceedance(rx_lat[i], rx_lon[i], distances[i].flatten(), fr, p, lambdas, forecast_len,
lambdas_perc, m_range, m_pdf, m_cdf, model, imt=imt, IMT_min = 0.0, lambdas_perc, m_range, m_pdf, m_cdf, model, imt=imt, IMT_min = 0.0,
IMT_max = 2.0, rx_label = i, rtol = 0.1, use_cython=True) IMT_max = IMT_max, rx_label = i, rtol = 0.1, use_cython=True)
iml.append(iml_i) iml.append(iml_i)
logger.info(f"Estimated {imt} at rx {i} is {iml_i}") logger.info(f"Estimated {imt} at rx {i} is {iml_i}")
iml_grid_raw.append(iml) iml_grid_raw.append(iml)