Source code for openquake.hmtk.plotting.seismicity.completeness.cumulative_rate_analysis
# -*- coding: utf-8 -*-# vim: tabstop=4 shiftwidth=4 softtabstop=4## LICENSE## Copyright (C) 2010-2023 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 openquake.hmtk.plotting.seismicity.completeness.simple_completeness is a graphicalfunction for estimating the completeness period of magnitude intervalsby plotting the cumulative rate of events with time in each interval"""importnumpyasnpimportpylabimportmatplotlib.pyplotaspltfromopenquake.hmtk.seismicity.completeness.baseimport(BaseCatalogueCompleteness,)
[docs]classSimpleCumulativeRate(BaseCatalogueCompleteness):""" Class to define the temporal variation in completess using simple changes in cumulative rates in individual completeness bins """
[docs]defcompleteness(self,catalogue,config,saveplot=False,filetype="png",timeout=120):""" :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 """ifsaveplotandnotisinstance(saveplot,str):raiseValueError("To save the figures enter a filename: ")# Get magntitude binsmagnitude_bins=self._get_magnitudes_from_spacing(catalogue.data["magnitude"],config["magnitude_bin"])dec_time=catalogue.get_decimal_time()completeness_table=np.zeros([len(magnitude_bins)-1,2],dtype=float)min_year=float(np.min(catalogue.data["year"]))max_year=float(np.max(catalogue.data["year"]))+1.0has_completeness=np.zeros(len(magnitude_bins)-1,dtype=bool)forilocinrange(0,len(magnitude_bins)-1):lower_mag=magnitude_bins[iloc]upper_mag=magnitude_bins[iloc+1]idx=np.logical_and(catalogue.data["magnitude"]>=lower_mag,catalogue.data["magnitude"]<upper_mag,)cumvals=np.cumsum(np.ones(np.sum(idx)))plt.plot(dec_time[idx],cumvals,".")plt.xlim(min_year,max_year+5)title_string="Magnitude %5.2f to %5.2f"%(lower_mag,upper_mag)plt.title(title_string)pts=pylab.ginput(1,timeout=timeout)[0]ifpts[0]<=max_year:# Magnitude bin has no completeness!has_completeness[iloc]=Truecompleteness_table[iloc,0]=np.floor(pts[0])completeness_table[iloc,1]=magnitude_bins[iloc]print(completeness_table[iloc,:],has_completeness[iloc])if(config["increment_lock"]and(iloc>0)and(completeness_table[iloc,0]>completeness_table[iloc-1,0])):completeness_table[iloc,0]=completeness_table[iloc-1,0]# Add marker line to indicate completeness pointmarker_line=np.array([[0.0,completeness_table[iloc,0]],[cumvals[-1],completeness_table[iloc,0]],])plt.plot(marker_line[:,0],marker_line[:,1],"r-")ifsaveplot:filename=(saveplot+"_"+("%5.2f"%lower_mag)+("%5.2f"%upper_mag)+"."+filetype)plt.savefig(filename,format=filetype)plt.close()returncompleteness_table[has_completeness,:]
def_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)# 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_bins