python tools for working with strainmeter data
# Configuration file for the Sphinx documentation builder.
#
# For the full list of built-in configuration values, see the documentation:
# https://www.sphinx-doc.org/en/master/usage/configuration.html

# -- Project information -----------------------------------------------------
# https://www.sphinx-doc.org/en/master/usage/configuration.html#project-information

project = 'Earthscope Strain Tools'
copyright = '2023, Earthscope Consortium'
author = 'Mike Gottlieb, Catherine Hanagan'
release = '1'

# If extensions (or modules to document with autodoc) are in another directory,
# add these directories to sys.path here.
# note that the absolute paths are referenced relative to the configuration directory (with conf.py)
import os
import sys
sys.path.insert(0, os.path.abspath('.'))
sys.path.insert(0, os.path.abspath('../../src'))
sys.path.insert(0, os.path.abspath('../../src/earthscopestraintools'))

# To keep same env as the main project, will allow notebooks in the docs
# by downloading appropriate packages
from sphinx_pyproject import SphinxConfig
config = SphinxConfig("../../pyproject.toml", globalns=globals())

# -- General configuration ---------------------------------------------------
# https://www.sphinx-doc.org/en/master/usage/configuration.html#general-configuration

extensions = ['sphinx.ext.autodoc',
    'sphinx.ext.autosummary',
    'sphinx.ext.duration',
    'sphinx.ext.doctest',
    'sphinxcontrib.apidoc',
    'sphinx.ext.napoleon',
    'nbsphinx',
    'myst_parser'] # myst will allow document building in markdown syntax
myst_enable_extensions = [ "amsmath",
    "dollarmath",
    "colon_fence"]

nbsphinx_kernel_name = 'python'
nbsphinx_allow_errors = True

# APIDOC configuration 
apidoc_module_dir = "../../src/earthscopestraintools"
apidoc_output_dir = "api"
apidoc_separate_modules = True
apidoc_toc_file = False
apidoc_module_first = True

templates_path = ['_templates']
exclude_patterns = []


# -- Options for HTML output -------------------------------------------------
# https://www.sphinx-doc.org/en/master/usage/configuration.html#options-for-html-output

html_theme = 'sphinx_rtd_theme'
html_logo = "./_static/EarthScope_Logo-color.png"



docs/source/conf.py
EarthScope/earthscopestraintools
# ETL from ascii level2 files to tiledb
# write tiledb arrays locally, then aws s3 sync
import tarfile
import tiledb
import pandas as pd
import datetime
import sys, os
import shutil

# from edid import find_station_edid
#from earthscopestraintools.datasources_api_interact import get_station_edid
import json
from io import BytesIO
import requests
import configparser

from earthscopestraintools.tiledbtools import (
    ProcessedStrainWriter,
    ProcessedStrainReader,
)


import logging

logging.basicConfig(
    level=logging.INFO,
    format="%(asctime)s %(levelname)s: %(message)s",
    datefmt="%Y-%m-%d %H:%M:%S",
)
logger = logging.getLogger(__name__)
workdir = ""


# def write_df_to_tiledb(df, array):
#     # print(df)
#     mode = "append"
#     tiledb.from_pandas(
#         uri=array.uri,
#         dataframe=df,
#         index_dims=["data_type", "timeseries", "time"],
#         mode=mode,
#         ctx=array.ctx,
#     )
#
#     # update the string dimension metadata
#     data_type = df["data_type"].unique()
#     timeseries = df["timeseries"].unique()
#     if type(data_type) == str:
#         data_type = [data_type]
#     if type(timeseries) == str:
#         timeseries = [timeseries]
#     with tiledb.open(array.uri, "r", ctx=array.ctx) as A:
#         try:
#             dimension_json = A.meta["dimensions"]
#         except KeyError:
#             dimension_json = '{"data_types":[], "timeseries":[]}'
#
#         dimension_dict = json.loads(dimension_json)
#         # print(dimension_dict)
#         for item in data_type:
#             if item not in dimension_dict["data_types"]:
#                 dimension_dict["data_types"].append(item)
#         for item in timeseries:
#             if item not in dimension_dict["timeseries"]:
#                 dimension_dict["timeseries"].append(item)
#
#         # print(dimension_dict)
#         with tiledb.open(array.uri, "w", ctx=array.ctx) as A:
#             A.meta["dimensions"] = json.dumps(dimension_dict)
#


def loop_through_ts(filebase, file, writer):
    df = pd.read_csv(filebase + "/" + file, delimiter="\s+")
    label = df["strain"][0] + "(mstrain)"
    df = df.rename(
        columns={
            "strain": "data_type",
            label: "microstrain",
            "s_offset": "offset_c",
            "detrend_c": "trend_c",
            "MJD": "mjd",
            "date": "time",
        }
    )
    df["data_type"] = df["data_type"].str.replace("gauge", "CH")
    # display(df)
    cols = ["data_type", "timeseries", "time", "data", "quality", "level", "version"]
    # convert time string to unix ms
    df["time"] = (pd.to_datetime(df["time"]).astype(int) / 10**6).astype(int)

    timeseries = "microstrain"
    tmp_df = df[["data_type", "time", timeseries, "strain_quality", "level", "version"]]
    tmp_df = tmp_df.rename(columns={timeseries: "data", "strain_quality": "quality"})
    tmp_df["timeseries"] = timeseries
    tmp_df = tmp_df[cols]
    logger.info(f"{timeseries}: {len(tmp_df)} samples")
    writer.write_df_to_tiledb(tmp_df)

    timeseries = "offset_c"
    tmp_df = df[["data_type", "time", timeseries, "level", "version"]]
    tmp_df = tmp_df.rename(columns={timeseries: "data"})
    tmp_df["timeseries"] = timeseries
    tmp_df["quality"] = " "
    tmp_df = tmp_df[cols]
    logger.info(f"{timeseries}: {len(tmp_df)} samples")
    writer.write_df_to_tiledb(tmp_df)

    timeseries = "tide_c"
    tmp_df = df[["data_type", "time", timeseries, "level", "version"]]
    tmp_df = tmp_df.rename(columns={timeseries: "data"})
    tmp_df["timeseries"] = timeseries
    tmp_df["quality"] = " "
    tmp_df = tmp_df[cols]
    logger.info(f"{timeseries}: {len(tmp_df)} samples")
    writer.write_df_to_tiledb(tmp_df)

    timeseries = "trend_c"
    tmp_df = df[["data_type", "time", timeseries, "level", "version"]]
    tmp_df = tmp_df.rename(columns={timeseries: "data"})
    tmp_df["timeseries"] = timeseries
    tmp_df["quality"] = " "
    tmp_df = tmp_df[cols]
    logger.info(f"{timeseries}: {len(tmp_df)} samples")
    writer.write_df_to_tiledb(tmp_df)

    timeseries = "atmp_c"
    tmp_df = df[["data_type", "time", timeseries, "atmp_c_quality", "level", "version"]]
    tmp_df = tmp_df.rename(columns={timeseries: "data", "atmp_c_quality": "quality"})
    tmp_df["timeseries"] = timeseries
    tmp_df = tmp_df[cols]
    logger.info(f"{timeseries}: {len(tmp_df)} samples")
    writer.write_df_to_tiledb(tmp_df)

    if df["data_type"][0] == "CH0":
        tmp_df = df[["data_type", "time", "atmp", "level", "version"]]
        tmp_df = tmp_df.rename(columns={"atmp": "data"})
        tmp_df = tmp_df.assign(data_type="atmp")
        tmp_df["timeseries"] = "hpa"
        tmp_df["quality"] = " "
        tmp_df = tmp_df[cols]
        logger.info(f"{timeseries}: {len(tmp_df)} samples")
        writer.write_df_to_tiledb(tmp_df)


def etl_yearly_ascii_file(
    network, station, year, delete_array=False, workdir="", print_it=False
):
    os.makedirs(workdir, exist_ok=True)
    # edid = get_station_edid(station)
    # uri = f"{workdir}/{edid}_level2.tdb"
    uri = f"{workdir}/{network}_{station}_l2_etl.tdb"
    logger.info(f"Array uri: {uri}")
    writer = ProcessedStrainWriter(uri)
    if delete_array:
        writer.array.delete()

    # create new array if needed.  note: array_exists only works locally not in s3.
    if not tiledb.array_exists(writer.array.uri):
        writer.array.create(schema_type="3d", schema_source="s3")
        writer.array.set_array_meta(network=network, station=station, period=300)

    filebase = station + "." + year + ".bsm.level2"
    url = "http://bsm.unavco.org/bsm/level2/" + station + "/" + filebase + ".tar"
    response = requests.get(url, stream=True)
    print(url)
    tar = tarfile.open(fileobj=BytesIO(response.raw.read()), mode="r")
    tar.extractall()
    files = os.listdir(filebase)
    for file in files:
        logger.info(file)
        loop_through_ts(filebase, file, writer)
    shutil.rmtree(filebase)

    writer.array.consolidate_array_meta()
    writer.array.vacuum_array_meta()
    writer.array.consolidate_fragment_meta()
    writer.array.vacuum_fragment_meta()
    writer.array.consolidate_fragments()
    writer.array.vacuum_fragments()
    if print_it:
        reader = ProcessedStrainReader(uri)
        logger.info(f"Network: {reader.array.get_network()}")
        logger.info(f"Station: {reader.array.get_station()}")
        logger.info(f"Period: {reader.array.get_period()}")
        start_str = f"{year}-01-01T00:00:00"
        end_str = f"{int(year)+1}-01-01T00:00:00"
        df = reader.to_df(
            data_types=["CH0", "CH1", "CH2", "CH3", "Eee+Enn", "Eee-Enn", "2Ene"],
            timeseries="microstrain",
            attrs="data",
            start_str=start_str,
            end_str=end_str,
        )
        logger.info(f"\n{df}")


# if __name__ == "__main__":
#     workdir = "arrays"
#     network = sys.argv[1]
#     station = sys.argv[2]
#     year = sys.argv[3]
#     etl_yearly_ascii_file(network, station, year, workdir=workdir)
# years = ["2005","2006","2007","2008","2009",
#          "2010","2011","2012","2013","2014","2015","2016","2017","2018","2019",
#          "2020","2021","2022"]
# for year in years:
#     etl_yearly_ascii_file(network, station, year)
#

src/earthscopestraintools/ascii2tdb.py
EarthScope/earthscopestraintools
import os
import datetime
import struct
import logging

logger = logging.getLogger(__name__)

BTL_EPOCH = datetime.datetime(1970, 1, 1, 0, 0, 0, 0, None)
BTL_TYPE = {0: 'INT2', 1: 'INT4', 2: 'FLOAT4'}
BTL_MISSING = {0: 32767, 1: 999999, 2: 999999}
BTL_TYPE_SIZE = {0: 2, 1: 4, 2: 4}
BTL_TYPE_TO_PACK = {0: 'h', 1: 'i', 2: 'f'}
BTL_TYPE_TO_PRINT = {0: 'd', 1: 'd', 2: 'g'}

class Bottle:
    def __init__(self, filepath, fileobj=None):
        """
        """
        self.file_metadata = {}
        self.file_metadata['filepath'] = filepath
        self.file_metadata['filename'] = filepath.split("/")[-1]
        if fileobj and filepath:
            self.file = fileobj
        elif filepath:
            if not os.path.isfile(self.file_metadata['filepath']):
                logger.error("Can not find file '%s'" % self.file_metadata['filepath'])

            self.file = open(self.file_metadata['filepath'], 'rb')

        magic = self.file.read(2)
        if chr(magic[0]) == chr(0x01) and chr(magic[1]) == chr(0x5D):
            self.file_metadata['endian'] = '>'
        elif chr(magic[0]) == chr(0x5D) and chr(magic[1]) == chr(0x01):
            self.file_metadata['endian'] = '<'
        else:
            print(magic)
            print("Not recognized as a bottle file: %s" % self.file_metadata['filename'])
            logger.error("Not a bottle", self.file_metadata['filename'])

    def parse_filename(self):
        """
        Add info encoded within file name of bottle file from GTSM
        Logger to file metadata.  hour and min may be None.  All values are strings, even
        if they contain only digits.
        """

        # FIXME: we are hardcoding here knowledge of how to pick apart
        # a GTSM file name.  no other option really, but just be wary
        # of trusting file names.  and if you feed in something which
        # isn't a GTSM bottle file name, you'll get junk back - or
        # raise an exception.

        f = self.file_metadata['filename']

        stn4char = f[0:4]
        year = "20" + f[4:6]
        dayofyear = f[6:9]

        if f[-3:] == '_20' and len(f) == 19:
            # this is a Min file
            self.file_metadata['session'] = "Min"
            hour = f[9:11]
            min = f[11:13]
            channel = f[13:16]

        elif f[-3:] == '_20' and len(f) == 17:
            # this is the concatenation of 60 Min files
            hour = f[9:11]
            min = None
            channel = f[11:14]

        elif f[-3:-1] == 'CH' and len(f) == 14:
            # this is an Hour file
            self.file_metadata['session'] = "Hour"
            hour = f[9:11]
            min = None
            channel = f[11:14]

        else:
            # well, we have to assume it is a Day file
            self.file_metadata['session'] = "Day"
            hour = None
            min = None
            channel = f[9:]

        self.file_metadata['fcid'] = stn4char
        self.file_metadata['channel'] = channel
        self.file_metadata['year'] = year
        self.file_metadata['dayofyear'] = dayofyear
        self.file_metadata['hour'] = hour
        self.file_metadata['min'] = min

        self.get_seed_codes()

        # return self.file_metadata #(stn4char, channel, year, dayofyear, hour, min)

    def get_seed_codes(self):
        codes = {"BatteryVolts": ["T0", "RE1"], "CH0": ["T0", "RS1"], "CH1": ["T0", "RS2"], "CH2": ["T0", "RS3"],
                 "CH3": ["T0", "RS4"], "CalOffsetCH0G1": ["T1", "RCA"], "CalOffsetCH0G2": ["T2", "RCA"],
                 "CalOffsetCH0G3": ["T3", "RCA"], "CalOffsetCH1G1": ["T1", "RCB"], "CalOffsetCH1G2": ["T2", "RCB"],
                 "CalOffsetCH1G3": ["T3", "RCB"], "CalOffsetCH2G1": ["T1", "RCC"], "CalOffsetCH2G2": ["T2", "RCC"],
                 "CalOffsetCH2G3": ["T3", "RCC"], "CalOffsetCH3G1": ["T1", "RCD"], "CalOffsetCH3G2": ["T2", "RCD"],
                 "CalOffsetCH3G3": ["T3", "RCD"], "CalStepCH0G1": ["T4", "RCA"], "CalStepCH0G2": ["T5", "RCA"],
                 "CalStepCH0G3": ["T6", "RCA"], "CalStepCH1G1": ["T4", "RCB"], "CalStepCH1G2": ["T5", "RCB"],
                 "CalStepCH1G3": ["T6", "RCB"], "CalStepCH2G1": ["T4", "RCC"], "CalStepCH2G2": ["T5", "RCC"],
                 "CalStepCH2G3": ["T6", "RCC"], "CalStepCH3G1": ["T4", "RCD"], "CalStepCH3G2": ["T5", "RCD"],
                 "CalStepCH3G3": ["T6", "RCD"], "DownholeDegC": ["T0", "RKD"], "LoggerDegC": ["T0", "RK1"],
                 "PowerBoxDegC": ["T0", "RK2"], "PressureKPa": ["TS", "RDO"], "RTSettingCH0": ["T0", "RCA"],
                 "RTSettingCH1": ["T0", "RCB"], "RTSettingCH2": ["T0", "RCC"], "RTSettingCH3": ["T0", "RCD"],
                 "Rainfallmm": ["TS", "RRO"], "SolarAmps": ["T0", "REO"], "SystemAmps": ["T0", "RE2"],
                 "CalOffsetCH0G0": ["T7", "RCA"], "CalOffsetCH1G0": ["T7", "RCB"], "CalOffsetCH2G0": ["T7", "RCC"],
                 "CalOffsetCH3G0": ["T7", "RCD"], "CalStepCH0G0": ["T8", "RCA"], "CalStepCH1G0": ["T8", "RCB"],
                 "CalStepCH2G0": ["T8", "RCC"], "CalStepCH3G0": ["T8", "RCD"]}

        self.file_metadata['seed_loc'] = codes[self.file_metadata['channel']][0]
        self.file_metadata['seed_ch'] = codes[self.file_metadata['channel']][1]
        
        if self.file_metadata['channel'].startswith("CH"):
            if self.file_metadata['session'] == "Min":
                self.file_metadata['seed_ch'] = self.file_metadata['seed_ch'].replace("R","B")
            elif self.file_metadata['session'] == "Hour":
                self.file_metadata['seed_ch'] = self.file_metadata['seed_ch'].replace("R","L")


    def read_header(self, print_it=False):
        """
        """
        self.parse_filename()
        self.file.seek(0)
        data = self.file.read(40)
        format = "%shhidfiiiii" % self.file_metadata['endian']
        (self.file_metadata['magic'],
         self.file_metadata['unused'],
         self.file_metadata['header_size'],
         self.file_metadata['start'],
         self.file_metadata['interval'],
         self.file_metadata['num_pts'],
         self.file_metadata['data_type'],
         self.file_metadata['missing'],
         self.file_metadata['usgs_lock'],
         self.file_metadata['id']
         ) = struct.unpack(format, data)
        self.file_metadata['start_timestamp'] = datetime.datetime.isoformat(
            BTL_EPOCH + datetime.timedelta(0, self.file_metadata['start']), ' ')
        if self.file_metadata['data_type'] == 2:
            self.file_metadata['seed_scale_factor'] = 10000
        else:
            self.file_metadata['seed_scale_factor'] = 1
        if print_it:
            print("")
            print("file:        %s" % self.file_metadata['filename'])
            print("magic:       %x" % self.file_metadata['magic'])
            print("unused:      %d" % self.file_metadata['unused'])
            print("header_size: %d" % self.file_metadata['header_size'])
            print("start:       %s" % self.file_metadata['start_timestamp'])
            print("interval:    %g" % self.file_metadata['interval'])
            print("num_pts:     %d" % self.file_metadata['num_pts'])
            if self.file_metadata['data_type'] in list(BTL_TYPE.keys()):
                print("data_type:   %s" % BTL_TYPE[self.file_metadata['data_type']])
            else:
                print("data_type:   %d" % self.file_metadata['data_type'])
            print("missing:     %d" % self.file_metadata['missing'])
            print("usgs_lock:   %d" % self.file_metadata['usgs_lock'])
            print("id:          %d" % self.file_metadata['id'])
            print("")

        return self.file_metadata

    def read_data(self, print_it=False, timestamps=True):
        """
        """
        if print_it:
            if timestamps:
                print_format = "%%s %%%s" % BTL_TYPE_TO_PRINT[self.file_metadata['data_type']]
                # We use the python timedelta class, since it can
                # better represent intervals such as 0.05.  this becomes
                # (0,0,500000) rather than the floating point approximation
                # of 0.05 which prints as 0.050000000000000003.
                seconds = int(self.file_metadata['interval'])
                microseconds = int((self.file_metadata['interval'] - seconds) * 1000000)
                interval = datetime.timedelta(0, seconds, microseconds)
                start_offset = datetime.timedelta(0, self.file_metadata['start'])
            else:
                print_format = "%%%s" % BTL_TYPE_TO_PRINT[self.file_metadata['data_type']]

        pack_format = "%s%s" % (self.file_metadata['endian'], BTL_TYPE_TO_PACK[self.file_metadata['data_type']])
        self.file.seek(self.file_metadata['header_size'])
        self.data = []
        for i in range(self.file_metadata['num_pts']):
            datum = self.file.read(BTL_TYPE_SIZE[self.file_metadata['data_type']])
            if datum == '':
                print("error: read() returned empty string?")
                continue
            if len(datum) != BTL_TYPE_SIZE[self.file_metadata['data_type']]:
                print("error: read() returned unexpected number of bytes")
                continue
            value = struct.unpack(pack_format, datum)
            self.data.append(value[0])
            if print_it:
                if timestamps:
                    print(print_format % (
                        datetime.datetime.isoformat(BTL_EPOCH + start_offset + (i * interval), ' '), value[0]))
                else:
                    print(print_format % (value[0]))

        return self.data

    def get_unix_ms_timestamps(self):
        interval = round(self.file_metadata['interval'], 3)
        timestamps = []
        for i in range(self.file_metadata['num_pts']):
            timestamps.append(int((self.file_metadata['start'] + (i * interval)) * 1000))
        return timestamps

    def get_timestamps(self):
        seconds = int(self.file_metadata['interval'])
        microseconds = int((self.file_metadata['interval'] - seconds) * 1000000)
        interval = datetime.timedelta(0, seconds, microseconds)
        start_offset = datetime.timedelta(0, self.file_metadata['start'])
        timestamps = []
        for i in range(self.file_metadata['num_pts']):
            timestamps.append(datetime.datetime.isoformat(BTL_EPOCH + start_offset + (i * interval)))
        return timestamps

    def get_datetime_timestamps(self):
        seconds = int(self.file_metadata['interval'])
        microseconds = int((self.file_metadata['interval'] - seconds) * 1000000)
        interval = datetime.timedelta(0, seconds, microseconds)
        start_offset = datetime.timedelta(0, self.file_metadata['start'])
        timestamps = []
        for i in range(self.file_metadata['num_pts']):
            timestamps.append(BTL_EPOCH + start_offset + (i * interval))
        return timestamps

src/earthscopestraintools/bottle.py
EarthScope/earthscopestraintools
from datetime import timedelta
import numpy as np
import obspy
from earthscopestraintools.bottletar import GtsmBottleTar
import logging

logger = logging.getLogger(__name__)
if logger.hasHandlers():
    logger.setLevel(logging.INFO)
else:
    logging.basicConfig(
        format="%(asctime)s - %(levelname)s - %(message)s", level=logging.INFO
    )

def build_stream(network: str, station: str, gbt: GtsmBottleTar, verbose=False, print_traces: bool = True):
    st = obspy.core.stream.Stream()
    try:
        gbt.load_bottles(verbose=verbose)
        try:
            for bottle in gbt.bottles:
                stats = obspy.core.trace.Stats()
                stats.station = station
                stats.network = network
                stats.location = bottle.file_metadata["seed_loc"]
                stats.channel = bottle.file_metadata["seed_ch"]
                stats.delta = bottle.file_metadata["interval"]
                bottle_num_pts = bottle.file_metadata["num_pts"]
                bottle_starttime = stats.starttime + timedelta(seconds=bottle.file_metadata["start"])
                #print(bottle.file_metadata["filename"])
                if bottle.file_metadata['seed_scale_factor'] == 10000:
                    #print("float", bottle.file_metadata['data_type'])
                    float_data = np.array(bottle.read_data(), dtype=np.float32) * bottle.file_metadata['seed_scale_factor']
                    int_data = np.array(bottle.read_data(), dtype=np.int32)
                    gap_indicies = np.where(int_data==999999)[0] 
                    np.put(float_data, gap_indicies, 999999)
                    data = np.int32(float_data)

                else:
                    #print("int", bottle.file_metadata['data_type'])
                    data = np.array(bottle.read_data(), dtype=np.int32)
                    gap_indicies = np.where(data==999999)[0] 
                
                
                #print(data)
                #print(f"{len(gap_indicies)} fill values at {gap_indicies}")
                start_indicies = []
                num_gaps = len(gap_indicies)
                if num_gaps:
                    for i, gap_index in enumerate(gap_indicies):
                        if i +1 < num_gaps:
                            if gap_index + 1 != gap_indicies[i+1]:
                                start_indicies.append(gap_index + 1)
                        else:
                            start_indicies.append(gap_index + 1)
                #print(f"{len(start_indicies)} start indicies at {start_indicies}")

                gap_times = {}
                for index in gap_indicies:
                    gap_times[index] = bottle_starttime + timedelta(seconds=index*stats.delta)
                #print(gap_times)

                start_times = [bottle_starttime]
                for index in start_indicies:
                    start_times.append(bottle_starttime + timedelta(seconds=index*stats.delta))
                #print(start_times)

                arrays = np.split(data, start_indicies)
                for i, array in enumerate(arrays):
                    tracedata = np.delete(array, np.where(array==999999)[0])
                    #print(i, tracedata)
                    stats.npts = len(tracedata)
                    if stats.npts > 0:
                        stats.starttime = start_times[i]
                        tr = obspy.core.trace.Trace(data=tracedata, header=stats)
                        st.append(tr)
            if print_traces:
                print(st.__str__(extended=True))
            return st
        except Exception as e:
            logger.error(f"{gbt.file_metadata['filename']}: {type(e)} unable to build stream from bottles")
            gbt.delete_bottles_from_disk()
            raise
    except Exception as e:
        logger.error(f"{gbt.file_metadata['filename']}: {type(e)} unable to load bottles")
        gbt.delete_bottles_from_disk()
        raise
    # gbt.delete_bottles_from_disk()
    # return None


