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

Merge pull request 'July2025updates' (!17) from July2025updates into master

Browse Source 
Reviewed-on: #17
Reviewed-by: asia <asia@noreply.example.org>
This commit is contained in:
asia 2025-07-08 14:49:31 +02:00
commit fc2e8b53c7
1 changed files with 46 additions and 34 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

80
src/seismic_hazard_forecasting.py
View File

@ -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