Removed unnecessary python snippets

This commit is contained in:
Mieszko Makuch 2024-08-29 12:59:04 +02:00
parent c81a456449
commit 450884838e
5 changed files with 0 additions and 702 deletions

View File

@ -1,120 +0,0 @@
import pandas as pd
from obspy.core.event import Catalog
from obspy.core.event import Event
from obspy.core.event import Pick
from obspy.core.event import Origin
from obspy.core.event import OriginQuality
from obspy import UTCDateTime
from obspy.core.event.base import WaveformStreamID, Comment
class CatalogConverter:
"""
Class for converting SeisBench and pyOcto detection results to ObsPy Catalog, which can be saved as QuakeML.
"""
def __init__(self, config, picks, catalog_df, assignments_df, name_of_algorithm):
self._catalog_df = catalog_df
self._assignments_df = assignments_df
self._picks = picks
self._config = config
self._name_of_algorithm = name_of_algorithm
self.catalog = None
def _convert_pick(self, p):
"""
Function converts picks from SeisBench to ObsPy format
:param p: SeisBench pick
:return: ObsPy pick
"""
pick = Pick()
pick.time = UTCDateTime(p.peak_time.datetime)
pick.waveform_id = WaveformStreamID(network_code=p.trace_id.split(".")[0],
station_code=p.trace_id.split(".")[1],
channel_code=p.trace_id.split(".")[2])
if p.phase == 'P':
pick.phase_hint = self._config["P_hint"]
elif p.phase == 'S':
pick.phase_hint = self._config["S_hint"]
pick.evaluation_mode = 'automatic'
pick.evaluation_status = 'preliminary'
return pick
def _convert_origin(self, origin_sb, list_of_picks_sb):
origin = Origin()
origin.time = UTCDateTime(pd.to_datetime(origin_sb.time, unit='s').to_pydatetime())
origin.latitude = origin_sb.latitude # float
origin.longitude = origin_sb.longitude # float
origin.depth = origin_sb.depth # float in kilometers (SWIP5 origin version) down the see level
origin.depth_type = 'operator assigned'
# TODO: make sure that region is not necessary
# origin.region = self._config["region"]
origin.evaluation_mode = "automatic"
origin.evaluation_status = 'preliminary'
origin.comments.append(Comment(text=f"Localized by: {self._name_of_algorithm}", force_resource_id=False))
origin.quality = OriginQuality(used_phase_count=len(list_of_picks_sb))
return origin
def _convert_event(self, origin_sb, list_of_picks_sb):
"""
Function convert GaMMa detection to ObsPy Event
:param origin_sb:
:param list_of_picks_sb:
:return:
"""
event = Event()
for p in list_of_picks_sb:
pick = self._convert_pick(p)
event.picks.append(pick)
origin = self._convert_origin(origin_sb, list_of_picks_sb)
event.origins.append(origin)
return event
@staticmethod
def _append_pick_trace_id(pick, stream):
"""
Function assigns channel to pick - it is useful for work with SWIP
:param pick:
:param stream:
:return:
"""
channel = stream[0].stats.channel
if pick.phase == "P":
pick.trace_id = pick.trace_id + channel[:-1] + "Z"
if pick.phase == "S":
pick.trace_id = pick.trace_id + channel[:-1] + "E"
return pick
def catalog2obspy(self):
"""
Function convert GaMMa catalog and SeisBench picks
:return: ObsPy Catalog object
"""
# TODO: make sure that resource id is necessary
#cat = Catalog(resource_id=self._config["resource_id"])
cat = Catalog()
for j, row in self._catalog_df.iterrows():
event = self._catalog_df.iloc[j]
event_picks = [self._picks[i] for i in
self._assignments_df[self._assignments_df["event_idx"] ==
event["idx"]]["pick_idx"]]
event_obspy = self._convert_event(event, event_picks)
cat.append(event_obspy)
self.catalog = cat
def save_catalog_to_file(self, file_path):
"""
Save ObsPy catalog to a file.
Args:
file_path (str): The file path where the catalog will be saved.
Returns:
None
"""
try:
self.catalog.write(file_path, format="QUAKEML")
print(f"Catalog saved successfully to {file_path}")
except Exception as e:
print(f"Error occurred while saving catalog: {e}")