# def build_stream(network: str, station: str, gbt: GtsmBottleTar):
#     st = obspy.core.stream.Stream()
#     gbt.load_bottles()
#     for bottle in gbt.bottles:
#         stats = obspy.core.trace.Stats()
#         stats.station = station
#         stats.network = network
#         stats.location = bottle.file_metadata["seed_loc"]
#         stats.channel = bottle.file_metadata["seed_ch"]
#         # (stats.station, channel, year, dayofyear, hour, min) = btl.parse_filename()
#         # metadata = btl.read_header()
#         stats.delta = bottle.file_metadata["interval"]
#         stats.npts = bottle.file_metadata["num_pts"]
#         stats.starttime += timedelta(seconds=bottle.file_metadata["start"])

#         data = np.array(bottle.read_data(), dtype=np.int32)
#         tr = obspy.core.trace.Trace(data=data, header=stats)
#         st.append(tr)
#     st.merge()
#     logger.info(st.__str__(extended=True))
#     gbt.delete_bottles_from_disk()
#     return st


def bottle2mseed(network, station, filename, session, verbose=False, print_traces=True, plot_traces=False):
    gbt = GtsmBottleTar(f"bottles/{filename}", session, verbose=verbose)
    try:
        st = build_stream(network, station, gbt, verbose=verbose, print_traces=print_traces)
        if plot_traces:
            st.plot()
        miniseed_filename = f"miniseed/{gbt.file_metadata['filebase']}.ms"
        st.write(miniseed_filename, format="MSEED", reclen=512)
        return st
    except Exception as e:
        logger.error(f"{gbt.file_metadata['filename']}: unable to translate to miniseed \n{e}")

src/earthscopestraintools/bottle2mseed.py
EarthScope/earthscopestraintools
from datetime import datetime, timedelta
import numpy as np
import pandas as pd
from earthscopestraintools.bottletar import GtsmBottleTar
import tiledb
from earthscopestraintools.tiledbtools import (
    StrainArray,
    RawStrainWriter,
    RawStrainReader,
)
import logging

logger = logging.getLogger(__name__)
if logger.hasHandlers():
    logger.setLevel(logging.INFO)
else:
    logging.basicConfig(
        format="%(asctime)s - %(levelname)s - %(message)s", level=logging.INFO
    )


def build_strain_buffer(gbt: GtsmBottleTar, session: str):
    '''
    Testing the doc?
    '''
    logger.info(
        f"{gbt.file_metadata['filename']}: loading strain bottles into dataframe"
    )
    bottle_dfs = []
    logger.info(
        f"{gbt.file_metadata['filename']}: contains {len(gbt.bottle_list)} bottle files"
    )
    for i, name in enumerate(gbt.bottle_list):
        bottle = gbt.load_bottle(name)
        bottle.parse_filename()
        if bottle.file_metadata["channel"] in ["CH0", "CH1", "CH2", "CH3"]:
            # logger.info(bottle.file_metadata["channel"])
            bottle.read_header()
            channel = bottle.file_metadata["channel"]
            timestamps = bottle.get_unix_ms_timestamps()
            data = bottle.read_data()
            bottle.file.close()
            d = {
                "channel": channel,
                "time": timestamps,
                "data": data,
            }
            bottle_df = pd.DataFrame(data=d)
            bottle_dfs.append(bottle_df)
    tiledb_buffer = pd.concat(bottle_dfs, axis=0).reset_index(drop=True)
    tiledb_buffer["data"] = tiledb_buffer["data"].astype(np.int32)
    return tiledb_buffer


def build_ancillary_buffer(gbt: GtsmBottleTar, session: str):
    logger.info(
        f"{gbt.file_metadata['filename']}: loading ancillary bottles into dataframe"
    )
    bottle_dfs = []
    for i, name in enumerate(gbt.bottle_list):
        bottle = gbt.load_bottle(name)
        bottle.parse_filename()
        if bottle.file_metadata["channel"] not in ["CH0", "CH1", "CH2", "CH3"]:
            logger.info(bottle.file_metadata["channel"])
            bottle.read_header()
            channel = bottle.file_metadata["channel"]
            timestamps = bottle.get_unix_ms_timestamps()
            data = bottle.read_data()
            bottle.file.close()
            d = {
                "channel": channel,
                "time": timestamps,
                "data": data,
            }
            bottle_df = pd.DataFrame(data=d)
            bottle_dfs.append(bottle_df)
    tiledb_buffer = pd.concat(bottle_dfs, axis=0).reset_index(drop=True)
    tiledb_buffer["data"] = tiledb_buffer["data"].astype(np.float64)
    return tiledb_buffer


def write_buffer(uri, buffer: pd.DataFrame, cleanup_meta=True):
    writer = RawStrainWriter(uri)
    writer.write_df_to_tiledb(buffer)
    if cleanup_meta:
        writer.array.cleanup_meta()


def check_result(buffer, uri):
    try:
        start = int(buffer.time.iloc[0])
        end = int(buffer.time.iloc[-1])
        channels = list(buffer["channel"].unique())
        period = (float(buffer.time.iloc[1]) - float(buffer.time.iloc[0])) / 1000
        reader = RawStrainReader(uri, period=period)
        result_df = reader.to_df(channels=channels, start_ts=start, end_ts=end)
        logger.info(f"\n{result_df}")
    except Exception as e:
        logger.exception(f"{type(e)}")


def bottle2tdb(
    filepath: str,
    strain_uri: str,
    ancillary_uri: str,
    session: str,
    write_it=True,
    print_it=False,
    check_it=False,
):
    gbt = GtsmBottleTar(filepath, session)
    strain_buffer = build_strain_buffer(gbt, session)
    if session.casefold() == "Day".casefold():
        ancillary_buffer = build_ancillary_buffer(gbt, session)
    gbt.delete_bottles_from_disk()

    if write_it:
        try:
            logger.info(f"{filepath}: Writing to {strain_uri}")
            write_buffer(strain_uri, strain_buffer)
        except tiledb.TileDBError as e:
            logger.error(e)
            array = StrainArray(uri=strain_uri)
            array.create(schema_type="2D_INT", schema_source="s3")
            if session.casefold() == "Min".casefold():
                period = 0.05
            elif session.casefold() == "Hour".casefold():
                period = 1
            elif session.casefold() == "Day".casefold():
                period = 600
            else:
                logger.error(
                    f"Invalid session {session}, cannot determine period for array metadata"
                )
            array.set_array_meta(
                station=gbt.file_metadata["fcid"], period=period
            )  # todo: set network meta?
            logger.info(f"{filepath}: Writing to {strain_uri}")
            write_buffer(strain_uri, strain_buffer)

        if session.casefold() == "Day".casefold():
            try:
                logger.info(f"{filepath}: Writing to {ancillary_uri}")
                write_buffer(ancillary_uri, ancillary_buffer)
            except tiledb.TileDBError as e:
                logger.error(e)
                ancillary_array = StrainArray(uri=ancillary_uri)
                ancillary_array.create(schema_type="2D_FLOAT", schema_source="s3")
                ancillary_array.set_array_meta(
                    station=gbt.file_metadata["fcid"], period=1800
                )  # todo: set network meta?
                logger.info(f"{filepath}: Writing to {ancillary_uri}")
                write_buffer(ancillary_uri, ancillary_buffer)

    if print_it:
        logger.info(f"Strain Buffer:\n{strain_buffer}")
        if session.casefold() == "Day".casefold():
            logger.info(f"Ancillary Buffer:\n{ancillary_buffer}")

    if check_it:
        logger.info(f"Reading data back from array")
        logger.info(f"Network: {array.get_network()}")
        logger.info(f"Station: {array.get_station()}")
        logger.info(f"Period: {array.get_period()}")
        check_result(strain_buffer, strain_uri)
        if session.casefold() == "Day".casefold():
            logger.info(f"Network: {ancillary_array.get_network()}")
            logger.info(f"Station: {ancillary_array.get_station()}")
            logger.info(f"Period: {ancillary_array.get_period()}")
            check_result(ancillary_buffer, ancillary_uri)

src/earthscopestraintools/bottle2tdb.py
EarthScope/earthscopestraintools
import os
import shutil
import tarfile
import datetime
from earthscopestraintools.bottle import Bottle
import logging

logger = logging.getLogger(__name__)
if logger.hasHandlers():
    logger.setLevel(logging.INFO)
else:
    logging.basicConfig(
        format="%(asctime)s - %(levelname)s - %(message)s", level=logging.INFO
    )


class GtsmBottleTar:
    """
    This class is designed to unpack the following 5 cases of GTSM tar file.
    1. 'Day' session.  24hrs, low rate data.  tgz of bottles.  Logger and DMC Archive format.
    2. 'Hour' session. 1hr, 1s data.  tgz of bottles.  Logger format
    3. 'Min' session.  1hr, 20hz data.  tar of 60 tgz of bottles.  Logger format
    4. 'Hour_Archive' session.  24hr, 1s data.  tar of 24 tgz of bottles.  Archive format
    5. 'Min_Archive' session.  24hr, 20hz data.  tar of 24 tar of 60 tgz of bottles.  Archive format

    Bottles are unpacked and temporarily written to filebase/bottles. 

    """

    def __init__(self, filename, session=None, fileobj=None, verbose=False):
        # requires filename, which is path to local object or original filename of fileobj
        self.file_metadata = {}
        self.file_metadata["filename"] = filename
        self.file_metadata["fcid"] = filename.split("/")[-1][0:4]
        self.file_metadata["filebase"] = filename.split("/")[-1][:-4]
        self.file_metadata["session"] = session
        if fileobj:
            self.fileobj = fileobj
        else:
            self.fileobj = None
        self.untar_files(verbose=verbose)
        self.bottle_list = self.list_bottle_dir()
        self.bottle_list.sort()

    def untar_files(self, verbose=False):
        # unpack tars as complex as  tar(tar(tgz(bottle)))
        logger.info(f"{self.file_metadata['filename']}: unpacking tar file")
        path1 = f"{self.file_metadata['filebase']}/level1"
        path2 = f"{self.file_metadata['filebase']}/level2"
        self.bottle_path = f"{self.file_metadata['filebase']}/bottles"
        bad_files = []
        if self.file_metadata["filename"].endswith(
            "tar"
        ):  # contains more tars or tgzs.  Min, Min_Archive, Hour_Archive session
            with tarfile.open(self.file_metadata["filename"], "r") as tar:
                names = tar.getnames()
                tar.extractall(path=path1)
                for name in names:
                    try:
                        if name.endswith(
                            "tar"
                        ):  # contains more tars or tgzs. Min_Archive session
                            with tarfile.open(f"{path1}/{name}", "r") as tar2:
                                names2 = tar2.getnames()
                                tar2.extractall(path=path2)
                                for name2 in names2:
                                    #print(name2)
                                    if name2.endswith(
                                        "tgz"
                                    ):  # only contains bottles. Min_Archive session
                                        with tarfile.open(
                                            f"{path2}/{name2}", "r:gz"
                                        ) as tgz2:
                                            tgz2.extractall(path=self.bottle_path)
                                    else:
                                        bad_files.append(name2)
                                        logger.error(
                                            f"{name2} expected to be a tgz but is not."
                                        )
                                shutil.rmtree(path2)
                        elif name.endswith(
                            "tgz"
                        ):  # only contains bottles.  Min or Hour_Archive session
                            #logger.info(name)
                            try:
                                with tarfile.open(f"{path1}/{name}", "r:gz") as tgz:
                                    tgz.extractall(path=self.bottle_path)
                            except Exception as e:
                                bad_files.append(name)
                                if verbose:
                                    logger.error(f"{self.file_metadata['filename']}: {type(e)} on {name}, skipping")
                    except Exception as e:
                        bad_files.append(name)
                        if verbose:
                            logger.error(f"{self.file_metadata['filename']}: {type(e)} on {name}, skipping")
                if len(bad_files):
                    logger.warning(f"{self.file_metadata['filename']}: Failed to open {len(bad_files)} tar/tgzs")
                shutil.rmtree(path1)

        elif self.file_metadata["filename"].endswith(
            "tgz"
        ):  # only contains bottles.  Day or Hour session
            try:
                with tarfile.open(self.file_metadata["filename"], "r:gz") as tgz:
                    tgz.extractall(path=self.bottle_path)
            except Exception as e:
                logger.error(f"{self.file_metadata['filename']}: {type(e)} cannot open tgz")
                raise

    def list_bottle_dir(self):
        return os.listdir(self.bottle_path)

    def load_bottles(self, verbose=False):
        # opens and adds bottlefiles to self.bottles
        self.bottles = []
        bottle_list = self.list_bottle_dir()
        bottle_list.sort()
        failed = []
        for bottlename in bottle_list:
            try:
                btl = Bottle(f"{self.bottle_path}/{bottlename}")
                btl.read_header()
                self.bottles.append(btl)
            except Exception as e:
                failed.append(bottlename)
                if verbose:
                    print(verbose)
                    logger.error(f"{self.file_metadata['filename']}: {type(e)} on {bottlename}, skipping")
        if len(failed):
            logger.warning(f"{self.file_metadata['filename']}: Failed to load {len(failed)} of {len(bottle_list)} bottles")
        if len(self.bottles) == 0:  
            raise 

    def load_bottle(self, bottlename):
        # open and returns a bottle
        return Bottle(f"{self.bottle_path}/{bottlename}")

    def delete_bottles_from_disk(self):
        # clean up everything extracted from original file
        shutil.rmtree(self.file_metadata["filebase"])

    def list_bottles(self):
        for bottle in self.bottles:
            print(bottle.file_metadata["filename"])

    def get_bottle_names(self):
        bottle_names = []
        for bottle in self.bottles:
            bottle_names.append(bottle.file_metadata["filename"])
        return bottle_names

    def test(self):
        logger.info(
            f"{self.file_metadata['filename']}: contains {len(self.bottle_list)} bottles"
        )
        self.load_bottles()
        logger.info(f"Successfully loaded {len(self.bottles)} bottles")


#
# if __name__ == "__main__":
#
#     # function for running local tests
#
#     t1 = datetime.datetime.now()
#
#     filename = "B001.2022001Day.tgz"  # 24 hr Day session (archive and logger format)
#     session = "Day"
#     # filename = 'B018.2016366_01.tar' #24 Hour, Hour Session (archive format)
#     # session = "Hour_Archive"
#     # filename = 'B018.2016366_20.tar' #24 Hour, Min session (archive format)
#     # session = "Min_Archive"
#     # filename = 'B0012200100.tgz'  #1 Hour, Hour Session (logger format)
#     # session = "Hour"
#     # filename = 'B0012200100_20.tar' #1 Hour, Min Session (logger format)
#     # session = "Min"
#     gbt = GtsmBottleTar(f"bottles/{filename}", session)
#     gbt.load_bottles()
#     for bottle in gbt.bottles:
#         print(bottle.file_metadata["filename"])
#     # print(bottle.get_unix_ms_timestamps())
#     gbt.delete_bottles_from_disk()
#
#     t2 = datetime.datetime.now()
#     elapsed_time = t2 - t1
#     logger.info(f'{gbt.file_metadata["filename"]}: Elapsed time {elapsed_time} seconds')
#     # print(gbt.file_metadata['filename'], gbt.file_metadata['fcid'], gbt.file_metadata['session'])

src/earthscopestraintools/bottletar.py
EarthScope/earthscopestraintools
# look up edids without dependecy on es_datasources_client

import os
import requests
import json
import logging

logger = logging.getLogger(__name__)

def get_network_edid(
    api_url: str,
    namespace: str,
    network: str,
    with_children_count=False,
    with_edid_only=True,
    with_parents=False,
    with_slugs=False,
    timeout: int = 10
):
    """Returns a single network edid"""

    parameters = {
        "network_name": f"{namespace}:{network}",
        "with_children_count": with_children_count,
        "with_edid_only": with_edid_only,
        "with_parents": with_parents,
        "with_slugs": with_slugs,
    }
    NETWORK_EDID_PATH = os.path.join(api_url, "network")
    try:
        # async with httpx.AsyncClient() as client:
        #     r = await client.get(NETWORK_EDID_PATH, 
        #                         params=parameters,
        #                         timeout=10)
        r = requests.get(NETWORK_EDID_PATH, params=parameters, timeout=timeout)
        # print(r.url)
        r.raise_for_status()

    except requests.exceptions.HTTPError as errh:
        logger.error(f'Http error: {network} {json.loads(r.content)["detail"]}')
        raise

    except requests.exceptions.ConnectionError as errc:
        logger.error(f"Error Connecting: {network} {errc}")
        raise

    except requests.exceptions.Timeout as errt:
        logger.error(f"Timeout Error: {network} {errt}")
        raise

    except requests.exceptions.RequestException as err:
        logger.error(f"Oops: Something Else: {network} {err}")
        raise

    except Exception as err:
        logger.error("Unknown error: {}".format(err))
        raise

    if len(r.json()):
        return r.json()["items"][0]
    else:
        logger.warning("No items returned")
        return None


def get_station_edid(
    api_url: str,
    namespace: str,
    station: str,
    with_children_count=False,
    with_edid_only=True,
    with_parents=False,
    with_slugs=False,
    timeout: int = 10
):
    """Returns a station edid in the specified station namespace"""

    parameters = {
        "station_name": f"{namespace}:{station}",
        "with_children_count": with_children_count,
        "with_edid_only": with_edid_only,
        "with_parents": with_parents,
        "with_slugs": with_slugs,
    }
    STATION_EDID_PATH = os.path.join(api_url, "station")
    try:
        # async with httpx.AsyncClient() as client:
        #     r = await client.get(STATION_EDID_PATH, 
        #                         params=parameters,
        #                         timeout=10)
        r = requests.get(STATION_EDID_PATH, params=parameters, timeout=timeout)
        # print(r.url)
        r.raise_for_status()

    except requests.exceptions.HTTPError as errh:
        logger.error(f'Http error: {station} {json.loads(r.content)["detail"]}')
        raise

    except requests.exceptions.ConnectionError as errc:
        logger.error(f"Error Connecting: {station} {errc}")
        raise

    except requests.exceptions.Timeout as errt:
        logger.error(f"Timeout Error: {station} {errt}")
        raise

    except requests.exceptions.RequestException as err:
        logger.error(f"Oops: Something Else: {station} {err}")
        raise

    except Exception as err:
        logger.error("Unknown error: {}".format(err))
        raise

    if len(r.json()):
        return r.json()["items"][0]
    else:
        logger.warning("No items returned")
        return None


def get_session_edid(
    api_url: str,
    namespace: str,
    station: str,
    session: str,
    with_children_count=False,
    with_edid_only=True,
    with_parents=False,
    with_slugs=False,
    timeout: int = 10
):
    """Returns a session edid in the specified station namespace"""

    parameters = {
        "station_name": f"{namespace}:{station}",
        "session_name": f"DATAFLOW:{session}",
        "with_children_count": with_children_count,
        "with_edid_only": with_edid_only,
        "with_parents": with_parents,
        "with_slugs": with_slugs,
    }
    SESSION_EDID_PATH = os.path.join(api_url, "session")
    try:
        # async with httpx.AsyncClient() as client:
        #     r = await client.get(SESSION_EDID_PATH, 
        #                         params=parameters,
        #                         timeout=10)
        r = requests.get(SESSION_EDID_PATH, params=parameters, timeout=timeout)
        r.raise_for_status()

    except requests.exceptions.HTTPError as errh:
        logger.error(f'Http error: {station} {json.loads(r.content)["detail"]}')
        raise

    except requests.exceptions.ConnectionError as errc:
        logger.error(f"Error Connecting: {station} {errc}")
        raise

    except requests.exceptions.Timeout as errt:
        logger.error(f"Timeout Error: {station} {errt}")
        raise

    except requests.exceptions.RequestException as err:
        logger.error(f"Oops: Something Else: {station} {err}")
        raise

    except Exception as err:
        logger.error("Unknown error: {}".format(err))
        raise

    if len(r.json()):
        return r.json()["items"][0]
    else:
        logger.warning("No items returned")
        return None


def get_all_stations_in_network(
    api_url: str,
    namespace: str,
    network: str,
    with_children_count=False,
    with_edid_only=False,
    with_parents=False,
    with_slugs=False,
    limit=1000,
    timeout: int = 10
):
    """Returns a dictionary of station:edid for a given namespace:network"""

    parameters = {
        "network_name": f"{namespace}:{network}",
        "with_children_count": with_children_count,
        "with_edid_only": with_edid_only,
        "with_parents": with_parents,
        "with_slugs": with_slugs,
        "limit": limit,
    }
    STATION_EDID_PATH = os.path.join(api_url, "station")
    try:
        # async with httpx.AsyncClient() as client:
        #     r = await client.get(STATION_EDID_PATH, 
        #                         params=parameters,
        #                         timeout=10)
        r = requests.get(STATION_EDID_PATH, params=parameters, timeout=timeout)
        # print(r.url)
        r.raise_for_status()

    except requests.exceptions.HTTPError as errh:
        logger.error(f'Http error: {network} {json.loads(r.content)["detail"]}')
        raise

    except requests.exceptions.ConnectionError as errc:
        logger.error(f"Error Connecting: {network} {errc}")
        raise

    except requests.exceptions.Timeout as errt:
        logger.error(f"Timeout Error: {network} {errt}")
        raise

    except requests.exceptions.RequestException as err:
        logger.error(f"Oops: Something Else: {network} {err}")
        raise

    except Exception as err:
        logger.error("Unknown error: {}".format(err))
        raise

    if len(r.json()):
        stations = {}
        more_pages = True
        while more_pages:
            more_pages = r.json()["has_next"]

            for item in r.json()["items"]:
                for name in item["names"]:
                    if namespace == "BSM":
                        if "BNUM" in name:
                            stations[name.split(":")[-1]] = item["edid"]
                    elif namespace == "FDSN":
                        if "FDSN" in name:
                            stations[name.split(":")[-1]] = item["edid"]
                    elif namespace == "PERM":
                        if "PNUM" in name:
                            stations[name.split(":")[-1]] = item["edid"]
            parameters["offset"] = r.json()["offset"] + r.json()["limit"]

            r = requests.get(STATION_EDID_PATH, params=parameters, timeout=10)
            r.raise_for_status()

        return dict(sorted(stations.items()))

    else:
        logger.warning("No items returned")
        return None


