Calculate PGA/PGV values using CampbellBozorgnia2014 class in Python

[amer...@gmail.com]

Hey folks,

So I have copied the source code for CampbellBozorgnia2014 which is shown HERE in my python IDE, but I do not know how I can calculate the PGA/PGV values for a certain site under a specific simulated earthquake! Usually there should be a command in which you define the input variables (magnitude, rupture distance, imt, etc.), but I cannot find it in the openquake documentation. For example, in pygmm we call: pygmm.CampbellBozorgnia2014(mag=6.9 ,dist_jb=5, dist_x=5, dist_rup=1, dip=90, v_s30=600, mechanism  = "SS") and this will give us the PGA/PGV values based on the earthquake magnitude, distance, shear velocity and fault rupture mechanism. Does anyone know how this works in openquake? 

Thanks! 

[Peter Pažák]

Hi,

if you have OQ installation and can use the OQ python, you could do something like this to get means:

(I just had a look into the code how things are used, but would be glad if a developer could review my amateur code whether I am using everything

correctly below, this works in OQ 3.3, slight problem is that it may not work in other versions if the objects/structure has changed)

from openquake.hazardlib.gsim.campbell_bozorgnia_2014 import CampbellBozorgnia2014
from openquake.hazardlib import (geo, site, source, imt, gsim, const)
import numpy as np

#which GMPE to use - see https://docs.openquake.org/oq-hazardlib/0.14/gsim/index.html#built-in-gsims
gmpei = CampbellBozorgnia2014()

#prepare sites,  and  used by CB2014,  does not have impact on results
nsites = 5
vs30 = [180,200,220,300,500] #vs30 values
z1pt0 = [60,80,100,120,150] #z1pt0 [m] - unused by CB2014
z2pt5 = [0.05,0.1,0.15,0.2,0.25] #z2pt5 [km] values
dist = [0.1,0.2,0.3,0.4,0.5] #in degrees
sites = []
for i in range(nsites):
  sites.append(site.Site(location=geo.Point(dist[i],0,0.0), vs30=vs30[i], vs30measured=False, z1pt0=z1pt0[i], z2pt5=z2pt5[i], backarc=True))

#prepare rupture - very simple planar
surface = geo.surface.planar.PlanarSurface.from_corner_points(geo.Point(0, -0.5, 5.0), geo.Point(0, 0.5, 5.0), geo.Point(0, 0.5, 15.0), geo.Point(0, -0.5, 15.0))
hypocenter = geo.Point(0, 0, 10.0) #depth 10km
#print('hypocenter: ',hypocenter)
print('strike, dip: ',surface.get_strike(), surface.get_dip())
source_typology = None #object() #does not need to be specified, all is manually defined here
magnitude = 8.0
rake = -90.0 #Rakes angles within 30 of horizontal are strike-slip, angles from 30 to 150 are reverse, and angles from -30 to -150 are normal.
rupture = source.rupture.BaseRupture(magnitude,rake,const.TRT.SUBDUCTION_INTERFACE,hypocenter,surface,source_typology)

#create contexts
sitectx, distctx = gsim.base.ContextMaker([gmpei]).make_contexts(site.SiteCollection(sites), rupture) #for defined sites and rupture
print(sitectx.vs30) #print vs30 for all sites
print(distctx.rjb) #original calculated distances
#distctx.rrup = np.linspace(0,200,nsites) #you can redefine distances to nice kilometers
stddev_types = set([const.StdDev.TOTAL]) #only total standard deviations will be calculated, StdDev not used here

acc = [] #calculated accelerations
#what intensity measure types to calculate
int_meas_type = [imt.PGA(), imt.PGV(),imt.SA(0.05,5.0),imt.SA(0.1,5.0),imt.SA(0.2,5.0),imt.SA(0.5,5.0),imt.SA(1.0,5.0),imt.SA(5.0,5.0)

for i in int_meas_type:
  [means,stddevs] = gmpei.get_mean_and_stddevs(sitectx, rupture, distctx, i, stddev_types) #use the selected GMPE to calculate means and StdDevs
  acc.append(gmpei.to_imt_unit_values(means)) #convert logarithms back and store means to array

print(acc) #for each of the 5 sites contains values of PGA, PGV, SA(0.05), SA(0.1), ... , SA(5.0)

[amer...@gmail.com]

Hi Peter,

Thank you very much for your response. Your code works like a charm. I have two questions though:

1) Is there a document/user manual that clearly explains the code you have written in details? Like I mean from where/how did you learn to codify these? It is pretty easy to understand your code, but I was just wondering if there is a source I can use? 

2) Is there a way that instead of longitude/latitude we use x/y values? I am working on a virtual community and there is nothing defined as longitude/latitude! Also, why are you using the lat/long distances in "degrees", but then for hypocenter the inputs are in "km"? This is another reason why I asked if there is a document that clearly explains these as I am not sure about the units.

Thank you very much, I really appreciate it.

Best,

Reza

[Peter Pažák]

Hi, as far as I know there is just user manual/documentation on how to use OpenQuake Engine, but these python internals do not have separate documentation (apart from the comments in code) - maybe developers can correct me, I would like to know too :-)

Peter

[Reza Ameri]

Thanks Peter, appreciate it!

Best,

Reza

-- 

Reza Ameri
Ph.D. Candidate in Structural Engineering
Center for Risk-Based Community Resilience Planning
Department of Civil and Environmental Engineering
Colorado State University
Fort Collins, CO 80523; USA
Tel: (970) 290-8134
Email: amer...@engr.colostate.edu

On Apr 10, 2019, at 9:34 AM, Peter Pažák <peter...@gmail.com> wrote:

Hi, as far as I know there is just user manual/documentation on how to use OpenQuake Engine, but these python internals do not have separate documentation (apart from the comments in code) - maybe developers can correct me, I would like to know too :-)

Peter

--
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/MO_VAfYsoI0/unsubscribe.
To unsubscribe from this group and all its topics, send an email to openquake-use...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Marco Pagani]

Hi Peter, thanks a lot for your help. You are right at the moment we do
not have specific documentation for the code like the one you shared.

> --
> 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