Python RuptureContext.rake方法代码示例

# 需要导入模块: from openquake.hazardlib.gsim.base import RuptureContext [as 别名]
# 或者: from openquake.hazardlib.gsim.base.RuptureContext import rake [as 别名]
    def test_mag_greater_8pt5(self):
        gmpe = SadighEtAl1997()

        sctx = SitesContext()
        rctx = RuptureContext()
        dctx = DistancesContext()

        rctx.rake =  0.0
        dctx.rrup = numpy.array([0., 1.])
        sctx.vs30 = numpy.array([800., 800.])

        rctx.mag = 9.0
        mean_rock_9, _ = gmpe.get_mean_and_stddevs(
            sctx, rctx, dctx, PGA(), [StdDev.TOTAL]
        )
        rctx.mag = 8.5
        mean_rock_8pt5, _ = gmpe.get_mean_and_stddevs(
            sctx, rctx, dctx, PGA(), [StdDev.TOTAL]
        )
        numpy.testing.assert_allclose(mean_rock_9, mean_rock_8pt5)

        sctx.vs30 = numpy.array([300., 300.])
        rctx.mag = 9.0
        mean_soil_9, _ = gmpe.get_mean_and_stddevs(
            sctx, rctx, dctx, PGA(), [StdDev.TOTAL]
        )
        rctx.mag = 8.5
        mean_soil_8pt5, _ = gmpe.get_mean_and_stddevs(
            sctx, rctx, dctx, PGA(), [StdDev.TOTAL]
        )
        numpy.testing.assert_allclose(mean_soil_9, mean_soil_8pt5)

# 需要导入模块: from openquake.hazardlib.gsim.base import RuptureContext [as 别名]
# 或者: from openquake.hazardlib.gsim.base.RuptureContext import rake [as 别名]
 def _get_event_context(self, idx, nodal_plane_index=1):
     """
     Returns the event contexts for a specific event
     """
     idx = idx[0]
     rctx = RuptureContext()
     rup = self.records[idx]
     setattr(rctx, 'mag', rup.event.magnitude.value)
     if nodal_plane_index == 2:
         setattr(rctx, 'dip',
             rup.event.mechanism.nodal_planes.nodal_plane_2['dip'])
         setattr(rctx, 'rake', 
             rup.event.mechanism.nodal_planes.nodal_plane_2['rake'])
     else:
         setattr(rctx, 'dip',
             rup.event.mechanism.nodal_planes.nodal_plane_1['dip'])
         setattr(rctx, 'rake',
             rup.event.mechanism.nodal_planes.nodal_plane_1['rake'])
     if not rctx.rake:
         rctx.rake = rup.event.mechanism.get_rake_from_mechanism_type()
     if rup.event.rupture:
         setattr(rctx, 'ztor', rup.event.rupture.depth)
         setattr(rctx, 'width', rup.event.rupture.width)
     setattr(rctx, 'hypo_depth', rup.event.depth)
     return rctx

# 需要导入模块: from openquake.hazardlib.gsim.base import RuptureContext [as 别名]
# 或者: from openquake.hazardlib.gsim.base.RuptureContext import rake [as 别名]
 def _get_event_context(self, idx, nodal_plane_index=1):
     """
     Returns the event contexts for a specific event
     """
     idx = idx[0]
     rctx = RuptureContext()
     rup = self.records[idx]
     setattr(rctx, 'mag', rup.event.magnitude.value)
     if nodal_plane_index == 2:
         setattr(rctx, 'strike',
             rup.event.mechanism.nodal_planes.nodal_plane_2['strike'])
         setattr(rctx, 'dip',
             rup.event.mechanism.nodal_planes.nodal_plane_2['dip'])
         setattr(rctx, 'rake',
             rup.event.mechanism.nodal_planes.nodal_plane_2['rake'])
     else:
         setattr(rctx, 'strike', 0.0)
         setattr(rctx, 'dip', 90.0)
         rctx.rake = rup.event.mechanism.get_rake_from_mechanism_type()
     if rup.event.rupture.surface:
         setattr(rctx, 'ztor', rup.event.rupture.surface.get_top_edge_depth())
         setattr(rctx, 'width', rup.event.rupture.surface.width)
         setattr(rctx, 'hypo_loc', rup.event.rupture.surface.get_hypo_location(1000))
     else:
         setattr(rctx, 'ztor', rup.event.depth)
         # Use the PeerMSR to define the area and assuming an aspect ratio
         # of 1 get the width
         setattr(rctx, 'width',
                 np.sqrt(DEFAULT_MSR.get_median_area(rctx.mag, 0)))
         # Default hypocentre location to the middle of the rupture
         setattr(rctx, 'hypo_loc', (0.5, 0.5))
     setattr(rctx, 'hypo_depth', rup.event.depth)
     setattr(rctx, 'hypo_lat', rup.event.latitude)
     setattr(rctx, 'hypo_lon', rup.event.longitude)
     return rctx