def get_all_sessions_in_network(
    api_url: str,
    namespace: str,
    network: str,
    with_children_count=False,
    with_edid_only=False,
    with_parents=True,
    with_slugs=False,
    limit=1000,
    timeout: int = 10
):
    """Returns a nested dictionary of each session and its parent station/network for a given namespace:network"""

    parameters = {
        "network_name": f"{namespace}:{network}",
        "with_children_count": with_children_count,
        "with_edid_only": with_edid_only,
        "with_parents": with_parents,
        "with_slugs": with_slugs,
        "limit": limit,
    }
    SESSION_EDID_PATH = os.path.join(api_url, "session")
    try:
        stations = {}
        more_pages = True

        while more_pages:
            # async with httpx.AsyncClient() as client:
            #     r = await client.get(SESSION_EDID_PATH, 
            #                     params=parameters,
            #                     timeout=10)
            r = requests.get(SESSION_EDID_PATH, params=parameters, timeout=timeout)
            #print(r.url)
            r.raise_for_status()
            more_pages = r.json()["has_next"]

            for item in r.json()["items"]:
                session_edid = item["edid"]
                session_name = item["names"][0].split(":")[-1]
                for name in item["station"]["names"]:
                    if namespace == "BSM":
                        if "BNUM" in name:
                            station_name = name.split(":")[-1]
                            station_edid = item["station"]["edid"]
                            if station_name not in stations:
                                stations[station_name] = {"edid": station_edid}

                    elif namespace == "FDSN":
                        if "FDSN" in name:
                            station_name = name.split(":")[-1].split("_")[-1]
                            station_edid = item["station"]["edid"]
                            if station_name not in stations:
                                stations[station_name] = {"edid": station_edid}

                    elif namespace == "PERM":
                        if "PNUM" in name:
                            station_name = name.split(":")[-1]
                            station_edid = item["station"]["edid"]
                            if station_name not in stations:
                                stations[station_name] = {"edid": station_edid}

                stations[station_name][session_name] = session_edid

            if more_pages:
                logger.info(
                    f"get_all_sessions_in_network request exceeded {limit} responses, "
                    f"returned only partial results, grabbing more.."
                )
            else:
                logger.info("Collected all results.")

            parameters["offset"] = r.json()["offset"] + r.json()["limit"]

        return dict(sorted(stations.items()))

    except requests.exceptions.HTTPError as errh:
        logger.error(f'Http error: {network} {json.loads(r.content)["detail"]}')
        raise

    except requests.exceptions.ConnectionError as errc:
        logger.error(f"Error Connecting: {network} {errc}")
        raise

    except requests.exceptions.Timeout as errt:
        logger.error(f"Timeout Error: {network} {errt}")
        raise

    except requests.exceptions.RequestException as err:
        logger.error(f"Oops: Something Else: {network} {err}")
        raise

    except Exception as err:
        logger.error("Unknown error: {}".format(err))
        raise


def get_bnum_list(api_url: str, 
                  namespace="BSM", 
                  network="NOTA",
                  timeout: int = 10):
    #returns all bnums with GTSM instruments (explicitly removes MC dilatometers)
    response = get_all_stations_in_network(api_url=api_url, namespace=namespace, network=network, timeout=timeout)
    try:
        bnums = sorted(list(response.keys()))
        mc = ["AIRS", "OLV1", "OLV2", "TRNT"]
        return [x for x in bnums if x not in mc]

    except Exception:
        logger.error("No stations found")
        return []


def lookup_fdsn_network_from_bnum(api_url: str,
                                  station: str,
                                  timeout: int = 10):
    parameters = {
        "station_name": f"BNUM:{station}",
        "with_parents": True,
    }
    STATION_EDID_PATH = os.path.join(api_url, "station")
    try:
        # async with httpx.AsyncClient() as client:
        #     r = await client.get(STATION_EDID_PATH, 
        #                         params=parameters,
        #                         timeout=10)
        r = requests.get(STATION_EDID_PATH, params=parameters, timeout=timeout)
        # print(r.url)
        r.raise_for_status()

    except requests.exceptions.HTTPError as errh:
        logger.error(f'Http error: {station} {json.loads(r.content)["detail"]}')
        raise

    except requests.exceptions.ConnectionError as errc:
        logger.error(f"Error Connecting: {station} {errc}")
        raise

    except requests.exceptions.Timeout as errt:
        logger.error(f"Timeout Error: {station} {errt}")
        raise

    except requests.exceptions.RequestException as err:
        logger.error(f"Oops: Something Else: {station} {err}")
        raise

    except Exception as err:
        logger.error("Unknown error: {}".format(err))
        raise

    if len(r.json()):
        networks = r.json()["items"][0]["networks"]
        for network in networks:
            for name in network["names"]:
                if "FDSN" in name:
                    fdsn_name = name.split(":")[-1]
                    return fdsn_name
        return None
    else:
        logger.warning("No items returned")
        return None
src/earthscopestraintools/datasources_api_interact.py
EarthScope/earthscopestraintools
import datetime
import urllib.request, json
import pytz


def version():
    print("Earthquake metadata module")
    print("Version 1.3.  M. Gottlieb 3-27-2023")


class Earthquake:
    """
    Class object containing parameters for a specific earthquake event
    """
    def __init__(self, event_id, focal_mechanism=False):
        self.event_id = event_id

        self.load_url_data()
        self.name = self.data["properties"]['title']
        self.mag = self.data["properties"]['mag']
        self.unix_time = self.data['properties']['time']
        self.time = datetime.datetime.fromtimestamp(self.unix_time / 1000.0).astimezone(pytz.utc)
        self.lat = self.data["geometry"]['coordinates'][1]
        self.long = self.data["geometry"]['coordinates'][0]
        self.depth = self.data["geometry"]['coordinates'][2]
        if focal_mechanism:
            self.get_focal_mechanism()


    def load_url_data(self):
        url = "https://earthquake.usgs.gov/fdsnws/event/1/query?format=geojson&eventid=" + self.event_id
        response = urllib.request.urlopen(url)
        self.data = json.loads(response.read())

    def get_focal_mechanism(self):
        focal_mechanism_props = self.data['properties']['products']['focal-mechanism'][0]['properties']
        np1_dip = focal_mechanism_props['nodal-plane-1-dip']
        np1_rake = focal_mechanism_props['nodal-plane-1-rake']
        np1_strike = focal_mechanism_props['nodal-plane-1-strike']
        #mag = self.data["properties"]['mag']
        self.strike = np1_strike
        self.dip = np1_dip
        self.rake = np1_rake

src/earthscopestraintools/event.py
EarthScope/earthscopestraintools
import numpy as np
import pandas as pd
from scipy import stats, signal
import matplotlib.pyplot as plt
from matplotlib import cm
from geopy.distance import distance
import datetime
import math
from earthscopestraintools.gtsm_metadata import get_metadata_df

import logging

logger = logging.getLogger(__name__)

def calc_hypocentral_dist(
        eq_latitude,
        eq_longitude,
        eq_depth,
        station_latitude,
        station_longitude):
    """
    Function calculates hypocentral distance (km) between lat,long and earthquake

    :param eq_latitude: latitude of earthquake
    :type eq_latitude: float
    :param eq_longitude: longitude of earthquake
    :type eq_longitude: float
    :param eq_depth: depth of earthquake
    :type eq_depth: float
    :param station_latitude: latitude of station
    :type station_latitude: float
    :param station_longitude: longitude of station
    :type station_longitude: float
    :return: hypocentral distance in km
    :rtype: int

    """

    ed = distance((eq_latitude, eq_longitude), (station_latitude, station_longitude)).km
    hypocentral_dist = int(math.sqrt(float(ed) ** 2 + float(eq_depth) ** 2))
    return hypocentral_dist


def calculate_p_s_arrival(eq_latitude,
                          eq_longitude,
                          eq_time,
                          station_latitude,
                          station_longitude):
    """
    Function calculates arrival times for P and S waves at a given lat and long

    :param eq_latitude: latitude of earthquake
    :type eq_latitude: float
    :param eq_longitude: longitude of earthquake
    :type eq_longitude: float
    :param eq_time: time of earthquake
    :type eq_time: datetime.datetime
    :param station_latitude: latitude of station
    :type station_latitude: float
    :param station_longitude: longitude of station
    :type station_longitude: float
    :return: p_arrival, s_arrival
    :rtype: datetime.datetime

    """

    event_loc = "[" + str(eq_latitude) + "," + str(eq_longitude) + "]"
    station_loc = "[" + str(station_latitude) + "," + str(station_longitude) + "]"
    url = "https://service.iris.edu/irisws/traveltime/1/query?evloc=" + event_loc + "&staloc=" + station_loc
    df = pd.read_table(url, sep="\s+", header=1, index_col=2, usecols=[2, 3])

    p_delta = datetime.timedelta(seconds=float(df.iloc[(df.index == 'P').argmax()].Travel))
    s_delta = datetime.timedelta(seconds=float(df.iloc[(df.index == 'S').argmax()].Travel))
    p_arrival = eq_time + p_delta
    s_arrival = eq_time + s_delta
    return p_arrival, s_arrival

def dynamic_strain(df, gauge_weights=[1, 1, 1, 1]):
    """Calculates dynamic strain as RMS of gauge strains

    :param df: dataframe containing gauge strains as columns and time as an index
    :type df: pandas.DataFrame
    :param gauge_weights: list of which gauges to include, defaults to [1, 1, 1, 1]
    :type gauge_weights: list, optional
    :return: dataframe containing a single column of dynamic strain, time as an index
    :rtype: pandas.DataFrame
    """
    logger.info(f"Calculating dynamic strain using gauge weights: {gauge_weights}")
    ser = np.sqrt(
        (
            np.square(df.CH0) * gauge_weights[0]
            + np.square(df.CH1) * gauge_weights[1]
            + np.square(df.CH2) * gauge_weights[2]
            + np.square(df.CH3) * gauge_weights[3]
        )
        / sum(gauge_weights)
    )
    return ser.to_frame(name="dynamic")


def pre_event_trend_correction(df, eq_time):
    """calculate a linear trend correction based on any data provided prior to event start time

    :param df: dataframe containing strain data
    :type df: pandas.DataFrame
    :param eq_time: time of event 
    :type eq_time: datetime.datetime
    :return: dataframe containing a linear trend correction
    :rtype: pandas.DataFrame
    """
    df_trend_c = pd.DataFrame(data=df.index)
    df_pre = pd.DataFrame(data=df.index[df.index < eq_time])
    # print(df_pre.index)
    for ch in df.columns:
        # print(df[ch][df.index < eq_time])
        slope, intercept, r_value, p_value, std_err = stats.linregress(
            df_pre.index, df[ch][df.index < eq_time].interpolate()
        )
        # print(slope)
        df_trend_c[ch] = df_trend_c.index * slope

    return df_trend_c[df.columns].set_index(df_trend_c["time"])

def calculate_magnitude(dynamic_strain_df, hypocentral_distance, site_term, longitude_term):
    """Calculates a magnitude estimate based on Barbour et al 2021.

    :param dynamic_strain_df: dataframe containing dynamic strain during an event
    :type dynamic_strain_df: pandas.DataFrame
    :param hypocentral_distance: distance from station to event hypocenter, in km
    :type hypocentral_distance: float
    :param site_term: site term from Barbour et al 2021
    :type site_term: float
    :param longitude_term: longitude term from Barbour et al 2021
    :type longitude_term: float
    :return: dataframe containing a dynamic strain based magnitude estimate as a function of time
    :rtype: pandas.DataFrame
    """
    logger.info(f"Calculating magnitude from dynamic strain using site term {site_term} "
                f"and longitude term {longitude_term}")
    df2 = pd.DataFrame(index=dynamic_strain_df.index)
    df2['magnitude'] = ((np.log10(dynamic_strain_df['dynamic'].cummax() / 1000000)) +
                             (0.00072 * hypocentral_distance) +
                             (1.45 * math.log10(hypocentral_distance)) +
                             8.52 - longitude_term - site_term
                             ) / 0.92
    return df2

def plot_coseismic_offset(
    df,
    title: str = "",
    remove_9s: bool = True,
    zero: bool = False,
    detrend: bool = None,
    ymin: float = None,
    ymax: float = None,
    plot_type: str = "scatter",
    units: str = None,
    eq_time: datetime.datetime = None,
    coseismic_offset: bool = False,
    color="black",
    save_as: str = None,
):
    """plot strain data from an event 

    :param df: strain data including an event signal
    :type df: pandas.DataFrame
    :param title: plot title, defaults to ""
    :type title: str, optional
    :param remove_9s: option to remove gap fill values, defaults to True
    :type remove_9s: bool, optional
    :param zero: option to zero the data against the first valid index, defaults to False
    :type zero: bool, optional
    :param detrend: option to linearly detrend data, will use pre-event data only if eq_time is provided, defaults to None
    :type detrend: bool, optional
    :param ymin: y-axis minimum for plot, defaults to None
    :type ymin: float, optional
    :param ymax: y-axis maximum for plot, defaults to None
    :type ymax: float, optional
    :param plot_type: matplotlib plot type. option of ['scatter','line'], defaults to "scatter"
    :type plot_type: str, optional
    :param units: units to display on y-axis, defaults to None
    :type units: str, optional
    :param eq_time: origin time of event, defaults to None
    :type eq_time: datetime.datetime, optional
    :param coseismic_offset: option to calculate and display the coseismic offset as a difference of the mean of the first quintile and last quintile of data, defaults to False
    :type coseismic_offset: bool, optional
    :param color: matplotlib color option, defaults to "black"
    :type color: str, optional
    :param save_as: filename to save plot, defaults to None
    :type save_as: str, optional
    """
    fig, axs = plt.subplots(len(df.columns), 1, figsize=(12, 10), squeeze=False)

    if remove_9s:
        df = df.replace(999999, np.nan)
    if zero:
        title += " Zeroed"
        df -= df.loc[df.first_valid_index()]
    if detrend:
        df = df.interpolate()
        if detrend == "linear":
            if eq_time:
                title += " Linearly Detrended on Pre-Event Data"
                trend_c = pre_event_trend_correction(
                    df - df.loc[df.first_valid_index()], eq_time
                )
                df -= trend_c
            else:
                title += " Linearly Detrended"
                for ch in df.columns:
                    df[ch] = signal.detrend(df[ch], type="linear")
        else:
            print("Only linear detrend implemented")
    if coseismic_offset:
        # use first and last quintile
        df_pre = df[
            (df.index > np.percentile(df.index, 0))
            & (df.index <= np.percentile(df.index, 20))
        ]
        df_post = df[
            (df.index > np.percentile(df.index, 20))
            & (df.index <= np.percentile(df.index, 100))
        ]
        coseismic_offsets = df_post.mean() - df_pre.mean()

    for i, ch in enumerate(df.columns):
        if plot_type == "line":
            axs[i][0].plot(df[ch], color=color, label=ch)
        elif plot_type == "scatter":
            axs[i][0].scatter(df.index, df[ch], color=color, s=2, label=ch)
        else:
            print("Plot type must be either 'line' or 'scatter'")
        #         if self.units:
        #             axs[i][0].set_ylabel(self.units)
        if coseismic_offset:
            label = f"Co-seismic offset: {round(coseismic_offsets[i], 3)} ustrain"
            axs[i][0].text(
                0.85,
                0.2,
                label,
                horizontalalignment="center",
                verticalalignment="center",
                transform=axs[i][0].transAxes,
            )
        if ymin or ymax:
            axs[i][0].set_ylim(ymin, ymax)
        if units:
            axs[i][0].set_ylabel(units)
        axs[i][0].ticklabel_format(axis="y", useOffset=False, style="plain")
        axs[i][0].grid()
        axs[i][0].legend()
    if title:
        fig.suptitle(title)
    fig.tight_layout()
    if save_as:
        print(f"Saving plot to {save_as}")
        plt.savefig(save_as, facecolor="white", transparent=False)



def magnitude_plot(dynamic_strain_df: pd.DataFrame,
                   magnitude_df:pd.DataFrame,
                   eq_time: datetime.datetime,
                   eq_mag: datetime.datetime,
                   title: str=None,
                   save_as: str = None,):
    """plot dynamic strain and associated magnitude estimate on the same plot

    :param dynamic_strain_df: dataframe containing dynamic strain during an event
    :type dynamic_strain_df: pd.DataFrame
    :param magnitude_df: dataframe containing strain-based magnitude estimate during an event 
    :type magnitude_df: pd.DataFrame
    :param eq_time: event origin time
    :type eq_time: datetime.datetime
    :param eq_mag: published event magnitude from COMCAT 
    :type eq_mag: datetime.datetime
    :param title: plot title, defaults to None
    :type title: str, optional
    :param save_as: filename to save plot, defaults to None
    :type save_as: str, optional
    """
    num_colors = 4
    fig, ax = plt.subplots(figsize=(12, 3))
    colors = [cm.gnuplot(x) for x in np.linspace(0, 0.8, num_colors)]
    if title:
        fig.suptitle(title)
    ax.plot(dynamic_strain_df['dynamic'], color=colors[0], label='microstrain')
    ax2 = ax.twinx()
    ax2.plot(magnitude_df['magnitude'], color=colors[1],
             label=f'est magnitude: {round(magnitude_df["magnitude"].iloc[-1], 2)}')
    ax.axvline(eq_time, color=colors[2], label='earthquake time')
    ax2.axhline(eq_mag, color=colors[3], label=f'usgs magnitude: {eq_mag}')
    ax.set_ylabel('microstrain')
    ax2.set_ylabel('magnitude')
    lines, labels = ax.get_legend_handles_labels()
    lines2, labels2 = ax2.get_legend_handles_labels()
    ax2.legend(lines + lines2, labels + labels2, loc='right')
    fig.tight_layout()
    if save_as:
        logger.info(f"Saving plot to {save_as}")
        plt.savefig(save_as)

def get_stations_in_radius(latitude, longitude, depth, radius, print_it=False):
    """determine list of stations within a given radius from an event, using lat/long in gtsm metadata summary table

    :param latitude: event latitude
    :type latitude: float
    :param longitude: event longitude
    :type longitude: float
    :param depth: event depth
    :type depth: float
    :param radius: radius from event, in km
    :type radius: float
    :param print_it: option to print list of stations within radius, and calculated distances from event, defaults to False
    :type print_it: bool, optional
    :return: list of four character station codes for stations within the given radius
    :rtype: list
    """
    meta_df = get_metadata_df()
    station_list = []
    for station in meta_df.index:
        dist = calc_hypocentral_dist(latitude,
                                     longitude,
                                     depth,
                                         meta_df.loc[station]["LAT"],
                                         meta_df.loc[station]["LONG"])
        if dist <= radius:
            if print_it:
                print(f"{station} at {dist} km")
            station_list.append(station)
    return station_list
src/earthscopestraintools/event_processing.py
EarthScope/earthscopestraintools
import numpy as np
import pandas as pd
import math
import urllib.request as request

from importlib.resources import files
from datetime import datetime
from earthscopestraintools.datasources_api_interact import lookup_fdsn_network_from_bnum
import requests

import logging

logger = logging.getLogger(__name__)


