Python Magnitude.evaluation_mode方法代码示例

# 需要导入模块: from obspy.core.event import Magnitude [as 别名]
# 或者: from obspy.core.event.Magnitude import evaluation_mode [as 别名]
    def _on_file_save(self):
        """
        Creates a new obspy.core.event.Magnitude object and writes the moment
        magnitude to it.
        """
        # Get the save filename.
        filename = QtGui.QFileDialog.getSaveFileName(caption="Save as...")
        filename = os.path.abspath(str(filename))
        mag = Magnitude()
        mag.mag = self.final_result["moment_magnitude"]
        mag.magnitude_type = "Mw"
        mag.station_count = self.final_result["station_count"]
        mag.evaluation_mode = "manual"
        # Link to the used origin.
        mag.origin_id = self.current_state["event"].origins[0].resource_id
        mag.method_id = "Magnitude picker Krischer"
        # XXX: Potentially change once this program gets more stable.
        mag.evaluation_status = "preliminary"
        # Write the other results as Comments.
        mag.comments.append( \
            Comment("Seismic moment in Nm: %g" % \
            self.final_result["seismic_moment"]))
        mag.comments.append( \
            Comment("Circular source radius in m: %.2f" % \
            self.final_result["source_radius"]))
        mag.comments.append( \
            Comment("Stress drop in Pa: %.2f" % \
            self.final_result["stress_drop"]))
        mag.comments.append( \
                Comment("Very rough Q estimation: %.1f" % \
            self.final_result["quality_factor"]))

        event = copy.deepcopy(self.current_state["event"])
        event.magnitudes.append(mag)
        cat = Catalog()
        cat.events.append(event)
        cat.write(filename, format="quakeml")

# 需要导入模块: from obspy.core.event import Magnitude [as 别名]
# 或者: from obspy.core.event.Magnitude import evaluation_mode [as 别名]
    def _map_netmag2magnitude(self, db):
        """
        Return an obspy Magnitude from an dict of CSS key/values
        corresponding to one record.
        
        Inputs
        ======
        db : dict of key/values of CSS fields from the 'netmag' table

        Returns
        =======
        obspy.core.event.Magnitude

        Notes
        =====
        Any object that supports the dict 'get' method can be passed as
        input, e.g. OrderedDict, custom classes, etc.
        """
        m = Magnitude()
        m.mag = db.get('magnitude')
        m.magnitude_type = db.get('magtype')
        m.mag_errors.uncertainty = db.get('uncertainty')
        m.station_count = db.get('nsta')
        
        posted_author = _str(db.get('auth'))
        mode, status = self.get_event_status(posted_author)
        m.evaluation_mode = mode
        m.evaluation_status = status
        
        m.creation_info = CreationInfo(
            creation_time = _utc(db.get('lddate')),
            agency_id = self.agency,
            version = db.get('magid'),
            author = posted_author,
            )
        m.resource_id = self._rid(m)
        return m

# 需要导入模块: from obspy.core.event import Magnitude [as 别名]
# 或者: from obspy.core.event.Magnitude import evaluation_mode [as 别名]
    def test_creating_minimal_quakeml_with_mt(self):
        """
        Tests the creation of a minimal QuakeML containing origin, magnitude
        and moment tensor.
        """
        # Rotate into physical domain
        lat, lon, depth, org_time = 10.0, -20.0, 12000, UTCDateTime(2012, 1, 1)
        mrr, mtt, mpp, mtr, mpr, mtp = 1E18, 2E18, 3E18, 3E18, 2E18, 1E18
        scalar_moment = math.sqrt(
            mrr ** 2 + mtt ** 2 + mpp ** 2 + mtr ** 2 + mpr ** 2 + mtp ** 2)
        moment_magnitude = 0.667 * (math.log10(scalar_moment) - 9.1)

        # Initialise event
        ev = Event(event_type="earthquake")

        ev_origin = Origin(time=org_time, latitude=lat, longitude=lon,
                           depth=depth, resource_id=ResourceIdentifier())
        ev.origins.append(ev_origin)

        # populate event moment tensor
        ev_tensor = Tensor(m_rr=mrr, m_tt=mtt, m_pp=mpp, m_rt=mtr, m_rp=mpr,
                           m_tp=mtp)

        ev_momenttensor = MomentTensor(tensor=ev_tensor)
        ev_momenttensor.scalar_moment = scalar_moment
        ev_momenttensor.derived_origin_id = ev_origin.resource_id

        ev_focalmechanism = FocalMechanism(moment_tensor=ev_momenttensor)
        ev.focal_mechanisms.append(ev_focalmechanism)

        # populate event magnitude
        ev_magnitude = Magnitude()
        ev_magnitude.mag = moment_magnitude
        ev_magnitude.magnitude_type = 'Mw'
        ev_magnitude.evaluation_mode = 'automatic'
        ev.magnitudes.append(ev_magnitude)

        # write QuakeML file
        cat = Catalog(events=[ev])
        memfile = io.BytesIO()
        cat.write(memfile, format="quakeml", validate=IS_RECENT_LXML)

        memfile.seek(0, 0)
        new_cat = _read_quakeml(memfile)
        self.assertEqual(len(new_cat), 1)
        event = new_cat[0]
        self.assertEqual(len(event.origins), 1)
        self.assertEqual(len(event.magnitudes), 1)
        self.assertEqual(len(event.focal_mechanisms), 1)
        org = event.origins[0]
        mag = event.magnitudes[0]
        fm = event.focal_mechanisms[0]
        self.assertEqual(org.latitude, lat)
        self.assertEqual(org.longitude, lon)
        self.assertEqual(org.depth, depth)
        self.assertEqual(org.time, org_time)
        # Moment tensor.
        mt = fm.moment_tensor.tensor
        self.assertTrue((fm.moment_tensor.scalar_moment - scalar_moment) /
                        scalar_moment < scalar_moment * 1E-10)
        self.assertEqual(mt.m_rr, mrr)
        self.assertEqual(mt.m_pp, mpp)
        self.assertEqual(mt.m_tt, mtt)
        self.assertEqual(mt.m_rt, mtr)
        self.assertEqual(mt.m_rp, mpr)
        self.assertEqual(mt.m_tp, mtp)
        # Mag
        self.assertAlmostEqual(mag.mag, moment_magnitude)
        self.assertEqual(mag.magnitude_type, "Mw")
        self.assertEqual(mag.evaluation_mode, "automatic")