View File

@ -1,62 +0,0 @@
import os
import argparse
from obspy import read_inventory
def convert_inventory(input_file,
output_format,
output_dir="."):
"""
Function to read Inventory from provided file and convert it to FDSN Station XML format.
The supported input file formats are: INVENTORYXML, RESP, SC3ML, SEED, STATIONTXT, STATIONXML, XSEED
:type input_file: str
:param input_file: File name or path to the input inventory file.
:type output_format: str
:param output_format: Format of the output inventory file.
Supported formats: CSS, KML, SACPZ, SHAPEFILE, STATIONTXT, STATIONXML.
:type output_dir: str, optional
:param output_dir: Directory to which output files are written.
Defaults to current directory.
"""
inv = read_inventory(input_file)
result_filename = os.path.splitext(os.path.basename(input_file))[0] + "." + _get_extension(output_format)
inv.write(output_dir + "/" + result_filename, format=output_format)
return result_filename
def _get_extension(output_format):
format = output_format.upper()
if format == 'STATIONXML':
return "xml"
elif format == 'STATIONTXT':
return "txt"
elif format == 'SHAPEFILE':
return "shp"
else:
return format.lower()
def main():
parser = argparse.ArgumentParser(description="Convert provided inventory file"
" to another inventory format.")
parser.add_argument("input_file", help="Provide inventory file to convert."
"The supported input file formats are: INVENTORYXML, RESP, SC3ML, "
"SEED, STATIONTXT, STATIONXML, XSEED")
parser.add_argument("--output_format",
help="Format of the output inventory file. "
"Supported formats: CSS, KML, SACPZ, SHAPEFILE, STATIONTXT, STATIONXML.",
type=str, default=None, required=True)
parser.add_argument("--output_dir",
help="Directory to which output files are written. "
"Defaults to current directory.",
type=str, default=".", required=False)
args = parser.parse_args()
filename = convert_inventory(**vars(args))
print('Created file:')
print(filename)
return
if __name__ == "__main__":
main()

View File

@ -1,98 +0,0 @@
from obspy.io.sac import SACTrace
def to_sac_trace(tr,
unit_type,
remove_response,
inv):
"""
Function converting the specified trace to SAC-specific format
:type tr: obspy.core.trace.Trace
:param tr: Trace to be converted
:type unit_type: str
:param unit_type: Unit type of the result
:type remove_response: bool
:param remove_response: Information whether the response was removed, used to determine scale
:type inv: obspy.core.inventory.inventory.Inventory
:param inv: Inventory read either from SEED file or Station XML
:rtype: obspy.core.trace.Trace
:return: Converted trace
"""
# Assure that tr has fields present in SAC-type trace
sac = SACTrace.from_obspy_trace(tr)
sac = _prepare_sac_for_writing(sac, remove_response, unit_type, inv, tr.id, tr.stats['starttime'])
return sac.to_obspy_trace()
def _prepare_sac_for_writing(sac, remove_response, unit_type, inv, tr_id, tr_starttime):
"""
Method that adds some additional information to ordinary SACTrace file.
It adds location info and information about unit type and scale.
When the response is removed, it scales down the data to nano units
nm/s or so.
:type sac: obspy.io.sac.SACTrace
:param sac: SACTrace object, might be converted out of ordinary trace
:type remove_response: bool
:param remove_response: Parameter that lets user decide whether the
response is going to be removed or not
:type unit_type: str
:param unit_type: Unit type determined by _get_unit_type
:type inv: obspy.core.inventory.inventory.Inventory
:param inv: inventory contents
:type tr_id: str
:param tr_id: Trace's id string
:type tr_starttime: UTCDateTime
:param tr_starttime: Trace's start time
:rtype: obspy.io.sac.SACTrace
:return: SACTrace object with added additional information
"""
# Copy details of station's location
sac.idep = _get_sac_idep_unit(unit_type)
if inv:
try:
coords = inv.get_coordinates(tr_id, tr_starttime)
sensitivity = inv.get_response(tr_id, tr_starttime).instrument_sensitivity.value
sac.stla = coords["latitude"]
sac.stlo = coords["longitude"]
sac.stel = coords["elevation"]
sac.user2 = sensitivity
if remove_response:
# Divide by nano multiplier to get a data with removed response
# scaled to nano unit (nm/s or so).
sac.data = sac.data / _get_nano_multiplier()
sac.scale = -12345
else:
sac.scale = (1 / sensitivity / _get_nano_multiplier())
except Exception as err:
raise Exception("Cannot read location and sensitivity information: " + str(err))
else:
sac.scale = -12345
return sac
def _get_nano_multiplier():
"""
Returns a multiplier to obtain nano value from a full value.
:rtype: float
:return: Multiplier value
"""
return 10**(-9)
def _get_sac_idep_unit(unit_abbrev):
if unit_abbrev == "DISP":
return "idisp"
elif unit_abbrev == "VEL":
return "ivel"
elif unit_abbrev == "ACC":
return "iacc"
else:
return "iunkn"