class GtsmMetadata:
    """
    Class for loading and storing GTSM processing metadata for a given station.  Metadata sources include 
    a summary table \n
    https://www.unavco.org/data/strain-seismic/bsm-data/lib/docs/bsm_metadata.txt \n
    and individual station pages ie \n
    http://bsm.unavco.org/bsm/level2/B001/B001.README.txt 
    
    :param network: 2 character FDSN network code
    :type network: str
    :param station: 4 character FDSN station code
    :type station: str
    :param gauge_weights: optional list of 'good' channels, where 1 is good and 0 is bad, defaults to [1, 1, 1, 1]
    :type gauge_weights: list, optional
    
    :ivar network: str
    :ivar station: str
    :ivar latitude: float
    :ivar longitude: float
    :ivar elevation: float
    :ivar gap: float, instrument gap in meters
    :ivar diameter: float, instrument sensing diameter in meters
    :ivar start_date: str, formatted as "%Y-%m-%d"
    :ivar orientation: float, degrees East of North for CH0
    :ivar reference_strains: dict, containing 'linear_date':'YYYY:DOY' and each channel ie 'CH0':reference counts 
    :ivar strain_matrices: dict, contains one or more calibration matrices, keyed to the name of the calibration.  
    :ivar atmp_response: dict, reponse coefficients for each channel
    :ivar tidal_params: dict, keys are tuple of (channel, tidal constituent, phz/amp/doodson)

    """
    def __init__(self, network, station, gauge_weights=None):
        """
        Loads metadata for given station.  
        """
        self.network = network
        self.station = station  
        self.meta_df = self.get_meta_table()
        self.latitude = self.get_latitude()
        self.longitude = self.get_longitude()
        self.elevation = self.get_elevation()
        self.gap = self.get_gap()
        self.diameter = 0.087
        self.orientation = self.get_orientation()
        self.reference_strains = self.get_reference_strains()
        self.start_date = self.get_start_date()
        matrices = {}
        if not gauge_weights:
            self.gauge_weights = [1, 1, 1, 1]
        # matrices['weighted_strain_matrix'] = self.make_weighted_strain_matrix(gauge_weights=gauge_weights)
        matrices["lab"] = self.get_lab_strain_matrix()
        matrices["ER2010"] = self.get_er2010_strain_matrix()
        matrices["CH_prelim"] = self.get_ch_prelim_strain_matrix()
        self.strain_matrices = {k: v for k, v in matrices.items()}
        self.atmp_response = self.get_atmp_response()
        self.tidal_params = self.get_tidal_params()

    #        self.xml = self.get_xml()
    #        self.detrend = self.get_detrend_xml()

    # def get_xml(self):
    #     metadir = "../xml/"
    #     xml_path = "ftp://bsm.unavco.org/pub/bsm/level2/" + self.station + "/"
    #     xml_file = self.station + ".xml"
    #     with closing(request.urlopen(xml_path + xml_file)) as r:
    #         with open(metadir + xml_file, "wb") as f:
    #             shutil.copyfileobj(r, f)
    #     return xmltodict.parse(open(metadir + xml_file).read(), process_namespaces=True)

    def get_meta_table(self):
        """loads summary metadata table into dataframe

        :return: basic metadata for all strainmeter stations
        :rtype: pandas.DataFrame
        """
        url = "https://www.unavco.org/data/strain-seismic/bsm-data/lib/docs/bsm_metadata.txt"
        meta_df = pd.read_csv(url, sep="\s+", index_col="BNUM")
        return meta_df

    def get_latitude(self):
        """parse latitude from metadata table

        :return: station latitude
        :rtype: float
        """
        try:
            lats = self.meta_df["LAT"]
            return float(lats[self.station])
        except:
            logger.info("no latitude found")
            exit(1)

    def get_longitude(self):
        """parse longitude from metadata table

        :return: station longitude
        :rtype: float
        """
        try:
            longs = self.meta_df["LONG"]
            return float(longs[self.station])
        except:
            logger.info("no longitude found")
            exit(1)

    def get_elevation(self):
        """parse elevation from metadata table

        :return: station elevation
        :rtype: float
        """
        try:
            elevations = self.meta_df["ELEV(m)"]
            return float(elevations[self.station])
        except:
            logger.info("no elevation found")
            exit(1)

    def get_gap(self):
        """parse instrument gap from metadata table

        :return: station instument gap in meters
        :rtype: float
        """
        try:
            gaps = self.meta_df["GAP(m)"]
            return float(gaps[self.station])
        except:
            logger.info("no gap found for %s, using .0001" % self.station)
            return 0.0001

    def get_orientation(self):
        """parse instrument orientation from metadata table

        :return: orientation of CH0 in degrees East of North
        :rtype: float
        """
        try:
            orientations = self.meta_df["CH0(EofN)"]
            return float(orientations[self.station])
        except Exception as e:
            logger.error(e)
            logger.error("No orientation found for %s, using 0 deg" % self.station)
            return 0

    def get_start_date(self):
        """parse data start date from metadata table

        :return: start date string as "%Y-%m-%d"
        :rtype: str
        """
        try:
            start_date = self.meta_df.loc[self.station]["DATA_START"]
            return datetime.strptime(start_date, "%Y:%j").strftime("%Y-%m-%d")

        except Exception as e:
            logger.error(e)
            return None

    # def get_orientation_xml(self):
    #     for i, dic in enumerate(
    #         self.xml["strain_xml"]["inst_info"]["sensor_information"]["sensor_response"]
    #     ):
    #         if dic["sensor_type"] == "Gladwin_BSM_component_1_":
    #             orientation = float(
    #                 self.xml["strain_xml"]["inst_info"]["sensor_information"][
    #                     "sensor_response"
    #                 ][i]["orientation"]["#text"]
    #             )
    #     try:
    #         return orientation
    #     except Exception as e:
    #         logger.error(e)
    #         logger.error("No orientation found for %s, using 0 deg" % self.station)
    #         return 0

    def make_weighted_strain_matrix(self, gauge_weights=[1, 1, 1, 1]):
        """calculates lab strain matrix, and allows exclusion of a single gauge

        :param gauge_weights: which gauges to use, set to 0 to exclude one gauge, defaults to [1, 1, 1, 1]
        :type gauge_weights: list, optional
        :return: lab calibration matrix
        :rtype: numpy.array
        """
        # make strain matrix from manufacturers coefficients
        # logger.info(gauge_weights)
        c = 1.5
        d = 3
        scale_factors = np.array([[c, 0, 0], [0, d, 0], [0, 0, d]])
        gage_weights = np.array(
            [
                [gauge_weights[0], 0, 0, 0],
                [0, gauge_weights[1], 0, 0],
                [0, 0, gauge_weights[2], 0],
                [0, 0, 0, gauge_weights[3]],
            ]
        )
        orientations = np.array(
            [
                [
                    0.5,
                    0.5 * math.cos(math.radians(2 * (90 - self.orientation))),
                    0.5 * math.sin(math.radians(2 * (90 - self.orientation))),
                ],
                [
                    0.5,
                    0.5 * math.cos(math.radians(2 * (-30 - self.orientation))),
                    0.5 * math.sin(math.radians(2 * (-30 - self.orientation))),
                ],
                [
                    0.5,
                    0.5 * math.cos(math.radians(2 * (-150 - self.orientation))),
                    0.5 * math.sin(math.radians(2 * (-150 - self.orientation))),
                ],
                [
                    0.5,
                    0.5 * math.cos(math.radians(2 * (-120 - self.orientation))),
                    0.5 * math.sin(math.radians(2 * (-120 - self.orientation))),
                ],
            ]
        )

        # remove row from orientation matrix corresponding to gage weight of 0
        orientations = np.matmul(gage_weights, orientations)
        orientations = orientations[~np.all(orientations == 0, axis=1)]
        # scale gage_weights down to 3x3 identity if only using 3 gages
        gage_weights = gage_weights[:, ~np.all(gage_weights == 0, axis=0)]
        gage_weights = gage_weights[~np.all(gage_weights == 0, axis=1)]

        # calculate the strain matrix
        strain_matrix = np.matmul(
            np.matmul(np.linalg.inv(scale_factors), np.linalg.pinv(orientations)),
            np.linalg.inv(gage_weights),
        )

        # insert a column of zeros back in if one gage was dropped, leaving a 4x3 matrix
        for i in [0, 1, 2, 3]:
            if gauge_weights[i] == 0:
                strain_matrix = np.insert(strain_matrix, i, 0, axis=1)
        # print("strain_matrix: \n",strain_matrix)
        return strain_matrix

    def get_lab_strain_matrix(self):
        """parse lab strain matrix from station metadata 

        :return: lab calibration matrix
        :rtype: np.array
        """
        try:
            url = f"http://bsm.unavco.org/bsm/level2/{self.station}/{self.station}.README.txt"

            with request.urlopen(url) as response:
                lines = response.readlines()
            for i, line in enumerate(lines):
                line = line.decode("utf-8").rstrip()
                if line.startswith("  Manufacturer's Isotropic Strain Matrix"):
                    lab = np.array(
                        [
                            lines[i + 1].decode("utf-8").rstrip().split()[1:],
                            lines[i + 2].decode("utf-8").rstrip().split()[1:],
                            lines[i + 3].decode("utf-8").rstrip().split()[1:],
                        ]
                    )
            return lab.astype(float)
        except Exception as e:
            logger.error("Could not load lab strain matrix")
            return None

    def get_er2010_strain_matrix(self):
        """parse ER2010 strain matrix from station metadata if available.  From
        Roeloffs, E. (2010), Tidal calibration of Plate Boundary Observatory borehole strainmeters: 
        Roles of vertical and shear coupling , J. Geophys. Res., 115, B06405, doi:10.1029/2009JB006407
        
        :return: ER2010 calibration matrix
        :rtype: np.array
        """
        url = f"http://bsm.unavco.org/bsm/level2/{self.station}/{self.station}.README.txt"
        try:
            with request.urlopen(url) as response:
                lines = response.readlines()
            for i, line in enumerate(lines):
                line = line.decode("utf-8").rstrip()
                if line.startswith("  Roeloffs 2010 Tidal Calibration"):
                    er2010 = np.array(
                        [
                            lines[i + 3].decode("utf-8").rstrip().split()[1:],
                            lines[i + 4].decode("utf-8").rstrip().split()[1:],
                            lines[i + 5].decode("utf-8").rstrip().split()[1:],
                        ]
                    )
                    return er2010.astype(float)
            return None
        except Exception as e:
            logger.error("Could not load ER2010 strain matrix")
            return None

    def get_ch_prelim_strain_matrix(self):
        """parse CH preliminary strain matrix from station metadata if available.  
        Calibrations of TABOO-STAR stations performed by Cassie Hanagan in 2022.
        
        :return: CH_PRELIM calibration matrix
        :rtype: np.array
        """
        url = f"http://bsm.unavco.org/bsm/level2/{self.station}/{self.station}.README.txt"
        try:
            with request.urlopen(url) as response:
                lines = response.readlines()
            for i, line in enumerate(lines):
                line = line.decode("utf-8").rstrip()
                if line.startswith("  CH Preliminary Tidal Calibration"):
                    ch_prelim = np.array(
                        [
                            lines[i + 3].decode("utf-8").rstrip().split()[1:],
                            lines[i + 4].decode("utf-8").rstrip().split()[1:],
                            lines[i + 5].decode("utf-8").rstrip().split()[1:],
                        ]
                    )
                    return ch_prelim.astype(float)
            return None
        except Exception as e:
            logger.exception("Could not load ch_prelim strain matrix")
            return None

    def get_reference_strains(self):
        """Parse reference strains from metadata 

        :return: reference strains
        :rtype: dict
        """
        reference_strains = {}
        reference_strains["linear_date"] = self.meta_df.loc[self.station]["L_DATE"]
        reference_strains["CH0"] = int(self.meta_df.loc[self.station]["L0(cnts)"])
        reference_strains["CH1"] = int(self.meta_df.loc[self.station]["L1(cnts)"])
        reference_strains["CH2"] = int(self.meta_df.loc[self.station]["L2(cnts)"])
        reference_strains["CH3"] = int(self.meta_df.loc[self.station]["L3(cnts)"])
        return reference_strains

    # def get_linearization_params_xml(self):
    #     # get reference strains from xml
    #     linear_dict = self.xml["strain_xml"]["inst_info"]["processing"][
    #         "bsm_processing_history"
    #     ][-1]["bsm_processing"]["linearization"]
    #     self.reference_strains = {}
    #     for key in linear_dict:
    #         # print(key, linear_dict[key])
    #         if key == "linear_date":
    #             self.reference_strains["linear_date"] = linear_dict[key]
    #         if key == "g0_value":
    #             self.reference_strains["CH0"] = float(linear_dict[key])
    #         if key == "g1_value":
    #             self.reference_strains["CH1"] = float(linear_dict[key])
    #         if key == "g2_value":
    #             self.reference_strains["CH2"] = float(linear_dict[key])
    #         if key == "g3_value":
    #             self.reference_strains["CH3"] = float(linear_dict[key])
    #
    # def get_gauge_weightings_xml(self):
    #     weight_dict = self.xml["strain_xml"]["inst_info"]["processing"][
    #         "bsm_processing_history"
    #     ][-1]["bsm_processing"]["gauge_weightings"]
    #     self.gauge_weightings = [
    #         int(weight_dict["gw0"]),
    #         int(weight_dict["gw1"]),
    #         int(weight_dict["gw2"]),
    #         int(weight_dict["gw3"]),
    #     ]
    #
    # def get_detrend_xml(self):
    #     detrend = {}
    #     for channel in [
    #         "detrend_start_date",
    #         "detrend_g0",
    #         "detrend_g1",
    #         "detrend_g2",
    #         "detrend_g3",
    #     ]:
    #         if channel == "detrend_start_date":
    #             detrend[channel] = self.xml["strain_xml"]["inst_info"]["processing"][
    #                 "bsm_processing_history"
    #             ][-1]["bsm_processing"]["timeseries_start_date"]
    #         else:
    #             detrend_dict = {}
    #             detrend_params = self.xml["strain_xml"]["inst_info"]["processing"][
    #                 "bsm_processing_history"
    #             ][-1]["bsm_processing"][channel]
    #             for key in detrend_params:
    #                 if key[0] != "@":
    #                     detrend_dict[key] = float(detrend_params[key])
    #             detrend[channel] = detrend_dict
    #     return detrend

    def get_atmp_response(self):
        """Parse atmospheric pressure response from station processing metadata

        :return: atmospheric response coefficients
        :rtype: dict
        """
        url = f"http://bsm.unavco.org/bsm/level2/{self.station}/{self.station}.README.txt"
        try:
            with request.urlopen(url) as response:
                lines = response.readlines()
            baro = False
            atmp_coefficients = {}
            for line in lines:
                line = line.decode("utf-8").rstrip()
                if line.startswith("Barometric Response"):
                    baro = True
                if baro:
                    if line.startswith("CH"):
                        line = line.split()
                        atmp_coefficients[line[0]] = (
                            float(line[1]) * 1e-3
                        )  # microstrain
                        if line[0] == "CH3":
                            baro = False
            return atmp_coefficients
        except Exception as e:
            logger.error("Could not load atmp response")
            return None

    def get_tidal_params(self):
        """Parse tidal params from station processing metadata

        :return: tidal coefficients
        :rtype: dict
        """
        url = f"http://bsm.unavco.org/bsm/level2/{self.station}/{self.station}.README.txt"
        try:
            with request.urlopen(url) as response:
                lines = response.readlines()
            tide = False
            tide_coefficients = {}
            for line in lines:
                line = line.decode("utf-8").rstrip()

                if line.startswith("Tidal Constituents"):
                    tide = True

                if tide:
                    if line.startswith("CH"):
                        line = line.split()
                        if line[1] == "M2":
                            doodson = "2 0 0 0 0 0"
                        elif line[1] == "O1":
                            doodson = "1-1 0 0 0 0"
                        elif line[1] == "P1":
                            doodson = "1 1-2 0 0 0"
                        elif line[1] == "K1":
                            doodson = "1 1 0 0 0 0"
                        elif line[1] == "N2":
                            doodson = "2-1 0 1 0 0"
                        elif line[1] == "S2":
                            doodson = "2 2-2 0 0 0"
                        tide_coefficients[(line[0], line[1], "phz")] = line[2]
                        tide_coefficients[(line[0], line[1], "amp")] = line[3]
                        tide_coefficients[(line[0], line[1], "doodson")] = doodson
            return tide_coefficients
        except Exception as e:
            logger.error("Could not load tidal parameters")
            return None

    def load_site_terms(self):
        """Read site terms from file"""
        terms_file = files("earthscopestraintools").joinpath("event_site_terms.txt")
        site_terms = pd.read_csv(terms_file, sep=" ", header=0).fillna(0)
        return site_terms

    def get_event_terms(self):
        """
        Loads site and longitude terms for a particular strainmeter using Barbour et al 2021.  
        These terms are used for event magnitude estimation.
        """
        site_terms = self.load_site_terms()
        if self.station in site_terms.Station.values:
            self.site_term = site_terms[site_terms.Station == self.station].delta_s.item()
        else:
            self.site_term = 0
        if self.longitude < -124:
            self.longitude_term = -0.41
        else:
            self.longitude_term = 0

    def show(self):
        """print metadata
        """
        # pp = pprint.PrettyPrinter()
        logger.info(f"network: {self.network}")
        logger.info(f"station: {self.station}")
        logger.info(f"latitude: {self.latitude}")
        logger.info(f"longitude: {self.longitude}")
        logger.info(f"gap: {self.gap}")
        logger.info(f"orientation (CH0EofN): {self.orientation}")
        logger.info(f"reference_strains:\n {self.reference_strains}")
        logger.info("strain_matrices:")
        if len(self.strain_matrices):
            for key in self.strain_matrices:
                logger.info(f"{key}:\n {self.strain_matrices[key]}")
                # logger.info(self.strain_matrices[key])
                # pp.pprint(self.strain_matrices[key])
        logger.info(f"atmp_response:\n {self.atmp_response}")
        # logger.info(self.atmp_response)
        logger.info(f"tidal_params:\n {self.tidal_params}")
        # logger.info(self.tidal_params)

        # print("reference strains:")
        # pp.pprint(self.linearization)
        # if len(self.strain_matrices):
        #     for key in self.strain_matrices:
        #         print(f"{key}:")
        #         print(self.strain_matrices[key])
        #         #pp.pprint(self.strain_matrices[key])
        # print("atmp coefficients:")
        # pp.pprint(self.atmp_response)
        # print("tidal params: ")
        # pp.pprint(self.tidal_params)
        # print("detrend params:")
        # pp.pprint(self.detrend)


def fdsn2bottlename(channel):
    """convert FDSN channel into bottlename

    :param channel: FDSN channel code
    :type channel: str
    :return: bottlename 
    :rtype: str
    """

    codes = {
        "RS1": "CH0",
        "LS1": "CH0",
        "BS1": "CH0",
        "RS2": "CH1",
        "LS2": "CH1",
        "BS2": "CH1",
        "RS3": "CH2",
        "LS3": "CH2",
        "BS3": "CH2",
        "RS4": "CH3",
        "LS4": "CH3",
        "BS4": "CH3",
        "RDO": "atmp",
        "LDO": "atmp",
    }

    return codes[channel]


def get_metadata_df():
    """Function loads strainmeter metadata into pandas dataframe

    :return: bsm metadata 
    :rtype: pandas.DataFrame
    """
        
    # """
    # 

    # Parameters
    # ----------
    # :return: metadata: pandas DataFrame

    # """
    url = (
        "https://www.unavco.org/data/strain-seismic/bsm-data/lib/docs/bsm_metadata.txt"
    )
    metadata = pd.read_csv(url, index_col="BNUM", sep="\s+", on_bad_lines="skip")
    return metadata


def get_fdsn_network(station):
    """Get FDSN network name from station fourcharid using es-datasources-api

    :param station: station four character id
    :type station: string
    :return: 2 char FDSN network code for station
    :rtype: string
    """
    try:
        DATASOURCES_API_URL = "https://datasources-api.prod.earthscope.org"
        return lookup_fdsn_network_from_bnum(DATASOURCES_API_URL,station)
    except requests.exceptions.ConnectionError:
        logger.error(
            "Error: Unable to connect to datasources-api to lookup FDSN network code"
        )
        return None

src/earthscopestraintools/gtsm_metadata.py
EarthScope/earthscopestraintools
import numpy as np
import pandas as pd
import datetime

from obspy import read, Stream
from obspy.clients.fdsn import Client
from earthscopestraintools.timeseries import Timeseries

inv_client = Client("IRIS")
from requests.exceptions import HTTPError


import logging


logger = logging.getLogger(__name__)


def load_mseed_to_df(
    net: str,
    sta: str,
    loc: str,
    cha: str,
    start: str,
    end: str,
    save_as: str = None,
    print_url: bool = False,
    print_traces: bool = True
):
    """
    Load miniseed data from fdsnws-dataselect a pandas dataframe, with time as the index and
    each channel as a column.  Uses obspy as an intermediate step.

    :param net: FDSN 2 character network code, accepts wildcards
    :type net: str
    :param sta: FDSN 4 character station code, accepts wildcards
    :type sta: str
    :param loc: FDSN 2 character location code, accepts wildcards
    :type loc: str
    :param cha: FDSN 3 character channel code, accepts wildcards
    :type cha: str
    :param start: start time of query as ISO formatted string
    :type start: str
    :param end: end time of query as ISO formatted string
    :type end: str
    :param save_as: filename to save miniseed data, defaults to None
    :type save_as: str, optional
    :param print_url: show dataselect url, defaults to False
    :type print_url: bool, optional
    :param print_traces: show traces loaded into obspy, defaults to True
    :type print_traces: bool, optional
    :return: data loaded into pandas dataframe
    :rtype: pd.DataFrame
    """
    st = download_mseed(
        net=net, sta=sta, loc=loc, cha=cha, start=start, end=end, print_url=print_url
    )
    if save_as:
        st.write(filename=save_as)
    df = mseed2pandas(st, print_traces=print_traces)
    return df


def load_mseed_file_to_df(filename: str):
    """read miniseed file from local disk into pandas dataframe, with time as the index and
    each channel as a column.  Uses obspy as an intermediate step.
 

    :param filename: miniseed filename
    :type filename: str
    :return: data loaded into pandas dataframe
    :rtype: pd.DataFrame
    """
    st = read(filename)
    df = mseed2pandas(st)
    return df


def save_mseed_file(
    net: str, sta: str, loc: str, cha: str, start: str, end: str, filename: str
):
    """Load miniseed data from fdsnws-dataselect into obspy, and save as miniseed.

    :param net: FDSN 2 character network code, accepts wildcards
    :type net: str
    :param sta: FDSN 4 character station code, accepts wildcards
    :type sta: str
    :param loc: FDSN 2 character location code, accepts wildcards
    :type loc: str
    :param cha: FDSN 3 character channel code, accepts wildcards
    :type cha: str
    :param start: start time of query as ISO formatted string
    :type start: str
    :param end: end time of query as ISO formatted string
    :type end: str
    :param filename: filename to save miniseed data
    :type filename: str
    """
    st = download_mseed(net=net, sta=sta, loc=loc, cha=cha, start=start, end=end)
    st.write(filename)


def download_mseed(
    net: str,
    sta: str,
    loc: str,
    cha: str,
    start: str,
    end: str,
    print_url: bool = False,
):
    """Load miniseed data from fdsnws-dataselect into obspy Stream object.

    :param net: FDSN 2 character network code, accepts wildcards
    :type net: str
    :param sta: FDSN 4 character station code, accepts wildcards
    :type sta: str
    :param loc: FDSN 2 character location code, accepts wildcards
    :type loc: str
    :param cha: FDSN 3 character channel code, accepts wildcards
    :type cha: str
    :param start: start time of query as ISO formatted string
    :type start: str
    :param end: end time of query as ISO formatted string
    :type end: str
    :param print_url: show dataselect url, defaults to False
    :type print_url: bool, optional
    :return: data as Stream object
    :rtype: obspy.Core.Stream
    """
    logger.info(
        f"{net} {sta} Loading {loc} {cha} from {start} to {end} from Earthscope DMC miniseed"
    )
    url = (
        f"https://service.iris.edu/fdsnws/dataselect/1/query?net={net}"
        f"&sta={sta}&loc={loc}&cha={cha}&starttime={start}&endtime={end}"
        f"&format=miniseed&nodata=404"
    )
    if print_url:
        logger.info(f"{net} {sta} Reading from {url}")
    try:
        st = read(url)
    except HTTPError as e:
        logger.error(f"{net} {sta} No data found for specified query {url}")
        st = Stream()
    return st


def mseed2pandas(st: Stream, print_traces=True):
    """Restructure data from obspy.Core.Stream to pandas.DataFrame with time as the index and
    each channel as a column.

    :param st: Stream object containing one or more Traces
    :type st: obspy.Core.Stream
    :param print_traces: show traces loaded into obspy, defaults to True
    :type print_traces: bool, optional
    :return: data loaded into pandas dataframe
    :rtype: pandas.DataFrame
    """
    df = pd.DataFrame()
    dfs = {}
    for i, tr in enumerate(st):
        channel = fdsn2bottlename(tr.stats.channel)
        start = tr.stats["starttime"]
        stop = tr.stats["endtime"]
        step = datetime.timedelta(seconds=tr.stats["delta"])
        if print_traces:
            logger.info(
                f"    Trace {i + 1}. {start}:{stop} mapping {tr.stats.channel} to {channel}"
            )
        time_buffer = np.arange(start, stop + step, step, dtype=np.datetime64)
        df2 = pd.DataFrame(index=time_buffer)
        df2[channel] = tr.data
        if channel not in dfs.keys():
            dfs[channel] = df2
        else:
            dfs[channel] = pd.concat([dfs[channel], df2])
    for i, key in enumerate(dfs.keys()):
        df[key] = dfs[key]
    df.index.name = "time"
    return df


def fdsn2bottlename(channel):
    """convert FDSN channel code into bottlename

    :param channel: FDSN channel
    :type channel: str
    :return: bottlename
    :rtype: str
    """

    codes = {
        "RS1": "CH0",
        "LS1": "CH0",
        "BS1": "CH0",
        "RS2": "CH1",
        "LS2": "CH1",
        "BS2": "CH1",
        "RS3": "CH2",
        "LS3": "CH2",
        "BS3": "CH2",
        "RS4": "CH3",
        "LS4": "CH3",
        "BS4": "CH3",
        "RDO": "atmp",
        "LDO": "atmp",
        "RRO": "rain",
        "RK1": "LoggerDegC",
        "RKD": "DownholeDegC",
        "RE1": "BatteryVolts",
        "REO": "SolarAmps",
        "RE2": "SystemAmps",
        "RK2": "PowerBoxDegC",
        "LDD": "PorePressure",
        "LKD": "PoreDegC",
        "VEP": "Q330Volts"
    }

    return codes[channel]


def ts_from_mseed(
    network: str,
    station: str,
    location: str,
    channel: str,
    start: str,
    end: str,
    name: str = None,
    period: float = None,
    series: str = "raw",
    units: str = "",
    to_nan: bool = True,
    scale_factor: float = None,
):
    """Load miniseed data from fdsnws-dataselect into Timeseries object.  

    :param network: FDSN 2 character network code, accepts wildcards
    :type network: str
    :param station: FDSN 4 character station code, accepts wildcards
    :type station: str
    :param location: FDSN 2 character location code, accepts wildcards
    :type location: str
    :param channel: FDSN 3 character channel code, accepts wildcards
    :type channel: str
    :param start: start time of query as ISO formatted string
    :type start: str
    :param end: end time of query as ISO formatted string
    :type end: str
    :param name: name of timeseries, including station name.  useful for showing stats. , defaults to None
    :type name: str, optional
    :param period: sample period of data, defaults to None
    :type period: float, optional
    :param series: description of timeseries, ie 'raw', 'microstrain', 'atmp_c', 'tide_c', 'offset_c', 'trend_c', defaults to "raw"
    :type series: str, optional
    :param units: units of timeseries, defaults to ""
    :type units: str, optional
    :param to_nan: option to convert 999999 gap fill values to numpy.nan, defaults to True
    :type to_nan: bool, optional
    :param scale_factor: scale factor to apply to miniseed data, defaults to None
    :type scale_factor: float, optional
    :return: Timeseries object containing data loaded from miniseed
    :rtype: earthscopestraintools.timeseries.Timeseries
    """
    df = load_mseed_to_df(
        net=network,
        sta=station,
        loc=location,
        cha=channel,
        start=start,
        end=end,
    )
    level = "0"
    if channel.startswith("RS"):
        period = 600
        series = series
        units = "counts"
    elif channel.startswith("LS"):
        period = 1
        series = series
        units = "counts"
    elif channel.startswith("BS"):
        period = 0.05
        series = series
        units = "counts"
    else:
        period = period
        series = series
    if name is None:
        name = f"{network}.{station}.{location}.{channel}"
    ts = Timeseries(
        data=df, series=series, units=units, level=level, period=period, name=name
    )
    if to_nan:
        logger.info("Converting missing data from 999999 to nan")
        ts = ts.remove_999999s()
    if scale_factor:
        ts.data = ts.data * scale_factor
    return ts

src/earthscopestraintools/mseed_tools.py
EarthScope/earthscopestraintools
import numpy as np
import pandas as pd
from scipy import signal, stats
import subprocess

#from earthscopestraintools.gtsm_metadata import GtsmMetadata
import logging

logger = logging.getLogger(__name__)


