Source code for openquake.hmtk.seismicity.declusterer.dec_gardner_knopoff
# -*- 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.declusterer.dec_gardner_knopoff`defines the Gardner and Knopoff declustering algorithm"""importnumpyasnpfromopenquake.hmtk.seismicity.declusterer.baseimport(BaseCatalogueDecluster,DECLUSTERER_METHODS,)fromopenquake.hmtk.seismicity.utilsimportdecimal_year,haversinefromopenquake.hmtk.seismicity.declusterer.distance_time_windowsimport(TIME_DISTANCE_WINDOW_FUNCTIONS,)
[docs]@DECLUSTERER_METHODS.add("decluster",time_distance_window=TIME_DISTANCE_WINDOW_FUNCTIONS,fs_time_prop=float,)classGardnerKnopoffType1(BaseCatalogueDecluster):""" This class implements the Gardner Knopoff algorithm as described in this paper: Gardner, J. K. and Knopoff, L. (1974). Is the sequence of aftershocks in Southern California, with aftershocks removed, poissonian?. Bull. Seism. Soc. Am., 64(5): 1363-1367. """
[docs]defdecluster(self,catalogue,config):""" The configuration of this declustering algorithm requires two objects: - A time-distance window object (key is 'time_distance_window') - A value in the interval [0,1] expressing the fraction of the time window used for aftershocks (key is 'fs_time_prop') :param catalogue: Catalogue of earthquakes :type catalogue: Dictionary :param config: Configuration parameters :type config: Dictionary :returns: **vcl vector** indicating cluster number, **flagvector** indicating which eq events belong to a cluster :rtype: numpy.ndarray """# Get relevant parametersneq=len(catalogue.data["magnitude"])# Number of earthquakes# Get decimal year (needed for time windows)year_dec=decimal_year(catalogue.data["year"],catalogue.data["month"],catalogue.data["day"],)# Get space and time windows corresponding to each event# Initial Position Identifiersw_space,sw_time=config["time_distance_window"].calc(catalogue.data["magnitude"],config.get("time_cutoff"))eqid=np.arange(0,neq,1)# Pre-allocate cluster index vectorsvcl=np.zeros(neq,dtype=int)# Sort magnitudes into descending orderid0=np.flipud(np.argsort(catalogue.data["magnitude"],kind="heapsort"))longitude=catalogue.data["longitude"][id0]latitude=catalogue.data["latitude"][id0]sw_space=sw_space[id0]sw_time=sw_time[id0]year_dec=year_dec[id0]eqid=eqid[id0]flagvector=np.zeros(neq,dtype=int)# Begin cluster identificationclust_index=0foriinrange(0,neq-1):ifvcl[i]==0:# Find Events inside both fore- and aftershock time windowsdt=year_dec-year_dec[i]vsel=np.logical_and(vcl==0,np.logical_and(dt>=(-sw_time[i]*config["fs_time_prop"]),dt<=sw_time[i],),)# Of those events inside time window,# find those inside distance windowvsel1=(haversine(longitude[vsel],latitude[vsel],longitude[i],latitude[i],)<=sw_space[i])vsel[vsel]=vsel1[:,0]temp_vsel=np.copy(vsel)temp_vsel[i]=Falseifany(temp_vsel):# Allocate a cluster numbervcl[vsel]=clust_index+1flagvector[vsel]=1# For those events in the cluster before the main event,# flagvector is equal to -1temp_vsel[dt>=0.0]=Falseflagvector[temp_vsel]=-1flagvector[i]=0clust_index+=1# Re-sort the catalog_matrix into original orderid1=np.argsort(eqid,kind="heapsort")eqid=eqid[id1]vcl=vcl[id1]flagvector=flagvector[id1]returnvcl,flagvector