# The Hazard Library# Copyright (C) 2012-2025 GEM Foundation## This program 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.## This program is distributed in the hope that it will be useful,# but WITHOUT ANY WARRANTY; without even the implied warranty of# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the# GNU Affero General Public License for more details.## You should have received a copy of the GNU Affero General Public License# along with this program. If not, see <http://www.gnu.org/licenses/>."""Module :mod:`openquake.hazardlib.source.base` defines a base class forseismic sources."""importabcimportzlibfromdataclassesimportdataclassimportnumpyfromopenquake.baselibimportgeneralfromopenquake.hazardlibimportmfdfromopenquake.hazardlib.pmfimportPMFfromopenquake.hazardlib.tomimportPoissonTOMfromopenquake.hazardlib.calc.filtersimportmagstr,split_sourcefromopenquake.hazardlib.geoimportPointfromopenquake.hazardlib.geo.surface.planarimportbuild_planar,PlanarSurfacefromopenquake.hazardlib.geo.surface.multiimportMultiSurfacefromopenquake.hazardlib.source.ruptureimport(ParametricProbabilisticRupture,NonParametricProbabilisticRupture,EBRupture)
[docs]defget_code2cls():""" :returns: a dictionary source code -> source class """dic={}forclsingeneral.gen_subclasses(BaseSeismicSource):ifhasattr(cls,'code'):dic[cls.code]=clsreturndic
[docs]defis_poissonian(src):""" :returns: True if the underlying source is poissonian, false otherwise """ifsrc.code==b'F':# multiFaultreturnsrc.infer_occur_rateselifsrc.code==b'N':# nonParametricreturnFalsereturnTrue
[docs]defpoisson_sample(src,eff_num_ses,seed):""" :param src: a poissonian source :param eff_num_ses: number of stochastic event sets * number of samples :param seed: stochastic seed :yields: triples (rupture, rup_id, num_occurrences) """rng=numpy.random.default_rng(seed)ifhasattr(src,'temporal_occurrence_model'):tom=src.temporal_occurrence_modelelse:# multifaulttom=PoissonTOM(src.investigation_time)rupids=src.offset+numpy.arange(src.num_ruptures)ifnothasattr(src,'nodal_plane_distribution'):ifsrc.code==b'F':# multifaults=src.get_sections()fori,rateinenumerate(src.occur_rates):# NB: rng.poisson called inside to save memorynum_occ=rng.poisson(rate*tom.time_span*eff_num_ses)ifnum_occ==0:# skipcontinueidxs=src.rupture_idxs[i]iflen(idxs)==1:sfc=s[idxs[0]]else:sfc=MultiSurface([s[idx]foridxinidxs])hypo=s[idxs[0]].get_middle_point()rup=ParametricProbabilisticRupture(src.mags[i],src.rakes[i],src.tectonic_region_type,hypo,sfc,src.occur_rates[i],tom)yieldrup,rupids[i],num_occelse:# simple or complex faultruptures=list(src.iter_ruptures())rates=numpy.array([rup.occurrence_rateforrupinruptures])occurs=rng.poisson(rates*tom.time_span*eff_num_ses)forrup,rupid,num_occinzip(ruptures,rupids,occurs):ifnum_occ:yieldrup,rupid,num_occreturn# else (multi)point sources and area sourcesusd=src.upper_seismogenic_depthlsd=src.lower_seismogenic_depthrup_args=[]rates=[]forpsinsplit_source(src):ifnothasattr(ps,'location'):# unsplit containing a single source[ps]=srclon,lat=ps.location.x,ps.location.yformag,mag_occ_rateinps.get_annual_occurrence_rates():fornp_prob,npinps.nodal_plane_distribution.data:forhc_prob,hc_depthinps.hypocenter_distribution.data:args=(mag_occ_rate,np_prob,hc_prob,mag,np,lon,lat,hc_depth,ps)rup_args.append(args)rates.append(mag_occ_rate*np_prob*hc_prob)eff_rates=numpy.array(rates)*tom.time_span*eff_num_sesoccurs=rng.poisson(eff_rates)fornum_occ,args,rupid,rateinzip(occurs,rup_args,rupids,rates):ifnum_occ:_,np_prob,hc_prob,mag,np,lon,lat,hc_depth,ps=argshc=Point(lon,lat,hc_depth)hdd=numpy.array([(1.,hc.depth)])[[[planar]]]=build_planar(ps.get_planin([(1.,mag)],[(1.,np)]),hdd,lon,lat,usd,lsd)rup=ParametricProbabilisticRupture(mag,np.rake,ps.tectonic_region_type,hc,PlanarSurface.from_(planar),rate,tom)yieldrup,rupid,num_occ
[docs]deftimedep_sample(src,eff_num_ses,seed):""" :param src: a time-dependent source :param eff_num_ses: number of stochastic event sets * number of samples :param seed: stochastic seed :yields: triples (rupture, rup_id, num_occurrences) """rng=numpy.random.default_rng(seed)rupids=src.offset+numpy.arange(src.num_ruptures)ifsrc.code==b'F':# time-dependent multifaults=src.get_sections()fori,probsinenumerate(src.probs_occur):cdf=numpy.cumsum(probs)num_occ=numpy.digitize(rng.random(eff_num_ses),cdf).sum()ifnum_occ==0:# ignore non-occurring rupturescontinueidxs=src.rupture_idxs[i]iflen(idxs)==1:sfc=s[idxs[0]]else:sfc=MultiSurface([s[idx]foridxinidxs])hypo=sfc.get_middle_point()pmf=PMF([(p,o)foro,pinenumerate(probs)])yield(NonParametricProbabilisticRupture(src.mags[i],src.rakes[i],src.tectonic_region_type,hypo,sfc,pmf),rupids[i],num_occ)else:# time-dependent nonparametricmutex_weight=getattr(src,'mutex_weight',1)forrup,rupidinzip(src.iter_ruptures(),rupids):occurs=rup.sample_number_of_occurrences(eff_num_ses,rng)ifmutex_weight<1:# consider only the occurrencies below the mutex_weightoccurs*=(rng.random(eff_num_ses)<mutex_weight)num_occ=occurs.sum()ifnum_occ:yieldrup,rupid,num_occ
[docs]classBaseSeismicSource(metaclass=abc.ABCMeta):""" Base class representing a seismic source, that is a structure generating earthquake ruptures. :param source_id: Some (numeric or literal) source identifier. Supposed to be unique within the source model. :param name: String, a human-readable name of the source. :param tectonic_region_type: Source's tectonic regime. See :class:`openquake.hazardlib.const.TRT`. """id=-1# to be settrt_smr=0# set by the enginensites=1# set when filtering the sourcesplittable=Truechecksum=0# set in source_readerweight=0.001# set in contextsesites=0# updated in estimate_weightoffset=0# set in fix_src_offsettrt_smr=-1# set by the enginenum_ruptures=0# set by the engineseed=None# set by the engine@abc.abstractpropertydefMODIFICATIONS(self):pass@propertydeftrt_smrs(self):""" :returns: a list of integers (usually of 1 element) """trt_smr=self.trt_smrreturn(trt_smr,)ifisinstance(trt_smr,int)elsetrt_smr
[docs]defserial(self,ses_seed):""" :returns: a random seed derived from source_id and ses_seed """returnzlib.crc32(self.source_id.encode('ascii'),ses_seed)
[docs]defis_gridded(self):""" :returns: True if the source contains only gridded ruptures """returnFalse
[docs]@abc.abstractmethoddefiter_ruptures(self,**kwargs):""" Get a generator object that yields probabilistic ruptures the source consists of. :returns: Generator of instances of sublclass of :class: `~openquake.hazardlib.source.rupture.BaseProbabilisticRupture`. """
[docs]defiter_meshes(self):""" Yields the meshes underlying the ruptures """forrupinself.iter_ruptures():ifisinstance(rup.surface,MultiSurface):forsurfinrup.surface.surfaces:yieldsurf.meshelse:yieldrup.surface.mesh
[docs]defsample_ruptures(self,eff_num_ses,ses_seed):""" :param eff_num_ses: number of stochastic event sets * number of samples :yields: triples (rupture, trt_smr, num_occurrences) """seed=self.serial(ses_seed)sample=poisson_sampleifis_poissonian(self)elsetimedep_sampleforrup,rupid,num_occinsample(self,eff_num_ses,seed):ifself.smweight<1andhasattr(rup,'occurrence_rate'):# defined only for poissonian sources# needed to get convergency of the frequency to the rate# tested only in oq-risk-tests etna0rup.occurrence_rate*=self.smweightebr=EBRupture(rup,self.id,self.trt_smr,num_occ,rupid)ebr.seed=ebr.id+ses_seedyieldebr
[docs]defget_mags(self):""" :returns: the magnitudes of the ruptures contained in the source """mags=set()ifhasattr(self,'get_annual_occurrence_rates'):formag,rateinself.get_annual_occurrence_rates():mags.add(mag)elifhasattr(self,'mags'):# MultiFaultSourcemags.update(self.mags)else:# nonparametricforrup,pmfinself.data:mags.add(rup.mag)returnsorted(mags)
[docs]defget_magstrs(self):""" :returns: the magnitudes of the ruptures contained as strings """ifhasattr(self,'mags'):# MultiFaultSourcemags={magstr(mag)formaginself.mags}elifhasattr(self,'data'):# nonparametricmags={magstr(item[0].mag)foriteminself.data}else:mags={magstr(item[0])foriteminself.get_annual_occurrence_rates()}returnsorted(mags)
def__iter__(self):""" Override to implement source splitting """yieldself
[docs]@abc.abstractmethoddefcount_ruptures(self):""" Return the number of ruptures that will be generated by the source. """
[docs]@abc.abstractmethoddefget_min_max_mag(self):""" Return minimum and maximum magnitudes of the ruptures generated by the source. """
[docs]defmodify(self,modification,parameters):""" Apply a single modificaton to the source parameters Reflects the modification method and calls it passing ``parameters`` as keyword arguments. Modifications can be applied one on top of another. The logic of stacking modifications is up to a specific source implementation. :param modification: String name representing the type of modification. :param parameters: Dictionary of parameters needed for modification. :raises ValueError: If ``modification`` is missing from the attribute `MODIFICATIONS`. """ifmodificationnotinself.MODIFICATIONS:raiseValueError('Modification %s is not supported by %s'%(modification,type(self).__name__))meth=getattr(self,'modify_%s'%modification)meth(**parameters)
[docs]defto_xml(self):""" Convert the source into an XML string, very useful for debugging """fromopenquake.hazardlibimportnrml,sourcewriterreturnnrml.to_string(sourcewriter.obj_to_node(self))
def__repr__(self):""" String representation of a source, displaying the source class name and the source id. """return'<%s%s, weight=%.2f>'%(self.__class__.__name__,self.source_id,self.weight)
[docs]classParametricSeismicSource(BaseSeismicSource,metaclass=abc.ABCMeta):""" Parametric Seismic Source generates earthquake ruptures from source parameters, and associated probabilities of occurrence are defined through a magnitude frequency distribution and a temporal occurrence model. :param mfd: Magnitude-Frequency distribution for the source. See :mod:`openquake.hazardlib.mfd`. :param rupture_mesh_spacing: The desired distance between two adjacent points in source's ruptures' mesh, in km. Mainly this parameter allows to balance the trade-off between time needed to compute the :meth:`distance <openquake.hazardlib.geo.surface.base.BaseSurface.get_min_distance>` between the rupture surface and a site and the precision of that computation. :param magnitude_scaling_relationship: Instance of subclass of :class:`openquake.hazardlib.scalerel.base.BaseMSR` to describe how does the area of the rupture depend on magnitude and rake. :param rupture_aspect_ratio: Float number representing how much source's ruptures are more wide than tall. Aspect ratio of 1 means ruptures have square shape, value below 1 means ruptures stretch vertically more than horizontally and vice versa. :param temporal_occurrence_model: Instance of :class:`openquake.hazardlib.tom.PoissonTOM` defining temporal occurrence model for calculating rupture occurrence probabilities :raises ValueError: If either rupture aspect ratio or rupture mesh spacing is not positive (if not None). """def__init__(self,source_id,name,tectonic_region_type,mfd,rupture_mesh_spacing,magnitude_scaling_relationship,rupture_aspect_ratio,temporal_occurrence_model):self.source_id=source_idself.name=nameself.tectonic_region_type=tectonic_region_typeifrupture_mesh_spacingisnotNoneandnotrupture_mesh_spacing>0:raiseValueError('rupture mesh spacing must be positive')ifrupture_aspect_ratioisnotNoneandnotrupture_aspect_ratio>0:raiseValueError('rupture aspect ratio must be positive')self.mfd=mfdself.rupture_mesh_spacing=rupture_mesh_spacingself.magnitude_scaling_relationship=magnitude_scaling_relationshipself.rupture_aspect_ratio=rupture_aspect_ratioself.temporal_occurrence_model=temporal_occurrence_model
[docs]defget_annual_occurrence_rates(self,min_rate=0):""" Get a list of pairs "magnitude -- annual occurrence rate". The list is taken from assigned MFD object (see :meth: `openquake.hazardlib.mfd.base.BaseMFD.get_annual_occurrence_rates`) with simple filtering by rate applied. :param min_rate: A non-negative value to filter magnitudes by minimum annual occurrence rate. Only magnitudes with rates greater than that are included in the result list. :returns: A list of two-item tuples -- magnitudes and occurrence rates. """scaling_rate=getattr(self,'scaling_rate',1)return[(mag,occ_rate*scaling_rate)formag,occ_rateinself.mfd.get_annual_occurrence_rates()ifmin_rateisNoneorocc_rate>min_rate]
[docs]defget_min_max_mag(self):""" Get the minimum and maximum magnitudes of the ruptures generated by the source from the underlying MFD. """returnself.mfd.get_min_max_mag()
[docs]defmodify_set_msr(self,new_msr):""" Updates the MSR originally assigned to the source :param new_msr: An instance of the :class:`openquake.hazardlib.scalerel.BaseMSR` """self.magnitude_scaling_relationship=new_msr
[docs]defmodify_set_slip_rate(self,slip_rate:float):""" Updates the slip rate assigned to the source :param slip_rate: The value of slip rate [mm/yr] """self.slip_rate=slip_rate
[docs]defmodify_set_mmax_truncatedGR(self,mmax:float):""" Updates the mmax assigned. This works on for parametric MFDs.s :param mmax: The value of the new maximum magnitude """# Check that the current src has a TruncatedGRMFD MFDmsg='This modification works only when the source MFD is a 'msg+='TruncatedGRMFD'assertself.mfd.__class__.__name__=='TruncatedGRMFD',msgself.mfd.max_mag
[docs]defmodify_recompute_mmax(self,epsilon:float=0):""" Updates the value of mmax using the msr and the area of the fault :param epsilon: Number of standard deviations to be added or substracted """msr=self.magnitude_scaling_relationshiparea=self.get_fault_surface_area()# area in km^2mag=msr.get_median_mag(area=area,rake=self.rake)std=msr.get_std_dev_mag(area=area,rake=self.rake)self.mfd.max_mag=mag+epsilon*std
[docs]defmodify_adjust_mfd_from_slip(self,slip_rate:float,rigidity:float,constant_term:float=9.1,recompute_mmax:float=None):""" :param slip_rate: A float defining slip rate [in mm] :param rigidity: A float defining material rigidity [in GPa] :param constant_term: Constant term of the equation used to compute log M0 from magnitude """# Check that the current src has a TruncatedGRMFD MFDmsg='This modification works only when the source MFD is a 'msg+='TruncatedGRMFD'assertself.mfd.__class__.__name__=='TruncatedGRMFD',msg# Compute momentarea=self.get_fault_surface_area()*1e6# area in m^2rigidity*=1e9# rigidity in Paslip_rate*=1e-3# slip rate in mmo=rigidity*area*slip_rate# Update the MFDmin_mag=self.mfd.min_magmax_mag=self.mfd.max_magbin_w=self.mfd.bin_widthb_val=self.mfd.b_valself.mfd=mfd.TruncatedGRMFD.from_moment(min_mag,max_mag,bin_w,b_val,mo,constant_term)