SeismoStats: How To#
In this notebook we will show how to:#
Make a catalog object
by downloading it
by converting a csv file
<li>Plot the seismicity</li> <li>Analyze the FMD</li> <ol> <li>Plot the FMD</li> <li>Estimate the magnitude of completeness</li> <li>Estimate b-values</li> </ol> <li>Generate synthetic magnitudes</li> <li>Bin magnitudes</li>
0. Import general packages#
[1]:
#%matplotlib widget
import datetime as dt
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
1. Make a catalog object#
The catalog object is a dataframe, with some additional methods and attributes. The columns are:
column header |
type |
importance |
|---|---|---|
latitude |
float |
required |
longitude |
float |
required |
depth |
float |
required |
mag_type |
string |
optional |
magnitude |
float |
required |
time |
pandas timestamp |
required |
event_type |
string |
optional |
1.1 Download catalog#
[2]:
from seismostats import FDSNWSEventClient
Start date and end date have to be defined as a datetime. In case something does not work out, the link to retrieve it manually is given back.
[3]:
start_time = pd.to_datetime('2020/01/01')
end_time = pd.to_datetime('2022/01/01')
min_longitude = 5
max_longitude = 11
min_latitude = 45
max_latitude = 48
min_magnitude = 0.5
url = 'http://eida.ethz.ch/fdsnws/event/1/query'
client = FDSNWSEventClient(url)
cat = client.get_events(
start_time=start_time,
end_time=end_time,
min_magnitude=min_magnitude,
min_longitude=min_longitude,
max_longitude=max_longitude,
min_latitude=min_latitude,
max_latitude=max_latitude)
The output is a catalog object
[4]:
cat.tail()
[4]:
| event_type | time | latitude | longitude | depth | evaluationmode | magnitude | magnitude_type | magnitude_MLhc | magnitude_MLh | |
|---|---|---|---|---|---|---|---|---|---|---|
| 2787 | quarry blast | 2020-01-03 14:43:49.025320 | 47.187560 | 7.185673 | -683.593750 | manual | 1.213024 | MLh | NaN | 1.213023609 |
| 2788 | earthquake | 2020-01-03 14:28:09.701876 | 46.444436 | 9.104820 | 2033.203125 | manual | 0.596428 | MLh | NaN | 0.5964283368 |
| 2789 | quarry blast | 2020-01-02 08:47:20.725352 | 47.673434 | 7.585486 | 12113.281250 | manual | 0.657227 | MLh | NaN | 0.6572274734 |
| 2790 | earthquake | 2020-01-01 17:42:48.508164 | 46.031975 | 6.892110 | 5295.898438 | manual | 0.826313 | MLh | NaN | 0.8263128507 |
| 2791 | earthquake | 2020-01-01 13:43:47.626410 | 45.704174 | 7.068708 | 3302.734375 | manual | 0.824352 | MLh | NaN | 0.8243519995 |
2. Seismicity Plots#
We have three basic plots of the seismicity
Seismicity in space
Cumulative count
Magnitudes in time
2.1. Plot in space#
It is possible to choose the resolutions 10, 50 and 110. Optionally you can choose the color scheme and the country of which you want to see the borders. For large areas/countries, the plotting can take some time.
[5]:
fig = plt.figure(figsize=(10, 10), linewidth=1)
ax = cat.plot_in_space(
resolution='10m',
include_map=True,
country='Switzerland',
colors='Greys_r',
dot_labels=[1,2,3])
You can also choose the interpolation power and the size of the smallest and largest dot.
[6]:
fig = plt.figure(figsize=(10, 10), linewidth=1)
ax = cat.plot_in_space(
resolution='10m',
include_map=True,
dot_smallest=1,
dot_largest=500,
dot_interpolation_power=3,
dot_labels=None)
The default is that the map is not included.
[7]:
fig = plt.figure(figsize=(10, 10), linewidth=1)
ax = cat.plot_in_space(
dot_smallest=0.1,
dot_largest=10,
dot_interpolation_power=0,
dot_labels=False)
2.2. Plot in time#
[8]:
ax = cat.plot_cum_count(mcs=np.arange(0.5, 4.0, 1), delta_m=0.1)
[9]:
ax = cat.plot_mags_in_time()
3. analysze the FMD#
3.1 Plot magnitude distributions#
[14]:
# bin the magnitudes (if the binning of the catalog is known, one can just set it instead with e.g. cat.delta_m = 0.01)
cat.bin_magnitudes(delta_m=0.01, inplace=True)
[15]:
ax = plt.subplots(figsize=(8, 6))[1]
cat.plot_cum_fmd(ax=ax, fmd_bin=0.01)
cat.plot_fmd(ax=ax, grid=True, fmd_bin=0.01)
# if you want to use a different binning for the plot_fmd, you can set it here
cat.plot_fmd(fmd_bin= 0.1, ax=ax, grid=True, legend='binned with delta_m=0.1')
plt.show()
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
Cell In[15], line 3
1 ax = plt.subplots(figsize=(8, 6))[1]
----> 3 cat.plot_cum_fmd(ax=ax, fmd_bin=0.01)
4 cat.plot_fmd(ax=ax, grid=True, fmd_bin=0.01)
6 # if you want to use a different binning for the plot_fmd, you can set it here
File ~/Documents/Codes_and_tools/SeismoStats/seismostats/utils/__init__.py:63, in require_cols.<locals>.decorator_require.<locals>.wrapper_require(self, *args, **kwargs)
57 if not _check_required_cols(self, require):
58 raise AttributeError(
59 "Catalog is missing the following columns "
60 f'for execution of the method "{func.__name__}": '
61 f"{set(require).difference(set(self.columns))}."
62 )
---> 63 value = func(self, *args, **kwargs)
64 return value
File ~/Documents/Codes_and_tools/SeismoStats/seismostats/catalogs/catalog.py:1265, in Catalog.plot_cum_fmd(self, mc, fmd_bin, b_value, ax, color, size, grid, bin_position, legend)
1263 if b_value is None:
1264 b_value = self.b_value
-> 1265 ax = plot_cum_fmd(self.magnitude,
1266 b_value=b_value,
1267 mc=mc,
1268 fmd_bin=fmd_bin,
1269 ax=ax,
1270 color=color,
1271 size=size,
1272 grid=grid,
1273 bin_position=bin_position,
1274 legend=legend)
1275 return ax
TypeError: plot_cum_fmd() got an unexpected keyword argument 'fmd_bin'
3.2 Estimate completeness magnitude#
[16]:
# estimate b-value by stability
mc_stab, _ = cat.estimate_mc_b_stability(stop_when_passed=True)
# the b-value is saved within the catalog
cat.mc
return_vals: {'best_b_value': np.float64(1.0296700460964268), 'mcs_tested': [0.45, 0.46, 0.47, 0.48, 0.49, 0.5, 0.51, 0.52, 0.53, 0.54, 0.55, 0.56, 0.57, 0.58, 0.59, 0.6, 0.61, 0.62, 0.63, 0.64, 0.65, 0.66, 0.67, 0.68, 0.69, 0.7, 0.71, 0.72, 0.73, 0.74, 0.75, 0.76, 0.77, 0.78, 0.79, 0.8, 0.81, 0.82, 0.83, 0.84, 0.85, 0.86, 0.87, 0.88, 0.89, 0.9, 0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99, 1.0, 1.01, 1.02, 1.03, 1.04, 1.05, 1.06, 1.07, 1.08, 1.09, 1.1, 1.11, 1.12, 1.13, 1.14, 1.15, 1.16, 1.17, 1.18, 1.19, 1.2, 1.21, 1.22, 1.23, 1.24, 1.25, 1.26, 1.27, 1.28, 1.29, 1.3, 1.31, 1.32, 1.33, 1.34, 1.35, 1.36, 1.37, 1.38, 1.39, 1.4, 1.41, 1.42, 1.43, 1.44, 1.45, 1.46], 'b_values_tested': [np.float64(0.5713514154969684), np.float64(0.5783420877137975), np.float64(0.5829741129130803), np.float64(0.5881182062220058), np.float64(0.5931481105494977), np.float64(0.5978346077286778), np.float64(0.6008150762022235), np.float64(0.6062973304255794), np.float64(0.6102844862308586), np.float64(0.6147864679936459), np.float64(0.6184120585365727), np.float64(0.622067779140773), np.float64(0.6262403145768893), np.float64(0.6307111585394682), np.float64(0.6372496096276045), np.float64(0.6421619327329338), np.float64(0.6456048469581808), np.float64(0.6522029070395186), np.float64(0.6562944898563844), np.float64(0.6593422197651272), np.float64(0.6651187821462178), np.float64(0.6707287970295714), np.float64(0.6724809969303653), np.float64(0.6770243433235135), np.float64(0.6813249424536868), np.float64(0.6871435360865622), np.float64(0.6924662950925543), np.float64(0.6966449929981451), np.float64(0.7027127678077805), np.float64(0.7082613565882054), np.float64(0.7113233439015858), np.float64(0.7146707531957751), np.float64(0.7163279988232324), np.float64(0.7202339727458241), np.float64(0.724136101085941), np.float64(0.7273334093371504), np.float64(0.730486254826891), np.float64(0.737564143305669), np.float64(0.7418536717941107), np.float64(0.7439035069023177), np.float64(0.7484998117521112), np.float64(0.750789906953446), np.float64(0.7525778685772294), np.float64(0.7558318748987053), np.float64(0.7614716672451004), np.float64(0.764685719785733), np.float64(0.7678306731999055), np.float64(0.7752001895096526), np.float64(0.7792190107246149), np.float64(0.7823108823543774), np.float64(0.7875830240482408), np.float64(0.791034772534308), np.float64(0.7962937755520441), np.float64(0.8025340967403194), np.float64(0.8059169395690938), np.float64(0.8077148639931954), np.float64(0.8138926634671096), np.float64(0.8196169085506544), np.float64(0.8222140703002485), np.float64(0.8273437527269127), np.float64(0.8291830867951706), np.float64(0.8308197822004573), np.float64(0.8339458062019462), np.float64(0.8346223434972214), np.float64(0.8409303573043546), np.float64(0.8436759302676176), np.float64(0.8456440899944452), np.float64(0.8511593150188479), np.float64(0.8534770325063848), np.float64(0.8542927277303405), np.float64(0.8627931839076718), np.float64(0.8653616861757002), np.float64(0.8670426176734574), np.float64(0.8748231496266673), np.float64(0.8784007647190608), np.float64(0.8833017972175069), np.float64(0.8873842030603308), np.float64(0.887513916075401), np.float64(0.8911773691904326), np.float64(0.896282138733421), np.float64(0.9053924150075393), np.float64(0.9080288949591073), np.float64(0.9155238631466034), np.float64(0.920484604707641), np.float64(0.9209117563614001), np.float64(0.9354286349562014), np.float64(0.9421352098699323), np.float64(0.9526189538049191), np.float64(0.9623747284735956), np.float64(0.9673667050965653), np.float64(0.9701843873508833), np.float64(0.9737310222958364), np.float64(0.9759602203021255), np.float64(0.9821631921952334), np.float64(0.9871749047361138), np.float64(0.9874825630502299), np.float64(0.9953936305262784), np.float64(1.0033265259565667), np.float64(1.0088469520138679), np.float64(1.0179290874668885), np.float64(1.0207332930916213), np.float64(1.0296700460964268)], 'diff_bs': [np.float64(14.517336942524617), np.float64(13.826726809225462), np.float64(13.55871840666206), np.float64(13.211827787501608), np.float64(12.895483604867094), np.float64(12.63836332694301), np.float64(12.650938523137713), np.float64(12.276967735547869), np.float64(12.14210514614908), np.float64(11.928242598666085), np.float64(11.849675684658797), np.float64(11.763099740110212), np.float64(11.598067939375975), np.float64(13.179485287381732), np.float64(10.894777104167098), np.float64(12.373503121115075), np.float64(10.565187918064002), np.float64(11.78557433765248), np.float64(9.935681788483993), np.float64(11.578723879169383), np.float64(9.553073140942661), np.float64(10.851831261538019), np.float64(9.356805765298965), np.float64(10.762961621553682), np.float64(8.994532368480371), np.float64(8.663467851934932), np.float64(8.394195826832961), np.float64(8.252024266942193), np.float64(7.906715069956867), np.float64(7.620631816162461), np.float64(7.5910818796128785), np.float64(7.542488434257369), np.float64(7.660074184192062), np.float64(7.563339471687974), np.float64(7.469748078457214), np.float64(7.438980146922742), np.float64(7.4293452174833545), np.float64(7.06109820802947), np.float64(8.37969764513728), np.float64(7.0680221691396685), np.float64(8.348375516558047), np.float64(7.035958304804902), np.float64(8.542098580816836), np.float64(7.161653550768327), np.float64(8.314158094300536), np.float64(6.964524025893549), np.float64(8.30185057499497), np.float64(6.6367834026815045), np.float64(7.882646816629248), np.float64(6.596133256293603), np.float64(6.450917146148853), np.float64(6.443066903536373), np.float64(6.305881064583034), np.float64(6.097886451549821), np.float64(6.100716058563534), np.float64(6.223123916330054), np.float64(6.032058564856592), np.float64(5.8743157723628405), np.float64(5.928426358551749), np.float64(5.821997788862371), np.float64(5.935065737646834), np.float64(6.058245048328768), np.float64(6.082118518844141), np.float64(6.263780252499514), np.float64(6.070044331227414), np.float64(6.092747067468511), np.float64(6.166715874975575), np.float64(6.016702225837636), np.float64(6.048527594146056), np.float64(6.162783476553231), np.float64(5.821034290131426), np.float64(5.825995182942661), np.float64(5.8760280317910425), np.float64(5.59049621069365), np.float64(5.546021503052255), np.float64(5.41897360785256), np.float64(5.3371370765254955), np.float64(5.46851904970325), np.float64(5.404109002276291), np.float64(5.265143206734002), np.float64(4.920434665934402), np.float64(4.882678960540737), np.float64(4.608066537116519), np.float64(4.448697659716064), np.float64(4.4994369876503475), np.float64(3.9205966979152214), np.float64(3.6922077534241944), np.float64(3.3195323910242807), np.float64(2.9887842953730965), np.float64(2.844003818164201), np.float64(2.7733040525390007), np.float64(2.679164290076714), np.float64(2.6209172478616165), np.float64(2.414373770122177), np.float64(2.2384357123092657), np.float64(2.2188303222114234), np.float64(1.9344586813250648), np.float64(1.6516379397097698), np.float64(1.4520476145123495), np.float64(1.1536097447268316), np.float64(1.0395166371473186), np.float64(0.7569908559232158)]}
[16]:
np.float64(1.46)
[17]:
# You have to chose a delta_m here. This is basicallhy the binning that is used for the maxc method,
# and might be different from the binning of the catalog. Originally, 0.1 was used, but the optimal
# value depends on the data. The best way is to plot_fmd and use the delta_m that makes the most sense for the user.
mc_maxc, _ = cat.estimate_mc_maxc(fmd_bin=0.1)
# now, the mc is set to the one that was estimated by the maxc method.
cat.mc
[17]:
np.float64(1.1)
This method takes longer, especially when the magnitude sample is large.
If Mc is known to be larger than a certain value, giving the Mc values that should be tested as an input can make the Mc estimation faster.
[18]:
mc_kstest, mc_info = cat.estimate_mc_ks(
mcs_test=np.arange(1.0, 3.0, 0.1),
p_value_pass=0.1,
)
print("Tested Mc values:", mc_info['mcs_tested'])
print("First Mc to pass the KS test:", mc_kstest)
print(f"Associated beta value: {mc_info['best_b_value']:.2f}")
/Users/vanilleritz/Documents/Codes_and_tools/SeismoStats/seismostats/analysis/estimate_mc.py:247: UserWarning: Mcs to test are not binned correctly. Test might fail because of this.
warnings.warn("Mcs to test are not binned correctly. "
Tested Mc values: [np.float64(1.0), np.float64(1.1), np.float64(1.2000000000000002), np.float64(1.3000000000000003), np.float64(1.4000000000000004), np.float64(1.5000000000000004)]
First Mc to pass the KS test: 1.5000000000000004
Associated beta value: 1.05
[20]:
ax = plt.subplots(figsize=(8, 6))[1]
cat.plot_cum_fmd(fmd_bin=0.1, ax=ax, color='blue')
cat.plot_fmd(fmd_bin=0.1, ax=ax, color='red')
plt.axvline(mc_maxc, color='black', linestyle='--', label='Maximum curvature $m_c$')
plt.axvline(mc_stab, color='grey', linestyle='--', label='Stability $m_c$')
plt.axvline(mc_kstest, color='lightgrey', linestyle='--', label='KS-test $m_c$')
plt.legend()
plt.show()
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
Cell In[20], line 3
1 ax = plt.subplots(figsize=(8, 6))[1]
----> 3 cat.plot_cum_fmd(fmd_bin=0.1, ax=ax, color='blue')
4 cat.plot_fmd(fmd_bin=0.1, ax=ax, color='red')
6 plt.axvline(mc_maxc, color='black', linestyle='--', label='Maximum curvature $m_c$')
File ~/Documents/Codes_and_tools/SeismoStats/seismostats/utils/__init__.py:63, in require_cols.<locals>.decorator_require.<locals>.wrapper_require(self, *args, **kwargs)
57 if not _check_required_cols(self, require):
58 raise AttributeError(
59 "Catalog is missing the following columns "
60 f'for execution of the method "{func.__name__}": '
61 f"{set(require).difference(set(self.columns))}."
62 )
---> 63 value = func(self, *args, **kwargs)
64 return value
File ~/Documents/Codes_and_tools/SeismoStats/seismostats/catalogs/catalog.py:1265, in Catalog.plot_cum_fmd(self, mc, fmd_bin, b_value, ax, color, size, grid, bin_position, legend)
1263 if b_value is None:
1264 b_value = self.b_value
-> 1265 ax = plot_cum_fmd(self.magnitude,
1266 b_value=b_value,
1267 mc=mc,
1268 fmd_bin=fmd_bin,
1269 ax=ax,
1270 color=color,
1271 size=size,
1272 grid=grid,
1273 bin_position=bin_position,
1274 legend=legend)
1275 return ax
TypeError: plot_cum_fmd() got an unexpected keyword argument 'fmd_bin'
3.3 Estimate the b-value#
[21]:
b_estimator = cat.estimate_b()
# take care:cat.delta_m and cat.mc are used here!
cat.b_value
[21]:
np.float64(1.0511521850648504)
[22]:
# the b-estimator that was given back above contains more information on the b-value estimation
print('b-value', b_estimator.b_value)
print('standard deviation', b_estimator.std)
print('number of events used', b_estimator.n)
b-value 1.0511521850648504
standard deviation 0.0403308551924757
number of events used 742
[24]:
cat.plot_cum_fmd(color=['blue', 'orange'], fmd_bin=0.01)
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
Cell In[24], line 1
----> 1 cat.plot_cum_fmd(color=['blue', 'orange'], fmd_bin=0.01)
File ~/Documents/Codes_and_tools/SeismoStats/seismostats/utils/__init__.py:63, in require_cols.<locals>.decorator_require.<locals>.wrapper_require(self, *args, **kwargs)
57 if not _check_required_cols(self, require):
58 raise AttributeError(
59 "Catalog is missing the following columns "
60 f'for execution of the method "{func.__name__}": '
61 f"{set(require).difference(set(self.columns))}."
62 )
---> 63 value = func(self, *args, **kwargs)
64 return value
File ~/Documents/Codes_and_tools/SeismoStats/seismostats/catalogs/catalog.py:1265, in Catalog.plot_cum_fmd(self, mc, fmd_bin, b_value, ax, color, size, grid, bin_position, legend)
1263 if b_value is None:
1264 b_value = self.b_value
-> 1265 ax = plot_cum_fmd(self.magnitude,
1266 b_value=b_value,
1267 mc=mc,
1268 fmd_bin=fmd_bin,
1269 ax=ax,
1270 color=color,
1271 size=size,
1272 grid=grid,
1273 bin_position=bin_position,
1274 legend=legend)
1275 return ax
TypeError: plot_cum_fmd() got an unexpected keyword argument 'fmd_bin'
We can also plot the b-value for different mc:
[25]:
from seismostats.analysis import ClassicBValueEstimator, BPositiveBValueEstimator
mcs = np.arange(0.5, 3, 0.1)
ax = cat.plot_mc_vs_b(mcs, b_method=BPositiveBValueEstimator, dmc=0.1, color='red', label='b positive, dmc=0.1')
cat.plot_mc_vs_b(mcs, b_method=BPositiveBValueEstimator, dmc=0.2, ax=ax, color='orange', label='b positive, dmc=0.2')
cat.plot_mc_vs_b(mcs, b_method=BPositiveBValueEstimator, dmc=0.4, ax=ax, color='yellow', label='b positive, dmc=0.4')
cat.plot_mc_vs_b(mcs, b_method=ClassicBValueEstimator, ax =ax, color='blue', label='classic $b$')
[25]:
<Axes: xlabel='Completeness magnitude $m_c$', ylabel='b-value'>
3.4 Check if the b-value changes significantly#
Once we picked a fitting completeness, we might also want to see how the b-value changes with time (or along any other dimension). Unfortunately, this one is not yet implemented for the catalog class, so we have to use the original function
[43]:
from seismostats.analysis import b_significant_1D
from seismostats.plots import plot_b_series_constant_nm, plot_b_significant_1D
from seismostats.analysis import BPositiveBValueEstimator
[44]:
mc = 1
delta_m = cat.delta_m
times = cat.time
mags = cat.magnitude
[45]:
n_m = 100 # number of magnitudes taken per estimate in the running window
idx = mags >= 0
ax = plot_b_series_constant_nm(mags, delta_m, mc, times, n_m=n_m,x_variable=times, color='#1f77b4', plot_technique='right', label='classical b-value')
ax = plot_b_series_constant_nm(mags[idx], delta_m, mc, times[idx], n_m=n_m,x_variable=times[idx], color='red', plot_technique='right', label='b-positive', ax=ax, b_method=BPositiveBValueEstimator)
_ = plt.xticks(rotation=45)
ax.set_xlabel('Time')
/Users/vanilleritz/Documents/Codes_and_tools/SeismoStats/seismostats/analysis/bvalue/base.py:110: UserWarning: No magnitudes in the lowest magnitude bin are present.Check if mc is chosen correctly.
warnings.warn(
[45]:
Text(0.5, 0, 'Time')
Looking at the time-series above, one could be interested if the variation of the b-value is larger than what one would expect just from random fluctuation of the estimate. In other words, we want to know if we can reject the null-hypothesis that the true b-value is constant. For this, we can apply the method of Mirwald et. al., 2024.
For this, we estimate the mean autocorrelation (MAC). The MAC can then be used to estimate a p-value. If the p-value is smaller than a threshold (which we have to choose, often 0.05 is used), then we can reject the null-hypothesis, and we are justified to believe that the b-vlue is in fact changing.
[46]:
n_m = 100
p, mac, mu_mac, std_mac = b_significant_1D(mags, mc, delta_m, times, n_m)
# b-positive
p, mac, mu_mac, std_mac = b_significant_1D(mags, mc, delta_m, times, n_m, method= BPositiveBValueEstimator)
/Users/vanilleritz/Documents/Codes_and_tools/SeismoStats/seismostats/analysis/b_significant.py:367: UserWarning: The number of subsamples is less than 25. The normality assumption of the autocorrelation might not be valid.
warnings.warn(
[47]:
p_threshold = 0.05
print('The p-value of a constant b-value hypothesis is {:.2f}'.format(p))
print('This is significantly larger than our threshold of {:.2f}. Therefore, we cannot reject the null-hypothesis'.format(p_threshold))
The p-value of a constant b-value hypothesis is 0.16
This is significantly larger than our threshold of 0.05. Therefore, we cannot reject the null-hypothesis
We found that the temporal variation was not significant, therefore further interpretation of how the b-value changes with time might not be reasonable to do. But this was specifically using a certain number of magnitudes per estimate. Maybe there is some other scale, where the b-value does change significantly?
We can test this easily by applying the same method with different n_m.
[48]:
ax = plot_b_significant_1D(
mags, times, mc, delta_m, x_variable = times, color = '#1f77b4', label='classical b-value')
plot_b_significant_1D(
mags, times, mc, delta_m, x_variable = times, color = 'red', b_method=BPositiveBValueEstimator, ax = ax, label='b-positive')
[48]:
<Axes: xlabel='$n_m$', ylabel='MAC'>
We found that in fact, the variation of the b-value is not significant within the scales that we can practically test.
4. Generate and bin synthetic earthquakes#
First we need to define the number of earthquakes, the b-value and the completeness magnitude. If binnning is applied, it is important to generate the magnitudes half a bin smaller than the smallest magnitude, otherwise the first bin will contain only half the events. For the b-value, note that beta is defined as the natural logarithm equivalent of the b-value.
[189]:
from seismostats import Catalog
from seismostats.analysis import estimate_b
from seismostats.utils import simulate_magnitudes_binned
from seismostats.plots import plot_cum_fmd
import matplotlib.pyplot as plt
[190]:
n = 200
b_value = 1.5
delta_m = 0.05
mc = 3
dmc = 0.3
Now we can generate a synthetic magnitude distribution:
[191]:
mags = simulate_magnitudes_binned(n, b_value, mc, delta_m)
cat = pd.DataFrame({'magnitude': mags})
cat = Catalog(cat)
cat.delta_m = delta_m
cat.mc = mc
Now we can estimate the b-value:
[193]:
# equivelently, you can use the Catalog method
cat.estimate_b()
print('b-value', cat.b_value)
b-value 1.57641970308759
We can plot the original and binned magnitudes and their respective b-value estimates now. Note that we choose the bin position to be left in order to align the binned and the original magnitudes.
[194]:
ax = plt.subplots(figsize=(8, 6))[1]
cat.plot_cum_fmd(ax=ax, b_value=b_value,
color=['blue', 'darkblue'], legend=['Magnitudes binned', 'True b-value'], size=3)
cat.plot_cum_fmd(ax=ax,
color=['blue', 'red'], legend=['_', 'estimated b-value'], size=3)
plt.show()
