Source code for openquake.hmtk.seismicity.completeness.comp_stepp_1971
# -*- coding: utf-8 -*-# vim: tabstop=4 shiftwidth=4 softtabstop=4## LICENSE## Copyright (C) 2010-2025 GEM Foundation, G. Weatherill, M. Pagani,# D. Monelli.## The Hazard Modeller's Toolkit is free software: you can redistribute# it and/or modify it under the terms of the GNU Affero General Public# License as published by the Free Software Foundation, either version# 3 of the License, or (at your option) any later version.## You should have received a copy of the GNU Affero General Public License# along with OpenQuake. If not, see <http://www.gnu.org/licenses/>## DISCLAIMER## The software Hazard Modeller's Toolkit (openquake.hmtk) provided herein# is released as a prototype implementation on behalf of# scientists and engineers working within the GEM Foundation (Global# Earthquake Model).## It is distributed for the purpose of open collaboration and in the# hope that it will be useful to the scientific, engineering, disaster# risk and software design communities.## The software is NOT distributed as part of GEM’s OpenQuake suite# (https://www.globalquakemodel.org/tools-products) and must be considered as a# separate entity. The software provided herein is designed and implemented# by scientific staff. It is not developed to the design standards, nor# subject to same level of critical review by professional software# developers, as GEM’s OpenQuake software suite.## Feedback and contribution to the software is welcome, and can be# directed to the hazard scientific staff of the GEM Model Facility# (hazard@globalquakemodel.org).## The Hazard Modeller's Toolkit (openquake.hmtk) is therefore distributed# WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License# for more details.## The GEM Foundation, and the authors of the software, assume no# liability for use of the software."""Module :mod:`openquake.hmtk.seismicity.completeness.comp_stepp_1972` definesthe openquake.hmtk implementation of the Stepp (1972) algorithm for analysingthe completeness of an earthquake catalogue"""importnumpyasnpfromscipy.optimizeimportfmin_l_bfgs_bfromopenquake.hmtk.seismicity.utilsimport(decimal_time,piecewise_linear_scalar,)fromopenquake.hmtk.seismicity.completeness.baseimport(BaseCatalogueCompleteness,COMPLETENESS_METHODS,)
[docs]defget_bilinear_residuals_stepp(input_params,xvals,yvals,slope1_fit):""" Returns the residual sum-of-squares value of a bilinear fit to a data set - with a segment - 1 gradient fixed by an input value (slope_1_fit) :param list input_params: Input parameters for the bilinear model [slope2, crossover_point, intercept] :param numpy.ndarray xvals: x-values of the data to be fit :param numpy.ndarray yvals: y-values of the data to be fit :param float slope1_fit: Gradient of the first slope :returns: Residual sum-of-squares of fit """params=np.hstack([slope1_fit,input_params])num_x=len(xvals)y_model=np.zeros(num_x,dtype=float)residuals=np.zeros(num_x,dtype=float)forilocinrange(0,num_x):y_model[iloc]=piecewise_linear_scalar(params,xvals[iloc])residuals[iloc]=(yvals[iloc]-y_model[iloc])**2.0returnnp.sum(residuals)
[docs]@COMPLETENESS_METHODS.add("completeness",magnitude_bin=float,time_bin=float,increment_lock=bool)classStepp1971(BaseCatalogueCompleteness):""" Implements the completeness analysis methodology of Stepp (1972) Stepp, J. C. (1972) Analysis of Completeness of the Earhquake Sample in the Puget Sound Area and Its Effect on Statistical Estimates of Earthquake Hazard, NOAA Environmental Research Laboratories. The original methodology of J. C. Stepp (1972) implements a graphical method in which the deviation of the observed rate from the expected Poisson rate is assessed by judgement. To implement the selection in an automated fashion this implementation uses optimisation of a 2-segment piecewise linear fit to each magnitude bin, using the segment intersection point to identify the completeness period. Adaptation implemented by Weatherill, G. A., GEM Model Facility, Pavia :attribute numpy.ndarray magnitude_bin: Edges of the magnitude bins :attribute numpy.ndarray sigma: Sigma lambda defined by Equation 4 in Stepp (1972) :attribute numpy.ndarray time_values: Duration values :attribute numpy.ndarray model_line: Expected Poisson rate for each magnitude bin :attribute numpy.ndarray completeness_table: Resulting completeness table """def__init__(self):BaseCatalogueCompleteness.__init__(self)self.magnitude_bin=Noneself.time_values=Noneself.sigma=Noneself.model_line=Noneself.completeness_table=Noneself.end_year=None
[docs]defcompleteness(self,catalogue,config):""" Gets the completeness table. :param catalogue: Earthquake catalogue as instance of :class:`openquake.hmtk.seismicity.catalogue.Catalogue` :param dict config: Configuration parameters of the algorithm, containing the following information: 'magnitude_bin' Size of magnitude bin (non-negative float) 'time_bin' Size (in dec. years) of the time window (non-negative float) 'increment_lock' Boolean to indicate whether to ensure completeness magnitudes always decrease with more recent bins :returns: 2-column table indicating year of completeness and corresponding magnitude numpy.ndarray """# If mag_bin is an array then bins are input, otherwise a single# parameter is inputdyear=decimal_time(catalogue.data["year"],catalogue.data["month"],catalogue.data["day"],catalogue.data["hour"],catalogue.data["minute"],catalogue.data["second"],)mag=catalogue.data["magnitude"]# Get magnitude binsself.magnitude_bin=self._get_magnitudes_from_spacing(catalogue.data["magnitude"],config["magnitude_bin"])# Get time bins_s_year,time_bin=self._get_time_limits_from_config(config,dyear)# Count magnitudes(self.sigma,_counter,n_mags,n_times,self.time_values,)=self._count_magnitudes(mag,dyear,time_bin)# Get completeness magnitudescomp_time,_gradient_2,self.model_line=self.get_completeness_points(self.time_values,self.sigma,n_mags,n_times)# If the increment lock is selected then ensure completeness time# does not decreaseifconfig["increment_lock"]:forilocinrange(0,len(comp_time)):cond=(iloc>0and(comp_time[iloc]<comp_time[iloc-1]))ornp.isnan(comp_time[iloc])ifcond:comp_time[iloc]=comp_time[iloc-1]self.completeness_table=np.column_stack([np.floor(self.end_year-comp_time),self.magnitude_bin[:-1]])returnself.completeness_table
[docs]defsimplify(self,deduplicate=True,mag_range=None,year_range=None):""" Simplify a completeness table result. Intended to work with 'increment_lock' enabled. """ifself.completeness_tableisNone:returnyears=self.completeness_table[:,0]mags=self.completeness_table[:,1]keep=np.array([True]*years.shape[0])ifdeduplicate:keep[1:]=years[1:]!=years[:-1]ifyear_rangeisnotNone:year_min,year_max=year_rangeifyear_minisnotNone:too_early=years<year_minkeep&=years>=years[too_early].max()self.completeness_table[too_early,0]=year_minifyear_maxisnotNone:keep&=years<=year_maxifmag_rangeisnotNone:mag_min,mag_max=mag_rangeifmag_minisnotNone:keep&=mags>=mag_minifmag_maxisnotNone:keep&=mags<=mag_maxself.completeness_table=self.completeness_table[keep,:]self.model_line=self.model_line[:,keep]self.sigma=self.sigma[:,keep]self.magnitude_bin=self.magnitude_bin[np.hstack((keep,True))]
def_get_time_limits_from_config(self,config,dec_year):""" Defines the time bins for consideration based on the config time_bin settings - also sets self.end_year (int) the latest year in catalogue :param dict config: Configuration for the Stepp (1971) algorithm :param numpy.ndarray dec_year: Time of the earthquake in decimal years :returns: * start_year: Earliest year found in the catalogue * time_bin: Bin edges of the time windows """cond=isinstance(config["time_bin"],list)orisinstance(config["time_bin"],np.ndarray)ifcond:# Check to make sure input years are in order from recent to oldestforivalinrange(1,len(config["time_bin"])):diff=config["time_bin"][ival]-config["time_bin"][ival-1]ifdiff>0.0:raiseValueError("Configuration time windows must be ""ordered from recent to oldest!")self.end_year=config["time_bin"][0]start_year=config["time_bin"][-1]time_bin=np.array(config["time_bin"])else:self.end_year=np.floor(np.max(dec_year))start_year=np.floor(np.min(dec_year))if(self.end_year-start_year)<config["time_bin"]:raiseValueError("Catalogue duration smaller than time bin"" width - change time window size!")time_bin=np.arange(self.end_year-config["time_bin"],start_year-config["time_bin"],-config["time_bin"],)returnstart_year,time_bindef_get_magnitudes_from_spacing(self,magnitudes,delta_m):"""If a single magnitude spacing is input then create the bins :param numpy.ndarray magnitudes: Vector of earthquake magnitudes :param float delta_m: Magnitude bin width :returns: Vector of magnitude bin edges (numpy.ndarray) """min_mag=np.min(magnitudes)max_mag=np.max(magnitudes)if(max_mag-min_mag)<delta_m:raiseValueError("Bin width greater than magnitude range!")mag_bins=np.arange(np.floor(min_mag),np.ceil(max_mag+delta_m),delta_m)# Check to see if there are magnitudes in lower and upper binsis_mag=np.logical_and(mag_bins-max_mag<delta_m,min_mag-mag_bins<delta_m)mag_bins=mag_bins[is_mag]returnmag_binsdef_count_magnitudes(self,mags,times,time_bin):""" For each completeness magnitude-year counts the number of events inside each magnitude bin. :param numpy.ndarray mags: Magnitude of earthquakes :param numpy.ndarray times: Vector of decimal event times :param numpy.ndarray time_bin: Vector of bin edges of the time windows :returns: * sigma - Poisson variance (numpy.ndarray) * counter - Number of earthquakes in each magnitude-time bin * n_mags - number of magnitude bins (Integer) * n_times - number of time bins (Integer) * n_years - effective duration of each time window (numpy.ndarray) """n_mags=len(self.magnitude_bin)-1n_times=len(time_bin)counter=np.zeros([n_times,n_mags],dtype=int)# Count all the magnitudes later than or equal to the reference timeforilocinrange(0,n_times):id0=times>time_bin[iloc]counter[iloc,:]=np.histogram(mags[id0],self.magnitude_bin)[0]# Get sigma_lambda = sqrt(n / nyears) / sqrt(n_years)sigma=np.zeros([n_times,n_mags],dtype=float)n_years=np.floor(np.max(times))-time_binforilocinrange(0,n_mags):id0=counter[:,iloc]>0ifany(id0):nvals=counter[id0,iloc].astype(float)sigma[id0,iloc]=np.sqrt((nvals/n_years[id0]))/np.sqrt(n_years[id0])returnsigma,counter,n_mags,n_times,n_years
[docs]defget_completeness_points(self,n_years,sigma,n_mags,n_time):"""Fits a bilinear model to each sigma-n_years combination in order to get the crossover point. The gradient of the first line must always be 1 / sqrt(T), but it is free for the other lines :param numpy.ndarray n_years: Duration of each completeness time window :param numpy.ndarray sigma: Poisson variances of each time-magnitude combination :param int n_mags: Number of magnitude bins :param int n_time: Number of time bins :returns: * comp_time (Completeness duration) * gradient_2 (Gradient of second slope of piecewise linear fit) * model_line (Expected Poisson rate for data (only used for plot) """time_vals=np.log10(n_years)sigma_vals=np.copy(sigma)valid_mapper=np.ones([n_time,n_mags],dtype=bool)valid_mapper[sigma_vals<1e-9]=Falsecomp_time=np.zeros(n_mags,dtype=float)gradient_2=np.zeros(n_mags,dtype=float)model_line=np.zeros([n_time,n_mags],dtype=float)forilocinrange(0,n_mags):id0=valid_mapper[:,iloc]ifnp.sum(id0)<3:# Too few events for fitting a bilinear model!comp_time[iloc]=np.nangradient_2[iloc]=np.nanmodel_line[:,iloc]=np.nanelse:(comp_time[iloc],gradient_2[iloc],model_line[id0,iloc],)=self._fit_bilinear_to_stepp(time_vals[id0],np.log10(sigma[id0,iloc]))returncomp_time,gradient_2,model_line
def_fit_bilinear_to_stepp(self,xdata,ydata,initial_values=None):""" Returns the residuals of a bilinear fit subject to the following constraints: 1) gradient of slope 1 = 1 / sqrt(T) 2) Crossover (x_c) < 0 3) gradient 2 is always < 0 :param numpy.ndarray xdata: x-value of the data set :param numpy.ndarray ydata: y-value of the data set :param list initial_values: For unit-testing allows the possibility to specify the initial values of the algorithm [slope_2, cross_over, intercept] :returns: * completeness_time: The duration of completeness of the bin * Gradient of the second slope * model_line: Expected Poisson model """fixed_slope=-0.5# f'(log10(T^-0.5)) === 0.5ifisinstance(initial_values,list)orisinstance(initial_values,np.ndarray):x_0=initial_valueselse:x_0=[-1.0,xdata[int(len(xdata)/2)],xdata[0]]bnds=((None,fixed_slope),(0.0,None),(None,None))result,_,convergence_info=fmin_l_bfgs_b(get_bilinear_residuals_stepp,x_0,args=(xdata,ydata,fixed_slope),approx_grad=True,bounds=bnds,disp=0,)ifconvergence_info["warnflag"]!=0:# Optimisation has failed to converge - print the reason whyprint(convergence_info["task"])returnnp.nan,np.nan,np.nan*np.ones(len(xdata))# Result contains three parameters = m_2, x_c, c_0# x_c is the crossover point (i.e. the completeness_time)# m_2 is the gradient of the latter slope# c_0 is the intercept - which helps locate the line at the datamodel_line=10.0**(fixed_slope*xdata+result[2])completeness_time=10.0**result[1]returncompleteness_time,result[0],model_line