def linearize(df: pd.DataFrame, reference_strains: dict, gap: float):
    """Linearize raw gauge strain.

    :param df: Dataframe with four columns corresponding to raw gauge data in units of counts.
    :type df: pd.DataFrame
    :param reference_strains: Dictionary with four entries noting the reference count value on each gauge to zero the timeseries against in the conversion to gauge strain (in microstrain). 
    :type reference_strains: dict
    :param gap: Instrument measurement gap. For most NOTA strainmeters, the reference gap is 0.0001 m.
    :type gap: float
    :return: DataFrame with with four columns of microstrain.
    :rtype: pd.DataFrame
    """

    logger.info(f"Converting raw counts to microstrain")
    # build a series of reference strains.  if using metadata from XML the /1e8 is already included.
    reference_strain_series = pd.Series(dtype="float64")
    for ch in df.columns:
        reference_strain_series[ch] = reference_strains[ch]
    linearized_data = (
        (
            ((df / 100000000) / (1 - (df / 100000000)))
            - (
                (reference_strain_series / 100000000)
                / (1 - (reference_strain_series / 100000000))
            )
        )
        * (gap / 0.087)
        * 1000000
    )
    return linearized_data


def apply_calibration_matrix(
    df: pd.DataFrame,
    calibration_matrix: np.array,
    calibration_matrix_name: str = None,
    use_channels: list = [1, 1, 1, 1],
):
    """Applies a calibration matrix to convert 4 gauges into areal, differential, and shear strains

    :param df: four channels of strain data or correction in microstrain
    :type df: pd.DataFrame
    :param calibration_matrix: calibration matrix 
    :type calibration_matrix: np.array
    :param calibration_matrix_name: name of calibration matrix used, defaults to None
    :type calibration_matrix_name: str, optional
    :param use_channels: not yet implemented, set to 0 to ignore a bad channel, defaults to [1, 1, 1, 1]
    :type use_channels: list, optional
    :return: areal, differential, and shear strains based on the given calibration
    :rtype: pd.DataFrame
    """
    # todo: implement use channels
    logger.info(f"Applying {calibration_matrix_name} matrix: {calibration_matrix}")
    areal = "Eee+Enn"
    differential = "Eee-Enn"
    shear = "2Ene"
    if calibration_matrix_name and calibration_matrix_name != "lab":
        areal += "." + calibration_matrix_name
        differential += "." + calibration_matrix_name
        shear += "." + calibration_matrix_name
    regional_strain_df = np.matmul(
        calibration_matrix, df[df.columns].transpose()
    ).transpose()
    regional_strain_df = regional_strain_df.rename(
        columns={0: areal, 1: differential, 2: shear}
    )
    return regional_strain_df


def butterworth_filter(
    df: pd.DataFrame,
    period: float,
    filter_type: str,
    filter_order: int,
    filter_cutoff_s: float,
):
    """Apply a butterworth filter to a DataFrame using scipy.signal.butter()

    :param df: data to filter
    :type df: pd.DataFrame
    :param period: sample period of data
    :type period: float
    :param filter_type: {'lowpass', 'highpass', 'bandpass', 'bandstop'}
    :type filter_type: str
    :param filter_order: the order of the filter
    :type filter_order: int
    :param filter_cutoff_s: the filter cutoff in seconds
    :type filter_cutoff_s: float
    :return: butterworth filtered data
    :rtype: pandas.DataFrame
    """
    logger.info(f"Applying Butterworth Filter")
    fc = 1 / filter_cutoff_s
    fs = 1 / period
    [bn, an] = signal.butter(filter_order, fc / (1 / 2 * fs), btype=filter_type)
    df2 = pd.DataFrame(index=df.index)
    for ch in df.columns:
        df2[ch] = signal.filtfilt(bn, an, df[ch])
    return df2


def interpolate(
    df: pd.DataFrame,
    replace: int = 999999,
    method: str = "linear",
    limit: int = 3600,
    limit_direction="both",
):
    """Interpolate across gaps in data using pd.DataFrame.interpolate()

    :param df: data to be interpolated
    :type df: pd.DataFrame
    :param replace: gap fill value to interpolate across, defaults to 999999
    :type replace: int, optional
    :param method: interpolation method, defaults to "linear"
    :type method: str, optional
    :param limit: max number of samples to interpolate, defaults to 3600
    :type limit: int, optional
    :param limit_direction: ['forward', 'backward', 'both'], defaults to "both"
    :type limit_direction: str, optional
    :return: interpolated data
    :rtype: pd.DataFrame
    """
    logger.info(f"Interpolating data using method={method} and limit={limit}")
    df2 = df.replace(replace, np.nan).interpolate(
        method=method, limit_direction=limit_direction, limit=limit
    )
    return df2


def decimate_to_hourly(df: pd.DataFrame):
    """decimates a timeseries to hourly by selecting the first minute of each hour

    :param df: time series data to decimate
    :type df: pd.DataFrame
    :return: decimated data
    :rtype: pd.DataFrame
    """
    logger.info(f"Decimating to hourly")
    return df[df.index.minute == 0]


def decimate_1s_to_300s(df: pd.DataFrame, method: str = "linear", limit: int = 3600):
    """decimate 1hz data to 5 min data using \n
    Agnew, Duncan Carr, and K. Hodgkinson (2007), Designing compact causal digital filters for 
    low-frequency strainmeter data , Bulletin Of The Seismological Society Of America, 97, No. 1B, 91-99

    :param df: 1 hz data
    :type df: pd.DataFrame
    :param method: method to interpolate across gaps, defaults to "linear"
    :type method: str, optional
    :param limit: largest gap to interpolate, defaults to 3600 samples
    :type limit: int, optional
    :return: 300s (5 min) data
    :rtype: pd.DataFrame
    """
    logger.info(f"Decimating to 300s")
    df2 = interpolate(df, method="linear", limit=3600)

    # zero the data to the first value prior to filtering
    initial_values = df2.iloc[0]
    df3 = df2 - initial_values

    # perform a 5 stage filter-and-decimate
    # 1s -> 2s -> 4s -> 12s -> 60s -> 300s
    # FIR filter coefficients from Agnew and Hodgkinson 2007

    wtsdby2d = [
        0.0983262,
        0.2977611,
        0.4086973,
        0.3138961,
        0.0494246,
        -0.1507778,
        -0.1123764,
        0.0376576,
        0.0996838,
        0.0154992,
        -0.0666489,
        -0.0346632,
        0.0322767,
        0.0399294,
        -0.0097461,
        -0.0341585,
        -0.0039241,
        0.0246776,
        0.0099725,
        -0.0157879,
        -0.0099098,
        0.0078510,
        0.0081126,
        -0.0026986,
        -0.0061424,
        0.0007108,
        0.0039659,
        -0.0006209,
        -0.0017117,
        0.0007240,
    ]

    wtsdby3c = [
        0.0373766,
        0.1165151,
        0.2385729,
        0.3083302,
        0.2887327,
        0.1597948,
        0.0058244,
        -0.0973639,
        -0.1051034,
        -0.0358455,
        0.0359044,
        0.0632477,
        0.0302351,
        -0.0168856,
        -0.0356758,
        -0.0190635,
        0.0126188,
        0.0159705,
        0.0082144,
        -0.0087978,
        -0.0037289,
        -0.0017068,
        0.0028335,
    ]

    wtsdby5b = [
        0.0218528,
        0.0458359,
        0.0908603,
        0.1359777,
        0.1830881,
        0.1993418,
        0.1957624,
        0.1561194,
        0.0994146,
        0.0346412,
        -0.0236544,
        -0.0580081,
        -0.0703257,
        -0.0555546,
        -0.0287709,
        0.0032613,
        0.0267938,
        0.0358952,
        0.0311186,
        0.0134283,
        -0.0028524,
        -0.0170042,
        -0.0176765,
        -0.0123123,
        -0.0036798,
        0.0057730,
        0.0059817,
        0.0083501,
        0.0000581,
        0.0005724,
        -0.0033127,
        0.0004411,
        -0.0030766,
        0.0016604,
    ]

    channels = df.columns

    stage1 = pd.DataFrame(index=df3.index)
    for ch in channels:
        data = df3[ch].values
        stage1[ch] = signal.lfilter(wtsdby2d, 1.0, data)
    stage1d = stage1.iloc[::2]

    stage2 = pd.DataFrame(index=stage1d.index)
    for ch in channels:
        data = stage1d[ch].values
        stage2[ch] = signal.lfilter(wtsdby2d, 1.0, data)
    stage2d = stage2.iloc[::2]

    stage3 = pd.DataFrame(index=stage2d.index)
    for ch in channels:
        data = stage2d[ch].values
        stage3[ch] = signal.lfilter(wtsdby3c, 1.0, data)
    stage3d = stage3.iloc[::3]

    stage4 = pd.DataFrame(index=stage3d.index)
    for ch in channels:
        data = stage3d[ch].values
        stage4[ch] = signal.lfilter(wtsdby5b, 1.0, data)
    stage4d = stage4.iloc[::5]

    stage5 = pd.DataFrame(index=stage4d.index)
    for ch in channels:
        data = stage4d[ch].values
        stage5[ch] = signal.lfilter(wtsdby5b, 1.0, data)
    stage5d = stage5.iloc[::5]

    # add back in the initial values
    decimated_data = (stage5d + initial_values).astype(int)
    return decimated_data


def calculate_offsets(df, limit_multiplier: int = 10, cutoff_percentile: float = 0.75):
    """Calculate offsets using first differencing method (add more details).  

    :param df: uncorrected data, as dataframe with datetime index and 1 channel per column
    :type df: pandas.DataFrame
    :param limit_multiplier: _description_, defaults to 10
    :type limit_multiplier: int, optional
    :param cutoff_percentile: _description_, defaults to 0.75
    :type cutoff_percentile: float, optional
    :return: _description_
    :rtype: _type_
    """
    logger.info(
        f"Calculating offsets using cutoff percentile of {cutoff_percentile} and limit multiplier of {limit_multiplier}."
    )
    first_diffs = df.diff()
    #  drop a percentage of 1st differences to estimate an offset cutoff
    drop = round(len(first_diffs) * (1 - cutoff_percentile))
    offset_limit = []
    df_offsets = pd.DataFrame(index=first_diffs.index)
    for ch in df.columns:
        # Offset limit is a multiplier x the average absolute value of first differences
        # within 2 st_dev of the first differences
        offset_limit.append(
            np.mean(abs(first_diffs[ch].sort_values().iloc[0:-drop])) * limit_multiplier
        )

        # CH edit. Calculate offsets from the detrended series
        # Justification: if the offset is calculated from the original series,
        # there may be an overcorrection that is noticeable, especially with large trends
        df_offsets[ch] = first_diffs[first_diffs[ch].abs() > offset_limit[-1]][ch]

    # make a dataframe of running total of offsets
    df_cumsum = df_offsets.fillna(0).cumsum()
    # us_limits =
    logger.info(f"Using offset limits of {[round(x,6) for x in offset_limit]}")
    return df_cumsum


def calculate_pressure_correction(df: pd.DataFrame, response_coefficients: dict):
    """Generate a pressure correction timeseries from pressure data and response coefficients

    :param df: atmospheric pressure data 
    :type df: pd.DataFrame
    :param response_coefficients: response coefficients for each channel loaded from metadata
    :type response_coefficients: dict
    :return: pressure corrections for each channel
    :rtype: pd.DataFrame
    """
    logger.info(f"Calculating pressure correction")
    data_df = pd.DataFrame(index=df.index)
    for key in response_coefficients:
        data_df[key] = df * float(response_coefficients[key])
    return data_df


def calculate_linear_trend_correction(df, trend_start=None, trend_end=None):
    """Generate a linear trend correction

    :param df: uncorrected data, as dataframe with datetime index and 1 channel per column
    :type df: pd.DataFrame
    :param trend_start: start of window to calculate trend, defaults to first_valid_index()
    :type trend_start: datetime.datetime, optional
    :param trend_end: end of window to calculate trend, defaults to last_valid_index()
    :type trend_end: datetime.datetime, optional
    :return: trend correction timeseries for each column/channel in input data
    :rtype: pd.DataFrame
    """
    logger.info(f"Calculating linear trend correction")
    if not trend_start:
        trend_start = df.first_valid_index()
    if not trend_end:
        trend_end = df.last_valid_index()
    logger.info(f"    Trend Start: {trend_start}")
    logger.info(f"    Trend Start: {trend_end}")
    df_trend_c = pd.DataFrame(data=df.index)
    windowed_df = df.copy()[trend_start:trend_end].interpolate().reset_index()
    for ch in df.columns:
        slope, intercept, r_value, p_value, std_err = stats.linregress(
            windowed_df.index, windowed_df[ch]
        )
        df_trend_c[ch] = df_trend_c.index * slope
    # print("df_trend_c", df_trend_c)
    return df_trend_c[df.columns].set_index(df_trend_c["time"])


def calculate_tide_correction(df, period, tidal_parameters, longitude):
    """Generate tidal correction timeseries using SPOTL hartid

    :param df: uncorrected data, as dataframe with datetime index and 1 channel per column
    :type df: pd.DataFrame
    :param period: sample period of data, must be >= 1
    :type period: int
    :param tidal_parameters: tidal parameters loaded from station metadata
    :type tidal_parameters: dict
    :param longitude: station longitude
    :type longitude: float
    :return: tidal correction timeseries for each column/channel in input data
    :rtype: pd.DataFrame
    """
    logger.info(f"Calculating tide correction")
    # check if hartid in path, or else use spotl container
    result = subprocess.run("which hartid", shell=True)
    if result.returncode == 0:
        hartid = "hartid"
    else:
        hartid = "docker run -i --rm ghcr.io/earthscope/spotl hartid"
    start = df.first_valid_index()
    datestring = (
        str(start.year).zfill(4)
        + " "
        + str(start.month).zfill(2)
        + " "
        + str(start.day).zfill(2)
        + " "
        + str(start.hour).zfill(2)
        + " "
        + str(start.minute).zfill(2)
        + " "
        + str(start.second).zfill(2)
    )

    nterms = str(len(df))
    samp = int(period)
    cmd = f"{hartid} {datestring} {nterms} {samp}"
    # cmd = hartid + " " + datestring + " " + nterms + " " + samp

    channels = df.columns
    tides = set()
    for key in tidal_parameters:
        tides.add(key[1])
    cmds = {}
    # for i, ch in enumerate(gauges):
    for ch in channels:
        inputfile = f"printf 'l\n{longitude}\n"
        for tide in tides:
            inputfile += f" {tidal_parameters[(ch, tide, 'doodson')]} {tidal_parameters[(ch, tide, 'amp')].ljust(7, '0')} {tidal_parameters[(ch, tide, 'phz')].ljust(8, '0')}\n"
        inputfile += "-1'"
        cmds[ch] = inputfile + " | " + cmd

    df2 = pd.DataFrame(index=df.index)
    for ch in channels:
        output = subprocess.check_output(cmds[ch], shell=True).decode("utf-8")
        df2[ch] = np.fromstring(output, dtype=float, sep="\n")
    df2 = df2 * 1e-3
    return df2

src/earthscopestraintools/processing.py
EarthScope/earthscopestraintools
# write ascii files from tiledb arrays
# slice by start and end time

import tarfile
import tiledb
import numpy as np
import pandas as pd
import os
import json
import datetime
from earthscopestraintools.tiledbtools import ProcessedStrainReader
from earthscopestraintools.datasources_api_interact import get_station_edid

# from earthscopestraintools.straintiledbarray import StrainTiledbArray
# from edid import get_station_edid

import logging


logging.basicConfig(
    level=logging.INFO,
    format="%(asctime)s %(levelname)s: %(message)s",
    datefmt="%Y-%m-%d %H:%M:%S",
)
logger = logging.getLogger(__name__)

workdir = "arrays"
outputdir = "ascii_output"
DATASOURCES_API_URL = "https://datasources-api.prod.earthscope.org"

def to_date(datetime64):
    ts = pd.to_datetime(str(datetime64))
    d = ts.strftime("%Y%m%d")
    return d


def read_date_range(network, station, start_str, end_str):
    edid = get_station_edid(api_url=DATASOURCES_API_URL, namespace="BNUM", station=station)
    uri = f"{workdir}/{edid}_level2.tdb"
    logger.info(f"Array uri: {uri}")
    reader = ProcessedStrainReader(uri)
    data_types = reader.array.get_data_types()
    timeseries = reader.array.get_timeseries()
    attrs = ["data", "quality", "level", "version"]
    df = reader.to_df(
        data_types=data_types,
        timeseries=timeseries,
        attrs=attrs,
        start_str=start_str,
        end_str=end_str,
        reindex=False,
    )
    return df


def write_ascii(df, station, date_range_string, write_it=True, print_it=False):
    filenames = []
    for data_type in ["CH0", "CH1", "CH2", "CH3", "Eee+Enn", "Eee-Enn", "2Ene"]:
        logger.info(f"generating ascii file for {data_type}")
        for timeseries in ["microstrain", "offset_c", "tide_c", "trend_c", "atmp_c"]:
            tmp_df = (
                df.sort_index()
                .loc[data_type, timeseries]
                .rename(columns={"data": timeseries})
            )
            if timeseries == "microstrain":
                tmp_df = tmp_df.rename(
                    columns={
                        "microstrain": f"{data_type}(mstrain)",
                        "quality": "strain_quality",
                    }
                )
                new_df = tmp_df
                new_df["strain"] = data_type
            elif timeseries == "offset_c":
                new_df["s_offset"] = tmp_df["offset_c"]
            elif timeseries == "trend_c":
                new_df["detrend_c"] = tmp_df["trend_c"]
            elif timeseries == "tide_c":
                new_df["tide_c"] = tmp_df["tide_c"]
            elif timeseries == "atmp_c":
                new_df["atmp_c"] = tmp_df["atmp_c"]
                new_df["atmp_c_quality"] = tmp_df["quality"]

        new_df["atmp"] = df["data"].sort_index().loc["atmp", "hpa"]
        # new_df["MJD"] = df["data"].loc["time_index", "mjd"]

        new_df = new_df.reset_index().rename(columns={"time": "date"})
        new_df["doy"] = new_df["date"].dt.dayofyear
        new_df["MJD"] = pd.DatetimeIndex(new_df["date"]).to_julian_date() - 2400000.5
        # logger.info(f"\n{new_df}")

        ordered_columns = [
            "strain",
            "date",
            "doy",
            "MJD",
            f"{data_type}(mstrain)",
            "s_offset",
            "strain_quality",
            "tide_c",
            "detrend_c",
            "atmp_c",
            "atmp_c_quality",
            "level",
            "version",
            "atmp",
        ]
        new_df = new_df[ordered_columns]
        output = f"{outputdir}/{station}.{date_range_string}.{data_type}.txt.gz"
        filenames.append(output)
        if print_it:
            logger.info(f"\n{new_df}")
        if write_it:
            new_df.to_csv(output, sep="\t", index=False)

    if write_it:
        # Create Tar file
        tarfilename = f"{outputdir}/{station}.bsm.level2.{date_range_string}.tar"
        logger.info(f"writing to {tarfilename}")
        tar = tarfile.open(tarfilename, "w")

        for f in filenames:
            tar.add(f, arcname=f.split("/")[-1])

        tar.close()


def write_date_range(
    network, station, start_str, end_str, write_it=True, print_it=False
):
    df = read_date_range(network, station, start_str, end_str)
    logger.info(f"\n{df}")
    date_range_string = f"{start_str}-{end_str}"
    os.makedirs(outputdir, exist_ok=True)
    write_ascii(df, station, date_range_string, write_it=write_it, print_it=print_it)

src/earthscopestraintools/tdb2ascii.py
EarthScope/earthscopestraintools
# write tarred tiledb arrays from tiledb arrays
# slice by start and end time
# use in lambda with web service to deliver data in tiledb
# only works for processed strain data currently

import numpy as np
from earthscopestraintools.tiledbtools import (
    ProcessedStrainReader,
    ProcessedStrainWriter,
)

import logging

logging.basicConfig(
    level=logging.INFO,
    format="%(asctime)s %(levelname)s: %(message)s",
    datefmt="%Y-%m-%d %H:%M:%S",
)
logger = logging.getLogger(__name__)
workdir = "arrays"

#
# def to_date(datetime64):
#     ts = pd.to_datetime(str(datetime64))
#     d = ts.strftime("%Y%m%d")
#     return d


def read_date_range(uri, start_str, end_str):
    logger.info(f"Array uri: {uri}")
    reader = ProcessedStrainReader(uri)
    data_types = reader.array.get_data_types()
    timeseries = reader.array.get_timeseries()
    attrs = ["data", "quality", "level", "version"]
    df = reader.to_df(
        data_types=data_types,
        timeseries=timeseries,
        attrs=attrs,
        start_str=start_str,
        end_str=end_str,
        reindex=False,
    )
    return df


def convert_time_to_unix_ms(df):
    df = df.reset_index()
    df["time"] = df["time"].astype(int) / 10**6
    df["time"] = df["time"].astype(np.int64)
    return df


def export_date_range(
    uri, start_str, end_str, save_as=None, write_it=True, print_it=False
):
    reader = ProcessedStrainReader(uri=uri)
    network = reader.array.get_network()
    station = reader.array.get_station()
    period = reader.array.get_period()

    df = read_date_range(uri, start_str, end_str)
    df = convert_time_to_unix_ms(df)

    if write_it:
        if save_as:
            uri2 = save_as
        else:
            uri2 = f"{workdir}/{network}_{station}_level2_{start_str}_{end_str}.tdb"
        writer = ProcessedStrainWriter(uri=uri2)
        writer.array.delete()
        writer.array.create(schema_type="3D", schema_source="s3")
        writer.array.set_array_meta(network=network, station=station, period=period)
        writer.write_df_to_tiledb(df)

    if print_it:
        reader = ProcessedStrainReader(uri=uri2)
        logger.info(f"Network: {reader.array.get_network()}")
        logger.info(f"Station: {reader.array.get_station()}")
        logger.info(f"Period: {reader.array.get_period()}")
        data_types = reader.array.get_data_types()
        timeseries = reader.array.get_timeseries()
        attrs = ["data", "quality", "level", "version"]
        df = reader.to_df(
            data_types=data_types,
            timeseries=timeseries,
            attrs=attrs,
            start_str=start_str,
            end_str=end_str,
            reindex=False,
        )
        logger.info(f"\n{df}")

src/earthscopestraintools/tdb2tdb.py
EarthScope/earthscopestraintools
import datetime
import tiledb
import numpy as np
import pandas as pd
import json
import logging
from typing import Union
from earthscopestraintools.timeseries import Timeseries
from earthscopestraintools.datasources_api_interact import get_session_edid

logger = logging.getLogger(__name__)

DATASOURCES_API_URL = "https://datasources-api.prod.earthscope.org"

def str_to_unix_ms(time_string: str):
    # convert string into unix ms
    return np.datetime64(time_string, "ms").astype(int)