# 需要导入模块: from openquake.hazardlib.gsim.base import RuptureContext [as 别名]
# 或者: from openquake.hazardlib.gsim.base.RuptureContext import rake [as 别名]
 def test_get_mean_and_stddevs_good(self):
     """
     Tests the full execution of the GMPE tables for valid data
     """
     gsim = GMPETable(gmpe_table=self.TABLE_FILE)
     rctx = RuptureContext()
     rctx.mag = 6.0
     rctx.rake = 90.0
     dctx = DistancesContext()
     # Test values at the given distances and those outside range
     dctx.rjb = np.array([0.5, 1.0, 10.0, 100.0, 500.0])
     sctx = SitesContext()
     stddevs = [const.StdDev.TOTAL]
     expected_mean = np.array([20.0, 20.0, 10.0, 5.0, 1.0E-19])
     # PGA
     mean, sigma = gsim.get_mean_and_stddevs(sctx, rctx, dctx,
                                             imt_module.PGA(),
                                             stddevs)
     np.testing.assert_array_almost_equal(np.exp(mean), expected_mean, 5)
     np.testing.assert_array_almost_equal(sigma[0], 0.25 * np.ones(5), 5)
     # SA
     mean, sigma = gsim.get_mean_and_stddevs(sctx, rctx, dctx,
                                             imt_module.SA(1.0),
                                             stddevs)
     np.testing.assert_array_almost_equal(np.exp(mean), expected_mean, 5)
     np.testing.assert_array_almost_equal(sigma[0], 0.4 * np.ones(5), 5)

# 需要导入模块: from openquake.hazardlib.gsim.base import RuptureContext [as 别名]
# 或者: from openquake.hazardlib.gsim.base.RuptureContext import rake [as 别名]
 def test_get_amplification_factors(self):
     """
     Tests the amplification tables
     """
     rctx = RuptureContext()
     rctx.rake = 45.0
     rctx.mag = 6.0
     dctx = DistancesContext()
     # Takes distances at the values found in the table (not checking
     # distance interpolation)
     dctx.rjb = np.copy(self.amp_table.distances[:, 0, 0])
     # Test Vs30 is 700.0 m/s midpoint between the 400 m/s and 1000 m/s
     # specified in the table
     sctx = SitesContext()
     stddevs = [const.StdDev.TOTAL]
     expected_mean = np.ones_like(dctx.rjb)
     # Check PGA and PGV
     mean_amp, sigma_amp = self.amp_table.get_amplification_factors(
         imt_module.PGA(), sctx, rctx, dctx.rjb, stddevs)
     np.testing.assert_array_almost_equal(
         mean_amp,
         midpoint(1.0, 1.5) * expected_mean)
     np.testing.assert_array_almost_equal(
         sigma_amp[0],
         0.9 * expected_mean)
     mean_amp, sigma_amp = self.amp_table.get_amplification_factors(
         imt_module.PGV(), sctx, rctx, dctx.rjb, stddevs)
     np.testing.assert_array_almost_equal(
         mean_amp,
         midpoint(1.0, 0.5) * expected_mean)
     np.testing.assert_array_almost_equal(
         sigma_amp[0],
         0.9 * expected_mean)
     # Sa (0.5)
     mean_amp, sigma_amp = self.amp_table.get_amplification_factors(
         imt_module.SA(0.5), sctx, rctx, dctx.rjb, stddevs)
     np.testing.assert_array_almost_equal(
         mean_amp,
         midpoint(1.0, 2.0) * expected_mean)
     np.testing.assert_array_almost_equal(
         sigma_amp[0],
         0.9 * expected_mean)

