diff --git a/python/phaseassoc/catalogConverter.py b/python/phaseassoc/catalogConverter.py deleted file mode 100644 index 8ec288b..0000000 --- a/python/phaseassoc/catalogConverter.py +++ /dev/null @@ -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}") diff --git a/python/seedconverter/inventoryconverter.py b/python/seedconverter/inventoryconverter.py deleted file mode 100644 index c48c6b6..0000000 --- a/python/seedconverter/inventoryconverter.py +++ /dev/null @@ -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() \ No newline at end of file diff --git a/python/seedconverter/sacutil.py b/python/seedconverter/sacutil.py deleted file mode 100644 index 57dabad..0000000 --- a/python/seedconverter/sacutil.py +++ /dev/null @@ -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" diff --git a/python/seedconverter/seedconverter.py b/python/seedconverter/seedconverter.py deleted file mode 100644 index 4feff78..0000000 --- a/python/seedconverter/seedconverter.py +++ /dev/null @@ -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) diff --git a/python/seedconverter/unit.py b/python/seedconverter/unit.py deleted file mode 100644 index c38ce96..0000000 --- a/python/seedconverter/unit.py +++ /dev/null @@ -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)