class StrainArray:
    def __init__(self, uri: str, period: float = None):
        self.set_default_ctx()
        self.uri = uri
        self.period = period

    def default_config(self):
        # return a tiledb config
        config = tiledb.Config()
        config["vfs.s3.region"] = "us-east-2"
        config["vfs.s3.scheme"] = "https"
        config["vfs.s3.endpoint_override"] = ""
        config["vfs.s3.use_virtual_addressing"] = "true"
        config["sm.consolidation.mode"] = "fragment_meta"
        config["sm.vacuum.mode"] = "fragment_meta"
        return config

    def set_default_ctx(self):
        config = self.default_config()
        self.ctx = tiledb.Ctx(config=config)

    def get_schema_from_s3(self, schema_type):
        bucket = "earthscope-tiledb-schema-dev-us-east-2-ebamji"
        s3_schema_uri = f"s3://{bucket}/STRAIN_SCHEMA_{schema_type}.tdb"
        logger.info(f"Using schema {s3_schema_uri}")
        config = tiledb.Config()
        config["vfs.s3.region"] = "us-east-2"
        config["vfs.s3.scheme"] = "https"
        config["vfs.s3.endpoint_override"] = ""
        config["vfs.s3.use_virtual_addressing"] = "true"
        config["sm.consolidation.mode"] = "fragment_meta"
        config["sm.vacuum.mode"] = "fragment_meta"

        with tiledb.open(s3_schema_uri, "r", config=config) as A:
            schema = A.schema
        return schema

    def get_schema(self, schema_type):
        filters1 = tiledb.FilterList([tiledb.ZstdFilter(level=7)])
        filters2 = tiledb.FilterList(
            [tiledb.ByteShuffleFilter(), tiledb.ZstdFilter(level=7)]
        )
        filters3 = tiledb.FilterList(
            [tiledb.BitWidthReductionFilter(), tiledb.ZstdFilter(level=7)]
        )
        filters4 = tiledb.FilterList(
            [tiledb.DoubleDeltaFilter(), tiledb.ZstdFilter(level=7)]
        )
        filters5 = tiledb.FilterList(
            [tiledb.FloatScaleFilter(1e-6, 0, bytewidth=8), tiledb.ZstdFilter(level=7)]
        )
        filters6 = tiledb.FilterList(
            [
                tiledb.PositiveDeltaFilter(),
                tiledb.BitWidthReductionFilter(),
                tiledb.ZstdFilter(level=7),
            ]
        )

        if schema_type == "3D":
            ## time dimension with micro-second precision and 24 hour tiles, domain 1970 to 2100
            d0 = tiledb.Dim(name="data_type", dtype="ascii", filters=filters1)
            d1 = tiledb.Dim(name="timeseries", dtype="ascii", filters=filters1)
            d2 = tiledb.Dim(
                name="time",
                domain=(0, 4102444800000),
                tile=86400000,
                dtype=np.int64,
                filters=filters4,
            )

            dom = tiledb.Domain(d0, d1, d2)

            a0 = tiledb.Attr(name="data", dtype=np.float64, filters=filters1)
            a1 = tiledb.Attr(name="quality", dtype="ascii", var=True, filters=filters1)
            a2 = tiledb.Attr(name="level", dtype="ascii", var=True, filters=filters1)
            a3 = tiledb.Attr(name="version", dtype=np.int64, filters=filters1)
            attrs = [a0, a1, a2, a3]

        elif schema_type == "2D_INT":
            ## time dimension with micro-second precision and 24 hour tiles, domain 1970 to 2100
            d0 = tiledb.Dim(name="channel", dtype="ascii", filters=filters1)
            d1 = tiledb.Dim(
                name="time",
                domain=(0, 4102444800000),
                tile=86400000,
                dtype=np.int64,
                filters=filters4,
            )
            dom = tiledb.Domain(d0, d1)

            a0 = tiledb.Attr(name="data", dtype=np.int32, filters=filters4)
            attrs = [a0]

        elif schema_type == "2D_FLOAT":
            ## time dimension with micro-second precision and 24 hour tiles, domain 1970 to 2100
            d0 = tiledb.Dim(name="channel", dtype="ascii", filters=filters1)
            d1 = tiledb.Dim(
                name="time",
                domain=(0, 4102444800000),
                tile=86400000,
                dtype=np.int64,
                filters=filters4,
            )
            dom = tiledb.Domain(d0, d1)

            a0 = tiledb.Attr(name="data", dtype=np.float64, filters=filters1)
            attrs = [a0]

        else:
            raise NameError

        schema = tiledb.ArraySchema(
            domain=dom,
            sparse=True,
            attrs=attrs,
            cell_order="row-major",
            tile_order="row-major",
            capacity=100000,
            offsets_filters=filters6,
        )

        return schema

    def exists(self):
        # may not work with S3 arrays
        if tiledb.array_exists(self.uri):
            return True
        else:
            return False

    def create(self, schema_type: str, schema_source: str = "s3"):
        # schema should be one of 2D_INT, 2D_FLOAT, 3D. case insensitive
        schema_type = schema_type.upper()
        if schema_type in ["2D_INT", "2D_FLOAT", "3D"]:
            if schema_source == "s3":
                self.schema = self.get_schema_from_s3(schema_type)
            else:
                self.schema = self.get_schema(schema_type)
        else:
            logger.error("Invalid schema.  Must be one of 2D_INT, 2D_FLOAT, 3D")
            raise NameError

        try:
            tiledb.Array.create(self.uri, self.schema, ctx=self.ctx)
            logger.info(f"Created array at {self.uri}")
        except tiledb.TileDBError as e:
            logger.warning(e)

    def set_array_meta(
        self,
        network: str = None,
        station: str = None,
        period: float = None,
    ):
        # don't overwrite existing metadata
        with tiledb.open(self.uri, "r", ctx=self.ctx) as A:
            keys = A.meta.keys()
            if "network" not in keys:
                if network:
                    logger.info(f"Setting array metadata network={network}")
                    with tiledb.open(self.uri, "w", ctx=self.ctx) as B:
                        B.meta["network"] = network
            if "station" not in keys:
                if station:
                    logger.info(f"Setting array metadata station={station}")
                    with tiledb.open(self.uri, "w", ctx=self.ctx) as B:
                        B.meta["station"] = station
            if "period" not in keys:
                if period:
                    logger.info(f"Setting array metadata period={period}")
                    with tiledb.open(self.uri, "w", ctx=self.ctx) as B:
                        B.meta["period"] = period

    def update_array_meta(
        self,
        network: str = None,
        station: str = None,
        period: float = None,
    ):
        # overwrite existing metadata
        if network:
            with tiledb.open(self.uri, "w", ctx=self.ctx) as A:
                A.meta["network"] = network
        if station:
            with tiledb.open(self.uri, "w", ctx=self.ctx) as A:
                A.meta["station"] = station
        if period:
            with tiledb.open(self.uri, "w", ctx=self.ctx) as A:
                A.meta["period"] = period

    def delete(self):
        try:
            tiledb.remove(self.uri, ctx=self.ctx)
            print("Deleted ", self.uri)
        except tiledb.TileDBError as e:
            print(e)

    def consolidate_fragment_meta(self, print_it=True):
        config = self.ctx.config()
        config["sm.consolidation.mode"] = "fragment_meta"
        ctx = tiledb.Ctx(config=config)
        tiledb.consolidate(self.uri, ctx=ctx)
        if print_it:
            logger.info("consolidated fragment_meta")

    def consolidate_array_meta(self, print_it=True):
        config = self.ctx.config()
        config["sm.consolidation.mode"] = "array_meta"
        ctx = tiledb.Ctx(config=config)
        tiledb.consolidate(self.uri, ctx=ctx)
        if print_it:
            logger.info("consolidated array_meta")

    def consolidate_fragments(self, print_it=True):
        config = self.ctx.config()
        config["sm.consolidation.mode"] = "fragments"
        ctx = tiledb.Ctx(config=config)
        tiledb.consolidate(self.uri, ctx=ctx)
        if print_it:
            logger.info("consolidated fragments")

    def vacuum_fragment_meta(self, print_it=True):
        config = self.ctx.config()
        config["sm.vacuum.mode"] = "fragment_meta"
        ctx = tiledb.Ctx(config=config)
        tiledb.vacuum(self.uri, ctx=ctx)
        if print_it:
            logger.info("vacuumed fragment_meta")

    def vacuum_array_meta(self, print_it=True):
        config = self.ctx.config()
        config["sm.vacuum.mode"] = "array_meta"
        ctx = tiledb.Ctx(config=config)
        tiledb.vacuum(self.uri, ctx=ctx)
        if print_it:
            logger.info("vacuumed array_meta")

    def vacuum_fragments(self, print_it=True):
        config = self.ctx.config()
        config["sm.vacuum.mode"] = "fragments"
        ctx = tiledb.Ctx(config=config)
        tiledb.vacuum(self.uri, ctx=ctx)
        if print_it:
            logger.info("vacuumed fragments")

    def cleanup_meta(self):
        self.consolidate_array_meta(print_it=False)
        self.vacuum_array_meta(print_it=False)
        self.consolidate_fragment_meta(print_it=False)
        self.vacuum_fragment_meta(print_it=False)
        logger.info("Consolidated and vacuumed metadata")

    def cleanup(self):
        self.consolidate_array_meta(print_it=False)
        self.vacuum_array_meta(print_it=False)
        self.consolidate_fragment_meta(print_it=False)
        self.vacuum_fragment_meta(print_it=False)
        self.consolidate_fragments(print_it=False)
        self.vacuum_fragments(print_it=False)
        logger.info("Consolidated and vacuumed fragments and metadata")

    def get_nonempty_domain(self):
        with tiledb.open(self.uri, "r", ctx=self.ctx) as A:
            return A.nonempty_domain()[-1][0], A.nonempty_domain()[-1][1]

    def get_data_types(self):
        with tiledb.open(self.uri, "r", ctx=self.ctx) as A:
            return json.loads(A.meta["dimensions"])["data_types"]

    def get_timeseries(self):
        with tiledb.open(self.uri, "r", ctx=self.ctx) as A:
            return json.loads(A.meta["dimensions"])["timeseries"]

    def get_channels(self):
        with tiledb.open(self.uri, "r", ctx=self.ctx) as A:
            return json.loads(A.meta["channels"])["channels"]

    def get_network(self):
        with tiledb.open(self.uri, "r", ctx=self.ctx) as A:
            try:
                return A.meta["network"]
            except KeyError:
                return None

    def get_station(self):
        with tiledb.open(self.uri, "r", ctx=self.ctx) as A:
            try:
                return A.meta["station"]
            except KeyError:
                return None

    def get_period(self):
        with tiledb.open(self.uri, "r", ctx=self.ctx) as A:
            try:
                return A.meta["period"]
            except KeyError:
                return None

    def print_schema(self):
        with tiledb.open(self.uri, "r", ctx=self.ctx) as A:
            print(A.schema)


class ProcessedStrainReader:
    def __init__(self, uri: str, period=300):
        self.array = StrainArray(uri=uri, period=period)

    def to_df(
        self,
        data_types: Union[list, str],
        timeseries: Union[list, str],
        attrs: Union[list, str],
        start_ts: int = None,
        end_ts: int = None,
        start_str: str = None,
        end_str: str = None,
        start_dt: datetime.datetime = None,
        end_dt: datetime.datetime = None,
        reindex=True,
        print_array_range=False,
    ):
        # generic read function. requires all query parameters, accepts strings or lists
        # reindex handles one or more data_types
        # reindex does not handle more than one timeseries or attr, so should be set to False

        try:
            if not start_ts:
                if start_str:
                    start_ts = str_to_unix_ms(start_str)
                elif start_dt:
                    start_ts = int(
                        start_dt.replace(tzinfo=datetime.timezone.utc).timestamp()
                        * 1000
                    )
                else:
                    logger.error("No start time provided for read")
            if not end_ts:
                if end_str:
                    end_ts = str_to_unix_ms(end_str)
                elif start_dt:
                    end_ts = int(
                        end_dt.replace(tzinfo=datetime.timezone.utc).timestamp() * 1000
                    )
                else:
                    logger.error("No start time provided for read")
            with tiledb.open(self.array.uri, "r", ctx=self.array.ctx) as A:
                if print_array_range:
                    logger.info(
                        f"Array date range: {A.nonempty_domain()[2][0]} to {A.nonempty_domain()[2][1]}"
                    )
                    logger.info(f"Query range: {start_ts} to {end_ts}")
                index_col = ["data_type", "timeseries", "time"]
                data_types = [data_types] if isinstance(data_types, str) else data_types
                timeseries = [timeseries] if isinstance(timeseries, str) else timeseries
                attrs = [attrs] if isinstance(attrs, str) else attrs
                df = A.query(index_col=index_col, attrs=attrs).df[
                    data_types, timeseries, start_ts:end_ts
                ]
                df.index = df.index.set_levels(
                    pd.to_datetime(df.index.levels[2], unit="ms"), level=2
                )
            if not isinstance(df, pd.DataFrame):
                df = df.to_frame()
            epochs = len(df) / len(data_types)
            self.check_query_result(epochs, start_ts, end_ts)
            if reindex:
                df = self.reindex_df(df, columns=data_types, attr=attrs[0])

            return df

        except (IndexError, KeyError) as e:  # tiledb.TileDBError
            logger.error(e)
            logger.error("No data found matching query parameters")
            return pd.DataFrame()

    def reindex_df(self, df: pd.DataFrame, columns: list = [], attr="data"):
        # removes a multi-index and makes each data_type a column
        data_types = columns
        for data_type in data_types:
            df_data_type = df.xs(data_type, level="data_type")[attr].droplevel(level=0)
            df_data_type.name = data_type
            if data_type == data_types[0]:
                df2 = df_data_type
            else:
                df2 = pd.concat([df2, df_data_type], axis=1)
        return df2

    def check_query_result(self, epochs, start, end):
        if self.array.period is not None:
            expected_epochs = (
                int((end - start) / 1000 / self.array.period) + 1
            )  # inclusive of last epoch
            if epochs >= expected_epochs:
                logger.info(
                    f"Query complete, expected {expected_epochs} epochs and returned {epochs}"
                )
            else:
                logger.info(
                    f"Query incomplete, expected {expected_epochs} epochs and returned {epochs}"
                )
        else:
            logger.info(f"Cannot check query completess without array period")

    def to_ts(
        self,
        data_types: Union[list, str],
        timeseries: str,
        attrs: list = None,
        name: str = "timeseries",
        units: str = None,
        start_ts: int = None,
        end_ts: int = None,
        start_str: str = None,
        end_str: str = None,
        start_dt: datetime.datetime = None,
        end_dt: datetime.datetime = None,
        print_array_range=False,
        to_na=True,
    ):
        if attrs is None:
            attrs = ["data", "quality", "level", "version"]
        if units is None:
            if timeseries == "counts":
                units = "counts"
            elif isinstance(data_types, str):
                if data_types == "atmp":
                    units = "hpa"
                else:
                    units = "microstrain"
            elif isinstance(data_types, list):
                if "atmp" in data_types:
                    units = "hpa"
                else:
                    units = "microstrain"

        # print(units)
        df = self.to_df(
            data_types=data_types,
            timeseries=timeseries,
            attrs=attrs,
            start_ts=start_ts,
            end_ts=end_ts,
            start_str=start_str,
            end_str=end_str,
            start_dt=start_dt,
            end_dt=end_dt,
            reindex=False,
            print_array_range=print_array_range,
        )
        # print(df)
        data = self.reindex_df(df=df, columns=data_types, attr="data")
        quality_df = self.reindex_df(df=df, columns=data_types, attr="quality")
        level = df.iloc[0]["level"]
        ts = Timeseries(
            data=data,
            quality_df=quality_df,
            series=timeseries,
            units=units,
            period=self.array.period,
            level=level,
            name=name,
        )
        if to_na:
            ts = ts.remove_999999s(interpolate=False)
        return ts


class ProcessedStrainWriter:
    def __init__(self, uri: str):
        self.array = StrainArray(uri=uri)

    def create(self, schema_type: str = "3D", schema_source: str = "s3"):
        if not self.array.exists():
            self.array.create(schema_type=schema_type, schema_source=schema_source)
        else:
            logger.info(f"Array exists at {self.array.uri}")

    def write_df_to_tiledb(self, df: pd.DataFrame):
        mode = "append"
        # make sure there aren't any nans, replace with 999999 if so
        df = df.replace(np.nan, 999999)
        logger.debug("writing buffer")
        tiledb.from_pandas(
            uri=self.array.uri,
            dataframe=df,
            index_dims=["data_type", "timeseries", "time"],
            mode=mode,
            ctx=self.array.ctx,
        )
        logger.debug("buffer written")
        # update the string dimension metadata
        data_type = df["data_type"].unique()
        timeseries = df["timeseries"].unique()
        if type(data_type) == str:
            data_type = [data_type]
        if type(timeseries) == str:
            timeseries = [timeseries]
        with tiledb.open(self.array.uri, "r", ctx=self.array.ctx) as A:
            try:
                dimension_json = A.meta["dimensions"]
            except KeyError:
                dimension_json = '{"data_types":[], "timeseries":[]}'

            dimension_dict = json.loads(dimension_json)
            # print(dimension_dict)
            for item in data_type:
                if item not in dimension_dict["data_types"]:
                    dimension_dict["data_types"].append(item)
            for item in timeseries:
                if item not in dimension_dict["timeseries"]:
                    dimension_dict["timeseries"].append(item)

            # print(dimension_dict)
            with tiledb.open(self.array.uri, "w", ctx=self.array.ctx) as A:
                A.meta["dimensions"] = json.dumps(dimension_dict)

    def df_2_tiledb(
        self,
        df: pd.DataFrame,
        data_types: list,
        timeseries: str,
        level: str,
        quality_df: pd.DataFrame = None,
        print_it: bool = False,
    ):
        """
        prepares a dataframe to write to array schema
        df - dataframe with time index, columns are timeseries data (one column per data_type)
        uri - string. which tiledb array to write to
        data_types - list of strings. data_types to map columns into.  ie CH0, 2Ene, pressure, time_index
        timeseries - string. name of timeseries.  counts, microstrain, offset_c, tide_c, atmp_c, atmp, pore, mjd, doy
        level - 2char string.  '0a', '1a', '2a', '2b'
        quality_df- dataframe, optional.  if not included, quality flags will be set to 'g'
        print_it - bool, optional.  show the constructed dataframe as it is being written to tiledb.
        """
        logger.debug("creating buffer")
        df_buffer = pd.DataFrame(
            columns=[
                "data_type",
                "timeseries",
                "time",
                "data",
                "level",
                "quality",
                "version",
            ]
        )
        if quality_df is None:
            quality_df = pd.DataFrame(index=df.index)
            for ch in data_types:
                quality_df[ch] = "g"
        for ch in data_types:
            data = df[ch].values
            # convert datetimeindex to unix ms
            timestamps = df.index.astype(int) / 10**6
            version = int(datetime.datetime.now().strftime("%Y%j%H%M%S"))
            quality = quality_df[ch].values

            d = {
                "data_type": ch,
                "timeseries": timeseries,
                "time": timestamps,
                "data": data,
                "level": level,
                "quality": quality,
                "version": version,
            }
            ch_df = pd.DataFrame(data=d)
            # ch_df.loc[ch_df['data'] == 999999, 'quality'] = 'm'
            df_buffer = pd.concat([df_buffer, ch_df], axis=0).reset_index(drop=True)
            df_buffer["time"] = df_buffer["time"].astype(np.int64)
            df_buffer["data"] = df_buffer["data"].astype(np.float64)
            df_buffer["version"] = df_buffer["version"].astype(np.int64)

        if print_it:
            print(df_buffer)
            print(df_buffer.isna().sum())
        self.write_df_to_tiledb(df_buffer)

    def ts_2_tiledb(self, ts: Timeseries, cleanup: bool = False):
        self.df_2_tiledb(
            df=ts.data,
            data_types=ts.columns,
            timeseries=ts.series,
            level=ts.level,
            quality_df=ts.quality_df,
            print_it=False,
        )
        logger.info(f"Wrote {ts.columns} {ts.series} to {self.array.uri}")
        if cleanup:
            self.array.cleanup()


class RawStrainReader:
    def __init__(self, uri: str, period=None):
        self.array = StrainArray(uri=uri, period=period)

    def to_df(
        self,
        channels: Union[list, str],
        start_ts: int = None,
        end_ts: int = None,
        start_str: str = None,
        end_str: str = None,
        start_dt: datetime.datetime = None,
        end_dt: datetime.datetime = None,
        print_array_range=False,
    ):
        # generic read function. requires all query parameters, accepts strings or lists
        try:
            if not start_ts:
                if start_str:
                    logger.info(f"Query start: {start_str}")
                    start_ts = str_to_unix_ms(start_str)
                elif start_dt:
                    start_ts = int(
                        start_dt.replace(tzinfo=datetime.timezone.utc).timestamp()
                        * 1000
                    )
                else:
                    logger.error("No start time provided for read")
            if not end_ts:
                if end_str:
                    logger.info(f"Query end: {end_str}")
                    end_ts = str_to_unix_ms(end_str)
                elif start_dt:
                    end_ts = int(
                        end_dt.replace(tzinfo=datetime.timezone.utc).timestamp() * 1000
                    )
                else:
                    logger.error("No start time provided for read")
            logger.info(f"Channel query range: {channels}")
            logger.info(f"Time query range: {start_ts} to {end_ts}")
            with tiledb.open(self.array.uri, "r", ctx=self.array.ctx) as A:
                if print_array_range:
                    logger.info(
                        f"Array date range: {A.nonempty_domain()[1][0]} to {A.nonempty_domain()[1][1]}"
                    )
                index_col = ["channel", "time"]
                channels = [channels] if isinstance(channels, str) else channels
                attrs = ["data"]
                # logger.info("query start")
                df = A.query(index_col=index_col, attrs=attrs).df[
                    channels, start_ts:end_ts
                ]
                # logger.info("query stop")
                df.index = df.index.set_levels(
                    pd.to_datetime(df.index.levels[1], unit="ms"), level=1
                )
            df = self.reindex_df(df, columns=channels, attr="data")
            if not isinstance(df, pd.DataFrame):
                df = df.to_frame()
            self.check_query_result(df, start_ts, end_ts)
            return df

        except (IndexError, KeyError, tiledb.TileDBError) as e:  #
            logger.error(f"{type(e)}: {e}")
            logger.error("No data found matching query parameters")
            return pd.DataFrame()

    def reindex_df(self, df: pd.DataFrame, columns: list = [], attr="data"):
        # removes a multi-index and makes each data_type a column
        channels = columns
        for channel in channels:
            df_channel = df.xs(channel, level="channel")[attr]  # .droplevel(level=0)
            df_channel.name = channel
            if channel == channels[0]:
                df2 = df_channel
            else:
                df2 = pd.concat([df2, df_channel], axis=1)
        return df2

    def check_query_result(self, df, start, end):
        if self.array.period is not None:
            expected_samples = (
                int((end - start) / 1000 / self.array.period) + 1
            )  # inclusive of last epoch
            # expected_samples = int((str_to_unix_ms(end) - str_to_unix_ms(start)) / 1000 / self.array.period)
            # print("expected samples:", expected_samples)
            if len(df) >= expected_samples:
                logger.info(
                    f"Query complete, expected {expected_samples} and returned {len(df)}"
                )
            else:
                logger.info(
                    f"Query incomplete, expected {expected_samples} and returned {len(df)}"
                )
        else:
            logger.info(f"Cannot check query completess without array period")

    def to_ts(
        self,
        channels: list,
        units: str,
        start_ts: int = None,
        end_ts: int = None,
        start_str: str = None,
        end_str: str = None,
        start_dt: datetime.datetime = None,
        end_dt: datetime.datetime = None,
        to_nan: bool = True,
        name: str = "",
    ):
        # if not uri:
        #     uri = lookup_s3_uri(network, station, period)
        # reader = RawStrainReader(uri, period)
        df = self.to_df(
            channels=channels,
            start_str=start_str,
            end_str=end_str,
            start_ts=start_ts,
            end_ts=end_ts,
            start_dt=start_dt,
            end_dt=end_dt,
        )
        series = "raw"
        level = "0"
        # if name == None:
        #    name = f"{network}.{station}.{series}.{units}"
        ts = Timeseries(
            data=df,
            series=series,
            units=units,
            level=level,
            period=self.array.period,
            name=name,
        )
        if to_nan:
            logger.info("Converting missing data from 999999 to nan")
            return ts.remove_999999s()
        else:
            return ts


class RawStrainWriter:
    def __init__(self, uri: str, array_type="int"):
        self.array = StrainArray(uri=uri)
        self.array_type = array_type

    def create(self, schema_type: str, schema_source: str = "s3"):
        if not self.array.exists():
            self.array.create(schema_type=schema_type, schema_source=schema_source)
        else:
            logger.info(f"Array exists at {self.array.uri}")

    def write_df_to_tiledb(self, df: pd.DataFrame):
        mode = "append"
        # make sure there arent any nans, replace with 999999 if so
        df = df.replace(np.nan, 999999)
        tiledb.from_pandas(
            uri=self.array.uri,
            dataframe=df,
            index_dims=["channel", "time"],
            mode=mode,
            ctx=self.array.ctx,
        )
        # update the string dimension metadata
        channels = df["channel"].unique()

        if type(channels) == str:
            channels = [channels]

        with tiledb.open(self.array.uri, "r", ctx=self.array.ctx) as A:
            try:
                channels_json = A.meta["channels"]
            except KeyError:
                channels_json = '{"channels":[]}'

            channels_dict = json.loads(channels_json)
            for item in channels:
                if item not in channels_dict["channels"]:
                    channels_dict["channels"].append(item)

            with tiledb.open(self.array.uri, "w", ctx=self.array.ctx) as A:
                A.meta["channels"] = json.dumps(channels_dict)
        logger.info("Write complete")

    def df_2_tiledb(
        self,
        df: pd.DataFrame,
        print_it: bool = False,
    ):
        """
        prepares a dataframe to write to array schema
        df - dataframe with time index, columns are timeseries data (one column per channel)
        print_it - bool, optional.  show the constructed dataframe as it is being written to tiledb.
        """
        df_buffer = pd.DataFrame(columns=["channel", "time", "data"])
        for ch in df.columns:
            data = df[ch].values
            # convert datetimeindex to unix ms
            timestamps = df.index.astype(int) / 10**6

            d = {"channel": ch, "time": timestamps, "data": data}
            ch_df = pd.DataFrame(data=d)
            df_buffer = pd.concat([df_buffer, ch_df], axis=0).reset_index(drop=True)
            if self.array_type == "int":
                df_buffer["data"] = df_buffer["data"].astype(np.int32)
            elif self.array_type == "float":
                df_buffer["data"] = df_buffer["data"].astype(np.float64)
            df_buffer["time"] = df_buffer["time"].astype(np.int64)
        if print_it:
            print(df_buffer)
        self.write_df_to_tiledb(df_buffer)
        self.array.cleanup_meta()

    def ts_2_tiledb(self, ts: Timeseries, cleanup: bool = False):
        self.df_2_tiledb(
            df=ts.data,
            print_it=False,
        )
        logger.info(f"Wrote {ts.columns} to {self.array.uri}")
        if cleanup:
            self.array.cleanup()


# def ts_from_tiledb_processed(
#     network: str,
#     station: str,
#     channels: list,
#     series: str,
#     units: str = "",
#     period: int = 300,
#     uri: str = None,
#     start_ts: int = None,
#     end_ts: int = None,
#     start_str: str = None,
#     end_str: str = None,
#     start_dt: datetime.datetime = None,
#     end_dt: datetime.datetime = None,
#     name: str = None,
# ):
#     reader = ProcessedStrainReader(uri, period)
#     df = reader.to_df(
#         data_types=channels,
#         timeseries=series,
#         start_str=start_str,
#         end_str=end_str,
#         start_ts=start_ts,
#         end_ts=end_ts,
#         start_dt=start_dt,
#         end_dt=end_dt,
#     )
#     if name == None:
#         name = f"{network}.{station}.{series}.{units}"
#     ts = Timeseries(
#         data=df, series=series, units=units, level=level, period=period, name=name
#     )
#     logger.info("Converting missing data from 999999 to nan")
#     return ts.remove_999999s()
#


def lookup_s3_uri(network, station, period):
    bucket = "tiledb-strain"
    if period == 600:
        session = "Day"
    elif period == 1:
        session = "Hour"
    elif period == 0.05:
        session = "Min"
    edid = get_session_edid(api_url=DATASOURCES_API_URL, 
                            namespace="BNUM", 
                            station=station, 
                            session=session)
    uri = f"s3://{bucket}/{edid}.tdb"
    return uri

src/earthscopestraintools/tiledbtools.py
EarthScope/earthscopestraintools
import os
import numpy as np
import pandas as pd
import datetime

from earthscopestraintools.processing import (
    linearize,
    interpolate,
    decimate_1s_to_300s,
    butterworth_filter,
    apply_calibration_matrix,
    calculate_offsets,
    calculate_pressure_correction,
    calculate_tide_correction,
    calculate_linear_trend_correction,
)
from earthscopestraintools.event_processing import dynamic_strain, calculate_magnitude
from scipy import signal, stats
import matplotlib.pyplot as plt
from matplotlib import cm
import logging

logger = logging.getLogger(__name__)


def test():
    print("test1")


class Timeseries:
    """
    Class for storing strainmeter data or correction timeseries.  Each Timeseries object 
    contains the following attributes: \n
    - data: (pd.DataFrame) with datetime index and one or more columns of timeseries data \n
    - quality_df: (pd.DataFrame) with same shape as data, but with a character mapped to each
      data point.  flags include "g"=good, "m"=missing, "i"=interpolated, "b"=bad \n
    - series: (str) description of timeseries, ie 'raw', 'microstrain', 'atmp_c', 'tide_c', 'offset_c', 'trend_c' \n
    - units: (str) units of timeseries \n
    - level: (str) level of data. ie. '0','1','2a','2b' \n
    - period: (float) sample period of data \n
    - name: (str) name of timeseries, including station name.  useful for showing stats.  
      defaults to network.station \n
    - network: (str) FDSN two character network code\n
    - station: (str) FDSN four character station code
    """
    def __init__(
        self,
        data: pd.DataFrame = None,
        quality_df: pd.DataFrame = None,
        series: str = "",
        units: str = "",
        level: str = "",
        period: float = 0,
        name: str = None,
        network: str = "",
        station: str = "",
    ):
        if data is not None:
            self.data = data
            self.columns = list(data.columns)
        else:
            self.data = pd.DataFrame()
            self.columns = []

        if quality_df is not None:
            self.quality_df = quality_df
        else:
            self.quality_df = self.set_initial_quality_flags()

        self.series = series
        self.units = units
        self.level = level
        self.period = period
        self.network = network
        self.station = station
        if name:
            self.name = name
        else:
            self.name = f"{self.network}.{self.station}"
        self.check_for_gaps()

    def set_data(self, df):
        self.data = df

    def set_period(self, period):
        self.period = period

    def set_units(self, units):
        self.units = units

    def set_local_tdb_uri(self, local_tdb_uri):
        if self.local_archive:
            local_tdb_uri = f"{self.local_archive}/{local_tdb_uri}"
        self.local_tdb_uri = local_tdb_uri

    def set_s3_tdb_uri(self, s3_tdb_uri):
        self.s3_tdb_uri = s3_tdb_uri

    def set_initial_quality_flags(self, missing_data=999999):
        """used to flag any missing data

        :param missing_data: value used to represent gaps, defaults to 999999
        :type missing_data: int, optional
        :return: dataframe of quality flags
        :rtype: pandas.DataFrame
        """
        qual_df = pd.DataFrame(index=self.data.index)
        for ch in self.columns:
            qual_df[ch] = "g"
            qual_df[ch][self.data[ch] == missing_data] = "m"

        # print(qual_df.value_counts())
        return qual_df

    def set_initial_level_flags(self, level):
        """used to set level flag 

        :param level: level of data ['0','1','2a','2b']
        :type level: str, optional
        :return: dataframe of level flags
        :rtype: pandas.DataFrame
        """
        level_df = pd.DataFrame(index=self.data.index)
        for ch in self.columns:
            level_df[ch] = level
        return level_df

    def set_initial_version(self):
        """adds a version timestamp to each data point

        :return: dataframe of version flags
        :rtype: pandas.DataFrame
        """
        version_df = pd.DataFrame(index=self.data.index)
        version = int(datetime.datetime.now().strftime("%Y%j%H%M%S"))
        for ch in self.columns:
            version_df[ch] = version
        return version_df

    def check_for_gaps(self):
        """generates some statistics around nans, fill values, and missing epochs.
        """
        if self.period != 0 and self.data is not None:
            days = (self.data.index[-1] - self.data.index[0]).days
            seconds = (self.data.index[-1] - self.data.index[0]).seconds
            self.nans = int(self.data.isna().sum().sum() / len(self.data.columns))
            self.nines = round(
                self.data[self.data == 999999].count().sum() / len(self.data.columns), 2
            )
            expected_epochs = int((days * 86400 + seconds) / self.period) + 1
            self.epochs = len(self.data) - self.nans - self.nines
            # print(self.nans, self.nines, self.epochs, expected_epochs)
            self.gap_percentage = round((1 - (self.epochs / expected_epochs)) * 100, 2)

            logger.info(
                f"    Found {self.nans} epochs with nans, {self.nines} epochs with 999999s, and "
                f"{expected_epochs - len(self.data)} missing epochs.\n    Total missing data is "
                f"{self.gap_percentage}%"
            )
        else:
            self.gaps = None
            self.gap_percentage = None

    def stats(self):
        """displays summary information describing the Timeseries object
        """
        if len(self.data):
            outputstring = f"{self.name}\n{'':4}| Channels: {str(self.columns):40} "
            outputstring += f"\n{'':4}| TimeRange: {self.data.index[0]} - {self.data.index[-1]}        | Period: {self.period:>13}s"
            outputstring += f"\n{'':4}| Series: {self.series:>11}| Units: {self.units:>13}| Level: {self.level:>10}| Gaps: {self.gap_percentage:>15}% "
            if len(self.quality_df):
                cols = len(self.quality_df.columns)
                outputstring += f"\n{'':4}| Epochs: {len(self.quality_df):11}"
                outputstring += (
                    f"| Good: {(self.quality_df == 'g').sum().sum() / cols:14}"
                )
                outputstring += (
                    f"| Missing: {(self.quality_df == 'm').sum().sum() / cols:8}"
                )
                outputstring += (
                    f"| Interpolated: {(self.quality_df == 'i').sum().sum() / cols:8}"
                )
                outputstring += f"\n{'':4}| Samples: {len(self.quality_df) * cols:10}"
                outputstring += f"| Good: {(self.quality_df == 'g').sum().sum():14}"
                outputstring += f"| Missing: {(self.quality_df == 'm').sum().sum():8}"
                outputstring += (
                    f"| Interpolated: {(self.quality_df == 'i').sum().sum():8}"
                )
            logger.info(f"{outputstring}")

    def show_flags(self):
        """returns a dataframe with all flags that are not 'g'

        :return: times and channels with flagged data within the timeseries
        :rtype: pandas.DataFrame
        """
        return self.quality_df[self.quality_df[self.quality_df != "g"].any(axis=1)]

    def show_flagged_data(self):
        """returns dataframe containing any data with a flag other than 'g'

        :return: data that has been flagged
        :rtype: pandas.DataFrame
        """
        return self.data[self.quality_df[self.quality_df != "g"].any(axis=1)]


    def save_csv(self, filename: str, datadir: str = "./", sep=",", compression=None):
        """save data attribute as csv.  flattens object, does not save quality flags, level, or version information

        :param filename: name of csv file to save
        :type filename: str
        :param datadir: path to local directory to save file, defaults to "./"
        :type datadir: str, optional
        :param sep: separator to use in csv, defaults to ","
        :type sep: str, optional
        :param compression: compression algorthim ['infer', 'gzip', 'bz2', 'zip', 'xz', 'zstd'], defaults to None
        :type compression: str, optional
        """
        filepath = os.path.join(datadir, filename)
        logger.info(f"saving to {filepath}")
        if compression:
            self.data.to_csv(filepath, sep=sep, compression={"method": compression})
        else:
            self.data.to_csv(filepath, sep=sep)

    # def save_tdb(self):
    # TODO: update, handle quality fields
    # if self.local_tdb_uri:
    #     logger.info(f"saving to {self.local_tdb_uri}")
    #     self.array = StrainArray(uri=self.local_tdb_uri, period=self.period)
    #     if not tiledb.array_exists(self.array.uri):
    #         logger.info(f"Creating array {self.array.uri}")
    #         self.array.create()
    #     self.array.df_2_tiledb(
    #         df=self.data,
    #         data_types=self.columns,
    #         timeseries=self.series,
    #         level=self.level,
    #         quality_df=self.quality_df,
    #     )
    #     self.array.consolidate_fragments()
    #     self.array.vacuum_fragments()
    # else:
    #     logger.error(
    #         "Error, no local array specified.  Set with Timeseries.set_local_tdb_uri()"
    #     )

    def remove_999999s(
        self,
        interpolate: bool = False,
        method: str = "linear",
        limit_direction: str = "both",
        limit: any = None,
    ):
        """remove 999999 gap fill values from data, options to either replace with nans or interpolate

        :param interpolate: boolean of whether to interpolate across gaps using pd.DataFrame.interpolate(), defaults to False
        :type interpolate: bool, optional
        :param method: interpolation method from pd.DataFrame.interpolate(), defaults to "linear"
        :type method: str, optional
        :param limit_direction: limit direction from pd.DataFrame.interpolate(), defaults to "both"
        :type limit_direction: str, optional
        :param limit: limit from pd.DataFrame.interpolate(), defaults to None
        :type limit: any, optional
        :return: Timeseries with 999999 gap fills removed, and appropriate flags set 
        :rtype: Timeseries
        """
        logger.info("  Converting 999999 values to nan")
        if interpolate:
            df = self.data.replace(999999, np.nan).interpolate(
                method=method, limit_direction=limit_direction, limit=limit
            )
            quality_df = self.quality_df.replace("m", "i")
        else:
            df = self.data.replace(999999, np.nan)
            quality_df = self.quality_df

        return Timeseries(
            data=df,
            quality_df=quality_df,
            series=self.series,
            period=self.period,
            units=self.units,
            level=self.level,
            name=self.name,
        )

    def decimate_1s_to_300s(
        self, method: str = "linear", limit: int = 3600, name: str = None
    ):
        """decimate 1hz data to 5 min data using \n
        Agnew, Duncan Carr, and K. Hodgkinson (2007), Designing compact causal digital filters for 
        low-frequency strainmeter data , Bulletin Of The Seismological Society Of America, 97, No. 1B, 91-99

        :param method: method to interpolate across gaps, defaults to "linear"
        :type method: str, optional
        :param limit: largest gap to interpolate, defaults to 3600 samples
        :type limit: int, optional
        :param name: name for new Timeseries, defaults to None
        :type name: str, optional
        :return: Timeseries containing 300s decimated data
        :rtype: Timeseries
        """
        
        data = decimate_1s_to_300s(self.data, method=method, limit=limit)
        quality_df = self.quality_df.copy().reindex(data.index)
        quality_df[quality_df.isna()] = "i"
        # find any differences using the original data index
        mask1 = (data.reindex(self.data.index) != self.data).any(axis=1)

        # any nans from the original index
        mask2 = self.data[mask1].isna()
        quality_df[mask2] = "i"

        # any 999999s from the original index
        mask3 = self.data[mask1] == 999999
        quality_df[mask3] = "i"
        # quality_df = self.quality_df.reindex(data.index)
        if not name:
            name = f"{self.name}.decimated"
        ts2 = Timeseries(
            data=data,
            quality_df=quality_df,
            series=self.series,
            units=self.units,
            level="1",
            period=300,
            name=name,
        )
        return ts2

    def linearize(self, reference_strains: dict, gap: float, name: str = None):
        """Processing step to convert digital counts to microstrain based on geometry of GTSM gauges

        :param reference_strains: dict containing keys of CHX and values of reference strains
        :type reference_strains: dict
        :param gap: instrument gap in meters
        :type gap: float
        :param name: name for new Timeseries, defaults to None
        :type name: str, optional
        :return: Timeseries of linearized data in microstrain
        :rtype: Timeseries
        """

        # remove any 999999 values in data, ok to leave as Nan rather than interpolate.
        if self.nines:
            logger.info(f"Found {self.nines} 999999s, replacing with nans")
            df = self.remove_999999s(interpolate=False)
        else:
            df = self.data
        linearized_data = linearize(df, reference_strains, gap)
        if not name:
            name = f"{self.name}.linearized"
        ts2 = Timeseries(
            data=linearized_data,
            quality_df=self.quality_df,
            series="microstrain",
            units="microstrain",
            level="1",
            period=self.period,
            name=name,
        )
        return ts2

    def interpolate(
        self,
        replace: int = 999999,
        method: str = "linear",
        limit_seconds: int = 3600,
        limit_direction="both",
        name: str = None,
        new_index: pd.DatetimeIndex = None,
        period=None,
        level=None,
        series=None,
    ):
        """Interpolate across gaps in data using pd.DataFrame.interpolate()

        :param replace: gap fill value to interpolate across, defaults to 999999
        :type replace: int, optional
        :param method: interpolation method, defaults to "linear"
        :type method: str, optional
        :param limit_seconds: max gap (in seconds) to interpolate , defaults to 3600
        :type limit_seconds: int, optional
        :param limit_direction: ['forward', 'backward', 'both'], defaults to "both"
        :type limit_direction: str, optional
        :param name: name for new Timeseries, defaults to None
        :type name: str, optional
        :param new_index: option to manually set the index of the interpolated data, defaults to None
        :type new_index: pd.DatetimeIndex, optional
        :param period: sample rate of data in seconds, defaults to None
        :type period: float, optional
        :param level: level of data, defaults to None
        :type level: str, optional
        :param series: series name, defaults to None
        :type series: str, optional
        :return: Timeseries containing interpolated data
        :rtype: Timeseries
        """
        
        if new_index is None:
            new_index = self.data.index
            period = self.period
        elif not period:
            period = (new_index[1] - new_index[0]).total_seconds()

        limit = int(limit_seconds / period)  # defaults to 1 hr
        data = interpolate(
            self.data.reindex(new_index),
            replace=replace,
            method=method,
            limit=limit,
            limit_direction=limit_direction,
        )
        quality_df = self.quality_df.copy().reindex(data.index)
        quality_df[quality_df.isna()] = "i"
        # find any differences using the original data index
        mask1 = (data.reindex(self.data.index) != self.data).any(axis=1)

        # any nans from the original index
        mask2 = self.data[mask1].isna()
        quality_df[mask2] = "i"

        # any 999999s from the original index
        mask3 = self.data[mask1] == 999999
        quality_df[mask3] = "i"
        if not name:
            name = f"{self.name}.interpolated"
        if not level:
            level = self.level
        if not series:
            series = self.series
        return Timeseries(
            data=data,
            quality_df=quality_df,
            series=series,
            units=self.units,
            level=level,
            period=period,
            name=name,
        )

    def butterworth_filter(
        self,
        filter_type: str,
        filter_order: int,
        filter_cutoff_s: float,
        series: str = "",
        name: str = None,
    ):
        """Apply a butterworth filter to a DataFrame using scipy.signal.butter()

        :param filter_type: {'lowpass', 'highpass', 'bandpass', 'bandstop'}
        :type filter_type: str
        :param filter_order: the order of the filter
        :type filter_order: int
        :param filter_cutoff_s: the filter cutoff in seconds
        :type filter_cutoff_s: float
        :param series: series name, defaults to ""
        :type series: str, optional
        :param name: name for new Timeseries, defaults to None
        :type name: str, optional
        :return: Timeseries containing filtered data
        :rtype: Timeseries
        """

        df2 = butterworth_filter(
            df=self.data,
            period=self.period,
            filter_type=filter_type,
            filter_order=filter_order,
            filter_cutoff_s=filter_cutoff_s,
        )
        if not name:
            name = f"{self.name}.filtered"
        return Timeseries(
            data=df2,
            quality_df=self.quality_df,
            series=series,
            units=self.units,
            period=self.period,
            level=self.level,
            name=name,
        )

    def calculate_offsets(
        self,
        limit_multiplier: int = 10,
        cutoff_percentile: float = 0.75,
        name: str = None,
    ):
        """Calculate offsets using first differencing method (add more details).  

        :param limit_multiplier: _description_, defaults to 10
        :type limit_multiplier: int, optional
        :param cutoff_percentile: _description_, defaults to 0.75
        :type cutoff_percentile: float, optional
        :param name: name for new Timeseries, defaults to None
        :type name: str, optional
        :return: _description_
        :rtype: _type_
        """
        data = calculate_offsets(self.data, limit_multiplier, cutoff_percentile)
        if not name:
            name = f"{self.name}.offset_c"
        ts = Timeseries(
            data=data,
            quality_df=self.quality_df,
            series="offset_c",
            units=self.units,
            period=self.period,
            level="2a",
            name=name,
        )
        return ts

    def apply_calibration_matrix(
        self,
        calibration_matrix: np.array,
        calibration_matrix_name: str = None,
        use_channels: list = [1, 1, 1, 1],
        name: str = None,
    ):
        """Applies a calibration matrix to convert 4 gauges into areal, differential, and shear strains

        :param calibration_matrix: calibration matrix 
        :type calibration_matrix: np.array
        :param calibration_matrix_name:  name of calibration matrix used, defaults to None
        :type calibration_matrix_name: str, optional
        :param use_channels: not yet implemented, set to 0 to ignore a bad channel, defaults to [1, 1, 1, 1]
        :type use_channels: list, optional
        :param name: name for new Timeseries, defaults to None
        :type name: str, optional
        :return: areal, differential, and shear strains based on the given calibration
        :rtype: Timeseries
        """
        # calculate areal and shear strains from gauge strains
        # todo: implement UseChannels to arbitrary matrices
        data = apply_calibration_matrix(
            self.data, calibration_matrix, calibration_matrix_name, use_channels
        )
        quality_df = pd.DataFrame(
            index=self.quality_df.index,
            columns=data.columns,
            data="g",
        )
        quality_df[self.quality_df[self.quality_df == "m"].any(axis=1)] = "m"
        if not name:
            name = f"{self.name}.calibrated"
        ts2 = Timeseries(
            data=data,
            quality_df=quality_df,
            series=self.series,
            units="microstrain",
            level="2a",
            period=self.period,
            name=name,
        )
        return ts2

    def calculate_pressure_correction(
        self, response_coefficients: dict, name: str = None
    ):
        """Generate a pressure correction timeseries from pressure data and response coefficients

        :param response_coefficients: response coefficients for each channel loaded from metadata
        :type response_coefficients: dict
        :param name: name for new Timeseries, defaults to None
        :type name: str, optional
        :return: pressure corrections for each channel
        :rtype: Timeseries
        """
        data = calculate_pressure_correction(self.data, response_coefficients)
        quality_df = pd.DataFrame(index=data.index)
        for key in response_coefficients:
            quality_df[key] = self.quality_df
        if not name:
            name = f"{self.name}.atmp_c"
        ts2 = Timeseries(
            data=data,
            quality_df=quality_df,
            series="atmp_c",
            units="microstrain",
            level="2a",
            period=self.period,
            name=name,
        )
        return ts2

    def calculate_tide_correction(
        self, tidal_parameters: dict, longitude: float, name: str = None
    ):
        """Generate tidal correction timeseries using SPOTL hartid

        :param tidal_parameters: tidal parameters loaded from station metadata
        :type tidal_parameters: dict
        :param longitude: station longitude
        :type longitude: float
        :param name: name for new Timeseries, defaults to None
        :type name: str, optional
        :return: tidal correction Timeseries calculated for each column/channel in input data
        :rtype: Timeseries
        """
        data = calculate_tide_correction(
            self.data, self.period, tidal_parameters, longitude
        )
        if not name:
            name = f"{self.name}.tide_c"
        ts = Timeseries(
            data=data,
            series="tide_c",
            units="microstrain",
            period=self.period,
            level="2a",
            name=name,
        )
        return ts

    def linear_trend_correction(
        self, trend_start=None, trend_end=None, name: str = None
    ):
        """Generate a linear trend correction

        :param trend_start: start of window to calculate trend, defaults to first_valid_index()
        :type trend_start: datetime.datetime, optional
        :param trend_end: end of window to calculate trend, defaults to last_valid_index()
        :type trend_end: datetime.datetime, optional
        :param name: name for new Timeseries, defaults to None
        :type name: str, optional
        :return: trend correction timeseries for each column/channel in input data
        :rtype: Timeseries
        """
        data = calculate_linear_trend_correction(self.data, trend_start, trend_end)
        if not name:
            name = f"{self.name}.trend_c"
        ts = Timeseries(
            data=data,
            series="trend_c",
            units="microstrain",
            period=self.period,
            level="2a",
            name=name,
        )
        return ts

    def apply_corrections(self, corrections: list = [], name: str = None):
        """applies one or more corrections to a Timeseries and returns the corrected Timeseries

        :param corrections: List of correction Timeseries to apply, defaults to []
        :type corrections: list, optional
        :param name: name for new Timeseries, defaults to None
        :type name: str, optional
        :return: corrected Timeseries
        :rtype: Timeseries
        """
        logger.info(f"Applying corrections")
        if len(corrections):
            data = self.data.copy()
            for ts in corrections:
                data -= ts.data
        if not name:
            name = f"{self.name}.corrected"
        ts2 = Timeseries(
            data=data,
            quality_df=self.quality_df,
            series="corrected",
            units="microstrain",
            level="2a",
            period=self.period,
            name=name,
        )
        return ts2

    def dynamic_strain(
        self,
        gauge_weights: list = [1, 1, 1, 1],
        series="dynamic",
        name=None,
    ):
        """calculates dynamic strain for a given Timeseries as RMS of gauge strains

        :param gauge_weights: list of which channels to use, defaults to [1, 1, 1, 1]
        :type gauge_weights: list, optional
        :param series: series name, defaults to "dynamic"
        :type series: str, optional
        :param name: name for new Timeseries, defaults to None
        :type name: str, optional
        :return: calculated dynamic strain as a Timeseries object
        :rtype: Timeseries
        """
        df2 = dynamic_strain(self.data, gauge_weights)
        quality_df = pd.DataFrame(index=self.data.index)
        quality_df["dynamic"] = "g"
        quality_df[self.quality_df[(self.quality_df == "m")].any(axis=1)] = "m"
        quality_df[self.quality_df[(self.quality_df == "i")].any(axis=1)] = "i"
        if not name:
            name = f"{self.name}.dynamic"
        return Timeseries(
            data=df2,
            quality_df=quality_df,
            series=series,
            units=self.units,
            period=self.period,
            level=self.level,
            name=name,
        )

    def calculate_magnitude(
        self,
        hypocentral_distance,
        site_term: float = 0,
        longitude_term=0,
        name: str = None,
    ):
        """Calculates a magnitude estimate based on Barbour et al 2021.
        :param hypocentral_distance: distance from station to event hypocenter, in km
        :type hypocentral_distance: float
        :param site_term: site term from Barbour et al 2021, defaults to 0
        :type site_term: float, optional
        :param longitude_term: longitude term from Barbour et al 2021, defaults to 0
        :type longitude_term: int, optional
        :param name: name for new Timeseries, defaults to None
        :type name: str, optional
        :return: Timeseries containing a dynamic strain based magnitude estimate as a function of time
        :rtype: Timeseries
        """
        data = calculate_magnitude(
            self.data, hypocentral_distance, site_term, longitude_term
        )
        if not name:
            name = f"{self.name}.magnitude"
        return Timeseries(
            data=data,
            quality_df=self.quality_df,
            series="magnitude",
            units="magnitude",
            period=self.period,
            level=self.level,
            name=name,
        )

    def plot(
        self,
        title: str = None,
        remove_9s: bool = False,
        zero: bool = False,
        detrend: str = None,
        ymin: float = None,
        ymax: float = None,
        type: str = "line",
        show_quality_flags: bool = False,
        atmp=None,
        rainfall=None,
        save_as: str = None,
    ):
        """Generic plotting function for Timeseries data

        :param title: plot title, defaults to None
        :type title: str, optional
        :param remove_9s: option to remove gap fill values, defaults to False
        :type remove_9s: bool, optional
        :param zero: option to zero against first_valid_index(), defaults to False
        :type zero: bool, optional
        :param detrend: signal.detrend type, only 'linear' implented currently, defaults to None
        :type detrend: str, optional
        :param ymin: y-axis minimum for plot, defaults to None
        :type ymin: float, optional
        :param ymax: y-axis maximum for plot, defaults to None
        :type ymax: float, optional
        :param type: matplotlib plot type. option of ['scatter','line'], defaults to "line"
        :type type: str, optional
        :param show_quality_flags: option to highlight missing data flags, defaults to False
        :type show_quality_flags: bool, optional
        :param atmp: optional Timeseries containing atmospheric pressure data to be plotted in an extra subplot, defaults to None
        :type atmp: Timeseries, optional
        :param rainfall: optional Timeseries containing rainfall data to be plotted in an extra subplot.  will also plot cumsum of rainfall during time window. defaults to None
        :type rainfall: Timeseries, optional
        :param save_as: filename to save as, defaults to None
        :type save_as: str, optional
        """
        num_plots = len(self.columns)
        num_colors = 1
        if atmp is not None:
            num_plots += 1
        if rainfall is not None:
            num_plots += 1
            num_colors += 1
        fig, axs = plt.subplots(
            num_plots, 1, figsize=(12, 3 * num_plots), squeeze=False
        )
        colors = [cm.gnuplot(x) for x in np.linspace(0, 0.8, num_colors)]
        if title:
            fig.suptitle(title)
        else:
            fig.suptitle(self.name)
        if remove_9s:
            df = self.remove_999999s(interpolate=False).data
        else:
            df = self.data.copy()
        if zero:
            df -= df.loc[df.first_valid_index()]
        if detrend:
            if detrend == "linear":
                for ch in self.columns:
                    df[ch] = signal.detrend(df[ch], type="linear")
            else:
                logger.error("Only linear detrend implemented")
        for i, ch in enumerate(self.columns):
            if type == "line":
                axs[i][0].plot(df[ch], color=colors[0], label=ch)
            elif type == "scatter":
                axs[i][0].scatter(df.index, df[ch], color=colors[0], s=2, label=ch)
            else:
                logger.error("Plot type must be either 'line' or 'scatter'")
            if self.units:
                axs[i][0].set_ylabel(self.units)
            if ymin or ymax:
                axs[i][0].set_ylim(ymin, ymax)
            axs[i][0].ticklabel_format(axis="y", useOffset=False, style="plain")
            if show_quality_flags:
                missing = self.quality_df[self.quality_df[ch] == "m"].index
                if len(missing) and len(missing) < 100:
                    axs[i][0].text(
                        0.02, 0.95, "missing", color="r", transform=axs[i][0].transAxes
                    )
                    for time in missing:
                        axs[i][0].axvline(time, color="r")
                elif len(missing) > 100:
                    logger.info("too many missing points to plot")
                interpolated = self.quality_df[self.quality_df[ch] == "i"].index
                if len(interpolated) and len(interpolated) < 100:
                    axs[i][0].text(
                        0.08,
                        0.95,
                        "interpolated",
                        color="b",
                        transform=axs[i][0].transAxes,
                    )
                    for time in interpolated:
                        axs[i][0].axvline(time, color="b")
                elif len(interpolated) > 100:
                    logger.info("too many interpolated points to plot")
            axs[i][0].legend()
        if atmp is not None:
            i += 1
            if type == "line":
                axs[i][0].plot(atmp.data, color=colors[0], label="atmp")
            elif type == "scatter":
                axs[i][0].scatter(
                    atmp.data.index, atmp.data, color=colors[0], s=2, label="atmp"
                )
            else:
                logger.error("Plot type must be either 'line' or 'scatter'")
            if atmp.units:
                axs[i][0].set_ylabel(atmp.units)
            axs[i][0].ticklabel_format(axis="y", useOffset=False, style="plain")
            axs[i][0].legend()
        if rainfall is not None:
            i += 1
            ax2 = axs[i][0].twinx()
            axs[i][0].bar(
                rainfall.data.index,
                rainfall.data.values.reshape(len(rainfall.data)),
                width=0.02,
                color=colors[0],
                label="rainfall",
            )

            # axs[i][0].scatter(
            #     rainfall.data.index,
            #     rainfall.data,
            #     color=colors[0],
            #     s=2,
            #     label="rainfall",
            # )
            ax2.plot(
                rainfall.data.cumsum(), color=colors[1], label="cumulative rainfall"
            )

            axs[i][0].set_ylabel("mm/30m")
            ax2.set_ylabel("mm")
            axs[i][0].ticklabel_format(axis="y", useOffset=False, style="plain")
            lines, labels = axs[i][0].get_legend_handles_labels()
            lines2, labels2 = ax2.get_legend_handles_labels()
            ax2.legend(lines + lines2, labels + labels2, loc=0)

        fig.tight_layout()
        if save_as:
            logger.info(f"Saving plot to {save_as}")
            plt.savefig(save_as)


