Python script to compute GMPE at a particular site

[hadi]

Hi OQ team,

I am trying to calculate a gmpe's mean and std values for a particular rupture plane at a particular site! I had done it before using older versions of OQ hazardlib, following the examples provided here: 

https://github.com/GEMScienceTools/notebooks/blob/master/hazardlib/2%20-%20GMPEs.ipynb

but using the latest version, I am getting the below error. I have attached my python script to this message.

I was wondering if you have a similar example but comatible with the latest version of OQ, if not if you can help me to fix the error in my python script.

Thank You!

Error message:

--------------------------, if so that would be really appreciated!-------------------------------------------------
TypeError                                 Traceback (most recent call last)
File ~/Projects/NSHA22_gmpe/scripts/2018_GMPE_selection_report/test_gmpe_mag_scaling.py:50, in <module>
     48 means = []
     49 mags = []
---> 50 context_maker = ContextMaker(gmpes)
     51 for rupture in src.iter_ruptures():
     52     mags.append(rupture.mag)

TypeError: __init__() missing 2 required positional arguments: 'gsims' and 'oq'

[Peter Pažák]

Hi,

I think the prefered way would be to use SMTK for plotting via smtk.trellis.trellis_plots.MagnitudeIMTTrellis as the code you have uses too much internals which do change a lot.
The "new way" might be something like this:

from openquake.hazardlib.gsim import get_available_gsims #building this takes time
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
import numpy as np

gs = get_available_gsims()
gmpes = [gs['AbrahamsonSilva2008'](), gs['ChiouYoungs2008'](), gs['CampbellBozorgnia2008']()] #avoids importing classes, G = 3

# explore magnitude scaling, by defining a Point source and calculating median ground shaking at the point source location
location = Point(9.1500, 45.1833)
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
)
ruptures = [r for r in src.iter_ruptures()] #create list of ruptures
mags = [r.mag for r in ruptures] #magnitude list - from 5.05 to 6.45 step 0.1 (15)

# this is the site for which we compute the median ground shaking
site_collection = SiteCollection([Site(location=location, vs30=760., vs30measured=True, z1pt0=40., z2pt5=1.0)])

# these are the intensity measure type for which we compute the median ground shaking
imtls = {s: [0] for s in ['PGA','SA(0.3)']} #required for context maker, M = 2 IMTs

context_maker = ContextMaker('*',gmpes,{'imtls': imtls}) #necessary contexts builder
ctxs = context_maker.get_ctxs(ruptures,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) #first 4 are median, std_total, std_intra, std_inter, then G=3 gsims, M=2 IMTs, 15 scenarios = magnitudes

print(np.exp(gms[0][2][0])) #median values, AbrahamsonSilva2008, PGA

outputs

(4, 3, 2, 15)
[0.12275133 0.13594099 0.15066125 0.16710007 0.18714439 0.20575711 0.22021613 0.23511299 0.24979486 0.26322846 0.28464689 0.31933286 0.36152975 0.40846104 0.45970467]

being medians corresponding to the mags at the site.

Would be grateful if someone can check this.

Peter

Dátum: streda 18. mája 2022, čas: 10:49:07 UTC+2, odosielateľ: hadi

[Peter Pažák]

Sorry  gms[0][2][0] is CampbellBozorgnia, of course...

Dátum: piatok 20. mája 2022, čas: 11:09:41 UTC+2, odosielateľ: Peter Pažák

[hadi]

Hi Peter,

Thanks a lot for providing the sample script! A follow-up question:

I tried to run the script for NGAEastGMPE by defining 'gmpes' in your example as:

gmpes = [gs['NGAEastGMPE'](gmpe_table='NGAEast_BOORE_A04_J15.hdf5')]

but got the below error (have also attached the script to this message)! any idea how to properly call/use NGAEastGMPE in the script? Thank You!

Error:

KeyError                                  Traceback (most recent call last)
File ~/Projects/NSHA22_gmpe/scripts/oq_resp_call_gmpe.py:45, in <module>
     43 context_maker = ContextMaker('*',gmpes,{'imtls': imtls}) #necessary contexts builder
     44 ctxs = context_maker.get_ctxs(ruptures,site_collection) #returns rupture contexts
---> 45 gms = context_maker.get_mean_stds(ctxs) #calculate ground motions and stds, returns array of shape (4, G, M, N)
     46 print(gms.shape) #first 4 are median, std_total, std_intra, std_inter, then G=3 gsims, M=2 IMTs, 15 scenarios = magnitudes
     48 print(np.exp(gms[0][2][0]))

File ~/miniconda/lib/python3.8/site-packages/openquake.engine-3.14.0-py3.8.egg/openquake/hazardlib/contexts.py:829, in ContextMaker.get_mean_stds(self, ctxs)
    827 for ctx in recarrays:
    828     slc = slice(start, start + len(ctx))
--> 829     adj = compute(gsim, ctx, self.imts, *out[:, g, :, slc])
    830     if adj is not None:
    831         self.adj[gsim].append(adj)

File ~/miniconda/lib/python3.8/site-packages/openquake.engine-3.14.0-py3.8.egg/openquake/hazardlib/gsim/nga_east.py:720, in NGAEastGMPE.compute(self, ctx, imts, mean, sig, tau, phi)
    718 [mag] = np.unique(np.round(ctx.mag, 6))  # by construction
    719 for m, imt in enumerate(imts):
--> 720     mean[m], _, _ = get_mean_amp(self, mag, ctx, imt)
    721     # Get standard deviation model
    722     sig[m], tau[m], phi[m] = get_stddevs(self, ctx.mag, imt)

File ~/miniconda/lib/python3.8/site-packages/openquake.engine-3.14.0-py3.8.egg/openquake/hazardlib/gsim/nga_east.py:557, in get_mean_amp(self, mag, ctx, imt)
    555 else:
    556     rock_imt = SA(0.01)
--> 557 pga_r = get_hard_rock_mean(self, mag, ctx, rock_imt)
    559 # Get the desired spectral acceleration on rock
    560 if imt.string != "PGA":
    561     # Calculate the ground motion at required spectral period for
    562     # the reference rock

File ~/miniconda/lib/python3.8/site-packages/openquake.engine-3.14.0-py3.8.egg/openquake/hazardlib/gsim/nga_east.py:459, in get_hard_rock_mean(self, mag, ctx, imt)
    454 """
    455 Returns the mean and standard deviations for the reference very hard
    456 rock condition (Vs30 = 3000 m/s)
    457 """
    458 # return Distance Tables
--> 459 imls = self.mean_table['%.2f' % mag, imt.string]
    460 # Get distance vector for the given magnitude
    461 idx = np.searchsorted(self.m_w, mag)

KeyError: ('5.05', 'PGA')

[Peter Pažák]

Hi, I do have OQ v3.13 and I can confirm for me your script works without a problem (just modified to print(np.exp(gms[0][0][0]))).

I thought it is because interpolation is not implemented, but even with magnitude 5.05 it returns numbers,

could body with 3.14 or 3.13 try too?

Thank you

Peter

Dátum: utorok 24. mája 2022, čas: 5:42:11 UTC+2, odosielateľ: hadi

[Michele Simionato]

On Friday, May 20, 2022 at 11:09:41 AM UTC+2 peter...@gmail.com wrote:

gs = get_available_gsims()
gmpes = [gs['AbrahamsonSilva2008'](), gs['ChiouYoungs2008'](), gs['CampbellBozorgnia2008']()] #avoids importing classes, G = 3

This is not the recommended way, Peter.  get_available_gsims() is just to get the list of valid names, but the way to instantiate a gsim given a string is via the validator:

>>> from openquake.hazardlib import valid

>>> valid.gsim('AbrahamsonSilva2008')

 [AbrahamsonSilva2008]

This also works in the case the GMPE has parameters, as in this example:

>>> valid.gsim("""    

                    [ModifiableGMPE]
                    gmpe.Campbell2003 = {}
                    add_between_within_stds.with_betw_ratio = 0.6

      """)

[ModifiableGMPE]
gmpe.Campbell2003 = {}
add_between_within_stds.with_betw_ratio = 0.6

HTH,

        Michele

[Michele Simionato]

On Tuesday, May 24, 2022 at 5:42:11 AM UTC+2 hadi wrote:

    458 # return Distance Tables
--> 459 imls = self.mean_table['%.2f' % mag, imt.string]
    460 # Get distance vector for the given magnitude
    461 idx = np.searchsorted(self.m_w, mag)

KeyError: ('5.05', 'PGA')

GMPETables are special: in order to build the tables by magnitude and IMT you must pass the magnitudes are 2-digit strings,

something like this:

context_maker = ContextMaker('*',gmpes,{'imtls': imtls, 'mags': ['%.2f' % mag for mag in mags]}) 

    Michele Simionato

 

[Peter Pažák]

Thank you Michele, it took me some minutes to figure out how to apply your suggestion for GMPE table:

from openquake.hazardlib.valid import gsim

gmpes = [gsim('AbrahamsonSilva2008'),gsim("""[NGAEastGMPE]
gmpe_table = 'NGAEast_BOORE_A04_J15.hdf5'""")]

that the new line is really required after each key and value/parameter..

Specifying magnitudes was not necessary in v3.13 it seems, but is now required for the tables, right?

context_maker = ContextMaker('*',gmpes,{'imtls': imtls, 'mags': ['%1.2f'%m for m in mags]}) #necessary contexts builder

Thank you

Peter

Dátum: štvrtok 26. mája 2022, čas: 8:44:34 UTC+2, odosielateľ: michele....@globalquakemodel.org

[Michele Simionato]

On Saturday, May 28, 2022 at 7:05:49 AM UTC+2 peter...@gmail.com wrote:

Thank you Michele, it took me some minutes to figure out how to apply your suggestion for GMPE table:

from openquake.hazardlib.valid import gsim

gmpes = [gsim('AbrahamsonSilva2008'),gsim("""[NGAEastGMPE]
gmpe_table = 'NGAEast_BOORE_A04_J15.hdf5'""")]

that the new line is really required after each key and value/parameter..

Yes and no, you can do

gmpes = [gsim('AbrahamsonSilva2008'),gsim("[NGAEastGMPE]\ngmpe_table = 'NGAEast_BOORE_A04_J15.hdf5'")]

Specifying magnitudes was not necessary in v3.13 it seems, but is now required for the tables, right?

Yes, but on the plus side, now GMPETables are 5+ times faster than before ;-)

             Michele 

[Indrayoga Putra]

Hi peter, sory for interupt your time. I want to ask something. I try the last code you made with  AbrahamsonSilva2008  gmpe but the code just return empty array in ctxs variable and gms variable (just like in the empty_array.jpg  file i attach). Do you know how to fix this problem?    

Thank you Peter. 

[Peter Pažák]

Hi, you probably are aware that using the internals of OQ can be tricky and things are changing frequently.

Please let us know which version of OQ you are using.

I found out that for 3.16, that I am currently using, I have to create ctxs differently:

#ctxs = context_maker.get_ctxs(ruptures,site_collection) #returns rupture contexts - older version
ctxs = context_maker.from_srcs([src],site_collection) #returns rupture contexts

Just changing this line produces reasonable result I believe:

