response spectrum

[song]

I use the following py code previously to calculate the response spectrum.  However, due to since of the  SitesContext. I don't think it is working no, what's best way to calcualte response spectrum from a GMM directly using latest version?

import numpy as np
from matplotlib import pyplot

# Import the needed objects ( distance, sites, rupture, coeffs table) to store inputs for the target GMPE.
from openquake.hazardlib.gsim.base import DistancesContext, SitesContext, RuptureContext, CoeffsTable

# Import the constant container which contains  
# (1) Tectonic Region Types (2) Ground motion components (3) Ground motion variability types
from openquake.hazardlib import const

# Import the target ground motion intensity types
from openquake.hazardlib.imt import PGA, PGV, SA


def storage():
    """
    Create the variables to store the information needed for the GMPE
    return:
    sctx, for store the site condition information
    rctx, for store the rupture information
    dctx, for store the distance information
    """
    # (1) sites (site information)
    sctx = SitesContext()
    # (2) rup (rupture property)
    rctx = RuptureContext()
    # (3) dists (distance matrix)
    dctx = DistancesContext()
    return sctx, rctx, dctx

def gmpe_info(GMPEName):
    '''
    See what are the required site, rupture, distance parameters and
    the supported tectonic region, intensity measure types, ground
    component, and standard deviations for the GMPE
   
    Input: Name of the GMPE
    '''
    GMPE = GMPEName()
    print('required site parameters: %s' % GMPE.REQUIRES_SITES_PARAMETERS)
    print('required rupture parameters: %s' % GMPE.REQUIRES_RUPTURE_PARAMETERS)
    print('required distance parameters: %s' % GMPE.REQUIRES_DISTANCES)
    print('supported tectonic region: %s'% GMPE.DEFINED_FOR_TECTONIC_REGION_TYPE)
    print('supported intensity measure types: %s' % ', '.join([imt.__name__ for imt in GMPE.DEFINED_FOR_INTENSITY_MEASURE_TYPES]))
    print('supported component: %s' % GMPE.DEFINED_FOR_INTENSITY_MEASURE_COMPONENT)
    print('supported standard deviations: %s' % ', '.join([std for std in GMPE.DEFINED_FOR_STANDARD_DEVIATION_TYPES]))


def computeGMPE(GMPEName, site_info, rup_info, dist_info, imts):
    '''
    Use a GMPE for computing the predicted median ground motion and total standard deviation
    for the target scenario
   
    Input:
    GMPEName, the target GMPE method, e.g.,AkkarEtAlRjb2014
    site_info, the site information
    rup_info, the rupture information
    dist_info, the distance information
    imts, the target ground motion intensity types

    Output:
    median ground motion prediction, and total standard deviation
    '''
    GMPE = GMPEName()
    stddevs = [const.StdDev.TOTAL]
    median = GMPE.get_mean_and_stddevs(site_info, rup_info, dist_info, imts,stddevs)[0]
    total_std = GMPE.get_mean_and_stddevs(site_info, rup_info, dist_info, imts,stddevs)[1][0]
    return np.exp(median), total_std


def response_spectrum(GMPEName, site_info, rup_info, dist_info):
    '''
    Use a GMPE for computing the response spectrum for the target scenario
   
    Input:
    GMPEName, the target GMPE method, e.g.,AkkarEtAlRjb2014
    site_info, the site information
    rup_info, the rupture information
    dist_info, the distance information

    Output:
    1. Figure of response spectrum
    2. Txt file of response spectrum contains 3 colomns which are period,
       median spectrum acceleration and total standard deviations.
    '''

    GMPE = GMPEName()
    periods = []
    median_value = []
    total_std_value = []
    median_plus_sigma = []
    median_minus_sigma = []

    for i in GMPE.COEFFS.sa_coeffs.keys():
        periods.append(i.period)
    periods = np.sort(periods)


    for t in periods:
        imt = SA(t)
        median, total_std = computeGMPE(GMPEName, site_info, rup_info, dist_info, imt)
        median_plus_sigma.append(np.exp(np.log(median) + total_std))
        median_minus_sigma.append(np.exp(np.log(median) - total_std))
        median_value.append(median)
        total_std_value.append(total_std)
    # Visualise the result
    pyplot.figure(figsize=(8,5))
    pyplot.plot(periods, median_value, linewidth=4)
    pyplot.xlabel('Period (sec)', fontsize=20)
    pyplot.ylabel('Spectral Acceleration (g)', fontsize=20)
    pyplot.title('Acceleration Response Spectrum', fontsize=24)
    pyplot.tick_params(labelsize=14)
    pyplot.fill_between((periods.flatten()),
                        np.array(median_minus_sigma).flatten(),  
                        np.array(median_plus_sigma).flatten(),
                        color="#3F5D7D", alpha=0.3)
    # Output the plot
    pyplot.savefig('response_spectrum.pdf')

    print("The plot files generated: ./response_spectrum.pdf")
    return(periods.flatten(), np.array(median_value).flatten(), np.array(total_std_value).flatten(), np.array(median_plus_sigma).flatten())

[Michele Simionato]

See https://docs.openquake.org/oq-engine/master/manual/contributing/developing-with-the-engine.html#working-with-gmpes-directly-the-contextmaker