Compute and return top of rupture depth term. See paragraph
        'Depth-to-Top of Rupture Model', page 1042.
323165
def _get_top_of_rupture_depth_term(self, C, imt, rup):
        
        if rup.ztor >= 20.0:
            return C['a15']
        else:
            return C['a15'] * rup.ztor / 20.0
208879
def get_ridx(self, iloc):
        """"""
        with h5py.File(self.source_file, "r") as hdf5:
            return hdf5[self.idx_set["geol"] + "/RuptureIndex"][iloc]
1050452
def _get_style_of_faulting_term(self, C, rake):
        
        if rake > -150.0 and rake <= -30.0:
            return C['fN']
        elif rake > 30.0 and rake <= 150.0:
            return C['fR']
        else:
            return C['fSS']
1264084
def build_rupture_node(rupt, probs_occur):
    
    s = sum(probs_occur)
    if abs(s - 1) > pmf.PRECISION:
        raise ValueError('The sum of %s is not 1: %s' % (probs_occur, s))
    h = rupt.hypocenter
    hp_dict = dict(lon=h.longitude, lat=h.latitude, depth=h.depth)
    rupt_nodes = [Node('magnitude', {}, rupt.mag),
                  Node('rake', {}, rupt.rake),
                  Node('hypocenter', hp_dict)]
    rupt_nodes.extend(rupt.surface.surface_nodes)
    geom = rupt.surface.surface_nodes[0].tag
    if len(rupt.surface.surface_nodes) > 1:
        name = 'multiPlanesRupture'
    elif geom == 'planarSurface':
        name = 'singlePlaneRupture'
    elif geom == 'simpleFaultGeometry':
        name = 'simpleFaultRupture'
    elif geom == 'complexFaultGeometry':
        name = 'complexFaultRupture'
    elif geom == 'griddedSurface':
        name = 'griddedRupture'
    return Node(name, {'probs_occur': probs_occur}, nodes=rupt_nodes)
1338102
def _get_soil_depth_term(self, C, z1pt0, vs30):
        
        # Get reference z1pt0
        z1ref = self._get_z1pt0ref(vs30)
        # Get z1pt0
        z10 = copy.deepcopy(z1pt0)
        # This is used for the calculation of the motion on reference rock
        idx = z1pt0 < 0
        z10[idx] = z1ref[idx]
        factor = np.log((z10 + 0.01) / (z1ref + 0.01))
        # Here we use a linear interpolation as suggested in the 'Application
        # guidelines' at page 1044
        # Above 700 m/s the trend is flat, but we extend the Vs30 range to
        # 6,000 m/s (basically the upper limit for mantle shear wave velocity
        # on earth) to allow extrapolation without throwing an error.
        f2 = interpolate.interp1d(
            [0.0, 150, 250, 400, 700, 1000, 6000],
            [C['a43'], C['a43'], C['a44'], C['a45'], C['a46'], C['a46'],
             C['a46']],
            kind='linear')
        return f2(vs30) * factor
832084
def bottom(self):
        """"""
        if self._has_real():
            return self._data.real_bottom
        return self._data.bottom
335082
def _compute_term2(self, C, mag, rrup):
        
        c78_factor = (C['c7'] * np.exp(C['c8'] * mag)) ** 2
        R = np.sqrt(rrup ** 2 + c78_factor)

        return C['c4'] * np.log(R) + (C['c5'] + C['c6'] * mag) * rrup
159423
def get_magnitude_term(self, C, rup):
        
        b0, stress_drop = self._get_sof_terms(C, rup.rake)
        if rup.mag <= C["m1"]:
            return b0
        else:
            # Calculate moment (equation 5)
            m_0 = 10.0 ** (1.5 * rup.mag + 16.05)
            # Get stress-drop scaling (equation 6)
            if rup.mag > C["m2"]:
                stress_drop += (C["b2"] * (C["m2"] - self.CONSTANTS["mstar"]) +
                                (C["b3"] * (rup.mag - C["m2"])))
            else:
                stress_drop += (C["b2"] * (rup.mag - self.CONSTANTS["mstar"]))
            stress_drop = np.exp(stress_drop)
            # Get corner frequency (equation 4)
            f0 = 4.9 * 1.0E6 * 3.2 * ((stress_drop / m_0) ** (1. / 3.))
            return 1. / f0
1146634
def _get_distance_term(self, C, rrup, backarc):
        
        # Geometric attenuation function
        distance_scale = -np.log10(np.sqrt(rrup ** 2 + 3600.0))
        # Anelastic attenuation in the backarc
        distance_scale[backarc] += (C["c2"] * rrup[backarc])
        # Anelastic Attenuation in the forearc
        idx = np.logical_not(backarc)
        distance_scale[idx] += (C["c1"] * rrup[idx])
        return distance_scale