# 需要导入模块: from openquake.hazardlib.gsim.base import RuptureContext [as 别名]
# 或者: from openquake.hazardlib.gsim.base.RuptureContext import rake [as 别名]
    def check_gmpe_adjustments(self, adj_gmpe_set, original_gmpe):
        """
        Takes a set of three adjusted GMPEs representing the "low", "middle"
        and "high" stress drop adjustments for Germany and compares them
        against the original "target" GMPE for a variety of magnitudes
        and styles of fauling.
        """
        low_gsim, mid_gsim, high_gsim = adj_gmpe_set
        tot_std = [const.StdDev.TOTAL]
        for imt in self.imts:
            for mag in self.mags:
                for rake in self.rakes:
                    rctx = RuptureContext()
                    rctx.mag = mag
                    rctx.rake = rake
                    rctx.hypo_depth = 10.
                    # Get "original" values
                    mean = original_gmpe.get_mean_and_stddevs(self.sctx, rctx,
                                                              self.dctx, imt,
                                                              tot_std)[0]
                    mean = np.exp(mean)
                    # Get "low" adjustments (0.75 times the original)
                    low_mean = low_gsim.get_mean_and_stddevs(self.sctx, rctx,
                                                             self.dctx, imt,
                                                             tot_std)[0]
                    np.testing.assert_array_almost_equal(
                        np.exp(low_mean) / mean, 0.75 * np.ones_like(low_mean))

                    # Get "middle" adjustments (1.25 times the original)
                    mid_mean = mid_gsim.get_mean_and_stddevs(self.sctx, rctx,
                                                             self.dctx, imt,
                                                             tot_std)[0]
                    np.testing.assert_array_almost_equal(
                        np.exp(mid_mean) / mean, 1.25 * np.ones_like(mid_mean))

                    # Get "high" adjustments (1.5 times the original)
                    high_mean = high_gsim.get_mean_and_stddevs(self.sctx, rctx,
                                                               self.dctx, imt,
                                                               tot_std)[0]
                    np.testing.assert_array_almost_equal(
                        np.exp(high_mean) / mean,
                        1.5 * np.ones_like(high_mean))

# 需要导入模块: from openquake.hazardlib.gsim.base import RuptureContext [as 别名]
# 或者: from openquake.hazardlib.gsim.base.RuptureContext import rake [as 别名]
 def test_rhypo_smaller_than_15(self):
     # test the calculation in case of rhypo distances less than 15 km
     # (for rhypo=0 the distance term has a singularity). In this case the
     # method should return values equal to the ones obtained by clipping
     # distances at 15 km.
     sctx = SitesContext()
     sctx.vs30 = numpy.array([800.0, 800.0, 800.0])
     rctx = RuptureContext()
     rctx.mag = 5.0
     rctx.rake = 0
     dctx = DistancesContext()
     dctx.rhypo = numpy.array([0.0, 10.0, 16.0])
     dctx.rhypo.flags.writeable = False
     mean_0, stds_0 = self.GSIM_CLASS().get_mean_and_stddevs(
         sctx, rctx, dctx, PGA(), [StdDev.TOTAL])
     setattr(dctx, 'rhypo', numpy.array([15.0, 15.0, 16.0]))
     mean_15, stds_15 = self.GSIM_CLASS().get_mean_and_stddevs(
         sctx, rctx, dctx, PGA(), [StdDev.TOTAL])
     numpy.testing.assert_array_equal(mean_0, mean_15)
     numpy.testing.assert_array_equal(stds_0, stds_15)

# 需要导入模块: from openquake.hazardlib.gsim.base import RuptureContext [as 别名]
# 或者: from openquake.hazardlib.gsim.base.RuptureContext import rake [as 别名]
 def get_response_spectrum(self, magnitude, distance, periods, rake=90, vs30=800, damping=0.05):
     """
     """
     responses = np.zeros((len(periods),))
     p_damping = damping * 100
     rup = RuptureContext()
     rup.mag = magnitude
     rup.rake = rake
     dists = DistancesContext()
     dists.rjb = np.array([distance])
     sites = SitesContext()
     sites.vs30 = np.array([vs30])
     stddev_types = [StdDev.TOTAL]
     for i, period in enumerate(periods):
         if period == 0:
             imt = _PGA()
         else:
             imt = _SA(period, p_damping)
         responses[i] = np.exp(self._gmpe.get_mean_and_stddevs(sites, rup, dists, imt, stddev_types)[0][0])
     return ResponseSpectrum(periods, responses, unit='g', damping=damping)

# 需要导入模块: from openquake.hazardlib.gsim.base import RuptureContext [as 别名]
# 或者: from openquake.hazardlib.gsim.base.RuptureContext import rake [as 别名]
import matplotlib.pyplot as plt

fig_dir = '/Users/vsahakian/anza/models/statistics/misc/oq_vs_matlab/'


## This all works..... ##

ASK14 = AbrahamsonEtAl2014()

