Source code for openquake.hmtk.faults.fault_models
# -*- 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.faults.fault_model implements the set of classesto allow for a calculation of the magnitude frequency distribution fromthe geological slip rate"""importnumpyasnpfrommathimportfabsfromopenquake.hazardlib.scalerelimportget_available_scalerelfromopenquake.hazardlib.mfd.evenly_discretizedimportEvenlyDiscretizedMFDfromopenquake.hmtk.modelsimportIncrementalMFDfromopenquake.hmtk.faults.fault_geometriesimport(SimpleFaultGeometry,ComplexFaultGeometry,)fromopenquake.hmtk.sources.simple_fault_sourceimportmtkSimpleFaultSourcefromopenquake.hmtk.sources.complex_fault_sourceimportmtkComplexFaultSourcefromopenquake.hmtk.faultsimportmfdMFD_MAP=mfd.get_available_mfds()SCALE_REL_MAP=get_available_scalerel()DEFAULT_MSR_SIGMA=[(0.0,1.0)]def_update_slip_rates_with_aseismic(slip_rate,aseismic):""" For all the slip rates in the slip rate tuple, multiply by the aseismic coefficient :param list slip_rate: List of tuples (Slip Value, Weight) defining the slip distribution :param float aseismic: Fractional proportion of slip release aseismically :returns: slip - List of tuples (Slip Value, Weight) for adjusted slip rates """return[(slip_val*(1.0-aseismic),weight)forslip_val,weightinslip_rate]
[docs]classRecurrenceBranch(object):""" :class:`openquake.hmtk.faults.fault_model.RecurrenceBranch` is an object to store a set of parameters for recurrence calculations and the corresponding total weight :param str branch_id: Unique branch id :param float area: Fault area (km ^ 2) :param float slip: Fault slip rate (mm / yr) :param msr: Magnitude scaling relation as instance of :class: openquake.hazardlib.scale_rel.base.BaseASR :param float rake: Rake of fault (degrees) :param float shear_modulus: Shear modulus of fault (GPa) :param float disp_length_ratio: Displacement to length ratio of the fault :param float weight: Weight of recurrence model branch :param recurrence: Magnitude frquency distribution as instance of :class: openquake.hmtk.models.IncrementalMFD :param float max_mag: Maximum magnitude from the magnitude frequency distribution :param numpy.ndarray magnitudes: Magnitudes of MFD """def__init__(self,area,slip,msr,rake,shear_modulus,disp_length_ratio=None,msr_sigma=0.0,weight=1.0,):self.branch_id=Noneself.area=areaself.slip=slipself.msr=msrself.msr_sigma=msr_sigmaself.rake=rakeself.shear_modulus=shear_modulusself.disp_length_ratio=disp_length_ratioself.weight=weightself.recurrence=Noneself.max_mag=Noneself.magnitudes=None
[docs]defupdate_weight(self,new_weight):""" Updates the weight by multiplying by the new weight :param float new_weight: Weight to be multiplied by existing weight """self.weight=self.weight*new_weight
[docs]defget_recurrence(self,config):""" Calculates the recurrence model for the given settings as an instance of the openquake.hmtk.models.IncrementalMFD :param dict config: Configuration settings of the magnitude frequency distribution. """model=MFD_MAP[config["Model_Name"]]()model.setUp(config)model.get_mmax(config,self.msr,self.rake,self.area)model.mmax=model.mmax+(self.msr_sigma*model.mmax_sigma)# As the Anderson & Luco arbitrary model requires the input of the# displacement to length ratioif"AndersonLucoAreaMmax"inconfig["Model_Name"]:ifnotself.disp_length_ratio:# If not defined then default to 1.25E-5self.disp_length_ratio=1.25e-5min_mag,bin_width,occur_rates=model.get_mfd(self.slip,self.area,self.shear_modulus,self.disp_length_ratio,)else:min_mag,bin_width,occur_rates=model.get_mfd(self.slip,self.area,self.shear_modulus)self.recurrence=IncrementalMFD(min_mag,bin_width,occur_rates)self.magnitudes=(min_mag+np.cumsum(bin_width*np.ones(len(occur_rates),dtype=float))-bin_width)self.max_mag=np.max(self.magnitudes)
[docs]classmtkActiveFault(object):""" Main class to represent fault source :param int identifier: Identifier Code :param str name: Fault Name :param geometry: Instance of :class:`openquake.hmtk.faults.fault_model.SimpleFaultGeometry` or :class:`openquake.hmtk.faults.fault_model.ComplexFaultGeometry` :param list slip_rate: Slip rate (mm/yr) as list of tuples [(Value, Weight)] :param float aseismic: Aseismic slip coefficient :param float rake: Rake of the fault slip (degrees) :param neotectonic_fault: Instance of :class:`openquake.hmtk.faults.faulted_earth.NeotectonicFault` :param str trt: Tectonic region type :param scale_rel: Scaling relation as list of tuples [(:class: openquake.hazardlib.scalerel.base.BaseASR, Weight)] :param float aspect_ratio: Aspect ratio on fault :param tuple mfd: Tuple ([MFD], [Weight], [Scale_Rel]) defining the magnitude frequency distribution :param list shear_modulus: Shear Modulus (GPa) as list of tuples [(Value, Weight)] :param list disp_length_ratio: Displacement to length ratio as list of tuples [(Value, Weight)] :param list mfd_models: Magnitude frequency distributions as list of instances of :class: openquake.hmtk.faults.fault_model.RecurrenceBranch :param list mfd_models: Magnitude frequency distributions as list of instances of :class: openquake.hmtk.models.IncrementalMFD :param float area: Area of fault (km ^ 2) :param dict config: Dictionary of configuration paramters for magnitude freuency distribution calculation """def__init__(self,identifier,name,geometry,slip_rate,rake,trt,aseismic=0.0,msr_sigma=None,neotectonic_fault=None,scale_rel=None,aspect_ratio=None,shear_modulus=None,disp_length_ratio=None,):""" """self.id=identifierself.name=namemsr_sigma=msr_sigmaorDEFAULT_MSR_SIGMAcond=notisinstance(geometry,SimpleFaultGeometry)andnotisinstance(geometry,ComplexFaultGeometry)ifcond:msg=("Geometry must be instance of "+"openquake.hmtk.faults.fault_geometries.BaseFaultGeometry")raiseIOError(msg)self.geometry=geometryself.aseismic=aseismic# Assert that the sum of the slip rates is 1.0iffabs(np.sum([val[1]forvalinslip_rate])-1.0)>1e-7:raiseValueError("Slip rate weightings must sum to 1.0")self.slip=_update_slip_rates_with_aseismic(slip_rate,self.aseismic)self.rake=rakeself.neotectonic_fault=neotectonic_faultself.trt=trtself.rupt_aspect_ratio=aspect_ratioself.mfd=([],[],[])self.shear_modulus=shear_modulusself.disp_length_ratio=disp_length_ratioself.mfd_models=[]self.msr=scale_relself.msr_sigma=msr_sigmaself.area=self.geometry.get_area()self.config=Noneself.regionalisation=Noneself.catalogue=None
[docs]defget_tectonic_regionalisation(self,regionalisation,region_type=None):""" Defines the tectonic region and updates the shear modulus, magnitude scaling relation and displacement to length ratio using the regional values, if not previously defined for the fault :param regionalistion: Instance of the :class: openquake.hmtk.faults.tectonic_regionalisaion.TectonicRegionalisation :param str region_type: Name of the region type - if not in regionalisation an error will be raised """ifregion_type:self.trt=region_typeifself.trtnotinregionalisation.key_list:raiseValueError("Tectonic region classification missing or ""not defined in regionalisation")foriloc,key_valinenumerate(regionalisation.key_list):ifself.trtinkey_val:self.regionalisation=regionalisation.regionalisation[iloc]# Update undefined shear modulus from tectonic regionalisationifnotself.shear_modulus:self.shear_modulus=self.regionalisation.shear_modulus# Update undefined scaling relation from tectonic# regionalisationifnotself.msr:self.msr=self.regionalisation.scaling_rel# Update undefined displacement to length ratio from tectonic# regionalisationifnotself.disp_length_ratio:self.disp_length_ratio=(self.regionalisation.disp_length_ratio)breakreturn
[docs]defselect_catalogue(self,selector,distance,distance_metric="rupture",upper_eq_depth=None,lower_eq_depth=None,):""" Select earthquakes within a specied distance of the fault """ifselector.catalogue.get_number_events()<1:raiseValueError("No events found in catalogue!")# rupture metric is selectedif"rupture"indistance_metric:# Use rupture distanceself.catalogue=selector.within_rupture_distance(self.geometry.surface,distance,upper_depth=upper_eq_depth,lower_depth=lower_eq_depth,)else:# Use Joyner-Boore distanceself.catalogue=selector.within_joyner_boore_distance(self.geometry.surface,distance,upper_depth=upper_eq_depth,lower_depth=lower_eq_depth,)
def_generate_branching_index(self):""" Generates a branching index (i.e. a list indicating the number of branches in each branching level. Current branching levels are: 1) Slip 2) MSR 3) Shear Modulus 4) DLR 5) MSR_Sigma 6) Config :returns: * branch_index - A 2-D numpy.ndarray where each row is a pointer to a particular combination of values * number_branches - Total number of branches (int) """branch_count=np.array([len(self.slip),len(self.msr),len(self.shear_modulus),len(self.disp_length_ratio),len(self.msr_sigma),len(self.config),])n_levels=len(branch_count)number_branches=np.prod(branch_count)branch_index=np.zeros([number_branches,n_levels],dtype=int)cumval=1dstep=1e-9forilocinrange(0,n_levels):idx=np.linspace(0.0,float(branch_count[iloc])-dstep,number_branches//cumval,)branch_index[:,iloc]=np.reshape(np.tile(idx,[cumval,1]),number_branches)cumval*=branch_count[iloc]returnbranch_index.tolist(),number_branches
[docs]defgenerate_config_set(self,config):""" Generates a list of magnitude frequency distributions and renders as a tuple :param dict/list config: Configuration paramters of magnitude frequency distribution """ifisinstance(config,dict):# Configuration list contains only one elementself.config=[(config,1.0)]elifisinstance(config,list):# Multiple configurations with correscponding weightstotal_weight=0.0self.config=[]forparamsinconfig:weight=params["Model_Weight"]total_weight+=params["Model_Weight"]self.config.append((params,weight))iffabs(total_weight-1.0)>1e-7:raiseValueError("MFD config weights do not sum to 1.0 for ""fault %s"%self.id)else:raiseValueError("MFD config must be input as dictionary or list!")
[docs]defgenerate_recurrence_models(self,collapse=False,bin_width=0.1,config=None,rendered_msr=None):""" Iterates over the lists of values defining epistemic uncertainty in the parameters and calculates the corresponding recurrence model At present epistemic uncertainty is supported for: 1) slip rate, 2) magnitude scaling relation, 3) shear modulus, 4) displacement to length ratio) and 5) recurrence model. :param list config: List of MFD model configurations :param bool collapse: Boolean flag indicating whether to collapse the logic tree branches :param float bin_width: If collapsing the logic tree branches the reference mfd must be defined. The minimum and maximum magnitudes are updated from the model, but the bin width must be specified here :param list/dict config: Configuration (or sets of configurations) of the recurrence calculations :param rendered_msr: If collapsing the logic tree branches a resulting magnitude scaling relation must be defined as instance of :class: openquake.hazardlib.scalerel.base.BaseASR """ifcollapseandnotrendered_msr:raiseValueError("Collapsing logic tree branches requires input ""of a single msr for rendering sources")# Generate a set of tuples with corresponding weightsifconfigisnotNone:self.generate_config_set(config)ifnotisinstance(self.config,list):raiseValueError("MFD configuration missing or incorrectly ""formatted")# Generate the branching indexbranch_index,_number_branches=self._generate_branching_index()mmin=np.infmmax=-np.infforidxinbranch_index:tuple_list=[]# Get sliptuple_list.append(self.slip[idx[0]])# Get msrtuple_list.append(self.msr[idx[1]])# Get shear modulustuple_list.append(self.shear_modulus[idx[2]])# Get displacement length ratiotuple_list.append(self.disp_length_ratio[idx[3]])# Get msr sigmatuple_list.append(self.msr_sigma[idx[4]])# Get configtuple_list.append(self.config[idx[5]])# Calculate branch weight as product of tuple weightsbranch_weight=np.prod(np.array([val[1]forvalintuple_list]))# Instantiate recurrence modelmodel=RecurrenceBranch(self.area,tuple_list[0][0],tuple_list[1][0],self.rake,tuple_list[2][0],tuple_list[3][0],tuple_list[4][0],weight=branch_weight,)model.get_recurrence(tuple_list[5][0])self.mfd_models.append(model)# Update the total minimum and maximum magnitudes for the faultifmodel.recurrence.min_mag<mmin:mmin=model.recurrence.min_magifnp.max(model.magnitudes)>mmax:mmax=np.max(model.magnitudes)ifcollapse:self.mfd=([self.collapse_branches(mmin,bin_width,mmax)],[1.0],[rendered_msr],)else:mfd_mods=[]mfd_wgts=[]mfd_msr=[]formodelinself.mfd_models:mfd_mods.append(IncrementalMFD(model.recurrence.min_mag,model.recurrence.bin_width,model.recurrence.occur_rates,))mfd_wgts.append(model.weight)mfd_msr.append(model.msr)self.mfd=(mfd_mods,mfd_wgts,mfd_msr)
[docs]defcollapse_branches(self,mmin,bin_width,mmax):""" Collapse the logic tree branches into a single IncrementalMFD :param float mmin: Minimum magnitude of reference mfd :param float bin_width: Bin width of reference mfd :param float mmax: Maximum magnitude of reference mfd :returns: :class: openquake.hmtk.models.IncrementalMFD """master_mags=np.arange(mmin,mmax+(bin_width/2.0),bin_width)master_rates=np.zeros(len(master_mags),dtype=float)formodelinself.mfd_models:id0=np.logical_and(master_mags>=np.min(model.magnitudes)-1e-9,master_mags<=np.max(model.magnitudes)+1e-9,)# Use interpolation in log10-y valuesyvals=np.log10(model.recurrence.occur_rates)interp_y=np.interp(master_mags[id0],model.magnitudes,yvals)master_rates[id0]=master_rates[id0]+(model.weight*10.0**interp_y)returnIncrementalMFD(mmin,bin_width,master_rates)
[docs]defgenerate_fault_source_model(self):""" Creates a resulting `openquake.hmtk` fault source set. :returns: source_model - list of instances of either the :class: `openquake.hmtk.sources.simple_fault_source.mtkSimpleFaultSource` or :class: `openquake.hmtk.sources.complex_fault_source.mtkComplexFaultSource` model_weight - Corresponding weights for each source model """source_model=[]model_weight=[]forilocinrange(0,self.get_number_mfd_models()):model_mfd=EvenlyDiscretizedMFD(self.mfd[0][iloc].min_mag,self.mfd[0][iloc].bin_width,self.mfd[0][iloc].occur_rates.tolist(),)ifisinstance(self.geometry,ComplexFaultGeometry):# Complex fault classsource=mtkComplexFaultSource(self.id,self.name,self.trt,self.geometry.surface,self.mfd[2][iloc],self.rupt_aspect_ratio,model_mfd,self.rake,)source.fault_edges=self.geometry.traceelse:# Simple Fault sourcesource=mtkSimpleFaultSource(self.id,self.name,self.trt,self.geometry.surface,self.geometry.dip,self.geometry.upper_depth,self.geometry.lower_depth,self.mfd[2][iloc],self.rupt_aspect_ratio,model_mfd,self.rake,)source.fault_trace=self.geometry.tracesource_model.append(source)model_weight.append(self.mfd[1][iloc])returnsource_model,model_weight
[docs]defget_number_mfd_models(self):""" Returns the number of mfd models for a given fault model """returnlen(self.mfd[0])