From 3fd952501857932bbbd8f09fe2af28c2607e8a17 Mon Sep 17 00:00:00 2001 From: Kirsty Bayliss Date: Fri, 24 May 2024 13:09:53 +0200 Subject: [PATCH 01/23] add descriptions of different norm options --- openquake/cat/completeness/norms.py | 205 ++++++++++++++++++++++------ 1 file changed, 160 insertions(+), 45 deletions(-) diff --git a/openquake/cat/completeness/norms.py b/openquake/cat/completeness/norms.py index 9c1d18fe6..8f714313f 100644 --- a/openquake/cat/completeness/norms.py +++ b/openquake/cat/completeness/norms.py @@ -101,6 +101,10 @@ def get_completeness_matrix(tcat, ctab, mbinw, ybinw): def get_norm_optimize(tcat, aval, bval, ctab, cmag, n_obs, t_per, last_year, info=False): """ + This calculates the difference between the number of observed vs computed events in each completeness bin. + The differences are normalised by the number of completeness bins in order to minimise the effect of the number + of bins. + :param aval: GR a-value :param bval: @@ -109,13 +113,19 @@ def get_norm_optimize(tcat, aval, bval, ctab, cmag, n_obs, t_per, last_year, completeness table :param cmag: An array with the magnitude values at the center of each occurrence - bins + bins. Output from hmtk.seismicity.occurrence.utils.get_completeness_counts :param t_per: + array indicating total duration (in years) of completeness + Output from hmtk.seismicity.occurrence.utils.get_completeness_counts :param n_obs: - Number of observations + number of events in completeness period + Output from hmtk.seismicity.occurrence.utils.get_completeness_counts :param last_year: + last year to consider in analysis :param info: boolean controlling whether to print information as the function proceeds + :returns: + calculated norm for input completeness. Smaller norm is better """ occ = np.zeros((ctab.shape[0])) @@ -177,9 +187,11 @@ def get_norm_optimize(tcat, aval, bval, ctab, cmag, n_obs, t_per, last_year, def get_norm_optimize_a(aval, bval, ctab, binw, cmag, n_obs, t_per, info=False): """ - Computes a norm using a slightly different strategy than the one used in - `get_norm_optimize` - based on the probability of observing n events in each - magnitude bin relative to an exponential (GR) frequency-magnitude distribution + Computes a norm based on the probability of observing n events in each + magnitude bin relative to an exponential (GR) frequency-magnitude distribution. + The norm is the log-likelihood of observing the number of events in each magnitude + bin given the expected rate from the GR parameters calculated for these + completeness window. Larger is better. :param aval: GR a-value @@ -187,15 +199,22 @@ def get_norm_optimize_a(aval, bval, ctab, binw, cmag, n_obs, t_per, info=False) GR b-value :param ctab: Completeness table + :param binw: + binwidth for completeness analysis, specified in toml :param cmag: An array with the magnitude values at the center of each occurrence - bins - :param t_per: + bin :param n_obs: - :param last_year: - :param info: + number of events in completeness period + Output from hmtk.seismicity.occurrence.utils.get_completeness_counts + :param t_per: + array indicating total duration (in years) of completeness + Output from hmtk.seismicity.occurrence.utils.get_completeness_counts + :returns: + norm - log-likelihood of observing the number of events in + each magnitude bin given GR params from completeness bin. + Larger is better. """ - #breakpoint() # Rates of occurrence in each magnitude bin in the completeness interval rates = (10**(-bval * (cmag - binw/2) + aval) - 10**(-bval * (cmag + binw/2) + aval))*t_per @@ -215,21 +234,21 @@ def get_norm_optimize_a(aval, bval, ctab, binw, cmag, n_obs, t_per, info=False) for i, obs in enumerate(n_obs): log_prob[i] = (-rates[i]) + (n_obs[i]*np.math.log(rates[i])) - np.math.log(np.math.factorial(n_obs[i])) - #log_prob = (-rates) + (n_obs*np.math.log(rates)) - np.math.log(np.math.factorial(n_obs)) - - #norm = 1. - np.prod(occ_prob) + #print("from log prob: ", np.exp(log_prob)) #print("log likelihood = ", np.sum(log_prob)) - #norm = np.sum(log_prob) - norm = 1 - np.prod(np.exp(log_prob)) + norm = np.sum(log_prob) + #norm = 1 - np.prod(np.exp(log_prob)) return norm def get_norm_optimize_b(aval, bval, ctab, tcat, mbinw, ybinw, back=5, mmin=-1, mmax=10): """ - Computes a norm using a slightly different strategy than the one used in - `get_norm_optimize` + Computes a norm based on occurrences inside and outside of the completeness windows. + Difference in expected vs observed occurrences in/outside of completeness are + normalised by the rates from a GR distribution given the completeness. + The norm is calculated from the ratio of events in/out the window, and should be maximised. :param aval: GR a-value @@ -249,6 +268,9 @@ def get_norm_optimize_b(aval, bval, ctab, tcat, mbinw, ybinw, back=5, mmin=-1, Minimum magnitude :param mmax: Maximum magnitude + :returns: + norm - a ratio of the difference in event numbers within/outwith + completeness relative to expected GR. Larger is better. """ # oin and out have shape (num mags bins) x (num years bins) @@ -303,8 +325,17 @@ def get_norm_optimize_b(aval, bval, ctab, tcat, mbinw, ybinw, back=5, mmin=-1, return norm - def get_idx_compl(mag, compl): + """ + Find completeness windows for a specified magnitude + + :param mag: + magnitude + :param compl: + completeness table + :returns: + + """ if mag < compl[0, 1]: return None for i, com in enumerate(compl[:-1, :]): @@ -312,10 +343,23 @@ def get_idx_compl(mag, compl): return i return len(compl)-1 -def poiss_prob_int_time(rate, n, t, log_out = False): +def poiss_prob_int_time(rate, n, t, log_out = True): """ Calculate poisson probability of observing n events in some time step t given rate - Should this be a log? Probably. Yes. factorials and powers make this kinda pesky + In most cases, a log probability is preferred, due to issues with powers and factorials + when large numbers of events are involved. + + :param rate: + Poisson rate + :param n: + Number of observed events + :param t: + time step, multiplied with rate to get poisson expected number + :param log_out: + boolean indicating if log probabilities are preferred. If n is large, this should + be set to true to avoid instabilities. + :returns: + Poisson probability of observing n events in time t given poisson rate """ # Should use log probabilities so this doesn't break at large n log_prob = -(rate*t) + n*(np.log(rate) + np.log(t)) - np.math.log(np.math.factorial(n)) @@ -329,8 +373,32 @@ def poiss_prob_int_time(rate, n, t, log_out = False): def get_norm_optimize_c(cat, agr, bgr, compl, last_year, ref_mag, mmax=None, binw=0.1): """ - Variation on Poisson optimization of completeness windows - + Variation on Poisson optimization of completeness windows that considers events + within and outwith the completeness windows. + Probability is calculated as total probability of observing events within the + completeness windows + the probability of observing the events outside of + the completeness, given GR from specified a, b values. + + :param cat: + catalogue + :param aval: + GR a-value + :param bval: + GR b-value + :param compl: + completeness table + :param last_year: + end year to consider + :param ref_mag: + reference magnitude at which to perform calculations + :param mmax: + maximum magnitude + :param binw: + binwidth + :returns: + total Poisson probability of observing the events given the FMD. Larger + is better + """ mags = cat.data['magnitude'] @@ -343,21 +411,17 @@ def get_norm_optimize_c(cat, agr, bgr, compl, last_year, ref_mag, mmax=None, bin pocc = rates / sum(rates) - #prob = 1 # If using log (and not multiplicative) set initial prob to 0 prob = 0 first_year = min(yeas) for imag, mag in enumerate(mvals[:-1]): tot_n = len(mags[(mags >= mag) & (mags < mvals[imag+1])]) - #print(tot_n) if tot_n < 1: continue idxco = get_idx_compl(mag, compl) - - #print("total events in bin: ", tot_n) # if this magnitude bin is < mc in this window, nocc will be zero # Rather this disgards events outwith the completeness window, as it should! if idxco is None: @@ -367,27 +431,21 @@ def get_norm_optimize_c(cat, agr, bgr, compl, last_year, ref_mag, mmax=None, bin elif mag >= compl[idxco, 1]: idx = (mags >= mag) & (mags < mvals[imag+1]) & (yeas >= compl[idxco, 0]) - #print(idx) nocc_in = sum(idx) - #elif mag < min(cat.data['magnitude']): - # nocc_in = 0 + else: - print("how did this get here?", compl[idxco, 0], mag ) + print("event is outside magnitude limits", compl[idxco, 0], mag ) nocc_in = 0 #print(nocc_in) delta = (last_year - compl[idxco, 0]) + # events in bin before completeness idx = ((mags >= mag) & (mags < mvals[imag+1]) & (yeas < compl[idxco, 0])) - #idx2 = ((mags >= mag) & (mags < mvals[imag+1]) & (yeas > (compl[idxco, 0] - delta))) - #idx = ((mags >= mag) & (mags < mvals[imag+1]) & (yeas < compl[idxco, 0]) & (yeas > (compl[idxco, 0] - delta))) + nocc_out = sum(idx) - #print("nocc_in: ", nocc_in, "nocc_out: ", nocc_out, "total: ", nocc_in + nocc_out, "total events in bin: ", tot_n) - - #if mag < compl[idxco, 0]: - # nocc_out = 0 # also is this right? should yeas[idx] extend to years before compl[idxco, 0] too? if np.any(idx): @@ -398,26 +456,44 @@ def get_norm_optimize_c(cat, agr, bgr, compl, last_year, ref_mag, mmax=None, bin # Compute the duration for the completeness interval and the time # interval outside of completeness dur_out_compl = (compl[idxco, 0] - ylow) - # I think I want to limit this to the time interval of the completeness window + # I want to limit this to the time interval of the completeness window dur_compl = last_year - compl[idxco, 0] - + + # pmf for events within completeness windows pmf = poiss_prob_int_time(rates[imag], nocc_in, dur_compl, log_out = True) std_in = poisson.std(dur_compl*rates[imag]) - # cdf = poisson.cdf(nocc_out, delta*rates[imag]) - # std_out = poisson.std(dur_out_compl*rates[imag]) - #pmf_out = poisson.pmf(nocc_out, delta*rates[imag]) + # pmf outside pmf_out = poiss_prob_int_time(rates[imag], nocc_out, dur_out_compl, log_out = True) prob += pmf + (np.log(1.0) - pmf_out) - return prob def get_norm_optimize_poisson(cat, agr, bgr, compl, last_year, mmax=None, binw=0.1): """ - loop over the time increments - this is where we should be checking for Poisson after all + Alternative to get_norm_optimize_c that loops over the time increments + Probability is calculated as total probability of observing events within the + completeness windows + the probability of observing the events outside of + the completeness, given GR from specified a, b values. + + :param cat: + catalogue + :param aval: + GR a-value + :param bval: + GR b-value + :param compl: + completeness table + :param last_year: + end year to consider + :param binw: + binwidth + :returns: + total Poisson probability of observing the events given the FMD. Larger + is better + """ mags = cat.data['magnitude'] @@ -485,7 +561,19 @@ def get_norm_optimize_poisson(cat, agr, bgr, compl, last_year, mmax=None, binw=0 return prob def get_norm_optimize_d(cat, agr, bgr, compl, last_year, mmax=None, binw=0.1): - + ''' + + :param cat: + catalogue + :param aval: + GR a-value + :param bval: + GR b-value + :param compl: + completeness table + :param last_year: + end year to consider + ''' mags = cat.data['magnitude'] yeas = cat.data['year'] @@ -560,7 +648,28 @@ def get_norm_optimize_weichert(cat, agr, bgr, compl, last_year, mmax=None, binw= completeness which keeps more events will result in a larger likelihood. So when we try to use this to condition, we will find that smaller Mc are preffered earlier because this increases the Weichert likelihood. + This is why we should not compare likelihoods when using a different + number of events within the calculation (as we do here with different + completeness windows). This can still be interesting, but the above needs to be considered! + + :param cat: + catalogue + :param aval: + GR a-value + :param bval: + GR b-value + :param compl: + completeness table + :param last_year: + end year to consider + :param mmax: + maximum magnitude + :param binw: + binwidth + :returns: + norm - Weichert likelihood of observing the number of events across + all bins. Larger is better. """ mags = cat.data['magnitude'] @@ -631,6 +740,8 @@ def get_norm_optimize_gft(tcat, aval, bval, ctab, cmag, n_obs, t_per, last_year) """ Optimize fit using a version of the goodness-of-fit completeness approach (Wiemer and Wyss, 2000), using a parameter R to compare the goodness of fit. + A returned norm of 100 implies a perfect fit between observed and expected + events. :param aval: GR a-value @@ -642,11 +753,15 @@ def get_norm_optimize_gft(tcat, aval, bval, ctab, cmag, n_obs, t_per, last_year) An array with the magnitude values at the center of each occurrence bins :param t_per: + time period of completeness :param n_obs: Number of observations - :param last_year: + :param last_year: + last year to consider :param info: boolean controlling whether to print information as the function proceeds + :returns: + norm calculated with the Wiemer and Wyss (2000) 'R' parameter. Larger is better. """ # Select only events within 'complete' part of the catalogue occ = np.zeros((ctab.shape[0])) From b9c07313936bcc1840f8750ed26e3b8472c12b0f Mon Sep 17 00:00:00 2001 From: Kirsty Bayliss Date: Fri, 24 May 2024 13:16:35 +0200 Subject: [PATCH 02/23] can supply bin_width for completeness separately --- openquake/cat/completeness/analysis.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/openquake/cat/completeness/analysis.py b/openquake/cat/completeness/analysis.py index 77b61e0ff..0b751e357 100644 --- a/openquake/cat/completeness/analysis.py +++ b/openquake/cat/completeness/analysis.py @@ -196,7 +196,6 @@ def check_criterion(criterion, rate, previous_norm, tvars): tmp_rate = -1 norm = get_norm_optimize_gft(tcat, aval, bval, ctab, cmag, n_obs, t_per, last_year) - elif criterion == 'weichert': tmp_rate = -1 norm = get_norm_optimize_weichert(tcat, aval, bval, ctab, last_year) @@ -466,7 +465,10 @@ def read_compl_params(config): key = 'completeness' ms = np.array(config[key]['mags'], dtype=float) yrs = np.array(config[key]['years']) - bw = config.get('bin_width', 0.1) + try: + bw = np.array(config[key]['bin_width'], dtype =float) + except: + bw = config.get('bin_width', 0.1) r_m = config[key].get('ref_mag', 5.0) r_up_m = config[key].get('ref_upp_mag', None) bmin = config[key].get('bmin', 0.8) From e048c276c77ece36912b5af6316dc056ab745ca3 Mon Sep 17 00:00:00 2001 From: Kirsty Bayliss Date: Wed, 29 May 2024 10:16:05 +0200 Subject: [PATCH 03/23] updating 'optimize' norm to use specified binwidths --- openquake/cat/completeness/norms.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/openquake/cat/completeness/norms.py b/openquake/cat/completeness/norms.py index 8f714313f..0f039e529 100644 --- a/openquake/cat/completeness/norms.py +++ b/openquake/cat/completeness/norms.py @@ -98,7 +98,7 @@ def get_completeness_matrix(tcat, ctab, mbinw, ybinw): return oin, out, cmags, cyeas -def get_norm_optimize(tcat, aval, bval, ctab, cmag, n_obs, t_per, last_year, +def get_norm_optimize(tcat, aval, bval, ctab, binw, cmag, n_obs, t_per, last_year, info=False): """ This calculates the difference between the number of observed vs computed events in each completeness bin. @@ -111,6 +111,9 @@ def get_norm_optimize(tcat, aval, bval, ctab, cmag, n_obs, t_per, last_year, GR b-value :param ctab: completeness table + :param binw: + binwidth for completeness analysis, specified in toml. This is the width of bins for checking fit + (previously for optimize this was = 1) :param cmag: An array with the magnitude values at the center of each occurrence bins. Output from hmtk.seismicity.occurrence.utils.get_completeness_counts @@ -130,7 +133,8 @@ def get_norm_optimize(tcat, aval, bval, ctab, cmag, n_obs, t_per, last_year, occ = np.zeros((ctab.shape[0])) dur = np.zeros((ctab.shape[0])) - mags = np.array(list(ctab[:, 1])+[10]) + #mags = np.array(list(ctab[:, 1])+[10]) + mags = np.array(np.arange(ctab[:, 1], 10, binw)) for i, mag in enumerate(mags[:-1]): idx = (cmag >= mags[i]) & (cmag < mags[i+1]) From bdd56f263e2bb86f39a9af0d617da1baaaacc841 Mon Sep 17 00:00:00 2001 From: Kirsty Bayliss Date: Tue, 23 Jul 2024 11:45:38 +0200 Subject: [PATCH 04/23] updates to take all files in a folder and manage point sources --- openquake/man/source_tests.py | 77 ++++++++++++++++++++++++++++++----- 1 file changed, 67 insertions(+), 10 deletions(-) diff --git a/openquake/man/source_tests.py b/openquake/man/source_tests.py index 201537880..f82a7f607 100644 --- a/openquake/man/source_tests.py +++ b/openquake/man/source_tests.py @@ -1,4 +1,5 @@ import os +import re import glob import numpy as np import pandas as pd @@ -50,9 +51,44 @@ def get_params(ssm, mmin = 0, total_for_zone = False): else: return np.array(data) +def get_params_points(ssm, mmin = 0, total_for_zone = False): + ''' + Retrieve source parameters from xml source file for point sources + Specify total_for_zone to retain only total a value summed over all points + + :param ssm: + seismic source model retreived from xml file with to_python + :param mmin: + minimum magnitude to consider + :param total_for_zone: + if True, returns the total a-value for the zone. + if False, returns a-values at point locations + ''' + total_rate_above_zero = 0 + data = [] + + import pdb + #pdb.set_trace() + for ps in ssm[0]: + agr = ps.mfd.a_val + bgr = ps.mfd.b_val + mmax = ps.mfd.max_mag + lo = ps.location.longitude + la = ps.location.latitude + log_rate_above = agr - bgr * mmin + total_rate_above_zero += 10**(agr) + + data.append([lo, la, log_rate_above, bgr, mmax]) + + if total_for_zone == True: + #total_rate_above += 10**(agr - bgr * m_min) + return np.log10(total_rate_above_zero), bgr, mmax + else: + return np.array(data) + def plot_sources(folder_name, region, mmin, fig_folder, rate_range = [-9, 3, 0.2], mmax_range = [6.0, 9.0, 0.2], - sconv = DEFAULT, poly_name = '', plot_poly = False): + sconv = DEFAULT, poly_name = '', plot_poly = False, pointsource = False): ''' Create plots of source rates and mmax for each source in folder_name @@ -74,6 +110,8 @@ def plot_sources(folder_name, region, mmin, fig_folder, location of file containing the model source polygon (optional) :param plot_poly: boolean specifying if polygon outline should be plotted + :param pointsource: + alternative where point sources are used instead of multipoint ''' path = pathlib.Path(folder_name) @@ -96,7 +134,7 @@ def plot_sources(folder_name, region, mmin, fig_folder, poly["x"] = poly.representative_point().x poly["y"] = poly.representative_point().y - for fname in sorted(path.glob('src*.xml')): + for fname in sorted(path.glob('src_*.xml')): ssm = to_python(fname, sconv) fig_a = pygmt.Figure() fig_a.basemap(region=region, projection="M15c", frame=True) @@ -108,12 +146,16 @@ def plot_sources(folder_name, region, mmin, fig_folder, vmin = +1e10 vmax = -1e10 - - for grp in ssm: - for src in grp: - name = src.name - data = get_params(ssm, mmin=mmin, total_for_zone = False) + if pointsource == False: + data = get_params(ssm, mmin=mmin, total_for_zone = False) + for grp in ssm: + for src in grp: + name = src.name + else: + data = get_params_points(ssm, mmin=mmin, total_for_zone = False) + name = ssm.name + vmin = np.min([vmin, min(data[:,2])]) vmax = np.max([vmin, max(data[:,2])]) @@ -162,7 +204,7 @@ def simulate_occurrence(agr, bgr, rate, minmag, mmax, time_span, N=2000): num_occ = np.random.poisson(rate*time_span, N) return(num_occ) -def occurence_table(path_oq_input, path_to_subcatalogues, minmag, minyear, maxyear, N, src_ids, sconv = DEFAULT): +def occurence_table(path_oq_input, path_to_subcatalogues, minmag, minyear, maxyear, N, src_ids = [], sconv = DEFAULT, pointsource = False): ''' Check number of events expected from the source model against the number of observations. Uses N samples from a Poisson distribution with rate from source a and b value. @@ -184,10 +226,21 @@ def occurence_table(path_oq_input, path_to_subcatalogues, minmag, minyear, maxye list of sources to use :param sconv: source converter object specifying model setup + :param pointsource: + specify True if using set point source model rather than ''' table = [] time_span = maxyear - minyear + if len(src_ids) == 0: + src_ids = [] + # if source ids are not specified, take them from names of files in the specified folder + dir_URL = Path(path_oq_input) + filename_list = [file.name for file in dir_URL.glob("*.xml")] + for f in filename_list: + m = re.search(r'(?<=_)\w+', f) + src_ids.append(m.group(0)) + for src_id in sorted(src_ids): fname_src = os.path.join(path_oq_input, "src_{:s}.xml".format(src_id)) @@ -199,7 +252,11 @@ def occurence_table(path_oq_input, path_to_subcatalogues, minmag, minyear, maxye df = df.loc[(df.year >= minyear) & (df.year <= maxyear)] obs = len(df) - agr, bgr, mmax = get_params(ssm, minmag, total_for_zone = True) + if pointsource == False: + agr, bgr, mmax = get_params(ssm, minmag, total_for_zone = True) + else: + agr, bgr, mmax = get_params_points(ssm, minmag, total_for_zone = True) + rate = 10.0**(agr-minmag*bgr)-10.0**(agr-mmax*bgr) num_occ_per_time_span = simulate_occurrence(agr, bgr, rate, minmag, mmax, time_span, N) @@ -261,4 +318,4 @@ def source_info_table(folder_name, sconv = DEFAULT): ] sdata.loc[len(sdata.index)] = row - print(tabulate(sdata, headers="keys", tablefmt="psql")) \ No newline at end of file + print(tabulate(sdata, headers="keys", tablefmt="psql")) From 5fed7e95efef477fa11f89844bbc811311029602 Mon Sep 17 00:00:00 2001 From: Kirsty Bayliss Date: Tue, 23 Jul 2024 14:38:06 +0200 Subject: [PATCH 05/23] adding option to plot all sources in one plot --- openquake/man/source_tests.py | 100 ++++++++++++++++++++++++++++++++++ 1 file changed, 100 insertions(+) diff --git a/openquake/man/source_tests.py b/openquake/man/source_tests.py index f82a7f607..f7b815d26 100644 --- a/openquake/man/source_tests.py +++ b/openquake/man/source_tests.py @@ -184,6 +184,106 @@ def plot_sources(folder_name, region, mmin, fig_folder, out = os.path.join(fig_folder, 'mmax', name+'_mmax.png') fig_b.savefig(out) +def plot_all_sources(folder_name, region, mmin, fig_folder, + rate_range = [-9, 3, 0.2], mmax_range = [6.0, 9.0, 0.2], + sconv = DEFAULT, poly_name = '', plot_poly = False, pointsource = False, plot_fname = ""): + ''' + Plots a single plot of source rates and mmax using each source in folder_name + + :param folder_name: + folder containing source xml files + :param region: + region for pygmt plotting + :param mmin: + minimum magnitude to be used in model + :param fig_folder: + folder in which to store plots + :param rate_range: + range of rate to use in colour scale, + specified as [min, max, interval] + :param mmax_range: + range of mmax for colour scale, + specified as [min, max, interval] + :param poly_name: + location of file containing the model source polygon (optional) + :param plot_poly: + boolean specifying if polygon outline should be plotted + :param pointsource: + alternative where point sources are used instead of multipoint + ''' + path = pathlib.Path(folder_name) + + # make rate and mmax folders if they do not exist + if os.path.exists(os.path.join(os.path.join(fig_folder, 'rate'))): + print("found folders at ", os.path.join(os.path.join(fig_folder, 'rate'))) + else: + os.mkdir(os.path.join(fig_folder, 'rate')) + os.mkdir(os.path.join(fig_folder, 'mmax')) + + # set up colour scales + cpt_rate = os.path.join(fig_folder, 'rate.cpt') + pygmt.makecpt(cmap="turbo", series=rate_range, output=cpt_rate) + cpt_mmax = os.path.join(fig_folder, 'mmax.cpt') + pygmt.makecpt(cmap="rainbow", series=mmax_range, output=cpt_mmax) + + if poly_name: + # set plotting polygon + poly = gpd.read_file(poly_name) + poly["x"] = poly.representative_point().x + poly["y"] = poly.representative_point().y + + dataset = np.empty([1,5]) + + for fname in sorted(path.glob('src_*.xml')): + ssm = to_python(fname, sconv) + fig_a = pygmt.Figure() + fig_a.basemap(region=region, projection="M15c", frame=True) + fig_a.coast(land="grey", water="white") + + fig_b = pygmt.Figure() + fig_b.basemap(region=region, projection="M15c", frame=True) + fig_b.coast(land="grey", water="white") + + vmin = +1e10 + vmax = -1e10 + + if pointsource == False: + data = get_params(ssm, mmin=mmin, total_for_zone = False) + dataset = np.concatenate((dataset, data)) + + else: + data = get_params_points(ssm, mmin=mmin, total_for_zone = False) + dataset = np.concatenate((dataset, data)) + + vmin = np.min([vmin, min(dataset[:,2])]) + vmax = np.max([vmin, max(dataset[:,2])]) + + fig_a.plot(x=dataset[:,0], + y=dataset[:,1], + style="h0.5", + color=dataset[:, 2], + cmap=cpt_rate) + + fig_b.plot(x=dataset[:,0], + y=dataset[:,1], + style="h0.5", + color=dataset[:, 4], + cmap=cpt_mmax) + + if plot_poly == True: + fig_a.plot(data=poly, pen=".5p,black") + fig_b.plot(data=poly, pen=".5p,black") + + fig_a.colorbar(frame=f'af+l"Log((N(m)>{mmin}))"', cmap=cpt_rate) + fig_b.colorbar(frame='af+l"Mmax"', cmap=cpt_mmax) + + out = os.path.join(fig_folder, 'rate', plot_fname+'_rate.png') + fig_a.savefig(out) + + out = os.path.join(fig_folder, 'mmax', plot_fname+'_mmax.png') + fig_b.savefig(out) + + def simulate_occurrence(agr, bgr, rate, minmag, mmax, time_span, N=2000): ''' Simulate number of occurrences from a Poisson distribution given the FMD parameters From 993c08fae82e74466c66be28c8e5098f458525cf Mon Sep 17 00:00:00 2001 From: Kirsty Bayliss Date: Wed, 24 Jul 2024 16:37:51 +0200 Subject: [PATCH 06/23] compare two source models --- openquake/man/source_tests.py | 88 +++++++++++++++++++++++++++++++++++ 1 file changed, 88 insertions(+) diff --git a/openquake/man/source_tests.py b/openquake/man/source_tests.py index f7b815d26..a98353abb 100644 --- a/openquake/man/source_tests.py +++ b/openquake/man/source_tests.py @@ -419,3 +419,91 @@ def source_info_table(folder_name, sconv = DEFAULT): sdata.loc[len(sdata.index)] = row print(tabulate(sdata, headers="keys", tablefmt="psql")) + + +def compare_sources_plot(src1_fname, src2_fname, region, mmin, fig_folder, + rate_range = [-9, 3, 0.2], mmax_range = [6.0, 9.0, 0.2], + sconv = DEFAULT, poly_name = '', plot_poly = False, pointsource = False, plot_fname = ""): + ''' + Creates a plot showing the differences between two source models. + + :param src1_fname: + location of a file describing a point or multipoint source model + :param src2_fname: + location of a file describing a point or multipoint source model + :param region: + region for pygmt plotting + :param mmin: + minimum magnitude to be used in model + :param fig_folder: + folder in which to store plots + :param rate_range: + range of rate to use in colour scale, + specified as [min, max, interval] + :param mmax_range: + range of mmax for colour scale, + specified as [min, max, interval] + :param poly_name: + location of file containing the model source polygon (optional) + :param plot_poly: + boolean specifying if polygon outline should be plotted + :param pointsource: + alternative where point sources are used instead of multipoint + ''' + + # set up colour scales + cpt_rate = os.path.join(fig_folder, 'rate.cpt') + pygmt.makecpt(cmap="turbo", series=rate_range, output=cpt_rate) + + if poly_name: + # set plotting polygon + poly = gpd.read_file(poly_name) + poly["x"] = poly.representative_point().x + poly["y"] = poly.representative_point().y + + ssm1 = to_python(src1_fname, sconv) + ssm2 = to_python(src2_fname, sconv) + + if pointsource == False: + data1 = get_params(ssm1, mmin=mmin, total_for_zone = False) + data2 = get_params(ssm2, mmin=mmin, total_for_zone = False) + + else: + data1 = get_params_points(ssm1, mmin=mmin, total_for_zone = False) + data2 = get_params_points(ssm2, mmin=mmin, total_for_zone = False) + + #assert data1[:,0] == data2[:,0] + rate_diff = data1[:,2] - data2[:,2] + fig = pygmt.Figure() + fig.basemap(region=region, projection="M15c", frame=True) + fig.coast(land="grey", water="white") + + fig.plot(x=data1[:,0], + y=data1[:,1], + style="h0.75", + color=rate_diff, + cmap=cpt_rate) + + fig.colorbar(frame=f'af+l"Difference in log((N(m)>{mmin}))"', cmap=cpt_rate) + + out = os.path.join(fig_folder, plot_fname+'_rate_diff.png') + fig.savefig(out) + ''' + path1 = pathlib.Path(src1_fname) + path2 = pathlib.Path(src2_fname) + + if os.path.isdir(src1_fname): + assert os.path.isdir(path2) + + dataset1 = np.empty([1,5]) + for fname in sorted(path.glob('src_*.xml')): + ssm = to_python(fname, sconv) + + if pointsource == False: + data = get_params(ssm, mmin=mmin, total_for_zone = False) + dataset = np.concatenate((dataset1, data)) + + else: + data = get_params_points(ssm, mmin=mmin, total_for_zone = False) + dataset1 = np.concatenate((dataset1, data)) + ''' \ No newline at end of file From 1c9cca3c93c494716f70602686b235c74437ccad Mon Sep 17 00:00:00 2001 From: Kirsty Bayliss Date: Mon, 26 Aug 2024 12:13:38 +0200 Subject: [PATCH 07/23] adding more tests for homogenisation as useful examples --- openquake/cat/tests/hmg_test.py | 74 +++++++++++++++++++++++++++++++++ 1 file changed, 74 insertions(+) diff --git a/openquake/cat/tests/hmg_test.py b/openquake/cat/tests/hmg_test.py index 4667a8e25..674526ba8 100644 --- a/openquake/cat/tests/hmg_test.py +++ b/openquake/cat/tests/hmg_test.py @@ -32,6 +32,8 @@ import toml from openquake.cat.hmg.hmg import process_dfs, process_magnitude +from openquake.cat.isf_catalogue import ISFCatalogue +from openquake.cat.parsers.isf_catalogue_reader import ISFReader aeq = np.testing.assert_equal aaeq = np.testing.assert_almost_equal @@ -47,7 +49,31 @@ sigma = [0.283, 0.283] """ +SETTINGS_ISF1 = """ +[magnitude.TIR.Ml] +low_mags = [0.0] +conv_eqs = ["0.96 * m + 0.23"] +sigma = [0.1] +""" + +SETTINGS_ISF2 = """ +[magnitude.AFAD.MW] +low_mags = [0.0] +conv_eqs = ["m"] +sigma = [0.1] + +[magnitude.ATH.ML] +low_mags = [0.0] +conv_eqs = ["1.1 * m -0.2"] +sigma = [0.1] + +[magnitude.TIR.Ml] +low_mags = [0.0] +conv_eqs = ["0.96 * m + 0.23"] +sigma = [0.1] + +""" class HomogeniseNEICmbTestCase(unittest.TestCase): @@ -77,3 +103,51 @@ def test_case01(self): results = save.magMw.values expected = [4.36, 4.76, 6.8 ] aaeq(results, expected, decimal=6) + +class HomogeniseIsfTestCase(unittest.TestCase): + + def setUp(self): + self.fname_isf6 = os.path.join(BASE_PATH, 'data', 'cat06.isf') + + def test_case01(self): + """ + tests homogenisation of an isf file with multiple agencies + only one magnitude conversion + """ + #parser = ISFReader(self.fname_isf2) + cat = ISFCatalogue(self.fname_isf6, name = "isf") + parser = ISFReader(self.fname_isf6) + catisf = parser.read_file("tisf", "Test isf") + + otab, mtab = catisf.get_origin_mag_tables() + work = pd.merge(pd.DataFrame(otab), pd.DataFrame(mtab), on=["eventID", "originID"]) + + rules = toml.loads(SETTINGS_ISF1) + save, work = process_magnitude(work, rules['magnitude']) + + results = save.magMw.values + expected = [3.878, 3.686, 3.686, 2.822, 3.302, 3.494, 4.07 ] + aaeq(results, expected, decimal=6) + + + def test_case02(self): + + """ + tests homogenisation of an isf file with multiple agencies + use multiple agencies + """ + #parser = ISFReader(self.fname_isf2) + cat = ISFCatalogue(self.fname_isf6, name = "isf") + parser = ISFReader(self.fname_isf6) + catisf = parser.read_file("tisf", "Test isf") + + otab, mtab = catisf.get_origin_mag_tables() + work = pd.merge(pd.DataFrame(otab), pd.DataFrame(mtab), on=["eventID", "originID"]) + + rules = toml.loads(SETTINGS_ISF2) + save, work = process_magnitude(work, rules['magnitude']) + + results = save.magMw.values + # Note magnitudes will be in the order selected! + expected = [4.0, 3.43, 3.32, 3.43, 2.66, 3.302, 3.494 ] + aaeq(results, expected, decimal=6) From ece90e8becdd52fed2ebb40a439cdf2821fc4d6d Mon Sep 17 00:00:00 2001 From: Kirsty Bayliss Date: Mon, 26 Aug 2024 12:14:13 +0200 Subject: [PATCH 08/23] tweak to stop depths = 0 returning depths < 0 warning --- openquake/cat/isf_catalogue.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/openquake/cat/isf_catalogue.py b/openquake/cat/isf_catalogue.py index b24a49f95..2c85e7911 100755 --- a/openquake/cat/isf_catalogue.py +++ b/openquake/cat/isf_catalogue.py @@ -1210,8 +1210,9 @@ def get_origin_mag_tables(self): orig.location.depth == '' or orig.location.depth is None): depth = np.NaN - elif orig.location.depth < 0.01: + elif orig.location.depth < -0.01: depth = orig.location.depth + print('Depth:', orig.location.depth) fmt = "Warning, depth <= 0.0 (id:{:s})" warnings.warn(fmt.format(eq.id)) elif orig.location.depth: From f1ea7b189135dae05a411288800e5ef21c71938d Mon Sep 17 00:00:00 2001 From: Kirsty Bayliss Date: Mon, 26 Aug 2024 12:15:42 +0200 Subject: [PATCH 09/23] explain preferences for homogenisation magnitudes --- docsrc/contents/cat.rst | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/docsrc/contents/cat.rst b/docsrc/contents/cat.rst index 495dfeff8..34772cc75 100644 --- a/docsrc/contents/cat.rst +++ b/docsrc/contents/cat.rst @@ -194,7 +194,7 @@ This would give us a different fit to our data and a different equation to suppl Where there are not enough events to allow for a direct regression or we are unhappy with the fit for our data, there are many conversions in the literature which may be useful. This process may take some revising and iterating - it is sometimes very difficult to identify a best fit, especially where we have few datapoints or highly uncertain data. Once we are happy with the fits to our data, we can add the regression equation to the homogenisation .toml file. This process should be repeated for every magnitude we wish to convert to Mw. -The final homogenisation step itself is also controlled by a toml file, where each observed magnitude is specified individually and the regression coefficients and uncertainty are included. It is also necessary to specify a hierarchy of catalogues so that a preferred catalogue is used for the magnitude where the event has multiple entries. In the example below, we merge the ISCGEM and a local catalogue, preferring ISCGEM magnitudes where available as specified in the ranking. Because the ISCGEM already provides magnitudes in Mw, we simply retain all Mw magnitudes from ISCGEM. In this example, our local catalogue has two different magnitude types for which we have derived a regression. We specify how to convert to the standardised Mw from the local.mw and the standard deviations, which are outputs of the fitting we carried out above. +The final homogenisation step itself is also controlled by a toml file, where each observed magnitude is specified individually and the regression coefficients and uncertainty are included. It is also necessary to specify a hierarchy of catalogues so that a preferred catalogue is used for the magnitude where the event has multiple entries. If you have an isf format catalogue, you can also specify here the hierarchy of individual agencies (or authors in the isf format) within the catalogue. In the example below, we merge the ISCGEM and a local catalogue, preferring ISCGEM magnitudes where available as specified in the ranking. Because the ISCGEM already provides magnitudes in Mw, we simply retain all Mw magnitudes from ISCGEM. In this example, our local catalogue has two different magnitude types for which we have derived a regression. We specify how to convert to the standardised Mw from the local.mw and the standard deviations, which are outputs of the fitting we carried out above. .. code-block:: ini @@ -229,6 +229,7 @@ The final homogenisation step itself is also controlled by a toml file, where ea conv_eqs = ["0.1928 + 0.9757 * m"] std_devs = [0.0091, 0.0016] +The order of conversions in the list will determine priority for conversion, so for the local events we will first convert all events with mw magnitudes and then use mww only where the mw magnitudes are not available, and the local conversions will not be used when we have an ISCGEM Mw. In this way we can specify hierarchies for both the agencies and the magnitudes. The actual homogenisation step is carried out by calling :code:`oqm cat homogenise $ARG1 $ARG2 $ARG3` as in the bash script example, where $ARG1 is the homogenisation toml file and and $ARG2 and $ARG3 are the hdf5 file outputs from the merge step, describing the origins and magnitude information for the merged catalogue respectively. From 7695aabe95bb99ccc28c0fab0655d15fe14650bb Mon Sep 17 00:00:00 2001 From: Kirsty Bayliss Date: Mon, 26 Aug 2024 12:19:45 +0200 Subject: [PATCH 10/23] fix positional argument issues --- openquake/mbi/wkf/create_smoothing_per_zone.py | 7 ++++--- openquake/wkf/seismicity/smoothing.py | 9 +++++---- 2 files changed, 9 insertions(+), 7 deletions(-) diff --git a/openquake/mbi/wkf/create_smoothing_per_zone.py b/openquake/mbi/wkf/create_smoothing_per_zone.py index ef9d79c9b..6f6b11ffe 100755 --- a/openquake/mbi/wkf/create_smoothing_per_zone.py +++ b/openquake/mbi/wkf/create_smoothing_per_zone.py @@ -5,15 +5,16 @@ from openquake.wkf.seismicity.smoothing import create_smoothing_per_zone -def main(fname_points: str, fname_polygons: str, folder_out: str='/tmp',*, - skip=[], use = []): +def main(fname_points: str, fname_polygons: str, folder_out: str='/tmp', + skip: str=[], use:str=[]): create_smoothing_per_zone(fname_points, fname_polygons, folder_out, skip, use) main.fname_points = ".csv file created by the smoothing code" main.fname_polygons = "Shapefile with the polygons" main.folder_out = "The name of the output folder where to save .csv files" -main.skip = "A string containing a list of source IDs" +main.skip = "A string containing a list of source IDs to skip" +main.use = "A string containing a list of source IDs to use" if __name__ == '__main__': sap.run(main) diff --git a/openquake/wkf/seismicity/smoothing.py b/openquake/wkf/seismicity/smoothing.py index bd15fccbe..2f8f71b2d 100755 --- a/openquake/wkf/seismicity/smoothing.py +++ b/openquake/wkf/seismicity/smoothing.py @@ -11,13 +11,13 @@ from shapely import wkt -def main(fname_points: str, fname_polygons: str, folder_out: str='/tmp', - skip=[]): - create_smoothing_per_zone(fname_points, fname_polygons, folder_out, skip) +def main(fname_points: str, fname_polygons: str, folder_out: str='/tmp', *, + skip: str=[], use: str=[]): + create_smoothing_per_zone(fname_points, fname_polygons, folder_out, skip, use) def create_smoothing_per_zone(fname_points: str, fname_polygons: str, - folder_out: str='/tmp', *, skip=[], use = []): + folder_out: str='/tmp', skip:str=[], use:str = []): """ Creates subsets of points, one for each of the polygons included in the `fname_polygons` shapefile. The attibute table must have an 'id' attribute. @@ -76,6 +76,7 @@ def create_smoothing_per_zone(fname_points: str, fname_polygons: str, main.fname_polygons = "Shapefile with the polygons" main.folder_out = "The name of the output folder where to save .csv files" main.skip = "A string containing a list of source IDs" +main.use = "A string containg a list of source IDs to use" if __name__ == '__main__': sap.run(main) From 83bccc75dee39b6a8f76bb31b19ba8cd566bd69d Mon Sep 17 00:00:00 2001 From: Kirsty Bayliss Date: Mon, 26 Aug 2024 12:21:38 +0200 Subject: [PATCH 11/23] add option to provide lower magnitude != ref mag --- openquake/cat/completeness/analysis.py | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/openquake/cat/completeness/analysis.py b/openquake/cat/completeness/analysis.py index d0dc80154..8bcd2e374 100644 --- a/openquake/cat/completeness/analysis.py +++ b/openquake/cat/completeness/analysis.py @@ -233,7 +233,7 @@ def _make_ctab(prm, years, mags): return 'skip' -def _completeness_analysis(fname, years, mags, binw, ref_mag, ref_upp_mag, +def _completeness_analysis(fname, years, mags, binw, ref_mag, mag_low, ref_upp_mag, bgrlim, criterion, compl_tables, src_id=None, folder_out_figs=None, rewrite=False, folder_out=None): @@ -247,6 +247,9 @@ def _completeness_analysis(fname, years, mags, binw, ref_mag, ref_upp_mag, :param ref_mag: The reference magnitude used to compute the rate and select the completeness table + :param mag_low: + The lowest magnitude to include in rate calculations. + Will default to ref_mag if not provided. :param ref_upp_mag: The reference upper magnitude limit used to compute the rate and select the completeness table @@ -264,7 +267,7 @@ def _completeness_analysis(fname, years, mags, binw, ref_mag, ref_upp_mag, # Checking input if criterion not in ['match_rate', 'largest_rate', 'optimize', 'weichert', - 'poisson', 'optimize_a', 'optimize_b', 'optimize_d']: + 'poisson', 'optimize_a', 'optimize_b', 'optimize_c','optimize_d']: raise ValueError('Unknown optimization criterion') tcat = _load_catalogue(fname) @@ -273,7 +276,7 @@ def _completeness_analysis(fname, years, mags, binw, ref_mag, ref_upp_mag, # Info # Should have option to specify a mag_low != ref_mag - mag_low = ref_mag + #mag_low = lower_mag idx = tcat.data["magnitude"] >= mag_low fmt = 'Catalogue contains {:d} events equal or above {:.1f}' print('\nSOURCE:', src_id) @@ -468,13 +471,14 @@ def read_compl_params(config): yrs = np.array(config[key]['years']) bw = config.get('bin_width', 0.1) r_m = config[key].get('ref_mag', 5.0) + m_low = config[key].get('mag_low', r_m) r_up_m = config[key].get('ref_upp_mag', None) bmin = config[key].get('bmin', 0.8) bmax = config[key].get('bmax', 1.2) # Options: 'largest_rate', 'match_rate', 'optimize' crit = config[key].get('optimization_criterion', 'optimize') - return ms, yrs, bw, r_m, r_up_m, bmin, bmax, crit + return ms, yrs, bw, r_m, m_low, r_up_m, bmin, bmax, crit def read_compl_data(folder_in): @@ -513,7 +517,7 @@ def completeness_analysis(fname_input_pattern, f_config, folder_out_figs, # Loading configuration config = toml.load(f_config) - ms, yrs, bw, r_m, r_up_m, bmin, bmax, crit = read_compl_params(config) + ms, yrs, bw, r_m, m_low, r_up_m, bmin, bmax, crit = read_compl_params(config) compl_tables, mags_chk, years_chk = read_compl_data(folder_in) @@ -544,7 +548,7 @@ def completeness_analysis(fname_input_pattern, f_config, folder_out_figs, else: var = {} - res = _completeness_analysis(fname, yrs, ms, bw, r_m, + res = _completeness_analysis(fname, yrs, ms, bw, r_m, m_low, r_up_m, [bmin, bmax], crit, compl_tables, src_id, folder_out_figs=folder_out_figs, From 1c41f38a80b9a601173a8b25b1a33bed0b93b0da Mon Sep 17 00:00:00 2001 From: Kirsty Bayliss Date: Mon, 26 Aug 2024 19:32:38 +0200 Subject: [PATCH 12/23] remove extra print statements --- openquake/cat/completeness/generate.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/openquake/cat/completeness/generate.py b/openquake/cat/completeness/generate.py index e98b82de2..2608a865f 100755 --- a/openquake/cat/completeness/generate.py +++ b/openquake/cat/completeness/generate.py @@ -108,7 +108,8 @@ def _get_completenesses(mags, years, folder_out=None, num_steps=0, min_mag = min(mags) else: min_mag = min_mag_compl - + + print("minimum magnitude is ", min_mag) if len(np.where(min_mag <= mags)) < 1: msg = 'None of the magnitude intervals above the min_mag_compl' raise ValueError(msg) From e35c2164797e9567428ea6e3b859b56f313eff45 Mon Sep 17 00:00:00 2001 From: Kirsty Bayliss Date: Mon, 26 Aug 2024 19:33:12 +0200 Subject: [PATCH 13/23] remove print statements --- openquake/cat/completeness/analysis.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/openquake/cat/completeness/analysis.py b/openquake/cat/completeness/analysis.py index 0b751e357..942ceb48c 100644 --- a/openquake/cat/completeness/analysis.py +++ b/openquake/cat/completeness/analysis.py @@ -309,7 +309,6 @@ def _completeness_analysis(fname, years, mags, binw, ref_mag, ref_upp_mag, print(f'Iteration: {iper:05d} norm: {norm:12.6e}', end="\r") ctab = _make_ctab(prm, years, mags) - print(ctab) if isinstance(ctab, str): continue @@ -552,7 +551,7 @@ def completeness_analysis(fname_input_pattern, f_config, folder_out_figs, folder_out_figs=folder_out_figs, folder_out=folder_out, rewrite=False) - #print(len(res)) + if len(res) == 0: continue From 7c9b06fd5d1468dd4d4e6940ff2a42288e229575 Mon Sep 17 00:00:00 2001 From: Kirsty Bayliss Date: Thu, 5 Sep 2024 14:17:32 +0200 Subject: [PATCH 14/23] adding missing importlib import --- openquake/mbi/wkf/create_nrml_sources.py | 1 + 1 file changed, 1 insertion(+) diff --git a/openquake/mbi/wkf/create_nrml_sources.py b/openquake/mbi/wkf/create_nrml_sources.py index 56491584c..557e53389 100755 --- a/openquake/mbi/wkf/create_nrml_sources.py +++ b/openquake/mbi/wkf/create_nrml_sources.py @@ -10,6 +10,7 @@ from glob import glob from openquake.wkf.utils import create_folder +import importlib from openquake.baselib import sap from openquake.hazardlib.sourcewriter import write_source_model from openquake.hazardlib.source import PointSource, MultiPointSource From f7f859db44ebb37f7b59555c4923c2c29f237da9 Mon Sep 17 00:00:00 2001 From: Kirsty Bayliss Date: Mon, 9 Sep 2024 18:15:04 +0200 Subject: [PATCH 15/23] updates to subduction documentation --- docsrc/contents/sub.rst | 167 ++++++++++++++++++++++++++++++---------- 1 file changed, 128 insertions(+), 39 deletions(-) diff --git a/docsrc/contents/sub.rst b/docsrc/contents/sub.rst index 60989f27e..daab5c466 100644 --- a/docsrc/contents/sub.rst +++ b/docsrc/contents/sub.rst @@ -65,7 +65,7 @@ Herein we provide a brief description of the various steps. Note that we use the dep_max = 700 -2. Create a pickled version of your hmtk formatted catalog:: +2. Create a pickled version of your hmtk formatted catalog (NB: the catalogue must be in a hmtk format):: > pickle_catalogue.py ./catalogues/cac.cat` @@ -96,15 +96,44 @@ Once launched, by clicking on the image it is possible to digitize a sequence of Second approach =============== -The second approach proposed is simpler than the first one. At the beginning, it requires to complete point 1 and point 3 described in the `first approach`_ section. Once we have a configuration file and a set of cross sections ready we can complete the construction of the set of profiles with the following command:: +The second approach proposed is simpler than the first one, but it uses directly subduction contours from Slab2.0 (Hayes et al). - > sub_create_sections_from_slab.py +1. Set-up configuration files. This approach requires an input toml file describing the locations of the Slab2.0 data files and some other input parameters as below -Where: +.. code-block:: toml + # Locations for strike and depth information from Slab2pt0 + fname_str ='./slab2pt0/Slab2Distribute_Mar2018/str_grd/kur_slab2_str_02.24.18.grd' + fname_dep ='./slab2pt0/Slab2Distribute_Mar2018/dep_grd/kur_slab2_dep_02.24.18.grd' + # Spacing to use for the profiles (in km) + spacing = 100.0 + # Folder to save profiles to + folder_out = './cs' + # Optional: filename to save a figure of the profiles + fname_fig = 'kur.png' + + +2. Create a pickled version of your hmtk formatted catalog and make sure this is located as in your toml file:: + + > pickle_catalogue.py ./catalogues/cac.cat` + +3. Make subduction profiles. In this approach we have two choices for creating profiles across the subduction zone. The first is an automatic procedure that defines cross-sections perpendicular to the average strike of the subduction zone with the spacing specified in the toml:: + + > sub_get_profiles_from_slab2pt0 + +Alternatively, we can manually define a set of profiles from a geojson file. This method is necessary in areas where a subduction zone is curved, because parallel profiles defined by the first approach will cross the slab at strange angles and the resulting 3D geometry will be misshapen. In this case, the geojson file should be specified in the toml with `fname_geojson` and the function to create the profiles is:: + + >sub get_profiles_from_slab2pt0 + +As with the first approach, you can plot maps and cross_sections with the same commands. +4. Check the new set of traces in a map with the command:: -- ```` is the name of the file -- ```` is the name of the folder where to write the profiles -- ```` is the name of the file (produced by ``create_multiple_cross_sections.py``) with information aboout the traces of the cross-sections. + > plot_multiple_cross_sections_map.py ./ini/central_america.ini cs_traces.cs + +5. Create one .pdf file for each cross-section with the available information: e.g., earthquake hypocentres, focal mechanism, slab 1.0 geometry, CRUST 1.0 Moho:: + + > plot_multiple_cross_sections.py cs_traces.cs + +This command will produce as many ``.pdf`` files as the number of cross-sections specified in the ``.cs`` file Building the top of the slab geometry ************************************* @@ -116,12 +145,13 @@ This part of the procedure can be completed by running the 1. Build the surface of the subduction interface using ``create_2pt5_model.py``. The input information in this case is: - The name of the folder ```` containing the ``cs_`` files created using either the procedure described in the `first approach`_ or `first approach`_ section; + - The output profile folder ````; - The maximum sampling distance along a trace [km]; - - The output folder ````; + Example:: - > create_2pt5_model.py + > create_2pt5_model.py The output is a set of interpolated profiles and edges that can be used to create a complex fault source for the OpenQuake engine. The results of the code ``create_2pt5_model.py`` can be plotted using ``plot_2pt5_model.py``. Example:: @@ -129,8 +159,20 @@ The output is a set of interpolated profiles and edges that can be used to creat where ```` is the configuration file used to build the cross-sections. +2. You can construct surfaces for both the interface and slab components using ``build_complex_surface``, which takes the depth limits of each component. The profiles for the interface and slab components are then stored in two seperate files as specified in the function call:: -Classifying an earthquake catalog using the top of the slab surface [incomplete] + > sub build_complex_surface 0 50 + > sub build_complex_surface 50 450 + +3. You can plot the 3D geometry created here with:: + + > sub plot_geometries {ini_fname} False False + +where the first flag controls whether to plot the catalogue (as specified in the ini file below) and the second specifies whether or not to plot the classification (see next step). + + + +Classifying an earthquake catalog using the top of the slab surface ******************************************************************************** The ``create_2pt5_model.py`` code produces a set of profiles and edges (i.e. .csv files with the 3D coordinates) describing the geometry of the top of the slab. With this information we can separate the seismicity in an earthquake catalog into a few subsets, each one representing a specific tectonic environment (e.g. `Abrahamson and Shedlock, 1997 `__ or `Chen et al., 2017 `__ ). The procedure required to complete this task includes the following steps. @@ -140,7 +182,7 @@ The ``create_2pt5_model.py`` code produces a set of profiles and edges (i.e. .cs The configuration file specifies the geometry of surfaces, along with buffer regions, that are used as references for each tectonic environment, and the catalogue to be classified. Additionally, the configuration includes a ``priority list`` that indicates how hypocenters that can occur in overlapping buffer regions should be labeled. An example configuration file is shown below. The format of the configuration is as follows. The ``[general]`` section, which includes: - - the directory ``distance_folder`` where the Euclidean distance between each hypocenter and surface will be stored (NB: this folder must be manually created by the user) + - the directory ``distance_folder`` where the Euclidean distance between each hypocenter and surface will be stored (NB: with the first method this folder must be manually created by the user, but using the second approach it will be automatically created when making the profiles) - an .hdf5 file ``treg_filename`` that will store the results of the classfication - the .pkl file ``catalogue_filename``, which is the pickeled catalogue in HMTK format to be classified. - an array ``priority`` lists the tectonic regions, sorting the labels in the order of increasing priority, and a later label overrides classification of a hypocenter to a previous label. For example, in the configuration file shown below, an earthquake that could be classified as both ``crustal`` and ``int_prt`` will be labeled as ``int_prt``. @@ -155,12 +197,13 @@ A geometry section for each labelled tectonic environment in the ``priority`` li .. code-block:: ini [general] + root_folder = /home/kbayliss/projects/geese/classification/kur distance_folder = ./model/catalogue/classification/distances/ treg_filename = ./model/catalogue/classification/classified.hdf5 catalogue_filename = ./model/catalogue/csv/catalogue.pkl - priority=[slab_A, slab_B, crustal, int_A] + priority=[slab, crustal, int] [crustal] @@ -170,55 +213,96 @@ A geometry section for each labelled tectonic environment in the ``priority`` li crust_filename = ./model/litho1pt0/litho_crust3bottom.xyz - [int_A] + [int] - label = int_A - folder = ./model/surfaces/edges_A-int + label = int + folder = ./sfc_in lower_depth = 60. distance_buffer_above = 10. distance_buffer_below = 10. - [slab_A] + [slab] - label = slab_A - folder = ./model/surfaces/edges_A-slab + label = slab + folder = ./sfc_sl distance_buffer_above = 30. distance_buffer_below = 30. - [slab_B] - - label = slab_B - folder = ./model/surfaces/edges_B-slab - distance_buffer_above = 30. - distance_buffer_below = 30. + 2. Run the classification The classification algorithm is run using the following command:: - > cat_classify.py + > cat_classify.py Where: - ``configuration_file`` is the name of the .ini configuration file - ``distance_flag`` is a flag indicating whether or not the distances to surfaces must be computed (i.e. *True* is used the first time a classification is run for a set of surfaces and tectonic environments, but *False* when only the buffer and delta distances are changed) - - ``root_folder`` is the root directory for all paths specified in the ``configuration_file`` -3. Separate the classified events into subcatalogues +3. Check event classifications with:: -The user must decide the exact way in which they would like to separate the classified events into subcatalogues for each tectonic environment. For example, one may want to decluster the entire catalogue before separating the events, or to decluster each tectonic environment separately. View the following link for an example of the latter case: + > sub plot_geometries {ini_fname} False True -.. toctree:: - sub_tutorials/make_trts +Which will plot the 3D geometry with the events coloured by their classification. +It may be necessary to manually classify some events (e.g. where the literature supports an event being in the slab despite it being very shallow etc.). Events can be manually re-classified using:: + > ccl change_class cat.pkl classified.hdf5 eqlist.csv + +where eqlist is a csv with eventIDs to be reassigned and the required classifications. -Creating inslab sources for the OpenQuake Engine [incomplete] +4. Separate the classified events into subcatalogues + +The user must decide the exact way in which they would like to separate the classified events into subcatalogues for each tectonic environment. For example, one may want to decluster the entire catalogue before separating the events, or to decluster each tectonic environment separately. To create subcatalogues based on the classifications, the function `create_sub_catalogues` can be used:: + + > ccl create_sub_catalogues cat.pkl classified.hdf5 -p subcats_folder + + which will create subcatalogues with each label found in `classified.hdf5` which are in turn those supplied in the classification ini. These catalogues will be stored in the specified subcats_folder. +Further, these catalogues can be declustered using the `decluster_multiple_TR` function. This function takes its own toml file (decluster.toml) that specifies the different subcatalogues that should be declustered together and the declustering algorithms from within the OQ engine to use, along with the necessary parameters. For example: + +.. code-block:: toml + + [main] + + catalogue = 'cat.pkl' + tr_file = 'classified_up.hdf5' + output = 'subcats_folder/' + + create_subcatalogues = 'true' + save_aftershocks = 'true' + catalogue_add_defaults = 'true' + + [method1] + name = 'GardnerKnopoffType1' + params = {'time_distance_window' = 'UhrhammerWindow', 'fs_time_prop' = 0.1} + label = 'UH' + + [method2] + name = 'GardnerKnopoffType1' + params = {'time_distance_window' = 'GardnerKnopoffWindow', 'fs_time_prop' = 0.1} + label = 'GK' + + [case1] + regions = ['int', 'crustal'] + label = 'int_cru' + + [case2] + regions = ['slab'] + label = 'slab' + +This toml file specifies that the window declustering with parameters of Gardner and Knopoff and Uhrhammer should be used to decluster the interface and crustal events jointly, and the slab separately. This will result in four output declustered catalogues (one for each case and each method) stored in the subcats_folder. To run this declustering, we can simply use:: + + > ccl decluster_multiple_TR decluster.toml + + +Creating inslab sources for the OpenQuake Engine ************************************************************* -The construction of subduction inslab sources involves the creation of `virtual faults` elongated along the stike of the slab surface and constrained within the slab volume. +The construction of subduction inslab sources involves the creation of `virtual faults` elongated along the stike of the slab surface and constrained within the slab volume. This requires a configuration file defining some parameters for generating the ruptures 1. Create a configuration file -.. code-block:: ini +.. code-block:: slab.ini [main] @@ -236,15 +320,15 @@ The construction of subduction inslab sources involves the creation of `virtual hspa = 20. vspa = 20. - #profile_folder contains: resampled profiles and edges - profile_folder = ./model/subduction/cs_profiles/kerton/edges_zone1_slab + #profile_folder contains: resampled profiles and edges for the slab + profile_folder = ./sfc_sl # the pickled catalogue has the hmtk format - catalogue_pickle_fname = ./data/catalogues/locations/PI_cat.p + catalogue_pickle_fname = ./cat.pkl # the file with labels identifying earthquakes belonging to a given class - treg_fname = ./model/catalogue/PI_class_segments.hdf5 - label = slab_kerton1 + treg_fname = ./classified.hdf5 + label = slab # output folder out_hdf5_fname = ./tmp/ruptures/ruptures_inslab_kerton_1.hdf5 @@ -252,7 +336,7 @@ The construction of subduction inslab sources involves the creation of `virtual # output smoothing folder out_hdf5_smoothing_fname = ./tmp/smoothing/smoothing_kerton_1.hdf5 - # this is a lists + # this is a lists dips = [45, 135] # this is a dictionary @@ -270,3 +354,8 @@ The construction of subduction inslab sources involves the creation of `virtual mmin = 6.5 mmax = 7.80 +kwargs = {'out_hdf5_fname': out_hdf5_fname, 'out_hdf5_smoothing_fname': out_hdf5_smoothing_fname} + +calculate_ruptures(ini_fname, **kwargs) + + From 9e608a0f67cbd272ee305ebc8b0f33604a5bc7d7 Mon Sep 17 00:00:00 2001 From: Kirsty Bayliss Date: Tue, 10 Sep 2024 16:25:09 +0200 Subject: [PATCH 16/23] further improvements to sub documentation --- docsrc/contents/sub.rst | 72 +++++++++++++++++++++++++++++++++++------ 1 file changed, 63 insertions(+), 9 deletions(-) diff --git a/docsrc/contents/sub.rst b/docsrc/contents/sub.rst index daab5c466..0f9b806c1 100644 --- a/docsrc/contents/sub.rst +++ b/docsrc/contents/sub.rst @@ -298,7 +298,7 @@ This toml file specifies that the window declustering with parameters of Gardner Creating inslab sources for the OpenQuake Engine ************************************************************* -The construction of subduction inslab sources involves the creation of `virtual faults` elongated along the stike of the slab surface and constrained within the slab volume. This requires a configuration file defining some parameters for generating the ruptures +The construction of subduction inslab sources involves the creation of `virtual faults` elongated along the stike of the slab surface and constrained within the slab volume. This requires a configuration file defining some parameters for generating the ruptures. 1. Create a configuration file @@ -310,6 +310,12 @@ The construction of subduction inslab sources involves the creation of `virtual profile_sd_topsl = 40. edge_sd_topsl = 40. + + # MFD + agr = 5.945 + bgr = 1.057 + mmin = 6.5 + mmax = 7.80 sampling = 10. @@ -348,14 +354,62 @@ The construction of subduction inslab sources involves the creation of `virtual # magnitude scaling relationship mag_scaling_relation = StrasserIntraslab - # MFD - agr = 5.945 - bgr = 1.057 - mmin = 6.5 - mmax = 7.80 + The MFD parameters should be set by the modeller and determined from some combination of the seismicity and tectonics. The `mmin` parameter defines the lower magnitude limit at which to generate slab ruptures. A lower `mmin` will result in many more smaller ruptures which will increase the size of the rupture object, so this parameter should be chosen carefully considering the size of ruptures at slab locations that might be relevant for hazard. + +The `sampling` parameter determines the spatial sampling to be used when simulating ruptures. The `float_strike` and `float_dip` parameters specify the strike and dip for floating ruptures, while the list of `dips` instead specifies dip angles used when creating virtual faults inside the rupture. -kwargs = {'out_hdf5_fname': out_hdf5_fname, 'out_hdf5_smoothing_fname': out_hdf5_smoothing_fname} - -calculate_ruptures(ini_fname, **kwargs) +The `uniform_fraction` determines the percentage of the ruptures to be uniformly distributed across the slab. A higher uniform fraction means that the distribution of ruptures will be randomly uniform, and a lower uniform fraction means that a larger percentage of ruptures will instead be distributed according to the smoothed distribution of seismicity in the slab. + +Ruptures can be created using:: + + > calculate_ruptures ini_fname out_hdf5_fname out_hdf5_smoothing_fname + +which will create two hdf5 files with the different ruptures. We then need to process these into xml files, which we can do using:: + + > create_inslab_nrml(model, out_hdf5_fname, out_path, investigation_t = 1) +You can plot the 3D model of the subduction zone with ruptures as below: + +.. code-block:plot_ruptures +xml_fname = "/home/kbayliss/projects/force/code/Pacific_Jan24/mariana_April24/xml/mar_newrup_unif/6.55.xml" +ssm = to_python(xml_fname) + +hy_lo, hy_la, hy_d = [],[],[] +mags = [] +rates = [] +for sg in ssm: + for src in sg: + for rup in src.data: + h = rup[0].hypocenter + hy_lo.append(h.longitude); hy_la.append(h.latitude); hy_d.append(-h.depth) + mags.append(rup[0].mag) + prob1 = rup[1].data[1][0] + rate = -1*np.log(1-prob1) + rates.append(rate) + +df = pd.DataFrame({'lons': hy_lo, 'lats': hy_la, 'depth': hy_d, 'mag': mags, 'rate': rates}) + +fig = plt.figure() +ax = fig.add_subplot(111, projection='3d') +ax.scatter(df.lons, df.lats, df.depth, c=df.rate) +plt.show() + +This process will create xml files for each magnitude bin, which is a little impractical. To combine these, we can do the following:: + +path1 = glob.glob('./xml/mar_newrup/*.xml') +data_all = [] +for p in path1: + ssm = to_python(p) + for sg in ssm: + for src in sg: + print(len(src.data)) + data_all.extend(src.data) +src.data = data_all +src.name = 'ruptures for src_slab' +src.source_id = 'src_slab' + +write_source_model('./xml/mar_newrup/slab.xml', [src], investigation_time=1.0) + +Which creates a single `slab.xml` file containing all the ruptures across all magnitude bins. This file can then be used directly in OQ as a non-parametric rupture source. +See the process in more detail in this notebook From d0449986e378b17f462e8eeeb2556078cfd7b925 Mon Sep 17 00:00:00 2001 From: Kirsty Bayliss Date: Tue, 17 Sep 2024 14:36:08 +0200 Subject: [PATCH 17/23] binwidths not implemented properly --- openquake/cat/completeness/analysis.py | 10 +++++----- openquake/cat/completeness/norms.py | 17 ++++++++++++----- 2 files changed, 17 insertions(+), 10 deletions(-) diff --git a/openquake/cat/completeness/analysis.py b/openquake/cat/completeness/analysis.py index c88e51c10..b0b9a6238 100644 --- a/openquake/cat/completeness/analysis.py +++ b/openquake/cat/completeness/analysis.py @@ -176,8 +176,8 @@ def check_criterion(criterion, rate, previous_norm, tvars): norm = abs(tmp_rate - rate_ma) elif criterion == 'optimize': - tmp_rate = -1 - norm = get_norm_optimize(tcat, aval, bval, ctab, cmag, n_obs, t_per, + tmp_rate = -1 + norm = get_norm_optimize(tcat, aval, bval, ctab, binw, cmag, n_obs, t_per, last_year, info=False) elif criterion == 'optimize_a': @@ -190,7 +190,7 @@ def check_criterion(criterion, rate, previous_norm, tvars): mmin=ref_mag, mmax=ref_upp_mag) elif criterion == 'optimize_c': tmp_rate = -1 - norm = get_norm_optimize_c(tcat, aval, bval, ctab, last_year) + norm = get_norm_optimize_c(tcat, aval, bval, ctab, last_year, ref_mag, ref_upp_mag, binw) elif criterion == 'gft': tmp_rate = -1 @@ -201,8 +201,8 @@ def check_criterion(criterion, rate, previous_norm, tvars): norm = get_norm_optimize_weichert(tcat, aval, bval, ctab, last_year) elif criterion == 'poisson': - tmp_rate = -1 - norm = get_norm_optimize_c(tcat, aval, bval, ctab, last_year, ref_mag) + tmp_rate = -1 + norm = get_norm_optimize_c(tcat, aval, bval, ctab, last_year, ref_upp_mag, binw) if norm is None or np.isnan(norm): return False, -1, previous_norm diff --git a/openquake/cat/completeness/norms.py b/openquake/cat/completeness/norms.py index 0f039e529..14ee8ba97 100644 --- a/openquake/cat/completeness/norms.py +++ b/openquake/cat/completeness/norms.py @@ -131,11 +131,18 @@ def get_norm_optimize(tcat, aval, bval, ctab, binw, cmag, n_obs, t_per, last_yea calculated norm for input completeness. Smaller norm is better """ - occ = np.zeros((ctab.shape[0])) - dur = np.zeros((ctab.shape[0])) - #mags = np.array(list(ctab[:, 1])+[10]) - mags = np.array(np.arange(ctab[:, 1], 10, binw)) - + #occ = np.zeros((ctab.shape[0])) + #dur = np.zeros((ctab.shape[0])) + mags = np.array(list(ctab[:, 1])+[10]) + print("old mags: ", mags) + mags = np.array(np.arange(ctab[:, 1][0], 10, binw)) + print("new mags: ", mags) + occ = np.zeros(len(mags)) + dur = np.zeros(len(mags)) + ## Problem: we *should* be testing all the bins > completeness + ## This is currently binning based on the completeness windows, + ## which feels bad. + for i, mag in enumerate(mags[:-1]): idx = (cmag >= mags[i]) & (cmag < mags[i+1]) if np.any(idx): From 12d9a55b3d1657e3b7677e8ee7375ab8d94cfde2 Mon Sep 17 00:00:00 2001 From: Kirsty Bayliss Date: Mon, 30 Sep 2024 10:56:07 +0200 Subject: [PATCH 18/23] adding missing norm import --- openquake/cat/completeness/analysis.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openquake/cat/completeness/analysis.py b/openquake/cat/completeness/analysis.py index b0b9a6238..2f15587fe 100644 --- a/openquake/cat/completeness/analysis.py +++ b/openquake/cat/completeness/analysis.py @@ -266,7 +266,7 @@ def _completeness_analysis(fname, years, mags, binw, ref_mag, mag_low, ref_upp_m # Checking input if criterion not in ['match_rate', 'largest_rate', 'optimize', 'weichert', - 'poisson', 'optimize_a', 'optimize_b', 'optimize_c','optimize_d']: + 'poisson', 'optimize_a', 'optimize_b', 'optimize_c','optimize_d', 'gft']: raise ValueError('Unknown optimization criterion') tcat = _load_catalogue(fname) From c1edb4fa5afd1f1748994f85a32dc838d1e8bdfb Mon Sep 17 00:00:00 2001 From: Kirsty Bayliss Date: Mon, 30 Sep 2024 10:56:41 +0200 Subject: [PATCH 19/23] adding norm tests --- .../cat/tests/completeness/norms_test.py | 26 ++++++++++++++++--- 1 file changed, 22 insertions(+), 4 deletions(-) diff --git a/openquake/cat/tests/completeness/norms_test.py b/openquake/cat/tests/completeness/norms_test.py index 46472055e..59173e59a 100644 --- a/openquake/cat/tests/completeness/norms_test.py +++ b/openquake/cat/tests/completeness/norms_test.py @@ -31,7 +31,7 @@ from openquake.hmtk.seismicity.occurrence.utils import get_completeness_counts from openquake.mbt.tools.model_building.dclustering import _add_defaults from openquake.cat.completeness.norms import ( - get_norm_optimize_b, get_norm_optimize_c, get_norm_optimize_poisson) + get_norm_optimize_b, get_norm_optimize_c, get_norm_optimize, get_norm_optimize_poisson) DATA = os.path.join(os.path.dirname(__file__), 'data') @@ -60,8 +60,10 @@ def test_case01(self): bval = 1.0 cmag, t_per, n_obs = get_completeness_counts(self.cat, self.compl, mbinw) - _ = get_norm_optimize_b(aval, bval, self.compl, self.cat, mbinw, ybinw) - + norm = get_norm_optimize_b(aval, bval, self.compl, self.cat, mbinw, ybinw) + print(f'{norm:.5e}') + self.assertAlmostEqual(norm,8.60607e-01, msg='rmag_rate', places=4) + def test_case02(self): mbinw = 0.1 @@ -79,8 +81,22 @@ def test_case02(self): cmag, t_per, n_obs = get_completeness_counts(cat, compl, mbinw) norm = get_norm_optimize_c(cat, aval, bval, compl, 2022, ref_mag=4.4) + + print(f'{norm:.5e}') + self.assertAlmostEqual(norm,5.60922e-01, msg='rmag_rate', places=4) + + def test_optimize(self): + mbinw = 0.5 + ybinw = 10.0 + aval = 2.0 + bval = 1.0 + binw = 0.1 + last_year = 2020 + cmag, t_per, n_obs = get_completeness_counts(self.cat, self.compl, + mbinw) + + norm = get_norm_optimize(self.cat, aval, bval, self.compl,binw, cmag, n_obs, t_per, last_year) print(f'{norm:.5e}') - def test_poisson(self): @@ -101,3 +117,5 @@ def test_poisson(self): cmag, t_per, n_obs = get_completeness_counts(cat, compl, mbinw) norm = get_norm_optimize_poisson(cat, aval, bval, compl, 2022) print(f'{norm:.5e}') + + self.assertAlmostEqual(norm,-16.1132, msg='rmag_rate', places=4) From 8336ccea6045e59feefe28188acbabd210e41c39 Mon Sep 17 00:00:00 2001 From: Kirsty Bayliss Date: Mon, 30 Sep 2024 10:57:00 +0200 Subject: [PATCH 20/23] adding warnings for small datasets --- openquake/wkf/compute_gr_params.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/openquake/wkf/compute_gr_params.py b/openquake/wkf/compute_gr_params.py index f6d2ee174..9ff3e0c7a 100644 --- a/openquake/wkf/compute_gr_params.py +++ b/openquake/wkf/compute_gr_params.py @@ -439,6 +439,10 @@ def _weichert_analysis(tcat, ctab, binw, cmag, n_obs, t_per): weichert_config = {'magnitude_interval': binw, 'reference_magnitude': numpy.min(ctab[:, 1])} weichert = Weichert() + + nev = len(tcat.data['magnitude']) + if nev < 10: + print("Few events in this catalogue (only ", nev, " events above completeness)") # weichert.calculate returns bGR and its standard deviation + log10(rate) # for the reference magnitude and its standard deviation. In this case @@ -447,6 +451,9 @@ def _weichert_analysis(tcat, ctab, binw, cmag, n_obs, t_per): # bval, sigmab, aval, sigmaa = fun(tcat, weichert_config, ctab) bval, sigmab, rmag_rate, rmag_rate_sigma, aval, sigmaa = fun( tcat, weichert_config, ctab) + + if bval < 0.5 or bval > 2: + print("suspicious b-value, recheck your catalogue (b = ", bval, ")") # Computing confidence intervals gwci = get_weichert_confidence_intervals From 541a23e3db2743dcd1b9b65eaece6d443f8fd715 Mon Sep 17 00:00:00 2001 From: Kirsty Bayliss Date: Thu, 31 Oct 2024 16:13:17 +0100 Subject: [PATCH 21/23] require sourceConverter object --- openquake/wkf/distributed_seismicity.py | 32 ++++++++++++++++--------- 1 file changed, 21 insertions(+), 11 deletions(-) diff --git a/openquake/wkf/distributed_seismicity.py b/openquake/wkf/distributed_seismicity.py index 61bd6310e..62483a60e 100644 --- a/openquake/wkf/distributed_seismicity.py +++ b/openquake/wkf/distributed_seismicity.py @@ -63,6 +63,14 @@ PLOTTING = False +DEFAULT = SourceConverter( + investigation_time=1.0, + rupture_mesh_spacing=2.0, + complex_fault_mesh_spacing=5.0, + width_of_mfd_bin=binw, + area_source_discretization=5.0, + ) + def get_mfd_moment(mfd): return np.sum( @@ -205,7 +213,7 @@ def explode(srcs, check_moment_rates=True): if check_moment_rates: src_moment = get_mfd_moment(src.mfd) nsrc_moment = get_mfd_moment(nsrc.mfd) - np.testing.assert_almost_equal(src_moment * wei, nsrc_moment) + #np.testing.assert_almost_equal(src_moment * wei, nsrc_moment) exploded_srcs.append(nsrc) @@ -219,6 +227,7 @@ def remove_buffer_around_faults( dst: float, threshold_mag: float = 6.5, use: str = '', + source_conv: SourceConverter ): """ Remove the seismicity above a magnitude threshold for all the point @@ -231,11 +240,18 @@ def remove_buffer_around_faults( `./../m01_asc/oq/zones/src_*.xml` :param out_path: The path where to write the output .xml file - :param dst: - The distance in km of the buffer :param dst: The threshold distance used to separate seismicity on the fault and in the distributed seismicity sources + :param threshold_mag: + the magnitude to cut at faults (ie if threshold_mag = 6 + events M<6 are kept around faults) + :param use: + optional list of sources to apply this to + :param source_conv: + An instance of the class + :class:`openquake.hazardlib.sourceconverter.SourceConverter` + :returns: A .xml file with the ajusted point sources """ @@ -246,13 +262,7 @@ def remove_buffer_around_faults( # Create a source converter binw = 0.1 - sourceconv = SourceConverter( - investigation_time=1.0, - rupture_mesh_spacing=5.0, - complex_fault_mesh_spacing=5.0, - width_of_mfd_bin=binw, - area_source_discretization=5.0, - ) + # Get the surfaces representing the faults faults = _get_fault_surfaces(fname, sourceconv) @@ -395,7 +405,7 @@ def _get_fault_surfaces(fname: str, sourceconv: SourceConverter) -> list: # Read file the fault sources ssm_faults = to_python(fname, sourceconv) - + # Check content of the seismic source model. We want only one group. msg = 'The seismic source model for fault contains more than one group' assert len(ssm_faults) == 1 From 63f1ad9fedc1748c6966a75b9aaca2e0115f6f38 Mon Sep 17 00:00:00 2001 From: Kirsty Bayliss Date: Wed, 13 Nov 2024 17:27:37 +0100 Subject: [PATCH 22/23] properly update source names when removing fault buffers --- openquake/mbi/wkf/remove_buffer_around_faults.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/openquake/mbi/wkf/remove_buffer_around_faults.py b/openquake/mbi/wkf/remove_buffer_around_faults.py index 7bfd86383..c18942b0e 100755 --- a/openquake/mbi/wkf/remove_buffer_around_faults.py +++ b/openquake/mbi/wkf/remove_buffer_around_faults.py @@ -8,14 +8,14 @@ def main(fname: str, path_point_sources: str, out_path: str, dst: float, - threshold_mag: float=6.5, use: str=''): + threshold_mag: float=6.5): # Create the output folder (if needed) create_folder(out_path) # Process sources remove_buffer_around_faults(fname, path_point_sources, out_path, dst, - threshold_mag, use) + threshold_mag) main.fname = "Pattern for input .xml file with fault sources" @@ -24,6 +24,7 @@ def main(fname: str, path_point_sources: str, out_path: str, dst: float, main.dst = "Distance [km] of the buffer around the fault" main.threshold_mag = "Threshold magnitude" main.use = 'A list with the ID of sources that should be considered' +main.source_conv = 'SourceConverter object' if __name__ == '__main__': sap.run(main) From ba82dde640ae4faf661414c5a164d456bc8aafd8 Mon Sep 17 00:00:00 2001 From: Kirsty Bayliss Date: Fri, 15 Nov 2024 09:29:14 +0100 Subject: [PATCH 23/23] add b_sigma and a_sigma to config --- openquake/wkf/compute_gr_params.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/openquake/wkf/compute_gr_params.py b/openquake/wkf/compute_gr_params.py index 9ff3e0c7a..065efd0b2 100644 --- a/openquake/wkf/compute_gr_params.py +++ b/openquake/wkf/compute_gr_params.py @@ -461,7 +461,7 @@ def _weichert_analysis(tcat, ctab, binw, cmag, n_obs, t_per): rmag = weichert_config['reference_magnitude'] return (aval, bval, lcl, ucl, exrates, exrates_scaled, rmag, rmag_rate, - rmag_rate_sigma) + rmag_rate_sigma, sigmab, sigmaa) def _get_gr_double_trunc_exceedance_rates(agr, bgr, cmag, binw, mmax): @@ -693,7 +693,7 @@ def weichert_analysis(fname_input_pattern, fname_config, folder_out=None, # Compute aGR and bGR using Weichert out = _weichert_analysis(tcat, ctab, binw, cent_mag, n_obs, t_per) - aval, bval, lcl, ucl, ex_rat, ex_rts_scl, rmag, rm_rate, rm_sig = out + aval, bval, lcl, ucl, ex_rat, ex_rts_scl, rmag, rm_rate, rm_sig, sigmab, sigmaa = out # Plot _weichert_plot(cent_mag, n_obs, binw, t_per, ex_rts_scl, @@ -714,6 +714,10 @@ def weichert_analysis(fname_input_pattern, fname_config, folder_out=None, model['sources'][src_id]['rmag_rate'] = float(tmp) tmp = f"{rm_sig:.5e}" model['sources'][src_id]['rmag_rate_sig'] = float(tmp) + tmp = f"{sigmab:.5e}" + model['sources'][src_id]['bgr_sig'] = float(tmp) + tmp = f"{sigmaa:.5e}" + model['sources'][src_id]['agr_sig'] = float(tmp) # Save figures if folder_out_figs is not None: