How to Calculate PGA/SA in a Particular Exposure Site Given an Earthquake Event using openquake.engine GMPE in Python

[Randhy Pratama]

Hi all,

I am trying to compute a PGA/SA in a particular exposure site given an earthquake event using openquake.engine GMPE in Python. I have seen the document about GMPE that was given in OQ at https://docs.openquake.org/oq-engine/2.7/openquake.hazardlib.gsim.html. Since there are a lot of arguments to be input when calculating the PGA/SA using one of those GMPEs in Python, can someone help me to give an example of how to calculate the mean and str. dev of PGA/SA? Or where can I read about this topic?

Thank you for your help.

Regards,
Randhy Pratama

[Peter Pažák]

Hi, from another thread this is an example (works in OQ 3.16, as it uses internal OQ code which may change, use with caution:)

from openquake.hazardlib.source import PointSource
from openquake.hazardlib.mfd import TruncatedGRMFD
from openquake.hazardlib.scalerel import WC1994
from openquake.hazardlib.geo import Point, NodalPlane
from openquake.hazardlib.pmf import PMF
from openquake.hazardlib.tom import PoissonTOM
from openquake.hazardlib.site import Site, SiteCollection
from openquake.hazardlib.imt import PGA
from openquake.hazardlib.const import StdDev
from openquake.hazardlib.contexts import ContextMaker
from openquake.hazardlib.valid import gsim
from openquake.hazardlib import valid
import numpy as np

#import 2 GMPEs
gmpes = [gsim('AbrahamsonSilva2008'),gsim('CampbellBozorgnia2008')]

# this would be suitable to explore magnitude scaling, by defining a Point source
location = Point(9.1500, 45.1833,30)
src = PointSource(
    source_id='1',
    name='point',
    tectonic_region_type='Active Shallow Crust',
    mfd=TruncatedGRMFD(min_mag=5., max_mag=6.5, bin_width=0.1, a_val=0.01, b_val=0.98),
    rupture_mesh_spacing=2,
    magnitude_scaling_relationship=WC1994(),
    rupture_aspect_ratio=1,
    temporal_occurrence_model=PoissonTOM(50.),
    upper_seismogenic_depth=2.,
    lower_seismogenic_depth=12.,
    location=location,
    nodal_plane_distribution=PMF([(1., NodalPlane(strike=45, dip=50, rake=0))]),
    hypocenter_distribution=PMF([(1, 10.)]) #depth 10 km
)

#these magnitudes were generated
mags = [r.mag for r in src.iter_ruptures()] #mags from 5.05 to 6.45 step 0.1 (15 magnitudes)

# site for which we compute
site_collection = SiteCollection([Site(location=Point(9.200, 45.2,0), vs30=450., vs30measured=True, z1pt0=150., z2pt5=1.5)])

# intensity measure types for which we compute
imtls = {s: [0] for s in ['PGA','SA(0.3)']} #required for context maker, M = 2 IMTs


context_maker = ContextMaker('Active Shallow Crust',gmpes,{'imtls': imtls, 'mags': ['%1.2f'%m for m in mags]}) #necessary contexts builder
ctxs = context_maker.from_srcs([src],site_collection) #returns rupture contexts

gms = context_maker.get_mean_stds(ctxs) #calculate ground motions and stds, returns array of shape (4, G, M, N)

print(gms.shape) #4 values (0=median, 1=std_total, 2=std_intra, 3=std_inter), then G=2 gsims, M=2 IMTs, 15 scenarios = magnitudes

print(np.exp(gms[0][0][0][0])) #median value, for AbrahamsonSilva2008, PGA, Mw = 5.05

print(gms[1][0][0][0]) #standard deviation, for AbrahamsonSilva2008, PGA, Mw = 5.05

Dátum: štvrtok 9. februára 2023, čas: 13:37:32 UTC+1, odosielateľ: pratama...@gmail.com

[Angela Stallone]

Hi Peter,

I changed

ctx = context_maker.get_ctxs([source], Site_Collection)

to

ctx = context_maker.from_srcs([source], Site_Collection)

This, however, is not working with BaseRupture.

How to get contexts in this case?

Thank you,

Angela

[Marco Pagani]

You can use something like this:

imt_list = ['PGA', 'SA(0.2)']
oqp = {'imtls': {k: [] for k in imt_list}, 'mags': mags_str}
ctxm = ContextMaker('fake', gmms, oqp)
ctxs = list(ctxm.get_ctx_iter([rup], sites))

--
You received this message because you are subscribed to the Google Groups "OpenQuake Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to openquake-use...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/openquake-users/8ee3752d-4781-46d2-b46e-b666869d7be6n%40googlegroups.com.

[Angela Stallone]

Thank you Marco!

After implementing that edit, I am now getting this error:
File "/home/shake/miniconda/envs/shakemap/lib/python3.9/site-packages/openquake/hazardlib/calc/gmf.py", line 156, in __init__
    self.seed = rupture.seed
AttributeError: 'BaseRupture' object has no attribute 'seed'

Note: after updating OQ, I have also started to get this warning:

TusaLanger2016RepiBA08SE cannot contain the methods ['_get_stddevs', '_compute_distance', '_compute_magnitude', '_get_site_amplification', '_get_site_type_dummy_variables']
TusaLanger2016RepiSP87SE cannot contain the methods ['_compute_distance', '_compute_magnitude']
TusaLanger2016Rhypo cannot contain the methods ['_compute_distance']

Anything to worry about?

Thank you!

Best,

Angela

[Peter Pažák]

Why are you trying to use BaseRupture?

BaseRupture is a basic class that is extended by the subclasses to implement particular rupture class, you should probably use those, not the meta-class itself.

Maybe if you explain in more detail what you are trying to achieve, we can recommend a working example...

Peter

Dátum: štvrtok 21. decembra 2023, čas: 11:24:54 UTC+1, odosielateľ: angela....@ingv.it

[Angela Stallone]

Hi Peter,

thanks for your reply!

I am developing a code which implements OpenQuake modules to generate GMFs for a set of source scenarios.

With previous versions of OQ, these lines of code were working as expected:

planar_surface = surface.PlanarSurface.from_hypocenter(

        hypoc=Hypocenter,

        msr=msr(), 

        mag=Mag,

        aratio=rupture_aratio,

        strike=Strike,

        dip=Dip,

        rake=Rake,

        )

source = BaseRupture(

        mag=Mag,

        rake=Rake,

        tectonic_region_type=tectonicRegionType,

        hypocenter=Hypocenter,

        surface=planar_surface

        )

ctx = context_maker.get_ctxs([source], Site_Collection)

For the code purpose, this was ideal, as each source scenario is defined by these parameters:

magnitude, longitude, latitude, hypocenter depth, strike, dip, rake.

Now, changing the last code line to

ctx = list(context_maker.get_ctx_iter([source], Site_Collection))

is causing the seed error. 

Hope I can still reproduce the same outputs I used to have using the old OQ version (3.13.0) with the new (3.19.0).

Thanks a lot!

Best regards,

Angela

You received this message because you are subscribed to a topic in the Google Groups "OpenQuake Users" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/openquake-users/kOLGwM8f_6g/unsubscribe.
To unsubscribe from this group and all its topics, send an email to openquake-use...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/openquake-users/4fdf3dda-3cde-48c8-9887-e53f2877bb76n%40googlegroups.com.

[Peter Pažák]

I believe it is safer to create rupture inside of the PointSource in the way I had above, because source is expected not rupture itself...

(not sure about your exact settings, but hope you can adopt):

from openquake.hazardlib.source import PointSource
from openquake.hazardlib.mfd import ArbitraryMFD

from openquake.hazardlib.scalerel import WC1994
from openquake.hazardlib.geo import Point, NodalPlane
from openquake.hazardlib.pmf import PMF
from openquake.hazardlib.tom import PoissonTOM
from openquake.hazardlib.site import Site, SiteCollection
from openquake.hazardlib.imt import PGA
from openquake.hazardlib.const import StdDev
from openquake.hazardlib.contexts import ContextMaker
from openquake.hazardlib.valid import gsim
from openquake.hazardlib import valid
import numpy as np

#import 2 GMPEs
gmpes = [gsim('AbrahamsonSilva2008'),gsim('CampbellBozorgnia2008')]

#better to have ruptures created from a source, so it fills additional information

Mag = 6.5

src = PointSource(

    source_id='1',
    name='point',
    tectonic_region_type='Active Shallow Crust',

    mfd=ArbitraryMFD([Mag], [0.01]),
    rupture_mesh_spacing=1,
    magnitude_scaling_relationship=WC1994(),
    rupture_aspect_ratio=1.5,
    temporal_occurrence_model=PoissonTOM(50.),
    upper_seismogenic_depth=0.,
    lower_seismogenic_depth=50.,
    location=Point(9.1500, 45.1833,10),

    nodal_plane_distribution=PMF([(1., NodalPlane(strike=45, dip=50, rake=0))]),
    hypocenter_distribution=PMF([(1, 10.)]) #depth 10 km
)

# site for which we compute
site_collection = SiteCollection([Site(location=Point(9.200, 45.2,0), vs30=450., vs30measured=True, z1pt0=150., z2pt5=1.5)])

# intensity measure types for which we compute
imtls = {s: [0] for s in ['PGA','SA(0.3)']} #required for context maker, M = 2 IMTs

context_maker = ContextMaker('Active Shallow Crust',gmpes,{'imtls': imtls, 'mags': [Mag]}) #necessary contexts builder
ctxs = context_maker.from_srcs([src],site_collection) #returns rupture context, from the point source

#we can calculate GMFs now

gms = context_maker.get_mean_stds(ctxs) #calculate ground motions and stds, returns array of shape (4, G, M, N)

print(gms.shape) #4 values (0=median, 1=std_total, 2=std_intra, 3=std_inter), then G=2 gsims, M=2 IMTs, 1 scenario = magnitude

print(np.exp(gms[0][0][0][0])) #median value, for AbrahamsonSilva2008, PGA, Mw = 5.05

print(gms[1][0][0][0]) #standard deviation, for AbrahamsonSilva2008, PGA, Mw = 5.05
 

this works for me in 3.19 without a problem.

Peter

Dátum: nedeľa 24. decembra 2023, čas: 12:37:13 UTC+1, odosielateľ: angela....@ingv.it

[Michele Simionato]

On Thursday, December 21, 2023 at 11:24:54 AM UTC+1 angela....@ingv.it wrote:

Thank you Marco!

Note: after updating OQ, I have also started to get this warning:

TusaLanger2016RepiBA08SE cannot contain the methods ['_get_stddevs', '_compute_distance', '_compute_magnitude', '_get_site_amplification', '_get_site_type_dummy_variables']
TusaLanger2016RepiSP87SE cannot contain the methods ['_compute_distance', '_compute_magnitude']
TusaLanger2016Rhypo cannot contain the methods ['_compute_distance']

Anything to worry about?

Yes, it means you have GMPEs like they were 3/4 years ago. At that time we changed the way to implement GMPEs and now you should use functions and not methods (among other things). The warnings may very well become errors in the future so you must absolutely change your implementation. Just see how the Tusa Langer GMPEs are implemented in the current master of the engine

and follow that approach. The change of 3/4 years ago was made to implement vectorization by rupture which is the reason why now single site calculations are up to 200x faster than before.

                  Michele Simionato

 

[Angela Stallone]

Hi Peter,

the setup of my problem is actually more similar to this one: https://groups.google.com/g/openquake-users/c/MO_VAfYsoI0/m/naieQ6eWBAAJ, where the BaseRupture class is used as well.

Within my code, OQ functionalities are incorporated solely to generate GMFs for different ruptures at a set of sites using GMMs. 

Ruptures are defined as basic planar surfaces models. Inputs such as temporal_occurrence_model or mfd are missing, either because they are out of scope or because they have been pre-determined (e.g. the magnitude).

The combination surface.PlanarSurface.from_hypocenter +  BaseRupture worked fine with older versions of OQ.

Indeed, it is the best option for my code's purpose and I hope to ensure continuity in results after updating to the new version. 

Thanks a lot for your help!

Best regards,

Angela

[Angela Stallone]

Hi Michele,

thank you for your help.

Best,

Angela

[Peter Pažák]

Hi Angela, yes, the approach you are referring to was possible in the older 3.11 version, for the newer versions I believe the updated

script in this thread is more reasonable and robust - it leaves the calculation of one rupture surface fully to the standard OQ algorithm

taking just magnitude, depth, strike and dip as and input. So you would like to setup the rupture geometry manually yourself?

Thank you

Peter

Dátum: utorok 2. januára 2024, čas: 17:50:07 UTC+1, odosielateľ: angela....@ingv.it

[Angela Stallone]

Hi Peter,

it is perfectly fine for me to implement the script you posted in the thread.

The problem is that, in order to use your script, I also need to provide the mfd, the rupture mesh spacing, the temporal occurrence model, the upper/lower seismogenic depths. 

Unfortunately, I do not have this information or I cannot set it by myself, given the purpose of my code.

My code only implements OpenQuake modules to generate GMFs for a set of source scenarios defined by the parameters below:

magnitude, longitude, latitude, hypocenter depth, strike, dip, rake, fault area, fault length, slip 

I cannot change this, as it is the output of another code.

How would you suggest to proceed in order to get simple planar ruptures with the information above using the new OQ version?

Setting it manually, as you were suggesting?

Is there a script example I could follow?

Thank you in advance.

Best regards,

Angela Stallone