def plot_timeseries_comparison(
    timeseries: list = [],
    title: str = None,
    names: list = [],
    remove_9s: bool = False,
    zero: bool = False,
    detrend: str = None,
    type: str = "line",
    save_as: str = None,
):
    """plot multiple Timeseries in the same plot to compare values.
       useful for viewing uncorrected vs corrected data

    :param timeseries: list of Timeseries to plot, defaults to []
    :type timeseries: list, optional
    :param title: plot title, defaults to None
    :type title: str, optional
    :param names: list of names to use in legend, defaults to []
    :type names: list, optional
    :param remove_9s: option to remove gap fill values, defaults to False
    :type remove_9s: bool, optional
    :param zero: option to zero against first_valid_index(), defaults to False
    :type zero: bool, optional
    :param detrend: signal.detrend type, only 'linear' implented currently, defaults to None
    :type detrend: str, optional
    :param type: matplotlib plot type. option of ['scatter','line'], defaults to "line"
    :type type: str, optional
    :param save_as: filename to save as, defaults to None
    :type save_as: str, optional
    """
    if not isinstance(timeseries, list):
        timeseries = [timeseries]
    colors = [cm.gnuplot(x) for x in np.linspace(0, 0.8, len(timeseries))]

    fig, axs = plt.subplots(
        len(timeseries[0].columns),
        1,
        figsize=(12, 3 * len(timeseries[0].columns)),
        squeeze=False,
    )
    for j, ts in enumerate(timeseries):
        if title:
            fig.suptitle(title)
        if remove_9s:
            df = ts.remove_999999s(interpolate=False).data
        else:
            df = ts.data.copy()
        if zero:
            df -= df.loc[df.first_valid_index()]
        if detrend:
            if detrend == "linear":
                for ch in ts.columns:
                    df[ch] = signal.detrend(df[ch], type="linear")
            else:
                logger.error("Only linear detrend implemented")
        for i, ch in enumerate(ts.columns):
            if len(names) > j:
                label = f"{names[j]}"
            else:
                label = ch
            if type == "line":
                axs[i][0].plot(df[ch], color=colors[j], label=label)
            elif type == "scatter":
                axs[i][0].scatter(df.index, df[ch], color=colors[j], s=2, label=label)
            else:
                logger.error("Plot type must be either 'line' or 'scatter'")
            if ts.units:
                axs[i][0].set_ylabel(ts.units)
            axs[i][0].ticklabel_format(axis="y", useOffset=False, style="plain")
            axs[i][0].set_title(ch)
            axs[i][0].legend()

        fig.tight_layout()
        if save_as:
            logger.info(f"Saving plot to {save_as}")
            plt.savefig(save_as)

src/earthscopestraintools/timeseries.py
EarthScope/earthscopestraintools
from earthscopestraintools import ascii2tdb

if __name__ == '__main__':
    workdir = "arrays"
    network = 'PB'
    fcid = "B003"
    year = "2022"
    ascii2tdb.etl_yearly_ascii_file(network, fcid, year, workdir=workdir, print_it=True)


test/test_ascii2tdb.py
EarthScope/earthscopestraintools
from earthscopestraintools.bottle2mseed import bottle2mseed

from datetime import datetime
import logging

logger = logging.getLogger(__name__)
if logger.hasHandlers():
    logger.setLevel(logging.INFO)
else:
    logging.basicConfig(
        format="%(asctime)s - %(levelname)s - %(message)s", level=logging.INFO
    )

if __name__ == '__main__':

    t1 = datetime.now()

    
    print_traces=True
    
    # network = "GF"
    # station = "BUY1"
    # filename = 'BUY12322201.tgz' #Hour Session
    # session = "Hour"
    # bottle2mseed(network, station, filename, session, print_traces)

    # network = "PB"
    # station = "B900"

    # filename = 'B9002323517_20.tar'  
    # session = "Min"
    
    # try:
    #     bottle2mseed(network, station, filename, session, verbose=False, print_traces=print_traces, plot_traces=True)
    # except Exception as e:
    #     logger.error(e)

    network = "PB"
    station = "B072"

    filename = 'B07223196Day.tgz'  #24 hr Day session (archive and logger format)
    session = "Day"
    st = bottle2mseed(network, station, filename, session, print_traces, plot_traces=False)
    #for tr in st:
    #    print(tr.stats.channel, tr.data)
    # filename = 'B0012200100.tgz'  # 1 Hour, Hour Session (logger format)
    # session = "Hour"
    # bottle2mseed(network, station, filename, session, print_traces, plot_traces=True)

    # filename = 'B0012200100_20.tar'  # 1 Hour, Min Session (logger format)
    # session = "Min"
    # bottle2mseed(network, station, filename, session, print_traces)

    # filename = 'B001.2022001_01.tar' #24 Hour, Hour Session (archive format)
    # session = "Hour"
    # bottle2mseed(network, station, filename, session, print_traces)

    # filename = 'B001.2022001_20.tar' #24 Hour, Min session (archive format)
    # session = "Min"
    # bottle2mseed(network, station, filename, session, print_traces)

    # network = "PB"
    # station = "B012"
    # filename = 'B0122136319_20.tar' #Min Session
    # session = "Min"
    # bottle2mseed(network, station, filename, session, print_traces)

    # network = "IV"
    # station = "TSM2"
    # filename = 'TSM22130619_20.tar' #Min Session
    # session = "Min"
    # bottle2mseed(network, station, filename, session, print_traces, plot_traces=True)

    

    # t2 = datetime.now()
    # elapsed_time = t2 - t1
    # logger.info(f'{filename}: Elapsed time {elapsed_time} seconds')

test/test_bottle2mseed.py
EarthScope/earthscopestraintools
from earthscopestraintools.bottle2tdb import bottle2tdb
from earthscopestraintools.edid import get_station_edid, get_session_edid
from earthscopestraintools.tiledbtools import RawStrainReader
from datetime import datetime
import logging

logger = logging.getLogger(__name__)
if logger.hasHandlers():
    logger.setLevel(logging.INFO)
else:
    logging.basicConfig(
        format="%(asctime)s - %(levelname)s - %(message)s", level=logging.INFO
    )

if __name__ == '__main__':

    t1 = datetime.now()

    network = "PB"
    station = "B001"
    station_edid = get_station_edid(station)
    ancillary_uri = f"arrays/{station_edid}_ancillary.tdb"
    write_it = True
    print_it = True
    check_it = True


    filepath = 'bottles/B001.2022001Day.tgz'  #24 hr Day session (archive and logger format)
    session = "Day"
    session_edid = get_session_edid(station, session)
    strain_uri = f"arrays/{session_edid}.tdb"
    bottle2tdb(filepath, strain_uri, ancillary_uri, session, write_it=write_it, print_it=print_it, check_it=check_it)

    filepath = 'bottles/B0012200100.tgz'  # 1 Hour, Hour Session (logger format)
    session = "Hour"
    session_edid = get_session_edid(station, session)
    strain_uri = f"arrays/{session_edid}.tdb"
    bottle2tdb(filepath, strain_uri, ancillary_uri, session, write_it=write_it, print_it=print_it, check_it=check_it)

    filepath = 'bottles/B0012200100_20.tar'  # 1 Hour, Min Session (logger format)
    session = "Min"
    session_edid = get_session_edid(station, session)
    strain_uri = f"arrays/{session_edid}.tdb"
    bottle2tdb(filepath, strain_uri, ancillary_uri, session, write_it=write_it, print_it=print_it, check_it=check_it)

    # filepath = 'bottles/B001.2022001_01.tar' #24 Hour, Hour Session (archive format)
    # session = "Hour"
    # session_edid = get_session_edid(network, station, session)
    # strain_uri = f"arrays/{session_edid}.tdb"
    # bottle2tdb(filepath, strain_uri, ancillary_uri, session, write_it=write_it, print_it=print_it, check_it=check_it)

    # filepath = 'bottles/B001.2022001_20.tar' #24 Hour, Min session (archive format)
    # session = "Min"
    # session_edid = get_session_edid(network, station, session)
    # strain_uri = f"arrays/{session_edid}.tdb"
    # bottle2tdb(filepath, strain_uri, ancillary_uri, session, write_it=write_it, print_it=print_it, check_it=check_it)

    t2 = datetime.now()
    elapsed_time = t2 - t1
    logger.info(f'test_bottle2tdb.py: Elapsed time {elapsed_time} seconds')
test/test_bottle2tdb.py
EarthScope/earthscopestraintools
from earthscopestraintools.bottletar import GtsmBottleTar
from datetime import datetime
import logging

logger = logging.getLogger(__name__)
if logger.hasHandlers():
    logger.setLevel(logging.INFO)
else:
    logging.basicConfig(
        format="%(asctime)s - %(levelname)s - %(message)s", level=logging.INFO
    )

if __name__ == '__main__':

    t1 = datetime.now()

    filename = 'B001.2022001Day.tgz'  #24 hr Day session (archive and logger format)
    session = "Day"
    gbt = GtsmBottleTar(f"bottles/{filename}", session)
    gbt.test()

    filename = 'B0012200100.tgz'  # 1 Hour, Hour Session (logger format)
    session = "Hour"
    gbt = GtsmBottleTar(f"bottles/{filename}", session)
    gbt.test()

    filename = 'B0012200100_20.tar'  # 1 Hour, Min Session (logger format)
    session = "Min"
    gbt = GtsmBottleTar(f"bottles/{filename}", session)
    gbt.test()

    filename = 'B001.2022001_01.tar' #24 Hour, Hour Session (archive format)
    session = "Hour"
    gbt = GtsmBottleTar(f"bottles/{filename}", session)
    gbt.test()

    filename = 'B001.2022001_20.tar' #24 Hour, Min session (archive format)
    session = "Min"
    gbt = GtsmBottleTar(f"bottles/{filename}", session)
    gbt.test()

    t2 = datetime.now()
    elapsed_time = t2 - t1
    logger.info(f'{filename}: Elapsed time {elapsed_time} seconds')
test/test_bottletar.py
EarthScope/earthscopestraintools

from earthscopestraintools.edid import get_network_edid, get_station_edid, get_session_edid 
from earthscopestraintools.edid import get_network_name, get_all_stations_in_network, get_all_sessions_in_network, get_bnum_list
import logging
logger = logging.getLogger()
logging.basicConfig(
        format="%(message)s", level=logging.INFO
    )

if __name__ == '__main__':
    namespace = "FDSN"
    network = "PB"

    station_edids = get_all_stations_in_network(namespace, network)
    for station in station_edids:
        logger.info(f"{network} {station}: {station_edids[station]}")

    
    session_edids = get_all_sessions_in_network(namespace, network)
    for station in session_edids:
        logger.info(f"{network} {station}: {session_edids[station]}")
    
    station = "B001"
    network_edid = get_network_edid(namespace, network)
    logger.info(f"{network}: {network_edid}")

    station_edid = get_station_edid(station)
    logger.info(f"{station}: {station_edid}")
    session = "Day"
    session_edid = get_session_edid(station, session)
    logger.info(f"{station}.{session}: {session_edid}")
    session = "Hour"
    session_edid = get_session_edid(station, session)
    logger.info(f"{station}.{session}: {session_edid}")
    session = "Min"
    session_edid = get_session_edid(station, session)
    logger.info(f"{station}.{session}: {session_edid}")

    fdsn_network = get_network_name(station)
    logger.info(f"FDSN network name: {fdsn_network}")

    bnums = get_bnum_list(namespace="BSM", network="NOTA")
    logger.info(f"BSM:NOTA:{bnums}")
test/test_edid.py
EarthScope/earthscopestraintools
from earthscopestraintools.gtsm_metadata import GtsmMetadata

import logging
logging.basicConfig(
    level=logging.INFO,
    format="%(asctime)s %(levelname)s: %(message)s",
    datefmt="%Y-%m-%d %H:%M:%S",
)
logger = logging.getLogger(__name__)


if __name__ == '__main__':
    network = 'PB'
    station = 'B005'
    meta = GtsmMetadata(network, station)
    meta.show()
test/test_gtsm_metadata.py
EarthScope/earthscopestraintools
from earthscopestraintools.tdb2ascii import write_date_range

import logging
logging.basicConfig(
    level=logging.INFO,
    format="%(asctime)s %(levelname)s: %(message)s",
    datefmt="%Y-%m-%d %H:%M:%S",
)
logger = logging.getLogger(__name__)


if __name__ == '__main__':
    fcid = "B001"
    net = "PB"
    start_str = "2022-01-01T00:00:00"
    end_str = "2022-02-01T00:00:00"
    write_date_range(net, fcid, start_str, end_str, write_it=False, print_it=True)
    
test/test_tdb2ascii.py
EarthScope/earthscopestraintools
from earthscopestraintools.tdb2tdb import export_date_range
workdir = "arrays"

if __name__ == '__main__':
    network = "PB"
    station = "B003"
    start_str = "2022-01-01T00:00:00"
    end_str = "2022-02-01T00:00:00"
    uri = f"{workdir}/{network}_{station}_l2_etl.tdb"
    export_date_range(uri, start_str, end_str, write_it=True, print_it=True)
test/test_tdb2tdb.py
EarthScope/earthscopestraintools
from earthscopestraintools.tiledbtools import ProcessedStrainWriter, ProcessedStrainReader, RawStrainReader, RawStrainWriter
from earthscopestraintools.timeseries import Timeseries
from earthscopestraintools.edid import get_station_edid
import logging

logger = logging.getLogger(__name__)
if logger.hasHandlers():
    logger.setLevel(logging.INFO)
else:
    logging.basicConfig(
        format="%(message)s", level=logging.INFO
    )
    
if __name__ == '__main__':
    workdir = "arrays"
    network = 'PB'
    station = "B001"
    year = "2022"
    #station_edid = get_station_edid(station)
    #uri = f"{workdir}/{station_edid}_level2.tdb"
    uri = f"{workdir}/01GVDXYEEKY7T5XCCVX2BM51EM_level2.tdb"
    reader = ProcessedStrainReader(uri)
    start_str = f"{year}-01-01T00:00:00"
    end_str = f"{int(year) + 1}-01-01T00:00:00"
    ts = reader.to_ts(
        data_types=["CH0", "CH1", "CH2", "CH3"],#, "Eee+Enn", "Eee-Enn", "2Ene"],
        timeseries="microstrain",
        attrs=['data','quality','level','version'],
        start_str=start_str,
        end_str=end_str,
        units="microstrain",
        name=f"from {uri}"
    )
    ts.stats()
    ts2 = ts.remove_999999s()
    ts2.stats()
    uri2 = "arrays/test_write.tdb"
    writer = ProcessedStrainWriter(uri=uri2)
    writer.array.delete()
    writer.create()
    writer.ts_2_tiledb(ts2)

    reader2 = ProcessedStrainReader(uri2)
    start_str = f"{year}-01-01T00:00:00"
    end_str = f"{int(year) + 1}-01-01T00:00:00"
    ts3 = reader.to_ts(
        data_types=["CH0", "CH1", "CH2", "CH3"],  # , "Eee+Enn", "Eee-Enn", "2Ene"],
        timeseries="microstrain",
        attrs=['data', 'quality', 'level', 'version'],
        start_str=start_str,
        end_str=end_str,
        units="microstrain",
        name=f"from {uri2}"
    )
    ts3.stats()
test/test_ts2tdb.py
EarthScope/earthscopestraintools
from earthscopestraintools.tiledbtools import RawStrainWriter, RawStrainReader
from earthscopestraintools.mseed_tools import load_miniseed_to_df, mseed2pandas
from earthscopestraintools.gtsm_metadata import GtsmMetadata
from earthscopestraintools.timeseries import Timeseries, ts_from_mseed

import logging
logger = logging.getLogger()
logging.basicConfig(
        format="%(message)s", level=logging.INFO
    )

network = 'PB'
station = 'B944'
meta = GtsmMetadata(network, station)
start = "2021-05-01T00:00:00"
end = "2021-05-01T01:00:00"
ts = ts_from_mseed(net=network, sta=station, loc='T0', cha='LS*', start=start, end=end)
ts.stats()
ts.plot()

test_lib.py
EarthScope/earthscopestraintools