IMT = imt.PGA()
rctx = RuptureContext()
dctx = DistancesContext()
sctx = SitesContext()
sctx_rock = SitesContext()

rctx.rake = 0.0
rctx.dip = 90.0
rctx.ztor = 7.13
rctx.mag = 3.0
#rctx.mag = np.linspace(0.1,5.)
rctx.width = 10.0
rctx.hypo_depth = 8.0

#dctx.rrup = np.logspace(1,np.log10(200),100)
dctx.rrup = np.logspace(np.log10(10),np.log10(10.0),1)


# Assuming average ztor, get rjb:
dctx.rjb = np.sqrt(dctx.rrup**2 - rctx.ztor**2)
dctx.rhypo = dctx.rrup
dctx.rx = dctx.rjb

# 需要导入模块: from openquake.hazardlib.gsim.base import RuptureContext [as 别名]
# 或者: from openquake.hazardlib.gsim.base.RuptureContext import rake [as 别名]
def signal_end(st, event_time, event_lon, event_lat, event_mag,
               method=None, vmin=None, floor=None,
               model=None, epsilon=2.0):
    """
    Estimate end of signal by using a model of the 5-95% significant
    duration, and adding this value to the "signal_split" time. This probably
    only works well when the split is estimated with a p-wave picker since
    the velocity method often ends up with split times that are well before
    signal actually starts.

    Args:
        st (StationStream):
            Stream of data.
        event_time (UTCDateTime):
            Event origin time.
        event_mag (float):
            Event magnitude.
        event_lon (float):
            Event longitude.
        event_lat (float):
            Event latitude.
        method (str):
            Method for estimating signal end time. Either 'velocity'
            or 'model'.
        vmin (float):
            Velocity (km/s) for estimating end of signal. Only used if
            method="velocity".
        floor (float):
            Minimum duration (sec) applied along with vmin.
        model (str):
            Short name of duration model to use. Must be defined in the
            gmprocess/data/modules.yml file.
        epsilon (float):
            Number of standard deviations; if epsilon is 1.0, then the signal
            window duration is the mean Ds + 1 standard deviation. Only used
            for method="model".

    Returns:
        trace with stats dict updated to include a
        stats['processing_parameters']['signal_end'] dictionary.

    """
    # Load openquake stuff if method="model"
    if method == "model":
        mod_file = pkg_resources.resource_filename(
            'gmprocess', os.path.join('data', 'modules.yml'))
        with open(mod_file, 'r') as f:
            mods = yaml.load(f)

        # Import module
        cname, mpath = mods['modules'][model]
        dmodel = getattr(import_module(mpath), cname)()

        # Set some "conservative" inputs (in that they will tend to give
        # larger durations).
        sctx = SitesContext()
        sctx.vs30 = np.array([180.0])
        sctx.z1pt0 = np.array([0.51])
        rctx = RuptureContext()
        rctx.mag = event_mag
        rctx.rake = -90.0
        dur_imt = imt.from_string('RSD595')
        stddev_types = [const.StdDev.INTRA_EVENT]

    for tr in st:
        if not tr.hasParameter('signal_split'):
            continue
        if method == "velocity":
            if vmin is None:
                raise ValueError('Must specify vmin if method is "velocity".')
            if floor is None:
                raise ValueError('Must specify floor if method is "velocity".')
            epi_dist = gps2dist_azimuth(
                lat1=event_lat,
                lon1=event_lon,
                lat2=tr.stats['coordinates']['latitude'],
                lon2=tr.stats['coordinates']['longitude'])[0] / 1000.0
            end_time = event_time + max(floor, epi_dist / vmin)
        elif method == "model":
            if model is None:
                raise ValueError('Must specify model if method is "model".')
            epi_dist = gps2dist_azimuth(
                lat1=event_lat,
                lon1=event_lon,
                lat2=tr.stats['coordinates']['latitude'],
                lon2=tr.stats['coordinates']['longitude'])[0] / 1000.0
            dctx = DistancesContext()
            # Repi >= Rrup, so substitution here should be conservative
            # (leading to larger durations).
            dctx.rrup = np.array([epi_dist])
            lnmu, lnstd = dmodel.get_mean_and_stddevs(
                sctx, rctx, dctx, dur_imt, stddev_types)
            duration = np.exp(lnmu + epsilon * lnstd[0])
            # Get split time
            split_time = tr.getParameter('signal_split')['split_time']
            end_time = split_time + float(duration)
        else:
            raise ValueError('method must be either "velocity" or "model".')
        # Update trace params
        end_params = {
#.........这里部分代码省略.........