--
You received this message because you are subscribed to a topic in the Google Groups "OpenQuake Users" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/openquake-users/kOLGwM8f_6g/unsubscribe.
To unsubscribe from this group and all its topics, send an email to openquake-use...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/openquake-users/9d302eb8-b7cc-463b-bde2-e941140acfe4n%40googlegroups.com.

[Peter Pažák]

Hi Angela,

I tried to adopt it a bit to your needs, actually it does very similar things as if you used the point source,

but here in the characteristic source you can pass the planar surface explicitly:

from openquake.hazardlib.geo import Point
from openquake.hazardlib.geo.surface.planar import PlanarSurface
from openquake.hazardlib.source.characteristic import CharacteristicFaultSource
from openquake.hazardlib.mfd import ArbitraryMFD
from openquake.hazardlib.tom import PoissonTOM
from openquake.hazardlib.scalerel import WC1994

from openquake.hazardlib.site import Site, SiteCollection

from openquake.hazardlib.contexts import ContextMaker
from openquake.hazardlib.valid import gsim

import numpy as np

#import 2 GMPEs
gmpes = [gsim('AbrahamsonSilva2008'),gsim('CampbellBozorgnia2008')]

Mag = 6.5
rupture_aratio = 1.5
Strike = 45
Dip = 50
Rake = 0
Hypocenter = Point(9.1500, 45.1833,10)

#according to the magnitude and MSR calculate planar surface
planar_surface = PlanarSurface.from_hypocenter(
        hypoc=Hypocenter,
        msr=WC1994(),

        mag=Mag,
        aratio=rupture_aratio,
        strike=Strike,
        dip=Dip,
        rake=Rake,
        )

# site for which we compute
site_collection = SiteCollection([Site(location=Point(9.200, 45.2,0), vs30=450., vs30measured=True, z1pt0=150., z2pt5=1.5)])

# intensity measure types for which we compute
imtls = {s: [0] for s in ['PGA','SA(0.3)']} #required for context maker, M = 2 IMTs

context_maker = ContextMaker('Active Shallow Crust',gmpes,{'imtls': imtls, 'mags': [Mag]}) #necessary contexts builder

src = CharacteristicFaultSource(source_id = 1,
                                name = 'rup',
                                tectonic_region_type = 'Active Shallow Crust',
                                mfd = ArbitraryMFD([Mag], [0.01]), #this does not have any effect
                                temporal_occurrence_model = PoissonTOM(50.), #this is also not really used
                                surface = planar_surface,
                                rake = Rake)

ctxs = context_maker.from_srcs([src],site_collection) #returns one context from the source for one rupture

#we can calculate GMFs now

gms = context_maker.get_mean_stds(ctxs) #calculate ground motions and stds, returns array of shape (4, G, M, N)

print(gms.shape) #4 values (0=median, 1=std_total, 2=std_intra, 3=std_inter), then G=2 gsims, M=2 IMTs, 1 scenario = magnitude

print(np.exp(gms[0][0][0][0])) #median value, for AbrahamsonSilva2008, PGA, Mw = 5.05

print(gms[1][0][0][0]) #standard deviation, for AbrahamsonSilva2008, PGA, Mw = 5.05

Dátum: pondelok 8. januára 2024, čas: 11:59:53 UTC+1, odosielateľ: angela....@ingv.it

[Peter Pažák]

And if you would like to have even more control over the surface creation, you can also supply length and width:

fault_area = 270
fault_length = 27
fault_width = fault_area/fault_length
rupture_aratio = fault_length / fault_width

from openquake.hazardlib.scalerel.base import BaseMSR
class CustomMSR(BaseMSR):
    def get_median_area(self, mag, rake):
        return np.full_like(mag, fault_area)

Mag = 6.5

Strike = 45
Dip = 50
Rake = 0
Hypocenter = Point(9.1500, 45.1833,10)

planar_surface = PlanarSurface.from_hypocenter(
        hypoc=Hypocenter,
        msr=CustomMSR(),

        mag=Mag,
        aratio=rupture_aratio,
        strike=Strike,
        dip=Dip,
        rake=Rake,
        )

Dátum: štvrtok 11. januára 2024, čas: 18:55:10 UTC+1, odosielateľ: Peter Pažák

[Joel Hirales]

Hi, Angela 

At this moment I am doing an exercise to define the GMPE model according to the 5 Seismotectonic model of Mexico proposed by the GEM study,,, for the case of the Baja California Peninsula..WHAT MODEL IS IT APPROPRIATE??
Thanks..Joel

[Geo]

Hi, 

In this proposed solution for computing point source ground motion prediction, is there a way to introduce damping for spectral accelerations?

Thank you. 

Januka