View File

@ -1,327 +0,0 @@
import sys
import os
import warnings
import obspy
import argparse
import unit
import sacutil
import numpy as np
from obspy.io.xseed import Parser
from obspy.io.xseed.utils import SEEDParserException
from obspy.core.util import AttribDict
from obspy.core import Stream
def convert_seed(seed_file,
inv_file=None,
input_unit=None,
output_filetype=None,
remove_response=False,
zero_mean=False,
pre_filtering=None,
output_unit=None,
single_output_file=False,
output_dir="."):
"""
Function to convert provided file to SAC, ASCII or MSEED
:type seed_file: str
:param seed_file: File name or path to the seed (seed or msd) file.
:type inv_file: str
:param inv_file: File name or path to inventory file in any inventory format supported by ObsPy.
:type input_unit: str
:param input_unit: Unit of the input file: 'ACC', 'DISP' or 'VEL'. Specified only for miniSEED files (full SEED
files are always in 'COUNTS'). If the value is not set, no unit ('COUNTS') is assumed.
:type output_filetype: str
:param output_filetype: Output filetype. Choose either SAC, MSEED, SLIST or INTERNAL_ASCII.
When used SLIST, sample values are saved with `%+.10e` formatting.
INTERNAL_ASCII is an internal format used for some apps, that uses only one column of values
This will be default formatting since ObsPy 1.1.0
:type remove_response: bool
:param remove_response: Response removal from the waveform
Defaults to False.
Removal of the instrument's response from the waveform
By default it will save waveform without response removal procedure
This function uses obspy.core.stream.Stream.remove_response.
When set to True, user has to specify parameters zero_mean,
pre_filtering and output_unit. Otherwise the script will raise an
exception.
:type zero_mean: bool, optional
:param zero_mean: If true the mean of the data is subtracted
For details check the ObsPy's documentation for
obspy.core.stream.Stream.remove_response
:type pre_filtering: (float, float, float, float), optional
:param pre_filtering: Apply a bandpass filter to the data trace before
deconvolution.
The list or tuple defines the four corner frequencies (f1,f2,f3,f4)
of a cosine taper which is one between f2 and f3
and tapers to zero for f1 < f < f2 and f3 < f < f4.
For details check the ObsPy's documentation for
obspy.core.stream.Stream.remove_response
:type output_unit: str, optional
:param output_unit: Unit of the output waveform.
This parameter converts counts to real values during response removal,
or, if response_removal is false is just written to the converted file.
When used ``auto`` the script will determine the unit automatically
with use of second letter of a channel code which determines
instrument type.
For details see Appendix A to SEED Manual v2.4.
In case user need to convert all channels to common unit,
there can be chosen ``DISP`` for displacement,
``VEL`` for velocity and ``ACC`` for acceleration.
When this parameter is None and remove_response is set to True
script will throw an Exception.
:type single_output_file: bool
:param single_output_file: if True, a single file with all traces will be created
(available only for conversion to MSEED), otherwise, separate files with names in format
``network.station.location.channel.starttime_microsecond.output_filetype`` will be created
:type output_dir: str, optional
:param output_dir: Directory to which output files are written.
Defaults to current directory.
:return: Filenames of created files. Filenames follow pattern:
``network.station.location.channel.starttime_microsecond.SAC``
:raises
.. note::
For more information on parameters used in this function see
ObsPy's documentation in ``obspy.core.stream.Stream.remove_response``
"""
st = obspy.read(seed_file)
if inv_file:
inv = obspy.read_inventory(inv_file)
elif remove_response:
raise ValueError("Cannot remove response if inventory is not specified")
else:
inv = None
if output_unit:
if inv:
# assuming all the traces have the same unit
input_type = unit.determine_instrument_type_from_inventory(inv, st[0].id, st[0].stats.starttime)
else:
try:
# if we don't have the inventory and we have a full SEED file, we can attempt to read unit from it
resp = Parser(seed_file)
# assuming all the traces have the same unit
input_type = unit.determine_instrument_type_from_blockette(resp, st[0].id)
except (OSError, SEEDParserException):
warnings.warn("The provided file is not a SEED volume - cannot determine unit automatically")
if output_unit == "auto":
output_unit = input_type
if remove_response and output_unit is None:
raise ValueError("You have to provide output_unit parameter"
" for response removal procedure.\n"
"Available output_unit values are DISP, "
"VEL, ACC. Provide one or skip removing "
"response by passing "
"remove_response=False parameter")
if remove_response:
if pre_filtering is None:
tr_sampling_rate = st[0].stats.sampling_rate
pre_filt = [0.005, 0.006, 0.9 * 0.5 * tr_sampling_rate, 0.5 * tr_sampling_rate]
else:
pre_filt = pre_filtering
try:
_remove_response(st, inv, input_type, output_unit, pre_filt, zero_mean)
except (ValueError, SEEDParserException) as exc:
raise Exception("There was a problem with removing response from the file: " + str(exc))
elif output_unit and input_unit and input_unit != output_unit:
_convert_unit(st, input_unit, output_unit)
# we are not converting anything, but the input file already had a unit (might happen only with mseed files),
# so it should be written to the output file
elif input_unit:
output_unit = input_unit
if single_output_file:
if output_filetype != "MSEED":
raise ValueError("Writing all traces to one file is available only for conversion to MSEED format")
return _convert_to_single_mseed(seed_file, st, output_dir)
else:
return _convert_to_separate_files(st, inv, remove_response, output_unit, output_filetype, output_dir)
# remove response and correct numerical errors (depending on the input and output type)
def _remove_response(stream, inv, input_type, output_unit, pre_filt, zero_mean):
if input_type == "ACC" and output_unit == "DISP":
stream.remove_response(inventory=inv, pre_filt=pre_filt, zero_mean=zero_mean, output="VEL")
stream.filter("highpass", freq=1.0)
stream.integrate(method="cumtrapz")
stream.filter("highpass", freq=0.5)
else:
stream.remove_response(inventory=inv, pre_filt=pre_filt, zero_mean=zero_mean, output=output_unit)
if ((input_type == "ACC" and output_unit == "VEL") or
(input_type == "VEL" and (output_unit == "DISP" or output_unit == "VEL"))):
stream.filter("highpass", freq=1.0)
# unit conversion - applied when the input data already has response removed and is converted to one of the units
def _convert_unit(stream, input_unit, output_unit):
if input_unit == "DISP":
if output_unit == "ACC": # DIFF x2
stream.differentiate(method="gradient")
stream.differentiate(method="gradient")
elif output_unit == "VEL": # DIFF
stream.differentiate(method="gradient")
elif input_unit == "VEL":
if output_unit == "ACC": # DIFF
stream.differentiate(method="gradient")
elif output_unit == "DISP": # INTEG + FILTER
stream.integrate(method="cumtrapz")
stream.filter("highpass", freq=1.0)
elif input_unit == "ACC":
if output_unit == "VEL": # INTEG + FILTER
stream.integrate(method="cumtrapz")
stream.filter("highpass", freq=1.0)
elif output_unit == "DISP": # (INTEG + FILTER) x2
stream.integrate(method="cumtrapz")
stream.filter("highpass", freq=1.0)
stream.integrate(method="cumtrapz")
stream.filter("highpass", freq=0.5)
else:
raise TypeError("Cannot convert from ", input_unit)
def _convert_to_single_mseed(seed_filename, stream, output_dir):
result_filename = os.path.splitext(os.path.basename(seed_filename))[0] + ".msd"
stream.write(output_dir + "/" + result_filename, format="MSEED")
return [result_filename]
def _convert_to_separate_files(stream, inv, remove_response, output_unit, output_filetype, output_dir):
output_stream = Stream()
output_file_extension = output_filetype
for tr in stream:
if output_filetype == "SAC":
output_file_extension = "sac"
tr = sacutil.to_sac_trace(tr, output_unit, remove_response, inv)
elif output_filetype in ["SLIST", "TSPAIR", "SH_ASC"]:
output_file_extension = output_filetype.lower() + ".ascii"
slist_unit = unit.translate_instrument_type_to_unit(output_unit) if output_unit else "COUNTS"
tr.stats.ascii = AttribDict({"unit": slist_unit})
elif output_filetype == "MSEED":
output_file_extension = "msd"
elif output_filetype == "INTERNAL_ASCII":
output_file_extension = "ascii"
output_stream.append(tr)
return _write_separate_output_files(output_stream, output_filetype, output_file_extension, output_dir)
def _write_separate_output_files(stream, output_filetype, file_extension, output_dir):
"""
Writes all Trace objects present in Stream to separate files. Filetype
is set according to User's will.
:type stream: obspy.core.stream.Stream
:param stream: Obspy stream with traces to save separately.
:type output_filetype: str
:param output_filetype: Contains filetype to save. Currently supported
``SAC``, ``SLIST`` and ``MSEED``.
"""
result_filenames = list()
for tr in stream:
result_filename = ".".join([tr.id,
str(tr.stats.starttime.microsecond),
file_extension])
result_filepath = output_dir + "/" + result_filename
if output_filetype == "INTERNAL_ASCII":
_write_ascii_one_column(tr, result_filepath)
else:
tr.write(result_filepath, format=output_filetype)
result_filenames.append(result_filename)
return result_filenames
def _write_ascii_one_column(trace, filename):
"""
Writes trace data into an ASCII file in one-column format (i.e. no header, single column of values).
"""
with open(filename, 'wb') as fh:
data = trace.data.reshape((-1, 1))
dtype = data.dtype.name
fmt = '%f'
if dtype.startswith('int'):
fmt = '%d'
elif dtype.startswith('float64'):
fmt = '%+.16e'
np.savetxt(fh, data, delimiter=b'\t', fmt=fmt.encode('ascii', 'strict'))
def main(argv):
def str2bool(v):
if v.lower() in ("True", "TRUE", "yes", "true", "t", "y", "1"):
return True
elif v.lower() in ("False", "FALSE", "no", "false", "f", "n", "0"):
return False
else:
raise argparse.ArgumentTypeError("Boolean value expected.")
parser = argparse.ArgumentParser(description="Convert provided file"
" to SAC or SLIST (ASCII).")
parser.add_argument("seed_file", help="Provide SEED or mSEED file to convert")
parser.add_argument("--inv_file", help="Provide inventory file")
parser.add_argument("--input_unit",
help="Provide input unit. "
"ACC, VEL or DISP are available.",
type=str, default=None, required=False)
parser.add_argument("--output_filetype",
help="Provide output filetype. "
"MSEED, SAC, SLIST and INTERNAL_ASCII are available.",
type=str, default=None, required=True)
parser.add_argument("--remove_response",
help="Remove instrument's response from the waveform",
type=str2bool, default=False)
parser.add_argument("--zero_mean",
help="If true the mean of the data is subtracted",
type=str2bool, default=False, required=False)
parser.add_argument("--pre_filtering",
help="Apply a bandpass filter to the data trace "
"before deconvolution",
type=float, nargs="+", default=None,
required=False)
parser.add_argument("--output_unit",
help="Unit to which waveform should be converted "
"during response removal. When set to auto, "
"the unit will be automatically determined with "
"use of the second letter of channel code which "
"determines the instrument type."
"For details see Appendix A, SEED Manual v2.4.",
type=str, default=None, required=False)
parser.add_argument("--single_output_file",
help="Flag deciding whether all the traces will be written to "
"a single file, if set to False, each trace will be "
"written to a separate file. True is available only "
"for conversion to MSEED",
type=str2bool, default=False, required=False)
parser.add_argument("--output_dir",
help="Directory to which output files are written. "
"Defaults to current directory.",
type=str, default=".", required=False)
args = parser.parse_args()
filenames = convert_seed(**vars(args))
print('Created files:')
print(', '.join(filenames))
return
if __name__ == "__main__":
main(sys.argv)

