diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py index 4ebc6ba..41b3871 100644 --- a/src/seismic_hazard_forecasting.py +++ b/src/seismic_hazard_forecasting.py @@ -88,9 +88,7 @@ def main(catalog_file, mc_file, pdf_file, m_file, m_select, mag_label, mc, m_max else: logger.setLevel(logging.INFO) - # temporary hard-coded configuration - # exclude_low_fxy = False - exclude_low_fxy = True + exclude_low_fxy = True # skip low probability areas of the map 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 @@ -125,10 +123,6 @@ verbose: {verbose}") logger.info("No magnitude label of catalog specified, therefore try Mw by default") 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') # check for null magnitude values @@ -459,15 +453,19 @@ verbose: {verbose}") if use_pp: # use dask parallel computing pbar = ProgressBar() pbar.register() - # iter = range(0,len(distances)) iter = indices iml_grid_raw = [] # raw ground motion grids for imt in products: 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, 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] iml = dask.compute(*imls) iml_grid_raw.append(list(iml)) @@ -476,11 +474,17 @@ verbose: {verbose}") iml_grid_raw = [] iter = indices 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 = [] for i in iter: 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, - 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) logger.info(f"Estimated {imt} at rx {i} is {iml_i}") iml_grid_raw.append(iml)