# 需要导入模块: from obspy.core.event import Magnitude [as 别名]
# 或者: from obspy.core.event.Magnitude import evaluation_mode [as 别名]
    def build(self):
        """
        Build an obspy moment tensor focal mech event

        This makes the tensor output into an Event containing:
        1) a FocalMechanism with a MomentTensor, NodalPlanes, and PrincipalAxes
        2) a Magnitude of the Mw from the Tensor

        Which is what we want for outputting QuakeML using
        the (slightly modified) obspy code.

        Input
        -----
        filehandle => open file OR str from filehandle.read()

        Output
        ------
        event => instance of Event() class as described above
        """
        p = self.parser
        event         = Event(event_type='earthquake')
        origin        = Origin()
        focal_mech    = FocalMechanism()
        nodal_planes  = NodalPlanes()
        moment_tensor = MomentTensor()
        principal_ax  = PrincipalAxes()
        magnitude     = Magnitude()
        data_used     = DataUsed()
        creation_info = CreationInfo(agency_id='NN')
        ev_mode = 'automatic'
        ev_stat = 'preliminary'
        evid = None
        orid = None
        # Parse the entire file line by line.
        for n,l in enumerate(p.line):
            if 'REVIEWED BY NSL STAFF' in l:
                ev_mode = 'manual'
                ev_stat = 'reviewed'
            if 'Event ID' in l:
                evid = p._id(n)
            if 'Origin ID' in l:
                orid = p._id(n)
            if 'Ichinose' in l:
                moment_tensor.category = 'regional'
            if re.match(r'^\d{4}\/\d{2}\/\d{2}', l):
                ev = p._event_info(n)
            if 'Depth' in l:
                derived_depth = p._depth(n)
            if 'Mw' in l:
                magnitude.mag = p._mw(n) 
                magnitude.magnitude_type = 'Mw'
            if 'Mo' in l and 'dyne' in l:
                moment_tensor.scalar_moment = p._mo(n)
            if 'Percent Double Couple' in l:
                moment_tensor.double_couple = p._percent(n)
            if 'Percent CLVD' in l:
                moment_tensor.clvd = p._percent(n)
            if 'Epsilon' in l:
                moment_tensor.variance = p._epsilon(n)
            if 'Percent Variance Reduction' in l:
                moment_tensor.variance_reduction = p._percent(n)
            if 'Major Double Couple' in l and 'strike' in p.line[n+1]:
                np = p._double_couple(n)
                nodal_planes.nodal_plane_1 = NodalPlane(*np[0])
                nodal_planes.nodal_plane_2 = NodalPlane(*np[1])
                nodal_planes.preferred_plane = 1
            if 'Spherical Coordinates' in l:
                mt = p._mt_sphere(n)
                moment_tensor.tensor = Tensor(
                    m_rr = mt['Mrr'],
                    m_tt = mt['Mtt'],
                    m_pp = mt['Mff'],
                    m_rt = mt['Mrt'],
                    m_rp = mt['Mrf'],
                    m_tp = mt['Mtf'],
                    )
            if 'Eigenvalues and eigenvectors of the Major Double Couple' in l:
                ax = p._vectors(n)
                principal_ax.t_axis = Axis(ax['T']['trend'], ax['T']['plunge'], ax['T']['ev'])
                principal_ax.p_axis = Axis(ax['P']['trend'], ax['P']['plunge'], ax['P']['ev'])
                principal_ax.n_axis = Axis(ax['N']['trend'], ax['N']['plunge'], ax['N']['ev'])
            if 'Number of Stations' in l:
                data_used.station_count = p._number_of_stations(n)
            if 'Maximum' in l and 'Gap' in l:
                focal_mech.azimuthal_gap = p._gap(n)
            if re.match(r'^Date', l):
                creation_info.creation_time = p._creation_time(n)
        # Creation Time
        creation_info.version = orid
        # Fill in magnitude values
        magnitude.evaluation_mode = ev_mode
        magnitude.evaluation_status = ev_stat
        magnitude.creation_info = creation_info.copy()
        magnitude.resource_id = self._rid(magnitude)
        # Stub origin
        origin.time = ev.get('time')
        origin.latitude = ev.get('lat')
        origin.longitude = ev.get('lon')
        origin.depth = derived_depth * 1000.
        origin.depth_type = "from moment tensor inversion"