ctxs: [rec.array([(50., 5.05, 0.,  3.64937396, 8.60130195, 760.,  True, 40., 0., 18.63779027, 1.17303525, 0,  0, 0, 0., 2.60233254e-06, [], 9.1605874 , 45.1758361 ),
           (50., 5.15, 0.,  4.04779365, 8.44859934, 760.,  True, 40., 0., 18.49385182, 1.30102707, 0,  1, 0, 0., 2.07664754e-06, [], 9.16174311, 45.17502117),
           (50., 5.25, 0.,  4.48971072, 8.27922546, 760.,  True, 40., 0., 18.3353988 , 1.44297532, 0,  2, 0, 0., 1.65715370e-06, [], 9.16302496, 45.17411725),
           (50., 5.35, 0.,  4.97987399, 8.09136022, 760.,  True, 40., 0., 18.16116178, 1.60039975, 0,  3, 0, 0., 1.32239985e-06, [], 9.1644467 , 45.17311464),
           (50., 5.45, 0.,  5.52359261, 7.76596965, 760.,  True, 40., 0., 18.07695424, 1.87317249, 0,  4, 0, 0., 1.05526806e-06, [], 9.16513792, 45.17262694),
           (50., 5.55, 0.,  6.12672149, 7.30372017, 760.,  True, 40., 0., 18.07681708, 2.26072301, 0,  5, 0, 0., 8.42098302e-07, [], 9.16513759, 45.17262671),
           (50., 5.65, 0.,  6.79571794, 6.79100469, 760.,  True, 40., 0., 18.07664833, 2.69054515, 0,  6, 0, 0., 6.71989971e-07, [], 9.16513723, 45.17262645),
           (50., 5.75, 0.,  7.5377779 , 6.22231359, 760.,  True, 40., 0., 18.07644072, 3.16724486, 0,  7, 0, 0., 5.36244427e-07, [], 9.16513683, 45.17262616),
           (50., 5.85, 0.,  8.36088404, 5.59153574, 760.,  True, 40., 0., 18.07618531, 3.6959289 , 0,  8, 0, 0., 4.27920204e-07, [], 9.16513638, 45.17262585),
           (50., 5.95, 0.,  9.27389188, 4.89189284, 760.,  True, 40., 0., 18.07587109, 4.28225897, 0,  9, 0, 0., 3.41478049e-07, [], 9.16513588, 45.17262549),
           (50., 6.05, 0., 10.28662546, 4.11586655, 760.,  True, 40., 0., 18.07548451, 4.93251173, 0, 10, 0, 0., 2.72497669e-07, [], 9.16513533, 45.1726251 ),
           (50., 6.15, 0., 11.4099834 , 3.25511774, 760.,  True, 40., 0., 18.07500892, 5.6536452 , 0, 11, 0, 0., 2.17451692e-07, [], 9.16513472, 45.17262467),
           (50., 6.25, 0., 12.65605676, 2.30039684, 760.,  True, 40., 0., 18.07442382, 6.45337228, 0, 12, 0, 0., 1.73525295e-07, [], 9.16513404, 45.17262418),
           (50., 6.35, 0., 13.04813997, 2.        , 760.,  True, 40., 0., 18.07309838, 6.70287344, 0, 13, 0, 0., 1.38472264e-07, [], 9.16513383, 45.17262403),
           (50., 6.45, 0., 13.04813655, 2.        , 760.,  True, 40., 0., 18.07081224, 6.69878565, 0, 14, 0, 0., 1.10500131e-07, [], 9.16513384, 45.17262402)],
          dtype=[('dip', '<f8'), ('mag', '<f8'), ('rake', '<f8'), ('width', '<f8'), ('ztor', '<f8'), ('vs30', '<f8'), ('vs30measured', '?'), ('z1pt0', '<f8'), ('rjb', '<f8'), ('rrup', '<f8'), ('rx', '<f8'), ('src_id', '<u4'), ('rup_id', '<u4'), ('sids', '<u4'), ('weight', '<f8'), ('occurrence_rate', '<f8'), ('probs_occur', '<f8', (0,)), ('clon', '<f8'), ('clat', '<f8')])]
(4, 1, 2, 15)
[0.07363788 0.07948661 0.08566536 0.09216913 0.09724708 0.10069594
 0.10372123 0.10622624 0.1081108  0.10927434 0.11308793 0.11989745
 0.12615699 0.14054073 0.16057652]

Dátum: pondelok 13. februára 2023, čas: 3:26:21 UTC+1, odosielateľ: indrayoga...@gmail.com