diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py index e1ad7ac..2213c35 100644 --- a/src/seismic_hazard_forecasting.py +++ b/src/seismic_hazard_forecasting.py @@ -17,7 +17,7 @@ def main(catalog_file, mc_file, pdf_file, m_file, m_select, mag_label, mc, m_max then the program looks for a label of 'Mw' for magnitude in the catalog mc: The magnitude of completeness (Mc) of the catalog m_max:M_max. The magnitude distribution is estimated for the range from Mc to M_max. If no value is provided, - then the program sets M_max to be 3 magnitude units above the maximum magnitude value in the catalog. + then the program sets M_max to be 1 magnitude units above the maximum magnitude value in the catalog. m_kde_method: The kernel density estimator to use. xy_select: If True, perform an estimation of the magnitude distribution using KDE with the chosen KDE method grid_dim: The grid cell size (in metres) of the final ground motion product map. A smaller cell size will @@ -45,20 +45,17 @@ def main(catalog_file, mc_file, pdf_file, m_file, m_select, mag_label, mc, m_max """ import sys + from importlib.metadata import version import logging from base_logger import getDefaultLogger from timeit import default_timer as timer from math import ceil, floor, isnan import numpy as np - import scipy - import obspy import dask from dask.diagnostics import ProgressBar # use Dask progress bar import kalepy as kale import utm from skimage.transform import resize - import psutil - import openquake.engine import igfash from igfash.io import read_mat_cat, read_mat_m, read_mat_pdf, read_csv from igfash.window import win_CTL, win_CNE @@ -105,25 +102,13 @@ verbose: {verbose}") # print key package version numbers logger.debug(f"Python version {sys.version}") - logger.debug(f"Numpy version {np.__version__}") - logger.debug(f"Scipy version {scipy.__version__}") - logger.debug(f"Obspy version {obspy.__version__}") - logger.debug(f"Openquake version {openquake.engine.__version__}") - logger.debug(f"Igfash version {igfash.__version__}") - - # print number of cpu cores available - ncpu = psutil.cpu_count(logical=False) - logger.debug(f"Number of cpu cores available: {ncpu}") - for process in psutil.process_iter(): - with process.oneshot(): - - # cpu = process.cpu_percent() - cpu = process.cpu_percent() / ncpu - - if cpu > 1: - logger.debug(f"{process.name()}, {cpu}") - - logger.debug(f"BASELINE CPU LOAD% {psutil.cpu_percent(interval=None, percpu=True)}") + logger.debug(f"Numpy version {version('numpy')}") + logger.debug(f"Scipy version {version('scipy')}") + logger.debug(f"Obspy version {version('obspy')}") + logger.debug(f"Openquake version {version('openquake.engine')}") + logger.debug(f"Igfash version {version('igfash')}") + logger.debug(f"Rbeast version {version('rbeast')}") + logger.debug(f"Dask version {version('dask')}") dask.config.set(scheduler='processes') @@ -145,6 +130,12 @@ verbose: {verbose}") time, mag, lat, lon, depth = read_mat_cat(catalog_file, mag_label=mag_label, catalog_label='Catalog') + # check for null magnitude values + m_null_idx = np.where(np.isnan(mag))[0] + if len(m_null_idx) > 0: + msg = f"There are null values in the magnitude column of the catalog at indices {m_null_idx}" + logger.error(msg) + raise Exception(msg) if mc != None: logger.info("Mc value provided by user") trim_to_mc = True @@ -166,9 +157,10 @@ verbose: {verbose}") lat = np.delete(lat, indices) lon = np.delete(lon, indices) - # if user does not provide a m_max, set m_max to 3 magnitude units above max magnitude in catalog + # if user does not provide a m_max, set m_max to 1 magnitude unit above max magnitude in catalog if m_max == None: - m_max = mag.max() + 3.0 + m_max = mag.max() + 1.0 + logger.info(f"No m_max was given. Therefore m_max is automatically set to: {m_max}") start = timer() @@ -233,8 +225,6 @@ verbose: {verbose}") y_min = y.min() x_max = x.max() y_max = y.max() - z_min = depth.min() - z_max = depth.max() grid_x_max = int(ceil(x_max / grid_dim) * grid_dim) grid_x_min = int(floor(x_min / grid_dim) * grid_dim) @@ -349,6 +339,12 @@ verbose: {verbose}") elif time_unit == 'years': multiplicator = 1 / 365 + # Raise an exception when time_win_duration from the user is too large relative to the catalog + if time_win_duration/multiplicator > 0.5*(time[-1] - time[0]): + msg = "Activity rate estimation time window must be less than half the catalog length. Use a shorter time window." + logger.error(msg) + raise Exception(msg) + # Selects dates in datenum format and procceeds to forecast value start_date = datenum_data[-1] - (2 * time_win_duration / multiplicator) dates_calc = [date for date in datenum_data if start_date <= date <= datenum_data[-1]] @@ -365,7 +361,7 @@ verbose: {verbose}") act_rate, bin_counts, bin_edges, out, pprs, rt, idx, u_e = calc_bins(np.array(datenum_data), time_unit, time_win_duration, dates_calc, rate_forecast, rate_unc_high, rate_unc_low, - multiplicator, quiet=True) + multiplicator, quiet=True, figsize=(14,9)) # Assign probabilities lambdas, lambdas_perc = lambda_probs(act_rate, dates_calc, bin_edges) @@ -377,18 +373,26 @@ verbose: {verbose}") if forecast_select: products = products_string.split() - logger.info( - f"Ground motion forecasting selected with ground motion model {model} and IMT products {products_string}") + logger.info(f"Ground motion forecasting selected with ground motion model {model} and IMT products {products_string}") + + # validate m_max against the grond motion model + models_anthro_limited = ['Lasocki2013', 'Atkinson2015', 'ConvertitoEtAl2012Geysers'] # these models require that m_max<=4.5 + if m_max > 4.5 and model in models_anthro_limited: + msg = f"Selected ground motion model {model} is only valid up to a maximum magnitude of 4.5. Please try again with a lower maximum magnitude." + logger.error(msg) + raise Exception(msg) if not xy_select: msg = "Event location distribution modeling was not selected; cannot continue..." logger.error(msg) raise Exception(msg) - elif m_pdf[0] == None: + + if m_pdf[0] == None: msg = "Magnitude distribution modeling was not selected and magnitude PDF file was not provided; cannot continue..." logger.error(msg) raise Exception(msg) - elif lambdas[0] == None: + + if lambdas[0] == None: msg = "Activity rate modeling was not selected and custom activity rate was not provided; cannot continue..." logger.error(msg) raise Exception(msg) @@ -426,7 +430,10 @@ verbose: {verbose}") rx_lat[i], rx_lon[i] = utm.to_latlon(x_rx[i], y_rx[i], utm_zone_number, utm_zone_letter) # get receiver location as lat,lon - # experimental - compute ground motion only at grid points that have minimum probability density of thresh_fxy + # convert distances from m to km because openquake ground motion models take input distances in kilometres + distances = distances/1000.0 + + # compute ground motion only at grid points that have minimum probability density of thresh_fxy if exclude_low_fxy: indices = list(np.where(fxy.flatten() > thresh_fxy)[0]) else: @@ -475,6 +482,11 @@ verbose: {verbose}") end = timer() logger.info(f"Ground motion exceedance computation time: {round(end - start, 1)} seconds") + if np.isnan(iml_grid_raw).all(): + msg = "No valid ground motion intensity measures were forecasted. Try a different ground motion model." + logger.error(msg) + raise Exception(msg) + # create list of one empty list for each imt iml_grid = [[] for _ in range(len(products))] # final ground motion grids iml_grid_prep = iml_grid.copy() # temp ground motion grids