567807
def get_ruptures(self, srcfilter=calc.filters.nofilter):
        
        ebrs = []
        with datastore.read(self.filename) as dstore:
            rupgeoms = dstore['rupgeoms']
            for rec in self.rup_array:
                if srcfilter.integration_distance:
                    sids = srcfilter.close_sids(rec, self.trt, rec['mag'])
                    if len(sids) == 0:  # the rupture is far away
                        continue
                else:
                    sids = None
                mesh = numpy.zeros((3, rec['sy'], rec['sz']), F32)
                geom = rupgeoms[rec['gidx1']:rec['gidx2']].reshape(
                    rec['sy'], rec['sz'])
                mesh[0] = geom['lon']
                mesh[1] = geom['lat']
                mesh[2] = geom['depth']
                rupture_cls, surface_cls = self.code2cls[rec['code']]
                rupture = object.__new__(rupture_cls)
                rupture.serial = rec['serial']
                rupture.surface = object.__new__(surface_cls)
                rupture.mag = rec['mag']
                rupture.rake = rec['rake']
                rupture.hypocenter = geo.Point(*rec['hypo'])
                rupture.occurrence_rate = rec['occurrence_rate']
                rupture.tectonic_region_type = self.trt
                if surface_cls is geo.PlanarSurface:
                    rupture.surface = geo.PlanarSurface.from_array(
                        mesh[:, 0, :])
                elif surface_cls is geo.MultiSurface:
                    # mesh has shape (3, n, 4)
                    rupture.surface.__init__([
                        geo.PlanarSurface.from_array(mesh[:, i, :])
                        for i in range(mesh.shape[1])])
                elif surface_cls is geo.GriddedSurface:
                    # fault surface, strike and dip will be computed
                    rupture.surface.strike = rupture.surface.dip = None
                    rupture.surface.mesh = Mesh(*mesh)
                else:
                    # fault surface, strike and dip will be computed
                    rupture.surface.strike = rupture.surface.dip = None
                    rupture.surface.__init__(RectangularMesh(*mesh))
                grp_id = rec['grp_id']
                ebr = EBRupture(rupture, rec['srcidx'], grp_id,
                                rec['n_occ'], self.samples)
                # not implemented: rupture_slip_direction
                ebr.sids = sids
                ebrs.append(ebr)
        return ebrs
528794
def get_ruptures_within(dstore, bbox):
    
    minlon, minlat, maxlon, maxlat = map(float, bbox.split(','))
    hypo = dstore['ruptures']['hypo'].T  # shape (3, N)
    mask = ((minlon <= hypo[0]) * (minlat <= hypo[1]) *
            (maxlon >= hypo[0]) * (maxlat >= hypo[1]))
    return dstore['ruptures'][mask]
741976
def _get_basin_depth_term(self, C, sites, period):
        
        return np.zeros(len(sites.vs30), dtype=float)
948212
def get_top_edge_depth(self):
        
        areas = self._get_areas()
        depths = numpy.array(
            [surf.get_top_edge_depth() for surf in self.surfaces])
        return numpy.sum(areas * depths) / numpy.sum(areas)
1103822
def convert_simpleFaultRupture(self, node):
        
        mag, rake, hypocenter = self.get_mag_rake_hypo(node)
        with context(self.fname, node):
            surfaces = [node.simpleFaultGeometry]
        rupt = source.rupture.BaseRupture(
            mag=mag, rake=rake, tectonic_region_type=None,
            hypocenter=hypocenter,
            surface=self.convert_surfaces(surfaces))
        return rupt
150664
def _get_mean_deep_soil(self, mag, rake, rrup, is_reverse, imt):
        
        if mag <= self.NEAR_FIELD_SATURATION_MAG:
            c4 = self.COEFFS_SOIL_IMT_INDEPENDENT['c4lowmag']
            c5 = self.COEFFS_SOIL_IMT_INDEPENDENT['c5lowmag']
        else:
            c4 = self.COEFFS_SOIL_IMT_INDEPENDENT['c4himag']
            c5 = self.COEFFS_SOIL_IMT_INDEPENDENT['c5himag']
        c2 = self.COEFFS_SOIL_IMT_INDEPENDENT['c2']
        c3 = self.COEFFS_SOIL_IMT_INDEPENDENT['c3']
        C = self.COEFFS_SOIL[imt]
        if is_reverse:
            c1 = self.COEFFS_SOIL_IMT_INDEPENDENT['c1r']
            c6 = C['c6r']
        else:
            c1 = self.COEFFS_SOIL_IMT_INDEPENDENT['c1ss']
            c6 = C['c6ss']
        # clip mag if greater than 8.5. This is to avoid
        # ValueError: negative number cannot be raised to a fractional power
        mag = 8.5 if mag > 8.5 else mag
        return (c1 + c2 * mag + c6 + C['c7'] * ((8.5 - mag) ** 2.5)
                - c3 * numpy.log(rrup + c4 * numpy.exp(c5 * mag)))
643763
def depth(self):
        
        return len(self.path.rstrip(os.sep).split(os.sep))
787184
def _check_seismogenic_depths(self, upper_depth, lower_depth):
        '''
        
        '''
        # Simple check on depths
        if upper_depth:
            if upper_depth < 0.:
                raise ValueError('Upper seismogenic depth must be greater than'
                                 ' or equal to 0.0!')
            else:
                self.upper_depth = upper_depth
        else:
            self.upper_depth = 0.0

        if lower_depth:
            if lower_depth < self.upper_depth:
                raise ValueError('Lower seismogenic depth must take a greater'
                                 ' value than upper seismogenic depth')
            else:
                self.lower_depth = lower_depth
        else:
            self.lower_depth = np.inf