Update src/seismic_hazard_forecasting.py

use AOI from user interface

LICENCE.txt

@@ -0,0 +1,13 @@
+Copyright 2025 Institute of Geophysics, Polish Academy of Sciences
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+    http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.

gsim/lasocki_2013.py

@@ -0,0 +1,142 @@
+# -*- coding: utf-8 -*-
+# vim: tabstop=4 shiftwidth=4 softtabstop=4
+#
+# Copyright (C) 2013-2023 GEM Foundation
+#
+# OpenQuake is free software: you can redistribute it and/or modify it
+# under the terms of the GNU Affero General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# OpenQuake is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU Affero General Public License for more details.
+#
+# You should have received a copy of the GNU Affero General Public License
+# along with OpenQuake. If not, see <http://www.gnu.org/licenses/>.
+
+"""
+Module exports :class:`Lasocki2013`.
+"""
+import numpy as np
+from scipy.constants import g
+
+from openquake.hazardlib.gsim.base import GMPE, CoeffsTable
+from openquake.hazardlib import const
+from openquake.hazardlib.imt import PGA, PGV, SA
+
+def get_magnitude_energy(mag):
+    """
+    Converts magnitude to energy-based magnitude term.
+    """
+    return 1.15 + 1.96 * mag
+
+def get_distance_term(coeffs, repi):
+    """
+    Computes the distance term using the given GMPE coefficients.
+    """
+    R_h = np.sqrt(repi ** 2 + coeffs["c7"] ** 2)
+    return np.log10(R_h)
+
+def get_standard_deviation(coeffs, coeffs_cov, magE, repi):
+    """
+    Computes the standard deviation term.
+    """
+    Cb = np.array(list(coeffs_cov)).reshape(3,3)
+    R_h = np.sqrt(repi ** 2 + coeffs["c7"] ** 2)
+    X0 = np.array([1, magE[0], np.log10(R_h[0]**2)])
+    variance_term = np.sqrt(X0 @ Cb @ X0 + coeffs["sigma"]**2)
+    return variance_term
+
+class Lasocki2013(GMPE):
+    """
+    Implement equation developed by Lasocki in "REPORT ON THE ATTENUATION 
+    RELATIONS OF PEAK GROUND ACCELERATION AND SPECTRAL ORDINATES OF GROUND 
+    MOTION FOR MINING-INDUCED SEISMIC EVENTS IN THE REGION OF THE ŻELAZNY 
+    MOST REPOSITORY", 2009.
+    Equation coefficients provided for the random horizontal component
+    """
+    #: Supported tectonic region type is induced, given
+    #: that the equations have been derived for the LGCD mining area
+    DEFINED_FOR_TECTONIC_REGION_TYPE = const.TRT.INDUCED
+
+    #: Supported intensity measure types are spectral acceleration,
+    #: and peak ground acceleration
+    DEFINED_FOR_INTENSITY_MEASURE_TYPES = {PGA, SA}
+
+    #: Supported intensity measure component is random horizontal
+    #: :attr:`~openquake.hazardlib.const.IMC.RANDOM_HORIZONTAL`,
+    DEFINED_FOR_INTENSITY_MEASURE_COMPONENT = const.IMC.RANDOM_HORIZONTAL
+    # DEFINED_FOR_INTENSITY_MEASURE_COMPONENT = const.IMC.VERTICAL
+
+    #: Supported standard deviation type is total
+    DEFINED_FOR_STANDARD_DEVIATION_TYPES = {const.StdDev.TOTAL}
+
+    #: site params are not required
+    REQUIRES_SITES_PARAMETERS = set()
+
+    #: Required rupture parameter is magnitude
+    REQUIRES_RUPTURE_PARAMETERS = {'mag'}
+
+    #: Required distance measure is epicentral distance
+    #: see paragraph 'Predictor Variables', page 6.
+    REQUIRES_DISTANCES = {'repi'}
+
+    def compute(self, ctx: np.recarray, imts, mean, sig, tau, phi):
+        """
+        Computes mean ground motion values and standard deviations for Lasocki (2013) GMPE.
+        
+        Parameters:
+            repi (float): Epicentral distance (m)
+            mag (float): Earthquake magnitude
+            imts (list of str): List of intensity measure types (e.g., 'PHA', 'PVA')
+            mean (np.array): Array to store computed mean values
+            sig (np.array): Array to store computed standard deviations
+        """
+        # Loop through each IMT and compute values
+        for i, imt in enumerate(imts):            
+            C = self.COEFFS[imt]
+            C_cov = self.COEFFS_COV[imt]
+            mag = ctx.mag
+            repi = ctx.repi*1000.0 # Convert distance from km to m
+            # Compute magnitude energy term
+            magE = get_magnitude_energy(mag)
+
+            # Compute GMPE terms
+            mag_term = C['c1'] + C['c2'] * magE
+            dist_term = C['c5'] * get_distance_term(C, repi)
+            
+            # Compute mean ground motion
+            imean = mag_term + dist_term
+            mean_value = np.log((10 ** imean) / g)  # Convert to natural log scale and divide by g
+            mean[i] = mean_value
+
+            # Compute standard deviation
+            sigma = get_standard_deviation(C, C_cov, magE, repi)
+            sig[i] = sigma
+
+    #: coefficient table provided by report
+    COEFFS = CoeffsTable(sa_damping=5, table="""\
+    IMT  c1     c2      c5      c7    sigma     N        
+    pga  1.25   0.31    -1.34   558   0.196     1196
+    0.6  -3.86  0.55    -0.65   183   0.242     1194
+    1.0  -3.94  0.62    -0.66   308   0.262     1197
+    2.0  -1.14  0.47    -1.02   741   0.235     1195
+    5.0  0.99   0.31    -1.15   690   0.234     1206
+    10.0 3.06   0.27    -1.65   906   0.203     1192
+    20.0 2.62   0.27    -1.60   435   0.196     1191
+    50.0 2.09   0.27    -1.48   375   0.204     1191
+    """)
+    
+    COEFFS_COV = CoeffsTable(sa_damping=5, table="""\
+    IMT  Cb00      Cb01      Cb02      Cb10      Cb11      Cb12      Cb20      Cb21      Cb22        
+    pga  0.005586  -0.000376 -0.000752 -0.000376 0.000103  -0.000111 -0.000752 -0.000111 0.000440
+    0.6  0.007509  -0.000662 -0.000688 -0.000662 0.000161  -0.000154 -0.000688 -0.000154 0.000516
+    1.0  0.009119  -0.00075  -0.000948 -0.00075  0.000189  -0.000187 -0.000948 -0.000187 0.000657                 
+    2.0  0.008563  -0.000514 -0.001282 -0.000514 0.000147  -0.000164 -0.001282 -0.000164 0.000696         
+    5.0  0.008283  -0.000516 -0.001202 -0.000516 0.000145  -0.000161 -0.001202 -0.000161 0.000668
+   10.0  0.006904  -0.00036  -0.001145 -0.00036  0.000108  -0.000126 -0.001145 -0.000126 0.000578
+   20.0  0.005389  -0.000396 -0.000658 -0.000396 0.000104  -0.000107 -0.000658 -0.000107 0.000408
+   50.0  0.005874  -0.000449 -0.000678 -0.000449 0.000114  -0.000114 -0.000678 -0.000114 0.000428
+    """)