#.........这里部分代码省略.........

# 需要导入模块: from obspy.core.event import Magnitude [as 别名]
# 或者: from obspy.core.event.Magnitude import evaluation_mode [as 别名]

#.........这里部分代码省略.........
            corner_freqs = []
            for trace in stream:
                # Get the index of the pick.
                pick_index = int(round((pick.time - trace.stats.starttime) / \
                    trace.stats.delta))
                # Choose date window 0.5 seconds before and 1 second after pick.
                data_window = trace.data[pick_index - \
                    int(TIME_BEFORE_PICK * trace.stats.sampling_rate): \
                    pick_index + int(TIME_AFTER_PICK * trace.stats.sampling_rate)]
                # Calculate the spectrum.
                spec, freq = mtspec.mtspec(data_window, trace.stats.delta, 2)
                try:
                    fit = fit_spectrum(spec, freq, pick.time - origin_time,
                            spec.max(), 10.0)
                except:
                    continue
                if fit is None:
                    continue
                Omega_0, f_c, err, _ = fit
                Omega_0 = np.sqrt(Omega_0)
                omegas.append(Omega_0)
                corner_freqs.append(f_c)
            M_0 = 4.0 * np.pi * DENSITY * velocity ** 3 * distance * \
                np.sqrt(omegas[0] ** 2 + omegas[1] ** 2 + omegas[2] ** 2) / \
                radiation_pattern
            r = 3 * k * V_S / sum(corner_freqs)
            moments.append(M_0)
            source_radii.append(r)
            corner_frequencies.extend(corner_freqs)
        if not len(moments):
            print "No moments could be calculated for event %s" % \
                event.resource_id.resource_id
            continue

        # Calculate the seismic moment via basic statistics.
        moments = np.array(moments)
        moment = moments.mean()
        moment_std = moments.std()

        corner_frequencies = np.array(corner_frequencies)
        corner_frequency = corner_frequencies.mean()
        corner_frequency_std = corner_frequencies.std()

        # Calculate the source radius.
        source_radii = np.array(source_radii)
        source_radius = source_radii.mean()
        source_radius_std = source_radii.std()

        # Calculate the stress drop of the event based on the average moment and
        # source radii.
        stress_drop = (7 * moment) / (16 * source_radius ** 3)
        stress_drop_std = np.sqrt((stress_drop ** 2) * \
            (((moment_std ** 2) / (moment ** 2)) + \
            (9 * source_radius * source_radius_std ** 2)))
        if source_radius > 0 and source_radius_std < source_radius:
            print "Source radius:", source_radius, " Std:", source_radius_std
            print "Stress drop:", stress_drop / 1E5, " Std:", stress_drop_std / 1E5

        Mw = 2.0 / 3.0 * (np.log10(moment) - 9.1)
        Mw_std = 2.0 / 3.0 * moment_std / (moment * np.log(10))
        Mws_std.append(Mw_std)
        Mws.append(Mw)
        Mls.append(local_magnitude)
        calc_diff = abs(Mw - local_magnitude)
        Mw = ("%.3f" % Mw).rjust(7)
        Ml = ("%.3f" % local_magnitude).rjust(7)
        diff = ("%.3e" % calc_diff).rjust(7)

        ret_string = colorama.Fore.GREEN + \
            "For event %s: Ml=%s | Mw=%s | " % (event.resource_id.resource_id,
            Ml, Mw)
        if calc_diff >= 1.0:
            ret_string += colorama.Fore.RED
        ret_string += "Diff=%s" % diff
        ret_string += colorama.Fore.GREEN
        ret_string += " | Determined at %i stations" % len(moments)
        ret_string += colorama.Style.RESET_ALL
        print ret_string

        mag = Magnitude()
        mag.mag = Mw
        mag.mag_errors.uncertainty = Mw_std
        mag.magnitude_type = "Mw"
        mag.origin_id = event.origins[0].resource_id
        mag.method_id = "smi:com.github/krischer/moment_magnitude_calculator/automatic/1"
        mag.station_count = len(moments)
        mag.evaluation_mode = "automatic"
        mag.evaluation_status = "preliminary"
        mag.comments.append(Comment( \
            "Seismic Moment=%e Nm; standard deviation=%e" % (moment,
            moment_std)))
        mag.comments.append(Comment("Custom fit to Boatwright spectrum"))
        if source_radius > 0 and source_radius_std < source_radius:
            mag.comments.append(Comment( \
                "Source radius=%.2fm; standard deviation=%.2f" % (source_radius,
                source_radius_std)))
        event.magnitudes.append(mag)

    print "Writing output file..."
    cat.write(output_file, format="quakeml")