diff --git a/docs/source/reference/analysis.md b/docs/source/reference/analysis.md index 32a18bb..70eebe8 100644 --- a/docs/source/reference/analysis.md +++ b/docs/source/reference/analysis.md @@ -14,6 +14,7 @@ shi_bolt_confidence analysis.estimate_b_tinti analysis.estimate_b_positive + analysis.estimate_b_more_positive analysis.estimate_b_utsu analysis.estimate_b_laplace analysis.estimate_b_weichert @@ -27,4 +28,14 @@ analysis.mc_ks analysis.mc_max_curvature +``` + +## Other +```{eval-rst} +.. autosummary:: + :toctree: api/ + + analysis.make_more_incomplete + analysis.beta_to_b_value + analysis.b_value_to_beta ``` \ No newline at end of file diff --git a/docs/source/reference/catalog.md b/docs/source/reference/catalog.md index bb44810..1780939 100644 --- a/docs/source/reference/catalog.md +++ b/docs/source/reference/catalog.md @@ -13,11 +13,35 @@ Catalog ``` -## Transformation +## Modify Catalog ```{eval-rst} .. autosummary:: :toctree: api/ Catalog.bin_magnitudes + Catalog.strip + Catalog.drop_ids + Catalog.drop_uncertainties +``` + +## Estimate from Catalog + +```{eval-rst} +.. autosummary:: + :toctree: api/ + + Catalog.estimate_b + Catalog.estimate_mc +``` + +## Transform from or to other format + +```{eval-rst} +.. autosummary:: + :toctree: api/ + + Catalog.to_quakeml + Catalog.from_quakeml + Catalog.from_dict ``` \ No newline at end of file diff --git a/docs/source/reference/utils.md b/docs/source/reference/utils.md index 430d5ec..2a1b675 100644 --- a/docs/source/reference/utils.md +++ b/docs/source/reference/utils.md @@ -23,7 +23,7 @@ :toctree: api/ utils.simulate_magnitudes - utils.simulated_magnitudes_binned + utils.simulate_magnitudes_binned ``` diff --git a/docs/source/user/getting_started.md b/docs/source/user/getting_started.md index ec9769e..608edc3 100644 --- a/docs/source/user/getting_started.md +++ b/docs/source/user/getting_started.md @@ -8,11 +8,11 @@ We didn't reinvent the wheel and rely on existing libraries and packages to perf ### GEOS The plotting of the seismicity requires [GEOS](https://libgeos.org/), a C/C++ library for computational geometry. If `GEOS` is not installed on your machine, you will need to get it, for example on a linux machine with -```terminal +```console sudo apt-get libgeos-dev ``` or on a mac with -```terminal +```console brew install geos ``` @@ -20,16 +20,16 @@ brew install geos ### Install from source This way of installing `SeismoStats` in another environement allows you to use the static version. -```terminal +```console pip install git+https://github.com/swiss-seismological-service/SeismoStats.git ``` If you want to install a specific branch: -```terminal +```console pip install git+https://github.com/swiss-seismological-service/SeismoStats.git@feature/branch ``` To update your environment to the latest version of `SeismoStats`: -```terminal +```console pip install --force-reinstall git+https://github.com/swiss-seismological-service/SeismoStats.git ``` \ No newline at end of file diff --git a/seismostats/analysis/__init__.py b/seismostats/analysis/__init__.py index c8d6e2e..553b6a1 100644 --- a/seismostats/analysis/__init__.py +++ b/seismostats/analysis/__init__.py @@ -5,5 +5,11 @@ estimate_b_positive, estimate_b_tinti, estimate_b_utsu, - estimate_b_weichert) + estimate_b_weichert, + estimate_b_kijko_smit, + estimate_b_more_positive, + make_more_incomplete, + beta_to_b_value, + b_value_to_beta + ) from seismostats.analysis.estimate_mc import mc_ks, mc_max_curvature diff --git a/seismostats/analysis/estimate_beta.py b/seismostats/analysis/estimate_beta.py index 86ca52d..61583cc 100644 --- a/seismostats/analysis/estimate_beta.py +++ b/seismostats/analysis/estimate_beta.py @@ -253,6 +253,7 @@ def differences(magnitudes: np.ndarray) -> np.ndarray: def estimate_b_positive( magnitudes: np.ndarray, delta_m: float = 0, + dmc: float | None = None, b_parameter: str = "b_value", return_std: bool = False, return_n: bool = False, @@ -270,6 +271,8 @@ def estimate_b_positive( To achieve the effect of reduced STAI, the magnitudes must be ordered in time. delta_m: discretization of magnitudes. default is no discretization. + dmc: cutoff value for the differences (diffferences below this + value are not considered). If None, the cutoff is set to delta_m b_parameter:either 'b-value', then the corresponding value of the Gutenberg-Richter law is returned, otherwise 'beta' from the exponential distribution [p(M) = exp(-beta*(M-mc))]. @@ -286,13 +289,95 @@ def estimate_b_positive( n: number of events used for the estimation """ + if dmc is None: + dmc = delta_m + elif dmc < 0: + raise ValueError("dmc must be larger or equal to 0") + elif dmc < delta_m: + warnings.warn("dmc is smaller than delta_m, not recommended") + mag_diffs = np.diff(magnitudes) + # only take the values where the next earthquake is d_mc larger than the + # previous one. delta_m is added to avoid numerical errors + mag_diffs = abs(mag_diffs[mag_diffs > dmc - delta_m / 2]) + + out = estimate_b_tinti( + mag_diffs, + mc=dmc, + delta_m=delta_m, + b_parameter=b_parameter, + return_std=return_std, + ) + + if return_n: + if type(out) is tuple: + return out + tuple([len(mag_diffs)]) + else: + return out, len(mag_diffs) + else: + return out + + +def estimate_b_more_positive( + magnitudes: np.ndarray, + delta_m: float = 0, + dmc: float | None = None, + b_parameter: str = "b_value", + return_std: bool = False, + return_n: bool = False, +) -> float | tuple[float, float] | tuple[float, float, float]: + """Return the b-value estimate calculated using the + next positive differences (this means that almost every magnitude has a + difference, as opposed to the b-positive method which results in half the + data). + + Source: + E. Lippiello and G. Petrillo. Journal of Geophysical Research: Solid + Earth, 129(2):e2023JB027849, 2024. + + Args: + magnitudes: array of magnitudes, sorted in time (first + entry is the earliest earthquake). + To achieve the effect + of reduced STAI, the magnitudes must be ordered in time. + delta_m: discretization of magnitudes. default is no discretization. + b_parameter:either 'b-value', then the corresponding value of the + Gutenberg-Richter law is returned, otherwise 'beta' from the + exponential distribution [p(M) = exp(-beta*(M-mc))]. + return_std: if True the standard deviation of beta/b-value (see above) + is returned. + return_n: if True the number of events used for the estimation is + returned. + + Returns: + b: maximum likelihood beta or b-value, depending on value of + input variable 'b_parameter'. Note that the difference is just a + factor [b_value = beta * log10(e)] + std: Shi and Bolt estimate of the beta/b-value estimate + n: number of events used for the estimation + """ + + if dmc is None: + dmc = delta_m + elif dmc < 0: + raise ValueError("dmc must be larger or equal to 0") + elif dmc < delta_m: + warnings.warn("dmc is smaller than delta_m, not recommended") + + mag_diffs = - np.ones(len(magnitudes) - 1) * delta_m + for ii in range(len(magnitudes) - 1): + for jj in range(ii + 1, len(magnitudes)): + mag_diff_loop = magnitudes[jj] - magnitudes[ii] + if mag_diff_loop > dmc - delta_m / 2: + mag_diffs[ii] = mag_diff_loop + break + # only take the values where the next earthquake is larger - mag_diffs = abs(mag_diffs[mag_diffs > 0]) + mag_diffs = abs(mag_diffs[mag_diffs > - delta_m / 2]) out = estimate_b_tinti( mag_diffs, - mc=delta_m, + mc=dmc, delta_m=delta_m, b_parameter=b_parameter, return_std=return_std, @@ -307,6 +392,61 @@ def estimate_b_positive( return out +def make_more_incomplete( + magnitudes: np.ndarray, + times: np.array, + delta_t: np.timedelta64 = np.timedelta64(60, "s"), + return_idx: bool = False, +) -> tuple[np.ndarray, np.ndarray]: + """Return filtered magnitudes and times. Filter the magnitudes and times in + the following way: If an earthquake is smaller than the previous one and + less than ``delta_t`` away, the earthquake is removed. + + Source: + E. Lippiello and G. Petrillo. Journal of Geophysical Research: Solid + Earth, 129(2):e2023JB027849, 2024. + + Args: + magnitudes: array of magnitudes, sorted in time (first + entry is the earliest earthquake). + times: array of datetime objects of occurrence of each earthquake + delta_t: time window in seconds to filter out events. default is 60 + seconds. + return_idx: if True the indices of the events that were kept are + returned + + Returns: + magnitudes: filtered array of magnitudes + times: filtered array of datetime objects + idx: indices of the events that were kept + """ + + # sort magnitudes in time + idx_sort = np.argsort(times) + magnitudes = magnitudes[idx_sort] + times = times[idx_sort] + + idx = np.full(len(magnitudes), True) + for ii in range(1, len(magnitudes)): + # get all times that are closer than delta_t + idx_close = np.where(times[ii] - times[:ii] < delta_t)[0] + + # check if these events are larger than the current event + idx_loop = magnitudes[idx_close] > magnitudes[ii] + + # if there are any, remove the current event + if sum(idx_loop) > 0: + idx[ii] = False + + magnitudes = magnitudes[idx] + times = times[idx] + + if return_idx is True: + return magnitudes, times, idx + + return magnitudes, times + + def estimate_b_laplace( magnitudes: np.ndarray, delta_m: float = 0, @@ -344,7 +484,7 @@ def estimate_b_laplace( mag_diffs = abs(mag_diffs) mag_diffs = mag_diffs[mag_diffs > 0] - return estimate_b_tinti( + out = estimate_b_tinti( mag_diffs, mc=delta_m, delta_m=delta_m, @@ -352,6 +492,14 @@ def estimate_b_laplace( return_std=return_std, ) + if return_n: + if type(out) is tuple: + return out + tuple([len(mag_diffs)]) + else: + return out, len(mag_diffs) + else: + return out + def estimate_b_weichert( magnitudes: np.ndarray, diff --git a/seismostats/analysis/estimate_mc.py b/seismostats/analysis/estimate_mc.py index 4687a44..0a6a2ee 100644 --- a/seismostats/analysis/estimate_mc.py +++ b/seismostats/analysis/estimate_mc.py @@ -303,12 +303,13 @@ def mc_max_curvature( catalogues: Estimating the magnitude of completeness and its uncertainty. Bulletin of the Seismological Society of America, 95(2), pp.684-698. + Args: - sample: Magnitudes to test - delta_m: Magnitude bins (sample has to be rounded to bins - beforehand) - correction_factor: Correction factor for the maximum curvature method - (default 0.2 after Woessner & Wiemer 2005) + sample: Magnitudes to test + delta_m: Magnitude bins (sample has to be rounded to bins beforehand) + correction_factor: Correction factor for the maximum curvature + method (default 0.2 after Woessner & Wiemer 2005) + Returns: mc: estimated completeness magnitude """ diff --git a/seismostats/analysis/tests/test_estimate_beta.py b/seismostats/analysis/tests/test_estimate_beta.py index 2d9c3c7..e35af24 100644 --- a/seismostats/analysis/tests/test_estimate_beta.py +++ b/seismostats/analysis/tests/test_estimate_beta.py @@ -3,6 +3,7 @@ import numpy as np import pytest +import datetime as dt # import functions to be tested from seismostats.analysis.estimate_beta import ( @@ -14,6 +15,8 @@ estimate_b_utsu, estimate_b_weichert, shi_bolt_confidence, + estimate_b_more_positive, + make_more_incomplete, ) from seismostats.utils.simulate_distributions import simulate_magnitudes_binned @@ -147,10 +150,10 @@ def test_differences(magnitudes: np.ndarray, mag_diffs: np.ndarray): @pytest.mark.parametrize( - "n,b,mc,delta_m,b_parameter,precision", + "n,b,mc,delta_m,b_parameter,dmc,precision", [ - (1000000, 1.2 * np.log(10), 3, 0, "beta", 0.005), - (1000000, np.log(10), 3, 0.1, "beta", 0.01), + (1000000, 1.2, 3, 0, "b_value", None, 0.005), + (1000000, np.log(10), 3, 0.1, "beta", 1, 0.01), ], ) def test_estimate_b_positive( @@ -159,13 +162,14 @@ def test_estimate_b_positive( mc: float, delta_m: float, b_parameter: str, + dmc: float, precision: float, ): mags = simulate_magnitudes_binned( n, b, mc, delta_m, b_parameter=b_parameter ) b_estimate = estimate_b_positive( - mags, delta_m=delta_m, b_parameter=b_parameter + mags, delta_m=delta_m, dmc=dmc, b_parameter=b_parameter ) assert abs(b - b_estimate) / b <= precision @@ -194,6 +198,65 @@ def test_estimate_b_laplace( assert abs(b - b_estimate) / b <= precision +@pytest.mark.parametrize( + "n,b,mc,delta_m,b_parameter,dmc,precision", + [ + (100000, 1.2 * np.log(10), 3, 0, "beta", None, 0.01), + (100000, 1, 3, 0.1, "b_value", 1, 0.04), + ], +) +def test_estimate_b_more_positive( + n: int, + b: float, + mc: float, + delta_m: float, + b_parameter: str, + dmc: float, + precision: float, +): + mags = simulate_magnitudes_binned( + n, b, mc, delta_m, b_parameter=b_parameter + ) + b_estimate = estimate_b_more_positive( + mags, delta_m=delta_m, dmc=dmc, b_parameter=b_parameter + ) + assert abs(b - b_estimate) / b <= precision + + +def test_make_more_incomplete(): + magnitudes = np.array([1, 2, 20, 3, 4, 9, 3]) + times = np.array([ + dt.datetime(2020, 1, 1), + dt.datetime(2020, 1, 2), + dt.datetime(2020, 1, 3), + dt.datetime(2020, 1, 4), + dt.datetime(2020, 1, 5), + dt.datetime(2020, 1, 6), + dt.datetime(2020, 1, 7), + ]) + + mags_inc, times_inc = make_more_incomplete( + magnitudes, times, delta_t=np.timedelta64(49, "h") + ) + + assert (mags_inc == [1, 2, 20, 9]).all() + assert ( + times_inc + == [ + dt.datetime(2020, 1, 1), + dt.datetime(2020, 1, 2), + dt.datetime(2020, 1, 3), + dt.datetime(2020, 1, 6), + ] + ).all() + + mags_inc, times_inc, idx = make_more_incomplete( + magnitudes, times, delta_t=np.timedelta64(49, "h"), return_idx=True + ) + + assert (mags_inc == magnitudes[idx]).all() + + def _create_test_catalog_poisson(a_val_true: float, b_val_true: float): """ creates a synthetic catalog with magnitudes diff --git a/seismostats/analysis/tests/test_estimate_mc.py b/seismostats/analysis/tests/test_estimate_mc.py index b6626e8..db3bc61 100644 --- a/seismostats/analysis/tests/test_estimate_mc.py +++ b/seismostats/analysis/tests/test_estimate_mc.py @@ -1,16 +1,14 @@ -from numpy.testing import assert_allclose, assert_equal -import numpy as np import pickle +import numpy as np import pytest -from seismostats.io.client import FDSNWSEventClient import pandas as pd +from numpy.testing import assert_allclose, assert_almost_equal, assert_equal -from seismostats.analysis.estimate_mc import ( - empirical_cdf, - mc_max_curvature, - mc_ks, - mc_by_bvalue_stability -) +from seismostats.io.client import FDSNWSEventClient +from seismostats.analysis.estimate_mc import (empirical_cdf, + mc_ks, + mc_max_curvature, + mc_by_bvalue_stability) @pytest.fixture @@ -177,7 +175,7 @@ def test_estimate_mc_ks(mags, mcs): mags, delta_m=0.1, mcs_test=mcs, p_pass=0.1 ) assert_equal(1.1, best_mc) - assert_equal(2.242124985031149, best_beta) + assert_almost_equal(2.242124985031149, best_beta) assert_equal([0.8, 0.9, 1.0, 1.1], mcs_tested) assert_allclose( [ @@ -219,7 +217,7 @@ def test_estimate_mc_ks(mags, mcs): mags, delta_m=0.1, p_pass=0.1 ) assert_equal(1.1, best_mc) - assert_equal(2.242124985031149, best_beta) + assert_almost_equal(2.242124985031149, best_beta) assert_equal([1.0, 1.1], mcs_tested) # test when b-positive is used @@ -227,7 +225,7 @@ def test_estimate_mc_ks(mags, mcs): mags, delta_m=0.1, b_method="positive" ) assert_equal(1.5, best_mc) - assert_equal(3.2542240043462796, best_beta) + assert_almost_equal(3.2542240043462796, best_beta) assert_equal(len(mcs_tested), 6) assert_equal(len(betas), 6) assert_equal(len(ks_ds), 6)