src/seismic_hazard_forecasting.py

@@ -1,9 +1,259 @@
 # -*- coding: utf-8 -*-
+from eqdist.rate import datenum_to_datetime
+import Rbeast as rb;
+from scipy.stats import bootstrap
+from matplotlib.dates import DateFormatter, AutoDateLocator
+from matplotlib.ticker import MultipleLocator
+import matplotlib.pyplot as plt
+import numpy as np
+
+global ncp_choice, tcp_max, torder_min, torder_max
+ncp_choice = 'default'
+tcp_max = 5
+torder_min = 0
+torder_max = 1
+#AOI_lat = np.array([51.48, 51.54])
+#AOI_lon = np.array([16.15, 16.24])
+#AOI_lat = np.array([None, None])
+#AOI_lon = np.array([None, None])
+
+def plot_results(act_rate, bin_edges, bin_edges_dt, rt, boundaries,
+                  bin_dur, unit, multiplicator,
+                  rate_forecast, rate_unc_high, rate_unc_low,
+                  datenum_data, mag_data):
+
+    end_date = bin_edges[-1]
+
+    fig, ax = plt.subplots(figsize=(14, 6))
+    ax.plot(bin_edges_dt[1:], act_rate, '-o', linewidth=2.5, markersize=6.5, label='Activity rate')
+
+    if rate_forecast is not None:
+        next_date = end_date + (bin_dur / multiplicator)
+        ax.plot(datenum_to_datetime(next_date), rate_forecast, 
+                'ro', label='Forecasted Rate', markersize=6.5)
+        ax.plot([bin_edges_dt[-1], datenum_to_datetime(next_date)], [act_rate[-1], rate_forecast], 'r-', linewidth=2.5)
+        ax.vlines(datenum_to_datetime(next_date), rate_unc_low, rate_unc_high, colors='r',
+                  linewidth=2, label='Bootstrap uncertainty')
+
+    ax.xaxis.set_major_locator(AutoDateLocator())
+    ax.xaxis.set_major_formatter(DateFormatter('%d-%b-%Y'))
+    plt.xticks(rotation=45)
+    plt.title(f'Activity rate (Time Unit: {unit}, Bin Duration: {bin_dur} {unit})',fontsize=18)
+    # plt.title(f'Activity rate (Bin Duration: {bin_dur} {unit})',fontsize=18)
+    plt.xlabel('Time (Bin Center Date)', fontsize=16)
+    ax.set_ylabel('Activity rate per selected time period',fontsize=16)
+    plt.grid(True)
+
+    if len(rt) > 0:
+        for i in range(len(rt)):
+            ax.plot(bin_edges_dt[1:][boundaries[i]:boundaries[i+1]],
+                    [rt[i]] * (boundaries[i+1] - boundaries[i]),
+                    linewidth=2, label=f'Rate period {i+1}')
+
+    # ---- Magnitude scatter on right y-axis ----
+    ax2 = ax.twinx()
+    event_dates = [datenum_to_datetime(d) for d in datenum_data]
+
+    #-------------extract magnitude bins from data---------------------
+    mags = np.array(mag_data)
+    min_mag = mags.min()
+    max_mag = mags.max()
+
+    low_thresh = int(np.floor(min_mag))
+    high_thresh = int(np.floor(max_mag))
+
+    thresholds = list(range(low_thresh, high_thresh + 1))
+
+    base_size = 15
+    size_step = 35
+
+    bins_def = []
+    for idx, t in enumerate(thresholds):
+        low = t
+        if idx < len(thresholds) - 1:
+            high = thresholds[idx + 1]
+            label = f'{low:.1f} \u2264 M < {high:.1f}'
+        else:
+            high = np.inf
+            label = f'M \u2265 {low:.1f}'
+        size = base_size + idx * size_step
+        bins_def.append((low, high, size, label))
+
+    for low, high, size, label in bins_def:
+        mask = (mags >= low) & (mags < high)
+        if np.any(mask):
+            sel_dates = [d for d, m in zip(event_dates, mask) if m]
+            sel_mags = mags[mask]
+            ax2.scatter(sel_dates, sel_mags, s=size,
+                        facecolor='purple', edgecolor='black',
+                        alpha=0.15, linewidth=1, label=label)
+
+    ax2.set_ylabel('Magnitude', color='purple',fontsize=16)
+    ax2.yaxis.set_major_locator(MultipleLocator(0.5))
+    ax2.yaxis.set_minor_locator(MultipleLocator(0.1))
+    ax2.spines['right'].set_color('purple')
+    ax2.tick_params(axis='y', colors='purple')
+
+    h1, l1 = ax.get_legend_handles_labels()
+    h2, l2 = ax2.get_legend_handles_labels()
+    handles = h1 + h2
+    labels = l1 + l2
+    n_legend = len(handles)
+
+    ncols = max(1, int(np.ceil(n_legend / 5)))  # ~5 entries per column
+    
+    #-------add 20% headroom above the data to make space for legend------
+    ymin, ymax = ax.get_ylim()
+    ax.set_ylim(ymin, ymax * 1.20) 
+    
+    ax.legend(handles, labels, loc='best',
+          ncol=ncols, borderaxespad=0,framealpha=0.7)
+
+    ax.set_zorder(ax2.get_zorder() + 1) # put scatter plot behind the line plot
+    ax.patch.set_visible(False)
+    
+    fig.tight_layout()
+    plt.savefig("activity_rate.png", dpi=600)
+    plt.show()
+
+def bootstrap_forecast(data):
+
+    window_data=data
+
+    if len(window_data) >= 5:
+        res95 = bootstrap((window_data,), np.mean, confidence_level=0.95,
+                          method='BCa', n_resamples=1000)
+    else:
+        res95 = bootstrap((window_data,), np.mean, confidence_level=0.95,
+                          method='BCa', n_resamples=int(len(window_data) ** len(window_data)))
+
+    forecast = np.mean(res95.bootstrap_distribution)
+    bca_conf95 = res95.confidence_interval
+    return forecast, bca_conf95
+
+def calc_rates(act_rate, cps):
+    """
+    Calculates mean activity rates between changepoints.
+    cps : sorted array of changepoint indices into act_rate
+    Returns rt (list of rates) and segment boundaries
+    """
+    boundaries = [0] + list(cps.astype(int)) + [len(act_rate)]
+    rt = [np.mean(act_rate[boundaries[i]:boundaries[i+1]]) 
+          for i in range(len(boundaries)-1)]
+    return rt, boundaries
+
+def apply_beast(act_rate):
+    """
+    Applies BEAST to the smmothed rate data using different smoothing windows.
+    Input
+    act_rate : The activity rate data array to smooth and apply BEAST.
+    Output
+    out : A list of BEAST results for each smoothed rate array.
+    prob : A list of probabilities and change points extracted from BEAST results.
+    """
+
+    mirror_len = int(np.ceil(0.20 * len(act_rate)))
+    left_mirror = act_rate[:mirror_len][::-1]
+    right_mirror = act_rate[-mirror_len:][::-1]
+    act_rate_mirrored = np.concatenate([left_mirror, act_rate, right_mirror])
+
+    mcmc_th = int(np.clip(np.ceil(len(act_rate) / 100), 2, 15))
+    beast_result = rb.beast(act_rate_mirrored, period=0,
+                            tcp_minmax=[0, tcp_max],
+                            torder_minmax=[torder_min, torder_max],
+                            tseg_minlength=2, mcmc_chains=10,
+                            mcmc_thin=mcmc_th, mcmc_seed=10)
+
+    # User-driven ncp selection
+    if ncp_choice == 'median':
+        ncp = beast_result.trend.ncp_median
+        if np.isnan(ncp) or ncp == 0:
+            return beast_result, np.array([])
+    elif ncp_choice == 'mode':
+        ncp = beast_result.trend.ncp_mode
+        if np.isnan(ncp) or ncp == 0:
+            return beast_result, np.array([])   
+    elif ncp_choice == 'pct90':
+        ncp = beast_result.trend.ncp_pct90
+        if np.isnan(ncp) or ncp == 0:
+            return beast_result, np.array([])
+    else:  # default: median with mode and pct90 fallback
+        ncp = beast_result.trend.ncp_median
+        if np.isnan(ncp) or ncp == 0:
+            ncp = beast_result.trend.mode
+        if np.isnan(ncp) or ncp == 0:
+            ncp = beast_result.trend.ncp_pct90
+        if np.isnan(ncp) or ncp == 0:    
+            return beast_result, np.array([])
+
+    ncp = int(ncp)
+    cps = beast_result.trend.cp[:ncp]
+
+    # Discard mirrored zone changepoints and correct indices
+    valid_mask = (cps > mirror_len) & (cps <= mirror_len + len(act_rate))
+    cps = cps[valid_mask] - mirror_len
+
+    return beast_result, np.sort(cps)
+
+def bins_and_beast(dates, unit, bin_dur, multiplicator):
+    start_date = dates.min()
+    end_date = dates.max()
+
+    valid_units = ['hours', 'days']
+    if unit not in valid_units:
+        unit = 'days'
+        bin_dur = 15
+    if (end_date - start_date) < 15 and unit == 'days':
+        unit = 'hours'
+        bin_dur = 12
+
+    bin_edges = [end_date]
+    while bin_edges[-1] > start_date:
+        bin_edges.append(bin_edges[-1] - (bin_dur / multiplicator))
+    bin_edges = bin_edges[::-1]
+
+    #-------Drop first bin or keep it if >80% of set duration------
+    first_width_days = bin_edges[1] - start_date
+    first_width_units = first_width_days * multiplicator
+
+    if first_width_units >= 0.8 * bin_dur:
+        bin_edges[0] = start_date  # edge of first bin is at data start
+    else:
+        bin_edges = bin_edges[1:]  # drop bin 0 (and its events)
+    
+    #------------Error if remaining bins are fewer than 2------------
+    if len(bin_edges) < 2:
+        raise ValueError(
+            f"Not enough data to form at least one full bin of duration "
+            f"{bin_dur} {unit}(s) after dropping the partial first bin "
+            f"({first_width_units:.2f} {unit}(s), below the 80% threshold). "
+            f"Try a shorter bin_dur or check your input data range."
+        )
+    
+    bin_edges_dt = [datenum_to_datetime(d) for d in bin_edges]
+    bin_counts, _ = np.histogram(dates, bins=bin_edges)
+
+    act_rate = [count / ((bin_edges[i + 1] - bin_edges[i]) * multiplicator / bin_dur)
+                for i, count in enumerate(bin_counts)]
+
+    out, cps = apply_beast(act_rate)
+
+    if len(cps) > 0:
+        rt, boundaries = calc_rates(act_rate, cps)
+        print(f'Changepoints detected at bins: {cps}')
+    else:
+        rt = []
+        boundaries = []
+        print('-----------------------------------------------------')
+        print('No changepoints detected by BEAST (Zhao et al., 2019)')
+        print('-----------------------------------------------------')
+
+    return act_rate, bin_counts, bin_edges, bin_edges_dt, out, cps, rt, boundaries, bin_dur, unit
 
 
 def main(catalog_file, mc_file, pdf_file, m_file, m_select, mag_label, mc, m_max,
          m_kde_method, xy_select, grid_dim, xy_win_method, rate_select, time_win_duration,
-         forecast_select, custom_rate, forecast_len, time_unit, model, products_string, verbose):
+         forecast_select, custom_rate, forecast_len, time_unit, AOI_extent, model, products_string, verbose):
     """
     Python application that reads an earthquake catalog and performs seismic hazard forecasting.
     Arguments:
@@ -17,7 +267,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
@@ -33,6 +283,8 @@ def main(catalog_file, mc_file, pdf_file, m_file, m_select, mag_label, mc, m_max
                     forecasting.
         forecast_len: Length of the forecast for seismic hazard assessment.
         time_unit: Times units for the inputs Time Window Duration, Custom Activity Rate, and Forecast Length.
+        AOI_extent: The forecast geographical area of interest specified as a latitude and longitude range in decimal degrees
+                    in the form [lat_min, lat_max, lon_min, lon_max].
         model: Select from the following ground motion models available. Other models in the Openquake library are
                     available but have not yet been tested.
         products_string: The ground motion intensity types to output. Use a space between names to select more than
@@ -45,22 +297,18 @@ 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.io import read_mat_cat, read_mat_m, read_mat_mc, read_mat_pdf, read_csv
     from igfash.window import win_CTL, win_CNE
     import igfash.kde as kde
     from igfash.gm import compute_IMT_exceedance
@@ -72,6 +320,7 @@ def main(catalog_file, mc_file, pdf_file, m_file, m_select, mag_label, mc, m_max
     from matplotlib.contour import ContourSet
     import xml.etree.ElementTree as ET
     import json
+    import multiprocessing as mp
 
     logger = getDefaultLogger('igfash')
 
@@ -84,17 +333,19 @@ def main(catalog_file, mc_file, pdf_file, m_file, m_select, mag_label, mc, m_max
         m_range = [None]
     else:
         m_range = read_mat_m(m_file)
+        m_max = m_range[-1] # take m_max from the m_file
 
     if verbose:
         logger.setLevel(logging.DEBUG)
     else:
         logger.setLevel(logging.INFO)
 
-    # temporary hard-coded configuration
-    # exclude_low_fxy = False
-    exclude_low_fxy = True
+    exclude_low_fxy = False # 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
 
+    AOI_lat = np.array(AOI_extent[:2]) 
+    AOI_lon = np.array(AOI_extent[2:])
+    
     # log user selections
     logger.debug(f"User input files\n Catalog: {catalog_file}\n Mc: {mc_file}\n Mag_PDF: {pdf_file}\n Mag: {m_file}")
     logger.debug(
@@ -105,25 +356,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')
 
@@ -139,18 +378,20 @@ 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
+        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
         elif mc_file != None:
             logger.info("Mc estimation output file provided; selecting largest Mc from the list")
-            mc = read_mc(mc_file)
+            mc = read_mat_mc(mc_file)
             trim_to_mc = True
         else:
             logger.info("No Mc provided; using all magnitudes from the catalog")
@@ -166,9 +407,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()
 
@@ -228,21 +470,37 @@ verbose: {verbose}")
         utm_zone_letter = u[3]
         logger.debug(f"Latitude / Longitude coordinates correspond to UTM zone {utm_zone_number}{utm_zone_letter}")
 
-        # define corners of grid based on global dataset
-        x_min = x.min()
-        y_min = y.min()
-        x_max = x.max()
-        y_max = y.max()
-        z_min = depth.min()
-        z_max = depth.max()
+        if (None not in AOI_lat) and (None not in AOI_lon):
+            use_AOI = True
+            logger.info(f"Area of Interest (AOI) selected with latitutde range {AOI_lat} and longitude range {AOI_lon}")
+            #convert AOI to UTM
+            u_AOI = utm.from_latlon(AOI_lat, AOI_lon) 
+            x_AOI = u_AOI[0]
+            y_AOI = u_AOI[1]
+            
+            # make sure grid contains the user's AOI
+            x_min = np.concatenate((x, x_AOI)).min()
+            y_min = np.concatenate((y, y_AOI)).min()
+            x_max = np.concatenate((x, x_AOI)).max()
+            y_max = np.concatenate((y, y_AOI)).max()
+            
+            exclude_low_fxy = False # don't exclude any points because we need to analyze all grid points in the AOI
+            
+        else:
+            use_AOI = False
+        
+            # define corners of grid based on global dataset
+            x_min = x.min()
+            y_min = y.min()
+            x_max = x.max()
+            y_max = y.max()
 
         grid_x_max = int(ceil(x_max / grid_dim) * grid_dim)
         grid_x_min = int(floor(x_min / grid_dim) * grid_dim)
         grid_y_max = int(ceil(y_max / grid_dim) * grid_dim)
         grid_y_min = int(floor(y_min / grid_dim) * grid_dim)
 
-        grid_lat_max, grid_lon_max = utm.to_latlon(grid_x_max, grid_y_max, utm_zone_number, utm_zone_letter)
-        grid_lat_min, grid_lon_min = utm.to_latlon(grid_x_min, grid_y_min, utm_zone_number, utm_zone_letter)
+       
 
         # rectangular grid
         nx = int((grid_x_max - grid_x_min) / grid_dim) + 1
@@ -257,17 +515,23 @@ verbose: {verbose}")
             nx = ny
             grid_x_max = int(grid_x_min + (nx - 1) * grid_dim)
 
+        # update grid extent in lat/lon
+        grid_lat_max, grid_lon_max = utm.to_latlon(grid_x_max, grid_y_max, utm_zone_number, utm_zone_letter)
+        grid_lat_min, grid_lon_min = utm.to_latlon(grid_x_min, grid_y_min, utm_zone_number, utm_zone_letter)
+        
         # new x and y range
         x_range = np.linspace(grid_x_min, grid_x_max, nx)
         y_range = np.linspace(grid_y_min, grid_y_max, ny)
 
+        logger.debug(f"Grid X range: {x_range}, Y range: {y_range}")
+        
         t_windowed = time
         r_windowed = [[x, y]]
 
         # %% compute KDE and extract PDF
         start = timer()
 
-        if xy_win_method == "TW":
+        if xy_win_method:
             logger.info("Time weighting function selected")
 
             x_weights = np.linspace(0, 15, len(t_windowed))
@@ -328,18 +592,21 @@ verbose: {verbose}")
 
     # run activity rate modeling
     lambdas = [None]
-    if custom_rate != None and forecast_select:
+    if custom_rate != None and forecast_select and not rate_select:
         logger.info(f"Using activity rate specified by user: {custom_rate} per {time_unit}")
-        lambdas = [custom_rate]
-        lambdas_perc = [1]
+        lambdas = np.array([custom_rate], dtype='d')
+        lambdas_perc = np.array([1], dtype='d')
 
     elif rate_select:
         logger.info(f"Activity rate modeling selected")
+                
+        datenum_data, mag_data, lat_dummy, lon_dummy, depth_dummy = read_mat_cat(catalog_file, mag_label=mag_label, output_datenum=True)
 
-        time, mag_dummy, lat_dummy, lon_dummy, depth_dummy = read_mat_cat(catalog_file, output_datenum=True)
-
-        datenum_data = time  # REMEMBER THE DECIMAL DENOTES DAYS
-
+        if trim_to_mc:
+            indices = np.argwhere(mag_data < mc)
+            mag_data = np.delete(mag_data, indices)
+            datenum_data = np.delete(datenum_data, indices)
+        
         if time_unit == 'hours':
             multiplicator = 24
         elif time_unit == 'days':
@@ -349,46 +616,79 @@ verbose: {verbose}")
         elif time_unit == 'years':
             multiplicator = 1 / 365
 
-        # 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]]
-        forecasts, bca_conf95, rate_mean = bootstrap_forecast_rolling(dates_calc, multiplicator)
+        # 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)
 
-        # FINAL VALUES OF RATE AND ITS UNCERTAINTY IN THE 5-95 PERCENTILE
-        unc_bca05 = [ci.low for ci in bca_conf95];
-        unc_bca95 = [ci.high for ci in bca_conf95]
-        rate_unc_high = multiplicator / np.array(unc_bca05);
-        rate_unc_low = multiplicator / np.array(unc_bca95);
-        rate_forecast = multiplicator / np.median(forecasts)  # [per time unit]
+        #-----------data are sorted in case they were not-----------------
+        sorted_pairs = sorted(zip(datenum_data, mag_data), key=lambda x: x[0])
+        datenum_data, mag_data = map(list, zip(*sorted_pairs))
 
-        # Plot of forecasted activity rate with previous binned activity rate
-        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)
+        #-------split the data into bins and apply BEAST for changepoint detection--------------------
+        act_rate, bin_counts, bin_edges, bin_edges_dt, out, cps, rt, boundaries, bin_dur, time_unit = bins_and_beast(
+            np.array(datenum_data), time_unit, time_win_duration, multiplicator)
+        
+        #------Forecasted rate is taken from BEAST or is equal to last value if no changepoints detected-----
+        if len(cps) > 0:
+            rate_forecast = rt[-1]
+            last_cp_bin = int(cps[-1])
+        else:
+            rate_forecast = act_rate[-1]
+            last_cp_bin = len(act_rate) - 1
 
-        # Assign probabilities 
-        lambdas, lambdas_perc = lambda_probs(act_rate, dates_calc, bin_edges)
+        last_cp_datenum = bin_edges[last_cp_bin]
+        dates_calc = [date for date in datenum_data if last_cp_datenum <= date <= datenum_data[-1]]
+        interevent_times = np.diff(dates_calc)
 
-        # print("Forecasted activity rates: ", lambdas, "events per", time_unit[:-1])
-        logger.info(f"Forecasted activity rates: {lambdas} events per {time_unit} with percentages {lambdas_perc}")
-        np.savetxt('activity_rate.csv', np.vstack((lambdas, lambdas_perc)).T, header="lambda, percentage",
+        #------------Use BCa for uncertainty intervals-----------------
+        forecast, bca_conf95 = bootstrap_forecast(interevent_times)
+        rate_unc_high = bin_dur / (bca_conf95.low * multiplicator)
+        rate_unc_low = bin_dur / (bca_conf95.high * multiplicator)
+
+        #----------------------Plot------------------------------------
+        plot_results(act_rate, bin_edges, bin_edges_dt, rt, boundaries,
+                      bin_dur, time_unit, multiplicator,
+                      rate_forecast, rate_unc_high, rate_unc_low,
+                      datenum_data, mag_data)
+
+        logger.info("\n----------------- Forecast Summary -----------------")
+        logger.info(f"Forecasted activity rate (next {bin_dur} {time_unit}(s)): {rate_forecast:.4f}")
+        logger.info(f"95% BCa confidence interval: [{rate_unc_low:.4f}, {rate_unc_high:.4f}]")
+        logger.info("------------------------------------------------------")
+        
+        lambdas = np.array([rate_forecast/bin_dur], dtype='d')
+        lambdas_perc = np.array([1], dtype='d')
+        
+        np.savetxt('activity_rate.csv', lambdas, header=f"Activity Rate (Events per {time_unit[:-1]})",
                    delimiter=',', fmt='%1.4f')
 
     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:
+            if m_file is None: 
+                msg = f"The selected ground motion model {model} is only valid for magnitudes up to 4.5. Please select a lower maximum magnitude."
+            else:
+                msg = f"The selected ground motion model {model} is only valid for magnitudes up to 4.5, but the provided magnitude file includes values up to {m_max}. Please adjust the magnitude range in the file accordingly."
+            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 == None:
             msg = "Activity rate modeling was not selected and custom activity rate was not provided; cannot continue..."
             logger.error(msg)
             raise Exception(msg)
@@ -406,18 +706,22 @@ verbose: {verbose}")
         fxy = xy_kde[0]
         logger.debug(f"Normalization check; sum of all f(x,y) values = {np.sum(fxy)}")
 
-        xx, yy = np.meshgrid(x_range, y_range, indexing='ij')  # grid points
+        xx, yy = np.meshgrid(x_range, y_range)  # grid points
 
         # set every grid point to be a receiver
+        grid_shape = xx.shape
         x_rx = xx.flatten()
         y_rx = yy.flatten()
 
+        num_points = x_rx.size
+        distances = np.zeros(shape=(num_points, grid_shape[0], grid_shape[1]))
+        
         # compute distance matrix for each receiver
-        distances = np.zeros(shape=(nx * ny, nx, ny))
+        #distances = np.zeros(shape=(nx * ny, nx, ny))
         rx_lat = np.zeros(nx * ny)
         rx_lon = np.zeros(nx * ny)
 
-        for i in range(nx * ny):
+        for i in range(num_points):
             # Compute the squared distances directly using NumPy's vectorized operations
             squared_distances = (xx - x_rx[i]) ** 2 + (yy - y_rx[i]) ** 2
             distances[i] = np.sqrt(squared_distances)
@@ -426,48 +730,67 @@ 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:
-            indices = range(0, len(distances))
+            indices = np.arange(num_points)
 
+        if use_AOI:
+            # Filter out receivers outside the AOI; Find indices where values are OUTSIDE the AOI
+            indices_outside_x = np.where((x_rx < x_AOI[0]) | (x_rx > x_AOI[1]))[0]
+            indices_outside_y = np.where((y_rx < y_AOI[0]) | (y_rx > y_AOI[1]))[0]      
+            indices_outside_AOI = np.unique(np.concatenate((indices_outside_x, indices_outside_y)))
+            indices_filtered = np.setdiff1d(indices, indices_outside_AOI)
+        else:
+            indices_filtered = indices
+        
         fr = fxy.flatten()
 
         # For each receiver compute estimated ground motion values
-        logger.info(f"Estimating ground motion intensity at {len(indices)} grid points...")
-
-        PGA = np.zeros(shape=(nx * ny))
+        logger.info(f"Estimating ground motion intensity at {len(indices_filtered)} grid points...")
 
         start = timer()
 
-        use_pp = False
+        use_pp = True
 
         if use_pp:         # use dask parallel computing
-            pbar = ProgressBar()
-            pbar.register()
-            # iter = range(0,len(distances))
-            iter = indices
+            mp.set_start_method("fork", force=True)
+            iter = indices_filtered
             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,
-                                                            rtol=0.1, use_cython=False) for i in iter]
+                                                            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))
 
         else:
             iml_grid_raw = []
-            iter = indices
+            iter = indices_filtered
             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=False)
+                                                    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)
@@ -475,58 +798,91 @@ verbose: {verbose}")
         end = timer()
         logger.info(f"Ground motion exceedance computation time: {round(end - start, 1)} seconds")
 
+        logger.debug(f"IMT values: {iml_grid_raw[0]}")
+        
+        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
 
-        if exclude_low_fxy:
-            for i in range(0, len(distances)):
-                if i in indices:
-                    for j in range(0, len(products)):
-                        iml_grid_prep[j].append(iml_grid_raw[j].pop(0))
-                else:
-                    list(map(lambda lst: lst.append(np.nan),
-                             iml_grid_prep))  # use np.nan to indicate grid point excluded
-        else:
-            iml_grid_prep = iml_grid_raw
+        #if use_AOI or exclude_low_fxy:
 
         for j in range(0, len(products)):
-            vmin = min(x for x in iml_grid_prep[j] if x is not np.nan)
-            vmax = max(x for x in iml_grid_prep[j] if x is not np.nan)
-            iml_grid[j] = np.reshape(iml_grid_prep[j], (nx, ny)).astype(
-                dtype=np.float64)  # this reduces values to 8 decimal places
-            iml_grid_tmp = np.nan_to_num(iml_grid[j])  # change nans to zeroes
+            # Reassemble the grid cleanly using the original shape
+            # Initialize a flat array filled entirely with NaNs
+            iml_grid_flat = np.full(num_points, np.nan, dtype=np.float64)
+            
+            # Assign the computed values to their exact original 1D index positions
+            iml_grid_flat[indices_filtered] = iml_grid_raw[j]
+            
+            # Reshape back using the exact shape of your original xx/yy grids
+            iml_grid_prep[j] = iml_grid_flat.reshape(grid_shape)
 
-            # upscale the grid
-            up_factor = 4
-            iml_grid_hd = resize(iml_grid_tmp, (up_factor * len(iml_grid_tmp), up_factor * len(iml_grid_tmp)),
-                                 mode='reflect', anti_aliasing=False)
-            iml_grid_hd[iml_grid_hd == 0.0] = np.nan  # change zeroes back to nan
+            
+            #for i in indices:
+            #   if i in indices_filtered:
+            #        for j in range(0, len(products)):
+            #            iml_grid_prep[j].append(iml_grid_raw[j].pop(0))
+            #    else:
+            #        list(map(lambda lst: lst.append(np.nan),
+            #                 iml_grid_prep))  # use np.nan to indicate grid point excluded
+        #elif exclude_low_fxy:
+        #    for i in range(0, len(distances)):
+        #        if i in indices:
+        #            for j in range(0, len(products)):
+        #                iml_grid_prep[j].append(iml_grid_raw[j].pop(0))
+        #        else:
+        #            list(map(lambda lst: lst.append(np.nan),
+        #                     iml_grid_prep))  # use np.nan to indicate grid point excluded
+        #else:
+        #    iml_grid_prep = iml_grid_raw
 
-            # trim edges so the grid is not so blocky
-            vmin_hd = min(x for x in iml_grid_hd.flatten() if not isnan(x))
-            vmax_hd = max(x for x in iml_grid_hd.flatten() if not isnan(x))
-            trim_thresh = vmin
-            iml_grid_hd[iml_grid_hd < trim_thresh] = np.nan
+        for j in range(0, len(products)):
+            vmin = np.nanmin(iml_grid_prep[j])
+            vmax = np.nanmax(iml_grid_prep[j])
+
+            #iml_grid[j] = np.reshape(iml_grid_prep[j], (nx, ny)).astype(dtype=np.float64)  # this reduces values to 8 decimal places
+            #iml_grid_tmp = np.nan_to_num(iml_grid[j])  # change nans to zeroes
+
+            # upscale the grid, trim, and interpolate if there are at least 10 grid values with range greater than 0.1
+            #if np.count_nonzero(iml_grid_tmp) >= 10 and vmax-vmin > 0.1:
+            #    up_factor = 1
+            #    iml_grid_hd = resize(iml_grid_tmp, (up_factor * len(iml_grid_tmp), up_factor * len(iml_grid_tmp)),
+            #                        mode='reflect', anti_aliasing=False)
+            #    trim_thresh = vmin
+            #    iml_grid_hd[iml_grid_hd < trim_thresh] = np.nan
+            #else:
+            #iml_grid_hd = iml_grid_tmp
+
+            #iml_grid_hd[iml_grid_hd == 0.0] = np.nan  # change zeroes back to nan
+
+            iml_grid_hd = iml_grid_prep[j]   
+            #vmin_hd = min(x for x in iml_grid_hd.flatten() if not isnan(x))
+            vmax_hd = np.nanmax(iml_grid_hd)
 
             # generate image overlay
-            north, south = lat.max(), lat.min()  # Latitude range
-            east, west = lon.max(), lon.min()  # Longitude range
+            north, south = grid_lat_max, grid_lat_min  # Latitude range
+            east, west = grid_lon_max, grid_lon_min  # Longitude range
             bounds = [[south, west], [north, east]]
-
+            
             map_center = [np.mean([north, south]), np.mean([east, west])]
 
             # Create an image from the grid
             cmap_name = 'YlOrRd'
             cmap = plt.get_cmap(cmap_name)
             fig, ax = plt.subplots(figsize=(6, 6))
-            ax.imshow(iml_grid_hd, origin='lower', cmap=cmap, vmin=vmin, vmax=vmax)
+            fig.add_axes([0, 0, 1, 1])
+            ax.imshow(iml_grid_hd, origin='lower', cmap=cmap, vmin=vmin, vmax=vmax, interpolation='bilinear', aspect='auto')
             ax.axis('off')
 
             # Save the figure
             fig.canvas.draw()
             overlay_filename = f"overlay_{j}.svg"
-            plt.savefig(overlay_filename, bbox_inches="tight", pad_inches=0, transparent=True)
+            plt.savefig(overlay_filename, pad_inches=0, transparent=True)
             plt.close(fig)
 
             # Embed geographic bounding box into the SVG

src/shf_wrapper.py

@@ -31,6 +31,9 @@ def main(argv):
             return False
         else:
             raise argparse.ArgumentTypeError("Boolean value expected.")
+        
+    def float_or_none(v):
+        return None if v.lower() == "none" else float(v)
 
     parser = argparse.ArgumentParser()
 
@@ -55,6 +58,7 @@ def main(argv):
     parser.add_argument("--time_unit", type=str)
     parser.add_argument("--model", type=str)
     parser.add_argument("--products_string", type=str)
+    parser.add_argument("--AOI_extent", nargs=4, type=float_or_none, default=[None] * 4, required=False)
     parser.add_argument("--verbose", type=str2bool)
 
     args = parser.parse_args()