View File

@ -1,95 +0,0 @@
def determine_instrument_type_from_blockette(parser, channel_id):
"""
Determines type of instrument used to record a channel with provided
channel_id.
:type parser: obspy.io.xseed.parser.Parser
:param parser: Parser object parsed from Full SEED or dataless or
StationXML.
:type channel_id: str
:param channel_id: Channel_id object generated by obspy.trace.Trace.id
:rtype: str
:return: Returns str determining the type of instrument
"""
return translate_unit_to_instrument_type(determine_unit_from_blockette(parser, channel_id))
def determine_unit_from_blockette(parser, channel_id):
"""
Determines the unit used to record a channel with provided
channel_id.
:type parser: obspy.io.xseed.parser.Parser
:param parser: Parser object parsed from Full SEED or dataless or
StationXML.
:type channel_id: str
:param channel_id: Channel_id object generated by obspy.trace.Trace.id
:rtype: str
:return: Returns str containing a real unit name that can be passed to
translate_unit_to_instrument_type method
to obtain a str compatible with response removal procedure
"""
for blkt in parser._select(channel_id):
if not blkt.id == 52:
continue
for bl in parser.blockettes[34]:
if bl.unit_lookup_code == blkt.units_of_signal_response:
return bl.unit_name
def determine_instrument_type_from_inventory(inv, channel_id, time):
"""
Determines type of instrument used to record a channel with provided channel_id at the provided time.
:type inv: obspy.core.inventory.inventory.Inventory
:param inv: ObsPy Inventory object parsed from a file
:type channel_id: str
:param channel_id: Channel_id object generated by obspy.trace.Trace.id
:type time: obspy.core.utcdatetime.UTCDateTime
:param time: time for which the unit should be determined in the inventory (e.g. start time of a trace)
:rtype: str
:return: Returns str determining the type of instrument
"""
return translate_unit_to_instrument_type(determine_unit_from_inventory(inv, channel_id, time))
def determine_unit_from_inventory(inv, channel_id, time):
"""
Determines unit used to record a channel with provided channel_id at the provided time.
:type inv: obspy.core.inventory.inventory.Inventory
:param inv: ObsPy Inventory object parsed from a file
:type channel_id: str
:param channel_id: Channel_id object generated by obspy.trace.Trace.id
:type time: obspy.core.utcdatetime.UTCDateTime
:param time: time for which the unit should be determined in the inventory (e.g. start time of a trace)
:rtype: str
:return: Returns str containing a real unit name that can be passed to
translate_unit_to_instrument_type method
to obtain a str compatible with response removal procedure
"""
resp = inv.get_response(channel_id, time)
return resp.instrument_sensitivity.input_units
def translate_unit_to_instrument_type(unit_name):
if unit_name == "M":
return "DISP"
elif unit_name == "M/S":
return "VEL"
elif unit_name == "M/S**2":
return "ACC"
else:
raise TypeError("Unknown unit code ", unit_name)
def translate_instrument_type_to_unit(unit_type):
if unit_type == "DISP":
return "M"
elif unit_type == "VEL":
return "M/S"
elif unit_type == "ACC":
return "M/S**2"
else:
raise TypeError("Unknown unit code ", unit_type)