July2025updates #17
@ -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
|
||||
|
||||
Loading…
Reference in New Issue
Block a user