Skip to content

Tides

gravtools.tides.abstract_tide_data

Abstract base class for tide data.

Copyright (C) 2023 Andreas Hellerschmied andreas.hellerschmied@bev.gv.at

This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses/.

AbstractTideData

Bases: ABC

Abstract base class for tidal data.

Source code in gravtools/tides/abstract_tide_data.py
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
class AbstractTideData(ABC):
    """Abstract base class for tidal data."""

    def get_tidal_gravity_correction(self, epoch_dt):
        pass

    @property
    @abstractmethod
    def get_tide_corr_df(self):
        """Returns the complete tidal gravity time-series as pandas DataFrame"""
        pass

    @property
    @abstractmethod
    def filename(self) -> str:
        """Returns the name of the data file."""
        pass

    @property
    @abstractmethod
    def filetype(self) -> str:
        """Returns the type of the data file."""
        pass

    @property
    @abstractmethod
    def number_channels(self) -> str:
        """Returns the number of data channels."""
        pass

    @property
    @abstractmethod
    def number_records(self) -> str:
        """Returns the number of data records."""
        pass

    @property
    @abstractmethod
    def starttime(self) -> dt.datetime:
        """Returns the start time of the tidal timeseries."""
        pass

    @property
    def starttime_str(self) -> str:
        """Returns the start time of the tidal timeseries as string."""
        return self.starttime.strftime('%Y-%d-%d %H:%M:%S')

    @property
    @abstractmethod
    def endtime(self) -> dt.datetime:
        """Returns the end time of the tidal timeseries."""
        pass

    @abstractmethod
    def get_channel_df(self, channel: [int, str]) -> pd.DataFrame:
        """Return the specified channel and the time reference (DateTime) as pandas DataFrame.

        Parameters
        ----------
        channel : str or int or 'last'
            If `channel` is as string the argument is interpreted as channel name and the according channel data is
            returned. If `channel` is an integer, it is interpreted as the number of the channel to be returned. `last`
            indicates that the last channel in the TSF file should be returned.


        Returns
        -------
        pandas.DataFrame containing the data records of the selected channel and the time reference as DateTime.
        """
        pass

    @property
    def endtime_str(self) -> str:
        """Returns the end time of the tidal timeseries as string."""
        return self.endtime.strftime('%Y-%d-%d %H:%M:%S')

    @property
    def filename_without_path(self) -> str:
        """Returns the filename without path."""
        pass

    def __str__(self):
        return f'Tidal data loaded from {self.filename} ({self.filetype}) with {self.number_channels} channels ' \
               f'and {self.number_records} datasets ({self.starttime_str} to {self.endtime_str})'

endtime abstractmethod property

Returns the end time of the tidal timeseries.

endtime_str property

Returns the end time of the tidal timeseries as string.

filename abstractmethod property

Returns the name of the data file.

filename_without_path property

Returns the filename without path.

filetype abstractmethod property

Returns the type of the data file.

get_tide_corr_df abstractmethod property

Returns the complete tidal gravity time-series as pandas DataFrame

number_channels abstractmethod property

Returns the number of data channels.

number_records abstractmethod property

Returns the number of data records.

starttime abstractmethod property

Returns the start time of the tidal timeseries.

starttime_str property

Returns the start time of the tidal timeseries as string.

get_channel_df(channel) abstractmethod

Return the specified channel and the time reference (DateTime) as pandas DataFrame.

Parameters:

Name Type Description Default
channel str or int or last

If channel is as string the argument is interpreted as channel name and the according channel data is returned. If channel is an integer, it is interpreted as the number of the channel to be returned. last indicates that the last channel in the TSF file should be returned.

required

Returns:

Type Description
pandas.DataFrame containing the data records of the selected channel and the time reference as DateTime.
Source code in gravtools/tides/abstract_tide_data.py
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
@abstractmethod
def get_channel_df(self, channel: [int, str]) -> pd.DataFrame:
    """Return the specified channel and the time reference (DateTime) as pandas DataFrame.

    Parameters
    ----------
    channel : str or int or 'last'
        If `channel` is as string the argument is interpreted as channel name and the according channel data is
        returned. If `channel` is an integer, it is interpreted as the number of the channel to be returned. `last`
        indicates that the last channel in the TSF file should be returned.


    Returns
    -------
    pandas.DataFrame containing the data records of the selected channel and the time reference as DateTime.
    """
    pass

gravtools.tides.tide_data_tfs

Class for TSF-formatted tidal data.

TSF Files are created e.g. by TSoft

Copyright (C) 2023 Andreas Hellerschmied andreas.hellerschmied@bev.gv.at

This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses/.

ChannelMetadata dataclass

Class for channel metadata in TSF files.

Source code in gravtools/tides/tide_data_tfs.py
33
34
35
36
37
38
39
40
41
42
43
44
45
@dataclass
class ChannelMetadata:
    """Class for channel metadata in TSF files."""
    location: str
    instrument_name: str
    data_type: str
    unit: str
    channel_number: int

    @property
    def channel_name(self):
        """Returns the channels name as read from the TSF file"""
        return self.location + ':' + self.instrument_name + ':' + self.data_type

channel_name property

Returns the channels name as read from the TSF file

TSF

Bases: AbstractTideData

Class for tidal data loaded from TSF files.

TSF files are the native format for the TSoft which allows to calculate time-series of synthetic tide for a specific station defined by geographical coordinates and height.

TSoft: Download TSoft and its manual: http://seismologie.oma.be/en/downloads/tsoft Reference: Michel Van Camp, Paul Vauterin: Tsoft: graphical and interactive software for the analysis of time series and Earth tides. Comput. Geosci. 31(5): 631-640 (2005) DOI: https://doi.org/10.1016/j.cageo.2004.11.015

Attributes:

Name Type Description
_filename str

Path and name of the TSF file.

tfs_format str

Format version of the TSF file (from block: [TSF-file]). Only version v01.0 is supported!

timeformat str

Format of the time reference (from block: [TIMEFORMAT]). Has to be "DATETIME" or "DATETIMEFRAC"!

undetval float

Determines which number is used to denote undetermined or unknown values in the [DATA] block (from block: [UNDETVAL]).

increment float

Determines the time delay (in seconds) between two consecutive data points (from block: [INCREMENT]).

channel_metadata list of ChannelMetadata objects

List of ChannelMetadata objects with one item per channel. Contains metadata such as location, instrument name, data type und unit.

data_df DataFrame

Contains all columns from the [DATA] block. Additionally, the column epoch_dt describes the reference time and date in as DateTime object.

countinfo (int, optiona(default=-1))

Specifies the number of data points in the file (from block: [COUNTINFO]). If this block is missing the in the TSF file the default value -1 is set.

comment (str, optional(default=''))

Optional multiline comment (from block: [COMMENT]).

Source code in gravtools/tides/tide_data_tfs.py
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
class TSF(AbstractTideData):
    """Class for tidal data loaded from TSF files.

    TSF files are the native format for the TSoft which allows to calculate time-series of synthetic tide for a
    specific station defined by geographical coordinates and height.

    TSoft:
    Download TSoft and its manual: http://seismologie.oma.be/en/downloads/tsoft
    Reference: Michel Van Camp, Paul Vauterin: Tsoft: graphical and interactive software for the analysis of time series
    and Earth tides. Comput. Geosci. 31(5): 631-640 (2005)
    DOI: https://doi.org/10.1016/j.cageo.2004.11.015

    Attributes
    ----------
    _filename : str
        Path and name of the TSF file.
    tfs_format : str
        Format version of the TSF file (from block: [TSF-file]). Only version v01.0 is supported!
    timeformat : str
        Format of the time reference (from block: [TIMEFORMAT]). Has to be "DATETIME" or "DATETIMEFRAC"!
    undetval : float
        Determines which number is used to denote undetermined or unknown values in the [DATA] block (from block:
        [UNDETVAL]).
    increment : float
        Determines the time delay (in seconds) between two consecutive data points (from block: [INCREMENT]).
    channel_metadata : list of ChannelMetadata objects
        List of ChannelMetadata objects with one item per channel. Contains metadata such as location, instrument name,
        data type und unit.
    data_df : pandas.DataFrame
        Contains all columns from the [DATA] block. Additionally, the column `epoch_dt` describes the reference time and
        date in as DateTime object.
    countinfo : int, optiona (default=-1)
        Specifies the number of data points in the file (from block: [COUNTINFO]). If this block is missing the in the
        TSF file the default value -1 is set.
    comment : str, optional (default='')
        Optional multiline comment (from block: [COMMENT]).
    """
    def __init__(self, filename: str, tfs_format: str, timeformat: str, undetval: float, increment: float,
                 channel_metadata: list, data_df: pd.DataFrame, countinfo: int = -1, comment: str = ''):
        """Default constructor."""
        self._filename = filename
        self.tfs_format = tfs_format
        self.timeformat = timeformat
        self.undetval = undetval
        self.increment = increment
        self.countinfo = countinfo
        self.comment = comment
        self.channel_metadata = channel_metadata
        self.data_df = data_df


    @classmethod
    def from_tfs_file(cls, filename):
        """ Create class instance from TSF file.

        Parameters
        ----------
        filename : str
            Path and name of TSF file.

        Returns
        -------
        Class instance
        """
        with open(filename, 'r') as f:
            tsf_str = f.read()

            # [TSF-file]: single attribute, single occurrence
            expr = r'\[TSF-file\]\s*(?P<tsf_format>\S+)\s*\n'
            expr_count = 0
            for expr_result in re.finditer(expr, tsf_str):
                expr_result_dict = expr_result.groupdict()
                expr_count += 1
            if expr_count == 1:  # OK
                tfs_format = expr_result_dict['tsf_format']
                if tfs_format != 'v01.0':
                    raise RuntimeError(f'Invalid TSF file format: {tfs_format} (valid: v01.0)')
            else:
                raise RuntimeError(f'TSF-file format tag is missing at the begin of the file, or it occurs more than '
                                   f'once!')

            # [TIMEFORMAT]: single attribute, single occurrence
            expr = r'\[TIMEFORMAT\]\s*(?P<timeformat>\S+)\s*\n'
            expr_count = 0
            for expr_result in re.finditer(expr, tsf_str):
                expr_result_dict = expr_result.groupdict()
                expr_count += 1
            if expr_count == 1:  # OK
                timeformat = expr_result_dict['timeformat']
            elif expr_count == 0:
                raise RuntimeError(f'TSF block [TIMEFORMAT] not found in input TSF file!')
            else:
                raise RuntimeError(f'TSF block [TIMEFORMAT] occurs {expr_count} time! Only once allowed.')
            if timeformat not in ('DATETIME', 'DATETIMEFRAC'):
                raise RuntimeError(f'Invalif [TIMEFORMAT] in TSF File: {timeformat}')

            # [UNDETVAL]: single attribute, single occurrence
            expr = r'\[UNDETVAL\]\s*(?P<undetval>\S+)\s*\n'
            expr_count = 0
            for expr_result in re.finditer(expr, tsf_str):
                expr_result_dict = expr_result.groupdict()
                expr_count += 1
            if expr_count == 1:  # OK
                undetval = expr_result_dict['undetval']
            elif expr_count == 0:
                raise RuntimeError(f'TSF block [UNDETVAL] not found in input TSF file!')
            else:
                raise RuntimeError(f'TSF block [UNDETVAL] occurs {expr_count} time! Only once allowed.')
            undetval = float(undetval)

            # [INCREMENT]: single attribute, single occurrence
            expr = r'\[INCREMENT\]\s*(?P<increment>\S+)\s*\n'
            expr_count = 0
            for expr_result in re.finditer(expr, tsf_str):
                expr_result_dict = expr_result.groupdict()
                expr_count += 1
            if expr_count == 1:  # OK
                increment = expr_result_dict['increment']
            elif expr_count == 0:
                raise RuntimeError(f'TSF block [INCREMENT] not found in input TSF file!')
            else:
                raise RuntimeError(f'TSF block [INCREMENT] occurs {expr_count} time! Only once allowed.')
            increment = float(increment)

            # [COUNTINFO]: optional, single attribute, single occurrence
            expr = r'\[COUNTINFO\]\s*(?P<countinfo>\S+)\s*\n'
            expr_count = 0
            for expr_result in re.finditer(expr, tsf_str):
                expr_result_dict = expr_result.groupdict()
                expr_count += 1
            if expr_count == 1:  # OK
                countinfo = int(expr_result_dict['countinfo'])
            else:
                countinfo = -1  # Error code, if the block wasn't found!

            # [COMMENT]: optional, single attribute, single occurrence
            expr = r'\[COMMENT\]\n(?P<comment>[\s\S]*?(?=\n+\[))'
            expr_count = 0
            for expr_result in re.finditer(expr, tsf_str):
                expr_result_dict = expr_result.groupdict()
                expr_count += 1
            if expr_count == 1:  # OK
                comment = expr_result_dict['comment']
                comment = comment.strip()
            else:
                comment = ''

            # [CHANNELS]: multiple lines, single occurrence
            expr = r'\[CHANNELS\]\n(?P<channels>(?: *.+[\n\r]{0,1})+)\s*'
            expr_count = 0
            for expr_result in re.finditer(expr, tsf_str):
                expr_result_dict = expr_result.groupdict()
                expr_count += 1
            if expr_count == 1:  # OK
                channels = expr_result_dict['channels']
            elif expr_count == 0:
                raise RuntimeError(f'TSF block [CHANNELS] not found in input TSF file!')
            else:
                raise RuntimeError(f'TSF block [CHANNELS] occurs {expr_count} time! Only once allowed.')
            channels = channels.splitlines()
            channels = [s.strip() for s in channels]

            # [UNITS]: multiple lines, single occurrence
            expr = r'\[UNITS\]\n(?P<units>(?: *.+[\n\r]{0,1})+)\s*'
            expr_count = 0
            for expr_result in re.finditer(expr, tsf_str):
                expr_result_dict = expr_result.groupdict()
                expr_count += 1
            if expr_count == 1:  # OK
                units = expr_result_dict['units']
            elif expr_count == 0:
                raise RuntimeError(f'TSF block [UNITS] not found in input TSF file!')
            else:
                raise RuntimeError(f'TSF block [UNITS] occurs {expr_count} time! Only once allowed.')
            units = units.splitlines()
            units = [s.strip() for s in units]

            # Create channel metadata:
            channel_metadata = []
            for count, channel in enumerate(channels):
                channel_metadata.append(ChannelMetadata(location=channel.split(':')[0],
                                instrument_name=channel.split(':')[1],
                                data_type=channel.split(':')[2],
                                unit=units[count],
                                channel_number=count+1))

            # [DATA] multiple blocks possible, multiple lines per block
            data_df_list = []
            if timeformat == 'DATETIME':  # Multiple of seconds
                expr = r'\[DATA\]\n(?P<obs_data>(?:[0-9]{4} [0-9]{2} [0-9]{2}  [0-9]{2} [0-9]{2} [0-9]{2}.+[\n\r]{0,1})+)\s*'
                for expr_result in re.finditer(expr, tsf_str):
                    expr_result_dict = expr_result.groupdict()
                    obs_data_str = expr_result_dict['obs_data']
                    data_df_list.append(pd.read_csv(io.StringIO(obs_data_str), delim_whitespace=True, header=None))
                    columns = ['year', 'month', 'day', 'hour', 'minute', 'second'] + [ch.channel_name for ch in
                                                                                channel_metadata]
                    num_channels_in_data_block = data_df_list[0].shape[1] - 6

            elif timeformat == 'DATETIMEFRAC':
                expr = r'\[DATA\]\n(?P<obs_data>(?:[0-9]{4} [0-9]{2} [0-9]{2}  [0-9]{2} [0-9]{2} [0-9]{2} [0-9]{3}.+[\n\r]{0,1})+)\s*'
                for expr_result in re.finditer(expr, tsf_str):
                    expr_result_dict = expr_result.groupdict()
                    obs_data_str = expr_result_dict['obs_data']
                    data_df_list.append(pd.read_csv(io.StringIO(obs_data_str), delim_whitespace=True, header=None))
                    columns = ['year', 'month', 'day', 'hour', 'minute', 'second','ms'] + [ch.channel_name for ch in
                                                                                           channel_metadata]
                    num_channels_in_data_block = data_df_list[0].shape[1] - 7

            # Concat DATA blocks and add column names:
            if len(data_df_list) == 0:
                raise RuntimeError(f'No [DATA] block found in TSF file!')
            elif len(data_df_list) > 0:
                data_df = pd.concat(data_df_list, ignore_index=True)
                data_df.columns = columns

            # Check, if metadata is available for each channel:
            if num_channels_in_data_block != len(channel_metadata):
                raise RuntimeError(f'Number of data channels in [DATA] block ({num_channels_in_data_block}) not equal '
                                   f'to number of channel metadata '
                                   f'sets in the [CHANNELS] and [UNITS] blocks ({len(channel_metadata)})!')

            # Create datetime column:
            if timeformat == 'DATETIME':
                data_df['epoch_dt'] = pd.to_datetime(data_df[['year', 'month', 'day', 'hour', 'minute', 'second']])
            elif timeformat == 'DATETIMEFRAC':
                data_df['epoch_dt'] = pd.to_datetime(data_df[['year', 'month', 'day', 'hour', 'minute', 'second',
                                                              'ms']])

        return cls(filename=filename, tfs_format=tfs_format, timeformat=timeformat, undetval=undetval,
                   increment=increment, countinfo=countinfo, comment=comment, channel_metadata=channel_metadata,
                   data_df=data_df)

        # super().__init__(filename, filetype='tfs')

    @property
    def get_tide_corr_df(self):
        """Returns the complete tidal gravity time-series as pandas DataFrame"""
        return self.data_df


    @property
    def filename(self) -> str:
        """Returns the name of the data file."""
        return self._filename

    @property
    def filetype(self) -> str:
        """Returns the type of the data file."""
        return 'TSoft TSF file'

    @property
    def number_channels(self) -> str:
        """Returns the number of data channels."""
        return len(self.channel_metadata)

    @property
    def number_records(self) -> str:
        """Returns the number of data records."""
        return len(self.data_df)

    @property
    def starttime(self) -> dt.datetime:
        """Returns the start time of the tidal timeseries."""
        return self.data_df['epoch_dt'].min()

    @property
    def endtime(self) -> dt.datetime:
        """Returns the end time of the tidal timeseries."""
        return self.data_df['epoch_dt'].max()

    def get_channel_df(self, channel: [int, str]):
        """Return the specified channel and the time reference (DateTime) as pandas DataFrame.

        Parameters
        ----------
        channel : str or int or 'last'
            If `channel` is as string the argument is interpreted as channel name and the according channel data is
            returned. If `channel` is an integer, it is interpreted as the number of the channel to be returned. `last`
            indicates that the last channel in the TSF file should be returned.

        Returns
        -------
        pandas.DataFrame containing the data records of the selected channel and the time reference as DateTime.
        """
        if isinstance(channel, str):
            if channel == 'last':
                channel = max([num.channel_number for num in self.channel_metadata])
            elif channel in self.data_df.columns and channel in self.channel_names:
                return self.data_df.loc[:, ['epoch_dt', channel]]
        if isinstance(channel, int):
            for ch in self.channel_metadata:
                if ch.channel_number == channel:
                    return self.data_df.loc[:, ['epoch_dt', ch.channel_name]]
        raise RuntimeError(f'"{channel}" is not a valid TSF data record channel name or number!')

    def get_channel_np(self, channel: [int, str]):
        """Return the time reference and the data of the specified channel as numpy array.

        Parameters
        ----------
        channel : str or int or 'last'
            If `channel` is as string the argument is interpreted as channel name and the according channel data is
            returned. If `channel` is an integer, it is interpreted as the number of the channel to be returned. `last`
            indicates that the last channel in the TSF file should be returned.

        Returns
        -------
        numpy.array of reference epochs (DateTime), numpy.array of channel data.
        """
        if isinstance(channel, str):
            if channel == 'last':
                channel = max([num.channel_number for num in self.channel_metadata])
            elif channel in self.data_df.columns and channel in self.channel_names:
                return self.data_df.loc[:, 'epoch_dt'].to_numpy(), self.data_df.loc[:, channel].to_numpy()
        if isinstance(channel, int):
            for ch in self.channel_metadata:
                if ch.channel_number == channel:
                    return self.data_df.loc[:, 'epoch_dt'].to_numpy(), self.data_df.loc[:, ch.channel_name].to_numpy()
        raise RuntimeError(f'"{channel}" is not a valid TSF data record channel name or number!')


    @property
    def channel_names(self):
        """Returns a list of the names of all data channels."""
        return [ch.channel_name for ch in self.channel_metadata]

    @property
    def channel_names(self):
        """Returns a list of the names of all data channels."""
        return [ch.channel_name for ch in self.channel_metadata]

    @property
    def locations(self):
        """Returns a list of the locations of all stations."""
        return [ch.location for ch in self.channel_metadata]

    @property
    def instruments(self):
        """Returns a list of the instruments of all stations."""
        return [ch.instrument_name for ch in self.channel_metadata]

    @property
    def data_types(self):
        """Returns a list of the data types of all stations."""
        return [ch.data_type for ch in self.channel_metadata]

    @property
    def units(self):
        """Returns a list of the units of all stations."""
        return [ch.unit for ch in self.channel_metadata]

    @property
    def filename_without_path(self):
        """Returns the TSF filename without path."""
        return os.path.basename(self.filename)

channel_names property

Returns a list of the names of all data channels.

data_types property

Returns a list of the data types of all stations.

endtime property

Returns the end time of the tidal timeseries.

filename property

Returns the name of the data file.

filename_without_path property

Returns the TSF filename without path.

filetype property

Returns the type of the data file.

get_tide_corr_df property

Returns the complete tidal gravity time-series as pandas DataFrame

instruments property

Returns a list of the instruments of all stations.

locations property

Returns a list of the locations of all stations.

number_channels property

Returns the number of data channels.

number_records property

Returns the number of data records.

starttime property

Returns the start time of the tidal timeseries.

units property

Returns a list of the units of all stations.

__init__(filename, tfs_format, timeformat, undetval, increment, channel_metadata, data_df, countinfo=-1, comment='')

Default constructor.

Source code in gravtools/tides/tide_data_tfs.py
85
86
87
88
89
90
91
92
93
94
95
96
def __init__(self, filename: str, tfs_format: str, timeformat: str, undetval: float, increment: float,
             channel_metadata: list, data_df: pd.DataFrame, countinfo: int = -1, comment: str = ''):
    """Default constructor."""
    self._filename = filename
    self.tfs_format = tfs_format
    self.timeformat = timeformat
    self.undetval = undetval
    self.increment = increment
    self.countinfo = countinfo
    self.comment = comment
    self.channel_metadata = channel_metadata
    self.data_df = data_df

from_tfs_file(filename) classmethod

Create class instance from TSF file.

Parameters:

Name Type Description Default
filename str

Path and name of TSF file.

required

Returns:

Type Description
Class instance
Source code in gravtools/tides/tide_data_tfs.py
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
@classmethod
def from_tfs_file(cls, filename):
    """ Create class instance from TSF file.

    Parameters
    ----------
    filename : str
        Path and name of TSF file.

    Returns
    -------
    Class instance
    """
    with open(filename, 'r') as f:
        tsf_str = f.read()

        # [TSF-file]: single attribute, single occurrence
        expr = r'\[TSF-file\]\s*(?P<tsf_format>\S+)\s*\n'
        expr_count = 0
        for expr_result in re.finditer(expr, tsf_str):
            expr_result_dict = expr_result.groupdict()
            expr_count += 1
        if expr_count == 1:  # OK
            tfs_format = expr_result_dict['tsf_format']
            if tfs_format != 'v01.0':
                raise RuntimeError(f'Invalid TSF file format: {tfs_format} (valid: v01.0)')
        else:
            raise RuntimeError(f'TSF-file format tag is missing at the begin of the file, or it occurs more than '
                               f'once!')

        # [TIMEFORMAT]: single attribute, single occurrence
        expr = r'\[TIMEFORMAT\]\s*(?P<timeformat>\S+)\s*\n'
        expr_count = 0
        for expr_result in re.finditer(expr, tsf_str):
            expr_result_dict = expr_result.groupdict()
            expr_count += 1
        if expr_count == 1:  # OK
            timeformat = expr_result_dict['timeformat']
        elif expr_count == 0:
            raise RuntimeError(f'TSF block [TIMEFORMAT] not found in input TSF file!')
        else:
            raise RuntimeError(f'TSF block [TIMEFORMAT] occurs {expr_count} time! Only once allowed.')
        if timeformat not in ('DATETIME', 'DATETIMEFRAC'):
            raise RuntimeError(f'Invalif [TIMEFORMAT] in TSF File: {timeformat}')

        # [UNDETVAL]: single attribute, single occurrence
        expr = r'\[UNDETVAL\]\s*(?P<undetval>\S+)\s*\n'
        expr_count = 0
        for expr_result in re.finditer(expr, tsf_str):
            expr_result_dict = expr_result.groupdict()
            expr_count += 1
        if expr_count == 1:  # OK
            undetval = expr_result_dict['undetval']
        elif expr_count == 0:
            raise RuntimeError(f'TSF block [UNDETVAL] not found in input TSF file!')
        else:
            raise RuntimeError(f'TSF block [UNDETVAL] occurs {expr_count} time! Only once allowed.')
        undetval = float(undetval)

        # [INCREMENT]: single attribute, single occurrence
        expr = r'\[INCREMENT\]\s*(?P<increment>\S+)\s*\n'
        expr_count = 0
        for expr_result in re.finditer(expr, tsf_str):
            expr_result_dict = expr_result.groupdict()
            expr_count += 1
        if expr_count == 1:  # OK
            increment = expr_result_dict['increment']
        elif expr_count == 0:
            raise RuntimeError(f'TSF block [INCREMENT] not found in input TSF file!')
        else:
            raise RuntimeError(f'TSF block [INCREMENT] occurs {expr_count} time! Only once allowed.')
        increment = float(increment)

        # [COUNTINFO]: optional, single attribute, single occurrence
        expr = r'\[COUNTINFO\]\s*(?P<countinfo>\S+)\s*\n'
        expr_count = 0
        for expr_result in re.finditer(expr, tsf_str):
            expr_result_dict = expr_result.groupdict()
            expr_count += 1
        if expr_count == 1:  # OK
            countinfo = int(expr_result_dict['countinfo'])
        else:
            countinfo = -1  # Error code, if the block wasn't found!

        # [COMMENT]: optional, single attribute, single occurrence
        expr = r'\[COMMENT\]\n(?P<comment>[\s\S]*?(?=\n+\[))'
        expr_count = 0
        for expr_result in re.finditer(expr, tsf_str):
            expr_result_dict = expr_result.groupdict()
            expr_count += 1
        if expr_count == 1:  # OK
            comment = expr_result_dict['comment']
            comment = comment.strip()
        else:
            comment = ''

        # [CHANNELS]: multiple lines, single occurrence
        expr = r'\[CHANNELS\]\n(?P<channels>(?: *.+[\n\r]{0,1})+)\s*'
        expr_count = 0
        for expr_result in re.finditer(expr, tsf_str):
            expr_result_dict = expr_result.groupdict()
            expr_count += 1
        if expr_count == 1:  # OK
            channels = expr_result_dict['channels']
        elif expr_count == 0:
            raise RuntimeError(f'TSF block [CHANNELS] not found in input TSF file!')
        else:
            raise RuntimeError(f'TSF block [CHANNELS] occurs {expr_count} time! Only once allowed.')
        channels = channels.splitlines()
        channels = [s.strip() for s in channels]

        # [UNITS]: multiple lines, single occurrence
        expr = r'\[UNITS\]\n(?P<units>(?: *.+[\n\r]{0,1})+)\s*'
        expr_count = 0
        for expr_result in re.finditer(expr, tsf_str):
            expr_result_dict = expr_result.groupdict()
            expr_count += 1
        if expr_count == 1:  # OK
            units = expr_result_dict['units']
        elif expr_count == 0:
            raise RuntimeError(f'TSF block [UNITS] not found in input TSF file!')
        else:
            raise RuntimeError(f'TSF block [UNITS] occurs {expr_count} time! Only once allowed.')
        units = units.splitlines()
        units = [s.strip() for s in units]

        # Create channel metadata:
        channel_metadata = []
        for count, channel in enumerate(channels):
            channel_metadata.append(ChannelMetadata(location=channel.split(':')[0],
                            instrument_name=channel.split(':')[1],
                            data_type=channel.split(':')[2],
                            unit=units[count],
                            channel_number=count+1))

        # [DATA] multiple blocks possible, multiple lines per block
        data_df_list = []
        if timeformat == 'DATETIME':  # Multiple of seconds
            expr = r'\[DATA\]\n(?P<obs_data>(?:[0-9]{4} [0-9]{2} [0-9]{2}  [0-9]{2} [0-9]{2} [0-9]{2}.+[\n\r]{0,1})+)\s*'
            for expr_result in re.finditer(expr, tsf_str):
                expr_result_dict = expr_result.groupdict()
                obs_data_str = expr_result_dict['obs_data']
                data_df_list.append(pd.read_csv(io.StringIO(obs_data_str), delim_whitespace=True, header=None))
                columns = ['year', 'month', 'day', 'hour', 'minute', 'second'] + [ch.channel_name for ch in
                                                                            channel_metadata]
                num_channels_in_data_block = data_df_list[0].shape[1] - 6

        elif timeformat == 'DATETIMEFRAC':
            expr = r'\[DATA\]\n(?P<obs_data>(?:[0-9]{4} [0-9]{2} [0-9]{2}  [0-9]{2} [0-9]{2} [0-9]{2} [0-9]{3}.+[\n\r]{0,1})+)\s*'
            for expr_result in re.finditer(expr, tsf_str):
                expr_result_dict = expr_result.groupdict()
                obs_data_str = expr_result_dict['obs_data']
                data_df_list.append(pd.read_csv(io.StringIO(obs_data_str), delim_whitespace=True, header=None))
                columns = ['year', 'month', 'day', 'hour', 'minute', 'second','ms'] + [ch.channel_name for ch in
                                                                                       channel_metadata]
                num_channels_in_data_block = data_df_list[0].shape[1] - 7

        # Concat DATA blocks and add column names:
        if len(data_df_list) == 0:
            raise RuntimeError(f'No [DATA] block found in TSF file!')
        elif len(data_df_list) > 0:
            data_df = pd.concat(data_df_list, ignore_index=True)
            data_df.columns = columns

        # Check, if metadata is available for each channel:
        if num_channels_in_data_block != len(channel_metadata):
            raise RuntimeError(f'Number of data channels in [DATA] block ({num_channels_in_data_block}) not equal '
                               f'to number of channel metadata '
                               f'sets in the [CHANNELS] and [UNITS] blocks ({len(channel_metadata)})!')

        # Create datetime column:
        if timeformat == 'DATETIME':
            data_df['epoch_dt'] = pd.to_datetime(data_df[['year', 'month', 'day', 'hour', 'minute', 'second']])
        elif timeformat == 'DATETIMEFRAC':
            data_df['epoch_dt'] = pd.to_datetime(data_df[['year', 'month', 'day', 'hour', 'minute', 'second',
                                                          'ms']])

    return cls(filename=filename, tfs_format=tfs_format, timeformat=timeformat, undetval=undetval,
               increment=increment, countinfo=countinfo, comment=comment, channel_metadata=channel_metadata,
               data_df=data_df)

get_channel_df(channel)

Return the specified channel and the time reference (DateTime) as pandas DataFrame.

Parameters:

Name Type Description Default
channel str or int or last

If channel is as string the argument is interpreted as channel name and the according channel data is returned. If channel is an integer, it is interpreted as the number of the channel to be returned. last indicates that the last channel in the TSF file should be returned.

required

Returns:

Type Description
pandas.DataFrame containing the data records of the selected channel and the time reference as DateTime.
Source code in gravtools/tides/tide_data_tfs.py
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
def get_channel_df(self, channel: [int, str]):
    """Return the specified channel and the time reference (DateTime) as pandas DataFrame.

    Parameters
    ----------
    channel : str or int or 'last'
        If `channel` is as string the argument is interpreted as channel name and the according channel data is
        returned. If `channel` is an integer, it is interpreted as the number of the channel to be returned. `last`
        indicates that the last channel in the TSF file should be returned.

    Returns
    -------
    pandas.DataFrame containing the data records of the selected channel and the time reference as DateTime.
    """
    if isinstance(channel, str):
        if channel == 'last':
            channel = max([num.channel_number for num in self.channel_metadata])
        elif channel in self.data_df.columns and channel in self.channel_names:
            return self.data_df.loc[:, ['epoch_dt', channel]]
    if isinstance(channel, int):
        for ch in self.channel_metadata:
            if ch.channel_number == channel:
                return self.data_df.loc[:, ['epoch_dt', ch.channel_name]]
    raise RuntimeError(f'"{channel}" is not a valid TSF data record channel name or number!')

get_channel_np(channel)

Return the time reference and the data of the specified channel as numpy array.

Parameters:

Name Type Description Default
channel str or int or last

If channel is as string the argument is interpreted as channel name and the according channel data is returned. If channel is an integer, it is interpreted as the number of the channel to be returned. last indicates that the last channel in the TSF file should be returned.

required

Returns:

Type Description
numpy.array of reference epochs (DateTime), numpy.array of channel data.
Source code in gravtools/tides/tide_data_tfs.py
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
def get_channel_np(self, channel: [int, str]):
    """Return the time reference and the data of the specified channel as numpy array.

    Parameters
    ----------
    channel : str or int or 'last'
        If `channel` is as string the argument is interpreted as channel name and the according channel data is
        returned. If `channel` is an integer, it is interpreted as the number of the channel to be returned. `last`
        indicates that the last channel in the TSF file should be returned.

    Returns
    -------
    numpy.array of reference epochs (DateTime), numpy.array of channel data.
    """
    if isinstance(channel, str):
        if channel == 'last':
            channel = max([num.channel_number for num in self.channel_metadata])
        elif channel in self.data_df.columns and channel in self.channel_names:
            return self.data_df.loc[:, 'epoch_dt'].to_numpy(), self.data_df.loc[:, channel].to_numpy()
    if isinstance(channel, int):
        for ch in self.channel_metadata:
            if ch.channel_number == channel:
                return self.data_df.loc[:, 'epoch_dt'].to_numpy(), self.data_df.loc[:, ch.channel_name].to_numpy()
    raise RuntimeError(f'"{channel}" is not a valid TSF data record channel name or number!')

gravtools.tides.correction_time_series

Class for correction data provided as time series per station and survey.

Copyright (C) 2023 Andreas Hellerschmied andreas.hellerschmied@bev.gv.at

This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses/.

CorrectionTimeSeries

Class for correction data provided as time series per station and survey.

Attributes:

Name Type Description
surveys dict of `gravtools.tides.correction_time_series.SurveyCorrections`

Dict of SurveyCorrection objects with ths survey names as keys.

Source code in gravtools/tides/correction_time_series.py
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
class CorrectionTimeSeries:
    """Class for correction data provided as time series per station and survey.

    Attributes
    ----------
    surveys: dict of `gravtools.tides.correction_time_series.SurveyCorrections`
        Dict of `SurveyCorrection` objects with ths survey names as keys.
    """

    def __init__(self):
        """Default initializer."""
        self.surveys = {}

    def load_tfs_file(self, survey_name: str, filename_tsf: str, location: str ='', instrument: str ='',
                      data_type: str ='', overwrite_channel: bool=True, is_correction: bool='False'):
        """Load time series data from a TSoft TSF file for one survey.

        Parameters
        ----------
        survey_name: str
            Name of the survey for which the data is loaded.
        filename_tsf: str
            Name and path of the TSF file from which the correction time series data is loaded.
        location: str, optional (default='')
            If not empty, only channels with matching locations are loaded.
        instrument: str, optional (default='')
            If not empty, only channels with matching instruments are loaded.
        data_type: str, optional (default='')
            If not empty, only channels with matching data types are loaded.
        overwrite_channel: bool, optional (default=`True`)
            `True` implies that an existing StationCorrection object of a station will be overwritten.
        is_correction: bool, optional (default=`False`)
            `True` implies that the channel data loaded from the TSF file model the gravity effect of phenomena, rather
            than corrections for gravity observations. Both options have the opposite sign: while corrections have to be
            added to gravity observations, effects have to be subtracted, in order to reduce observations.
            See TSoft manual (version 2.2.4 Release date 2015-09-09), p. 15: "The loading calculation by Tsoft will
            result in the correction, not the effect.
            On the other hand the prediction of the solid Earth tides using the WDD parameter set
            provided by Tsoft directly will result in the effect. Therefore both must be treated with
            different sign in order to reduce a gravity time series correctly."
        """
        # Load TSF file:
        tsf_data = TSF.from_tfs_file(filename_tsf)

        # Loop over channel and retrieve required data:
        retrieved_channel_locations = []
        retrieved_channel_numbers = []
        stations_corrections_dict = {}
        for channel_metadata in tsf_data.channel_metadata:
            ch_channel_name = channel_metadata.channel_name
            ch_instrument = channel_metadata.instrument_name
            ch_data_type = channel_metadata.data_type
            ch_location = channel_metadata.location
            channel_number = channel_metadata.channel_number

            if location:
                if location != ch_location:
                    continue
            if instrument:
                if instrument != ch_instrument:
                    continue
            if data_type:
                if data_type != ch_data_type:
                    continue
            if ch_location in retrieved_channel_locations:
                raise RuntimeError(f'Ambiguous TSF channel matching! The channels '
                                   f'{retrieved_channel_numbers[retrieved_channel_locations.index(ch_location)]} and '
                                   f'{channel_number} share the same location ({ch_location}). '
                                   f'Please constrain the channel selection by providing additional arguments ('
                                   f'instrument, data type) or edit the TSF file.')
            retrieved_channel_locations.append(ch_location)
            retrieved_channel_numbers.append(channel_number)
            ch_epoch_dt, ch_data  = tsf_data.get_channel_np(channel=channel_number)
            time_series = TimeSeries(ref_time_dt=ch_epoch_dt, data=ch_data, unit=channel_metadata.unit,
                                     data_source=f'{tsf_data.filename_without_path} ({tsf_data.filetype})',
                                     description=ch_channel_name,
                                     is_correction=is_correction)


            stations_corrections_dict[ch_location] = StationCorrections(station_name=ch_location,
                                                                        tidal_correction=time_series)

        # Checks:
        if len(retrieved_channel_locations) == 0:
            tmp_str_list = []
            if location:
                tmp_str_list.append(f'station "{location}"')
            if instrument:
                tmp_str_list.append(f'instrument "{instrument}"')
            if data_type:
                tmp_str_list.append(f'data_type "{data_type}"')
            if tmp_str_list:
                raise RuntimeError(f'Could not retrieve any data from TSF file {filename_tsf} with the following '
                                   f'restrictions:' + ', '.join(tmp_str_list) + '.')
            else:
                raise RuntimeError(f'Could not retrieve any data from TSF file {filename_tsf}.')

        # Add loaded data:
        if survey_name in self.surveys.keys():
            # Add station corrections to existing survey correction object:
            for stat_name, stat_corr in stations_corrections_dict.items():
                _ = self.surveys[survey_name].add_station_correction(station_name=stat_name,
                                                                     station_correction=stat_corr,
                                                                     overwrite=overwrite_channel)
        else:
            # Create survey correction object and add the station correction data:
            self.surveys[survey_name] = SurveyCorrections(survey_name=survey_name, stations=stations_corrections_dict)

    def delete_survey_correction(self, survey_name):
        """Removes a survey correction object.

        Parameters
        ----------
        survey_name: str
            Survey name.
        """
        del self.surveys[survey_name]

__init__()

Default initializer.

Source code in gravtools/tides/correction_time_series.py
42
43
44
def __init__(self):
    """Default initializer."""
    self.surveys = {}

delete_survey_correction(survey_name)

Removes a survey correction object.

Parameters:

Name Type Description Default
survey_name

Survey name.

required
Source code in gravtools/tides/correction_time_series.py
141
142
143
144
145
146
147
148
149
def delete_survey_correction(self, survey_name):
    """Removes a survey correction object.

    Parameters
    ----------
    survey_name: str
        Survey name.
    """
    del self.surveys[survey_name]

load_tfs_file(survey_name, filename_tsf, location='', instrument='', data_type='', overwrite_channel=True, is_correction='False')

Load time series data from a TSoft TSF file for one survey.

Parameters:

Name Type Description Default
survey_name str

Name of the survey for which the data is loaded.

required
filename_tsf str

Name and path of the TSF file from which the correction time series data is loaded.

required
location str

If not empty, only channels with matching locations are loaded.

''
instrument str

If not empty, only channels with matching instruments are loaded.

''
data_type str

If not empty, only channels with matching data types are loaded.

''
overwrite_channel bool

True implies that an existing StationCorrection object of a station will be overwritten.

True
is_correction bool

True implies that the channel data loaded from the TSF file model the gravity effect of phenomena, rather than corrections for gravity observations. Both options have the opposite sign: while corrections have to be added to gravity observations, effects have to be subtracted, in order to reduce observations. See TSoft manual (version 2.2.4 Release date 2015-09-09), p. 15: "The loading calculation by Tsoft will result in the correction, not the effect. On the other hand the prediction of the solid Earth tides using the WDD parameter set provided by Tsoft directly will result in the effect. Therefore both must be treated with different sign in order to reduce a gravity time series correctly."

'False'
Source code in gravtools/tides/correction_time_series.py
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
def load_tfs_file(self, survey_name: str, filename_tsf: str, location: str ='', instrument: str ='',
                  data_type: str ='', overwrite_channel: bool=True, is_correction: bool='False'):
    """Load time series data from a TSoft TSF file for one survey.

    Parameters
    ----------
    survey_name: str
        Name of the survey for which the data is loaded.
    filename_tsf: str
        Name and path of the TSF file from which the correction time series data is loaded.
    location: str, optional (default='')
        If not empty, only channels with matching locations are loaded.
    instrument: str, optional (default='')
        If not empty, only channels with matching instruments are loaded.
    data_type: str, optional (default='')
        If not empty, only channels with matching data types are loaded.
    overwrite_channel: bool, optional (default=`True`)
        `True` implies that an existing StationCorrection object of a station will be overwritten.
    is_correction: bool, optional (default=`False`)
        `True` implies that the channel data loaded from the TSF file model the gravity effect of phenomena, rather
        than corrections for gravity observations. Both options have the opposite sign: while corrections have to be
        added to gravity observations, effects have to be subtracted, in order to reduce observations.
        See TSoft manual (version 2.2.4 Release date 2015-09-09), p. 15: "The loading calculation by Tsoft will
        result in the correction, not the effect.
        On the other hand the prediction of the solid Earth tides using the WDD parameter set
        provided by Tsoft directly will result in the effect. Therefore both must be treated with
        different sign in order to reduce a gravity time series correctly."
    """
    # Load TSF file:
    tsf_data = TSF.from_tfs_file(filename_tsf)

    # Loop over channel and retrieve required data:
    retrieved_channel_locations = []
    retrieved_channel_numbers = []
    stations_corrections_dict = {}
    for channel_metadata in tsf_data.channel_metadata:
        ch_channel_name = channel_metadata.channel_name
        ch_instrument = channel_metadata.instrument_name
        ch_data_type = channel_metadata.data_type
        ch_location = channel_metadata.location
        channel_number = channel_metadata.channel_number

        if location:
            if location != ch_location:
                continue
        if instrument:
            if instrument != ch_instrument:
                continue
        if data_type:
            if data_type != ch_data_type:
                continue
        if ch_location in retrieved_channel_locations:
            raise RuntimeError(f'Ambiguous TSF channel matching! The channels '
                               f'{retrieved_channel_numbers[retrieved_channel_locations.index(ch_location)]} and '
                               f'{channel_number} share the same location ({ch_location}). '
                               f'Please constrain the channel selection by providing additional arguments ('
                               f'instrument, data type) or edit the TSF file.')
        retrieved_channel_locations.append(ch_location)
        retrieved_channel_numbers.append(channel_number)
        ch_epoch_dt, ch_data  = tsf_data.get_channel_np(channel=channel_number)
        time_series = TimeSeries(ref_time_dt=ch_epoch_dt, data=ch_data, unit=channel_metadata.unit,
                                 data_source=f'{tsf_data.filename_without_path} ({tsf_data.filetype})',
                                 description=ch_channel_name,
                                 is_correction=is_correction)


        stations_corrections_dict[ch_location] = StationCorrections(station_name=ch_location,
                                                                    tidal_correction=time_series)

    # Checks:
    if len(retrieved_channel_locations) == 0:
        tmp_str_list = []
        if location:
            tmp_str_list.append(f'station "{location}"')
        if instrument:
            tmp_str_list.append(f'instrument "{instrument}"')
        if data_type:
            tmp_str_list.append(f'data_type "{data_type}"')
        if tmp_str_list:
            raise RuntimeError(f'Could not retrieve any data from TSF file {filename_tsf} with the following '
                               f'restrictions:' + ', '.join(tmp_str_list) + '.')
        else:
            raise RuntimeError(f'Could not retrieve any data from TSF file {filename_tsf}.')

    # Add loaded data:
    if survey_name in self.surveys.keys():
        # Add station corrections to existing survey correction object:
        for stat_name, stat_corr in stations_corrections_dict.items():
            _ = self.surveys[survey_name].add_station_correction(station_name=stat_name,
                                                                 station_correction=stat_corr,
                                                                 overwrite=overwrite_channel)
    else:
        # Create survey correction object and add the station correction data:
        self.surveys[survey_name] = SurveyCorrections(survey_name=survey_name, stations=stations_corrections_dict)

StationCorrections dataclass

Correction time series data for one station.

Notes

For each station and correction type (e.g. tidal correction) only ONE time series can be added!

Source code in gravtools/tides/correction_time_series.py
309
310
311
312
313
314
315
316
317
318
@dataclass
class StationCorrections:
    """Correction time series data for one station.

    Notes
    -----
    For each station and correction type (e.g. tidal correction) only ONE time series can be added!
    """
    station_name: str
    tidal_correction: TimeSeries

SurveyCorrections dataclass

Correction time series data for one survey with multiple stations.

Source code in gravtools/tides/correction_time_series.py
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
@dataclass
class SurveyCorrections:
    """Correction time series data for one survey with multiple stations."""
    survey_name: str
    stations: typing.Dict[str, StationCorrections]

    def add_station_correction(self, station_name: str, station_correction: StationCorrections, overwrite: bool=True,
                               warn: bool=True):
        """Add a StationCorrection object to `self.stations`.

        Parameters
        ----------
        station_name: str
            Station name.
        station_correction: `StationCorrections`
            StationCorrections object containing correction time series data for a specific station.
        overwrite: bool, optional (default=`True`)
            `True` implies that an existing StationCorrection object of a station will be overwritten.
        warn: bool, optional (default=`True`)
            `True` implies that a warning is issued if an existing StationCorrection will or would bne be overwritten.

        Returns
        -------
        bool: `True`, if the station correction data was added. Otherwise, `False`.
        """
        # Check, if data for this station already exists:
        if station_name in self.stations.keys():
            if not overwrite:
                if warn:
                    warnings.warn(f'Station correction data for station {station_name} already exists! Overwriting '
                                  f'data is not permitted!')
                return False
            else:
                if warn:
                    warnings.warn(f'Station correction data for station {station_name} already exists! Overwriting '
                                  f'data is permitted and the existing data will be overwritten!')
                self.stations[station_name] = station_correction
                return True
        else:
            self.stations[station_name] = station_correction
            return True

    def delete_station_correction(self, station_name):
        """Removes a station correction object.

        Parameters
        ----------
        station_name: str
            Station name.
        """
        del self.stations[station_name]

    @property
    def number_of_stations(self):
        """Returns the number of station correction objects."""
        return len(self.stations)

    @property
    def station_names(self):
        """Returns the names of all stations."""
        return list(self.stations.keys())

number_of_stations property

Returns the number of station correction objects.

station_names property

Returns the names of all stations.

add_station_correction(station_name, station_correction, overwrite=True, warn=True)

Add a StationCorrection object to self.stations.

Parameters:

Name Type Description Default
station_name str

Station name.

required
station_correction StationCorrections

StationCorrections object containing correction time series data for a specific station.

required
overwrite bool

True implies that an existing StationCorrection object of a station will be overwritten.

True
warn bool

True implies that a warning is issued if an existing StationCorrection will or would bne be overwritten.

True

Returns:

Name Type Description
bool `True`, if the station correction data was added. Otherwise, `False`.
Source code in gravtools/tides/correction_time_series.py
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
def add_station_correction(self, station_name: str, station_correction: StationCorrections, overwrite: bool=True,
                           warn: bool=True):
    """Add a StationCorrection object to `self.stations`.

    Parameters
    ----------
    station_name: str
        Station name.
    station_correction: `StationCorrections`
        StationCorrections object containing correction time series data for a specific station.
    overwrite: bool, optional (default=`True`)
        `True` implies that an existing StationCorrection object of a station will be overwritten.
    warn: bool, optional (default=`True`)
        `True` implies that a warning is issued if an existing StationCorrection will or would bne be overwritten.

    Returns
    -------
    bool: `True`, if the station correction data was added. Otherwise, `False`.
    """
    # Check, if data for this station already exists:
    if station_name in self.stations.keys():
        if not overwrite:
            if warn:
                warnings.warn(f'Station correction data for station {station_name} already exists! Overwriting '
                              f'data is not permitted!')
            return False
        else:
            if warn:
                warnings.warn(f'Station correction data for station {station_name} already exists! Overwriting '
                              f'data is permitted and the existing data will be overwritten!')
            self.stations[station_name] = station_correction
            return True
    else:
        self.stations[station_name] = station_correction
        return True

delete_station_correction(station_name)

Removes a station correction object.

Parameters:

Name Type Description Default
station_name

Station name.

required
Source code in gravtools/tides/correction_time_series.py
363
364
365
366
367
368
369
370
371
def delete_station_correction(self, station_name):
    """Removes a station correction object.

    Parameters
    ----------
    station_name: str
        Station name.
    """
    del self.stations[station_name]

TimeSeries dataclass

Uniformly or irregularly sampled correction or effect time series for a single channel.

This dataclass stores a paired array of reference epochs and scalar values together with metadata describing the source and physical meaning of the data. It is the lowest-level building block used by :class:StationCorrections and :class:SurveyCorrections to hold, e.g., tidal corrections or atmospheric loading effects that are later interpolated onto individual gravity observations.

The sign convention follows the TSoft convention (see TSoft manual v2.2.4, p. 15):

  • Correction (is_correction=True): the value is added to the raw gravity observation (e.g. the tidal loading correction as exported by TSoft).
  • Effect (is_correction=False): the value represents the gravitational effect of a phenomenon and must be subtracted from the observation to reduce it.

Attributes:

Name Type Description
ref_time_dt numpy.ndarray of datetime64

Reference epochs of the time series, in UTC.

data numpy.ndarray of float

Scalar values associated with each epoch. Must have the same length as ref_time_dt.

unit str

Physical unit of data (e.g. 'nm/s^2', 'µGal'). Used for unit conversion via :data:gravtools.settings.UNIT_CONVERSION_TO_MUGAL.

data_source str

Human-readable identification of the data origin, e.g. the file name and file type from which the series was loaded.

description (str, optional)

Short free-text label for the time series (e.g. the TSF channel name). Defaults to ''.

is_correction (bool, optional)

True if data contains corrections (to be added to observations); False if data contains effects (to be subtracted). Defaults to True.

created_utc_dt datetime

UTC timestamp set automatically in :meth:__post_init__ when the object is created. Not a constructor argument.

Source code in gravtools/tides/correction_time_series.py
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
@dataclass
class TimeSeries:
    """Uniformly or irregularly sampled correction or effect time series for a single channel.

    This dataclass stores a paired array of reference epochs and scalar values together with
    metadata describing the source and physical meaning of the data.  It is the lowest-level
    building block used by :class:`StationCorrections` and :class:`SurveyCorrections` to hold,
    e.g., tidal corrections or atmospheric loading effects that are later interpolated onto
    individual gravity observations.

    The sign convention follows the TSoft convention (see TSoft manual v2.2.4, p. 15):

    * **Correction** (``is_correction=True``): the value is *added* to the raw gravity
      observation (e.g. the tidal loading correction as exported by TSoft).
    * **Effect** (``is_correction=False``): the value represents the *gravitational effect*
      of a phenomenon and must be *subtracted* from the observation to reduce it.

    Attributes
    ----------
    ref_time_dt : numpy.ndarray of datetime64
        Reference epochs of the time series, in UTC.
    data : numpy.ndarray of float
        Scalar values associated with each epoch.  Must have the same length as
        ``ref_time_dt``.
    unit : str
        Physical unit of ``data`` (e.g. ``'nm/s^2'``, ``'µGal'``).  Used for unit
        conversion via :data:`gravtools.settings.UNIT_CONVERSION_TO_MUGAL`.
    data_source : str
        Human-readable identification of the data origin, e.g. the file name and
        file type from which the series was loaded.
    description : str, optional
        Short free-text label for the time series (e.g. the TSF channel name).
        Defaults to ``''``.
    is_correction : bool, optional
        ``True`` if ``data`` contains *corrections* (to be added to observations);
        ``False`` if ``data`` contains *effects* (to be subtracted).
        Defaults to ``True``.
    created_utc_dt : datetime
        UTC timestamp set automatically in :meth:`__post_init__` when the object is
        created.  Not a constructor argument.
    """
    ref_time_dt: np.array  # np.array of DateTime objects
    data: np.array
    unit: str
    data_source: str
    description: str = ''
    is_correction: bool = True
    # created_utc_dt: datetime = None  # Is initialized in __post_init__

    def __post_init__(self):
        """Executed at the very end of __init__"""
        self.created_utc_dt = datetime.now(tz=timezone.utc)

    @property
    def created_datetime_utc(self):
        """Returns the time and date (in UTC) when the times series was created as datetime object."""
        return self.created_utc_dt

    @property
    def created_datetime_utc_str(self):
        """Returns the time and date (in UTC) when the times series was created as string."""
        return self.created_utc_dt.strftime('%Y-%m-%d %H:%M:%S.%f')

    @property
    def start_datetime(self):
        """Returns the start date and time."""
        return min(self.ref_time_dt)

    @property
    def start_datetime_str(self):
        """Returns the start date and time as string."""
        return pd.to_datetime(self.start_datetime).strftime('%Y-%m-%d %H:%M:%S.%f')

    @property
    def end_datetime(self):
        """Returns the end date and time."""
        return max(self.ref_time_dt)

    @property
    def end_datetime_str(self):
        """Returns the end date and time as string."""
        return pd.to_datetime(self.end_datetime).strftime('%Y-%m-%d %H:%M:%S.%f')

    @property
    def duration_timedelta64_ns(self):
        """Returns the duration as `numpy.timedelta64[ns]` object."""
        return max(self.ref_time_dt) - min(self.ref_time_dt)

    def duration_dhms(self):
        """Returns the duration as days, hours, minutes and seconds."""
        dt = max(self.ref_time_dt) - min(self.ref_time_dt)
        days = dt.astype('timedelta64[D]').astype(int)
        hours = dt.astype('timedelta64[h]').astype(int) - days * 24
        minutes = dt.astype('timedelta64[m]').astype(int) - days*24*60 - hours*60
        seconds = dt.astype('timedelta64[s]').astype(float) - days*24*60*60 - hours*60*60 - minutes*60
        return days, hours, minutes, seconds

    @property
    def duration_dhms_str(self):
        """Returns the duration as string (days, hours, minutes and seconds)"""
        days, hours, minutes, seconds = self.duration_dhms()
        return f'{days} days, {hours} hours, {minutes} minutes, {seconds} seconds'

    @property
    def number_of_datapoints(self):
        """Returns the number of datapoints."""
        return len(self.data)

    @property
    def model_type(self):
        """Returns the model type, i.e. whether effects or corrections are modeled."""
        if self.is_correction:
            return 'Correction'
        else:
            return 'Effect'

    def interpolate(self, interp_times: [np.array, typing.List[datetime]], kind: str = 'quadratic',
                    return_correction: bool = False) -> np.ndarray :
        """Return an interpolated value for the given interpolation epoch.

        Parameters
        ----------
        interp_times: list(datetime) or np.array(datetime)
            Interpolation epoch given as `datetime` object w.r.t. UTC! The interpolation times have to bin within the
            time range of the data series.
        kind: str or int, optional (default='quadratic')
            Specifies, which interpolation method is used. This argument is directly passed to
            `scipy.interpolate.inter1`. For options see scipy reference.
        return_correction: bool, optional (default=`False`)
            If `True`, corrections are returned considering the `self.is_correction` attribute.

        Returns
        -------
        `numpy.ndarray`: Interpolated value for the given interpolation times.
        """
        # Convert datetimes to UNIX timestamps, because numerical values are needed for interpolation:
        interp_times_unix = np.fromiter((t.replace(tzinfo=timezone.utc).timestamp() for t in interp_times), float)
        x = self.ref_time_unix
        y = self.data
        interp_func = scipy.interpolate.interp1d(x, y, kind=kind)
        interp_values = interp_func(interp_times_unix)
        if return_correction:
            if not self.is_correction:
                interp_values = interp_values * -1  # Convert from effect to correction!
        return interp_values

    def to_df(self):
        """Returns the time series as pandas DataFrame with the time reference as index (sorted!)."""
        df = pd.DataFrame({'epoch_dt': self.ref_time_dt,'data': self.data})
        df.set_index('epoch_dt', inplace=True)
        return df

    @property
    def ref_time_unix(self):
        """Returns the reference times as UNIX timestamps (seconds since Jan 1, 1970)."""
        return to_unix_seconds(self.ref_time_dt)

created_datetime_utc property

Returns the time and date (in UTC) when the times series was created as datetime object.

created_datetime_utc_str property

Returns the time and date (in UTC) when the times series was created as string.

duration_dhms_str property

Returns the duration as string (days, hours, minutes and seconds)

duration_timedelta64_ns property

Returns the duration as numpy.timedelta64[ns] object.

end_datetime property

Returns the end date and time.

end_datetime_str property

Returns the end date and time as string.

model_type property

Returns the model type, i.e. whether effects or corrections are modeled.

number_of_datapoints property

Returns the number of datapoints.

ref_time_unix property

Returns the reference times as UNIX timestamps (seconds since Jan 1, 1970).

start_datetime property

Returns the start date and time.

start_datetime_str property

Returns the start date and time as string.

__post_init__()

Executed at the very end of init

Source code in gravtools/tides/correction_time_series.py
200
201
202
def __post_init__(self):
    """Executed at the very end of __init__"""
    self.created_utc_dt = datetime.now(tz=timezone.utc)

duration_dhms()

Returns the duration as days, hours, minutes and seconds.

Source code in gravtools/tides/correction_time_series.py
239
240
241
242
243
244
245
246
def duration_dhms(self):
    """Returns the duration as days, hours, minutes and seconds."""
    dt = max(self.ref_time_dt) - min(self.ref_time_dt)
    days = dt.astype('timedelta64[D]').astype(int)
    hours = dt.astype('timedelta64[h]').astype(int) - days * 24
    minutes = dt.astype('timedelta64[m]').astype(int) - days*24*60 - hours*60
    seconds = dt.astype('timedelta64[s]').astype(float) - days*24*60*60 - hours*60*60 - minutes*60
    return days, hours, minutes, seconds

interpolate(interp_times, kind='quadratic', return_correction=False)

Return an interpolated value for the given interpolation epoch.

Parameters:

Name Type Description Default
interp_times [array, List[datetime]]

Interpolation epoch given as datetime object w.r.t. UTC! The interpolation times have to bin within the time range of the data series.

required
kind str

Specifies, which interpolation method is used. This argument is directly passed to scipy.interpolate.inter1. For options see scipy reference.

'quadratic'
return_correction bool

If True, corrections are returned considering the self.is_correction attribute.

False

Returns:

Type Description
`numpy.ndarray`: Interpolated value for the given interpolation times.
Source code in gravtools/tides/correction_time_series.py
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
def interpolate(self, interp_times: [np.array, typing.List[datetime]], kind: str = 'quadratic',
                return_correction: bool = False) -> np.ndarray :
    """Return an interpolated value for the given interpolation epoch.

    Parameters
    ----------
    interp_times: list(datetime) or np.array(datetime)
        Interpolation epoch given as `datetime` object w.r.t. UTC! The interpolation times have to bin within the
        time range of the data series.
    kind: str or int, optional (default='quadratic')
        Specifies, which interpolation method is used. This argument is directly passed to
        `scipy.interpolate.inter1`. For options see scipy reference.
    return_correction: bool, optional (default=`False`)
        If `True`, corrections are returned considering the `self.is_correction` attribute.

    Returns
    -------
    `numpy.ndarray`: Interpolated value for the given interpolation times.
    """
    # Convert datetimes to UNIX timestamps, because numerical values are needed for interpolation:
    interp_times_unix = np.fromiter((t.replace(tzinfo=timezone.utc).timestamp() for t in interp_times), float)
    x = self.ref_time_unix
    y = self.data
    interp_func = scipy.interpolate.interp1d(x, y, kind=kind)
    interp_values = interp_func(interp_times_unix)
    if return_correction:
        if not self.is_correction:
            interp_values = interp_values * -1  # Convert from effect to correction!
    return interp_values

to_df()

Returns the time series as pandas DataFrame with the time reference as index (sorted!).

Source code in gravtools/tides/correction_time_series.py
297
298
299
300
301
def to_df(self):
    """Returns the time series as pandas DataFrame with the time reference as index (sorted!)."""
    df = pd.DataFrame({'epoch_dt': self.ref_time_dt,'data': self.data})
    df.set_index('epoch_dt', inplace=True)
    return df

convert_to_mugal(data, data_unit) staticmethod

Converts data to from a given unit to µGal.

Source code in gravtools/tides/correction_time_series.py
383
384
385
386
387
388
389
@staticmethod
def convert_to_mugal(data, data_unit: str):
    """Converts data to from a given unit to µGal."""
    if data_unit not in settings.UNIT_CONVERSION_TO_MUGAL:
        raise RuntimeError(f'Conversion factor from {data_unit} to µGal not defined. Please add the factor to'
                           f'gravtools.settings.UNIT_CONVERSION_TO_MUGAL.')
    return data * settings.UNIT_CONVERSION_TO_MUGAL[data_unit]

gravtools.tides.longman1959

Longman Earth Tide Calculator - Adapted from John Leeman's implementation of I. M. Longman's earth tide calculations (see references below)

Parts of this program (c) 2017 John Leeman Any other modifications from the base LongmanTide (c) 2018 Zachery Brady

Licensed under the MIT License, see LICENSE text below

References

I.M. Longman "Formulas for Computing the Tidal Accelerations Due to the Moon and the Sun" Journal of Geophysical Research, vol. 64, no. 12, 1959, pp. 2351-2355 P. Schureman "Manual of harmonic analysis and prediction of tides" U.S. Coast and Geodetic Survey, 1958 John Leeman's GitHub page for the original implementation of this library: https://github.com/jrleeman/LongmanTide

Notes

Unicode greek symbols are used to more clearly name variables based on the equations in Longman's paper. This is simply a style decision and may not reflect Python best practices. Python's Datetime.datetime objects are timezone-naive - e.g. when creating a datetime of

datetime(1899, 12, 31, 12, 0, 0) The object refers to exactly the time specified, with no knowledge of the time-zone - when doing calculations against such an object we need to make sure that whatever datetime being used is specified in the same time zone.

LICENSE

MIT License

Copyright (c) 2017 John Leeman Copyright (c) 2018-2020 Zachery Brady

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

calculate_julian_century(dates)

Calculate the decimal Julian century and floating point hour, as referenced from 1899, December 31 at 12:00:00 This function accepts either a numpy ndarray or pandas DatetimeIndex as input, and returns a 2-tuple of the corresponding Julian century decimals and floating point hours.

Parameters:

Name Type Description Default
dates Union[ndarray, DatetimeIndex]

1-dimensional array of DateTime objects to convert into century/hours arrays

required
Notes

All DateTimes should be supplied as UTC values, using a timezone-naive DateTime object. This can be accomplished either by using datetime.utcnow(), or by constructing datetime objects where the times are specified as UTC values Reference Date: 1899 December 31 12:00:00 Reference ordinal: 693961.5 (MATLAB Serial date from January 0, 0000) Delta 366 days Reference ordinal: 694327.5 (Python ordinal from January 1, 0001)

Returns:

Type Description
2-tuple of:

T : np.ndarray Number of Julian centuries (36525 days) from GMT Noon on December 31, 1899 t0 : np.ndarray Greenwich civil dates measured in hours

Source code in gravtools/tides/longman1959.py
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
def calculate_julian_century(dates: Union[np.ndarray, pd.DatetimeIndex]):
    """
    Calculate the decimal Julian century and floating point hour, as referenced from
    1899, December 31 at 12:00:00
    This function accepts either a numpy ndarray or pandas DatetimeIndex as input,
    and returns a 2-tuple of the corresponding Julian century decimals and floating
    point hours.

    Parameters
    ----------
    dates : Union[np.ndarray, pd.DatetimeIndex]
        1-dimensional array of DateTime objects to convert into
        century/hours arrays

    Notes
    -----
    All DateTimes should be supplied as UTC values, using a timezone-naive DateTime object.
    This can be accomplished either by using datetime.utcnow(), or by constructing datetime
    objects where the times are specified as UTC values
    Reference Date: 1899 December 31 12:00:00
    Reference ordinal: 693961.5 (MATLAB Serial date from January 0, 0000) Delta 366 days
    Reference ordinal: 694327.5 (Python ordinal from January 1, 0001)

    Returns
    -------
    2-tuple of:
        T : np.ndarray
            Number of Julian centuries (36525 days) from GMT Noon on
            December 31, 1899
        t0 : np.ndarray
            Greenwich civil dates measured in hours
    """
    if isinstance(dates, np.ndarray):
        delta = dates - origin_date
        days = np.array([x.days + x.seconds / 3600. / 24. for x in delta])
        t0 = np.array([x.hour + x.minute / 60. + x.second / 3600. for x in dates])
        return days / 36525, t0
    elif isinstance(dates, pd.DatetimeIndex):
        delta = dates - origin_date
        days = delta.days + delta.seconds / 3600. / 24.
        t0 = dates.hour + dates.minute / 60. + dates.second / 3600.
        return days / 36525, np.array(t0, dtype=float)

solve_longman_tide(lat, lon, alt, time)

Find the total gravity correction due to the Sun/Moon for the given latitude, longitude, altitude, and time - supplied as numpy 1-d arrays. Corrections are calculated for each corresponding set of data in the supplied arrays, all arrays must be of the same shape (length and dimension) This function returns the Lunar, Solar, and Total gravity corrections as a 3-tuple of numpy 1-d arrays.

Parameters:

Name Type Description Default
lat ndarray

1-dimensional array of float values denoting latitudes

required
lon ndarray

1-dimensional array of float values denoting longitudes

required
alt ndarray

1-dimensional array of float values denoting altitude in meters

required
time ndarray

1-dimensional array of DateTime objects denoting the time series

required

Returns:

Type Description
3 - Tuple

gMoon (gm): Vertical component of tidal acceleration due to the moon gSun (gs): Vertical component of tidal acceleration due to the sun gTotal (g0): Total vertical component of tidal acceleration (moon + sun)

Source code in gravtools/tides/longman1959.py
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
def solve_longman_tide(lat: np.ndarray, lon: np.ndarray, alt: np.ndarray, time: np.ndarray):
    """
    Find the total gravity correction due to the Sun/Moon for the given
    latitude, longitude, altitude, and time - supplied as numpy 1-d arrays.
    Corrections are calculated for each corresponding set of data in the
    supplied arrays, all arrays must be of the same shape (length and dimension)
    This function returns the Lunar, Solar, and Total gravity corrections as
    a 3-tuple of numpy 1-d arrays.

    Parameters
    ----------
    lat : np.ndarray
        1-dimensional array of float values denoting latitudes
    lon : np.ndarray
        1-dimensional array of float values denoting longitudes
    alt : np.ndarray
        1-dimensional array of float values denoting altitude in meters
    time : np.ndarray
        1-dimensional array of DateTime objects denoting the time series

    Returns
    -------
    3-Tuple
        gMoon (gm): Vertical component of tidal acceleration due to the moon
        gSun (gs): Vertical component of tidal acceleration due to the sun
        gTotal (g0): Total vertical component of tidal acceleration (moon + sun)
    """
    assert lat.shape == lon.shape == alt.shape == time.shape

    T, t0 = calculate_julian_century(time)
    T2 = T ** 2
    T3 = T ** 3

    # t0 must be an ndarray, not a pandas DatetimeIndex
    t0[t0 < 0] += 24
    t0[t0 >= 24] -= 24

    # longitude is defined with West positive (+) in the Longman paper
    φ = -lon
    λ = np.radians(lat)  # λ Latitude of point P
    cosλ = np.cos(λ)
    sinλ = np.sin(λ)
    H = alt * 100  # height above sea-level of point P in centimeters (cm)

    # Lunar Calculations

    # s Mean longitude of moon in its orbit reckoned from the referred equinox
    # Constants from Bartels [1957 pp. 747] eq (10')
    # 270°26'11.72" + (1336 rev. + 1,108,406.05")T + 7.128" * T2 + 0.0072" * T3
    # Where rev. is revolutions expressed as radians: 1336 rev. = 1336 * 2 * pi = 8394.335570392
    s = 4.72000889397 + 8399.70927456 * T + 3.45575191895e-05 * T2 + 3.49065850399e-08 * T3
    # p Mean longitude of lunar perigee
    # constants from Bartels [1957] eq (11')
    # p = 334° 19' 46.42" + (11 rev. + 392,522.51") T - 37.15" * T2 - 0.036" T3
    p = 5.83515162814 + 71.0180412089 * T + 0.000180108282532 * T2 + 1.74532925199e-07 * T3
    # (h) Mean longitude of the sun
    h = 4.88162798259 + 628.331950894 * T + 5.23598775598e-06 * T2
    # (N) Longitude of the moon's ascending node in its orbit reckoned from the referred equinox
    N = 4.52360161181 - 33.757146295 * T + 3.6264063347e-05 * T2 + 3.39369576777e-08 * T3
    cosN = np.cos(N)
    sinN = np.sin(N)

    # I (uppercase i) Inclination of the moon's orbit to the equator
    I = np.arccos(np.cos(ω) * np.cos(i) - np.sin(ω) * np.sin(i) * cosN)
    # ν (nu) Longitude in the celestial equator of its intersection A with the moon's orbit
    ν = np.arcsin(np.sin(i) * sinN / np.sin(I))
    # t Hour angle of mean sun measured west-ward from the place of observations
    t = np.radians(15. * (t0 - 12) - φ)

    # χ (chi) right ascension of meridian of place of observations reckoned from A
    χ = t + h - ν
    # cos α (alpha) where α is defined in eq. 15 and 16
    cos_α = cosN * np.cos(ν) + sinN * np.sin(ν) * np.cos(ω)
    # sin α (alpha) where α is defined in eq. 15 and 16
    sin_α = np.sin(ω) * sinN / np.sin(I)
    # (α) α is defined in eq. 15 and 16
    α = 2 * np.arctan(sin_α / (1 + cos_α))
    # ξ (xi) Longitude in the moon's orbit of its ascending intersection with the celestial equator
    ξ = N - α

    # σ (sigma) Mean longitude of moon in radians in its orbit reckoned from A
    σ = s - ξ
    # l (lowercase el) Longitude of moon in its orbit reckoned from its ascending intersection with the equator
    l = σ + 2 * e * np.sin(s - p) + (5. / 4) * e * e * np.sin(2 * (s - p)) + (15. / 4) * m * e * np.sin(s - 2 * h + p)\
        + (11. / 8) * m * m * np.sin(2 * (s - h))

    # Solar Calculations

    # p1 (p-one) Longitude of solar perigee (Schureman [1941, pp. 162])
    # p1 = 281° 13' 15.0" + 6189.03" T + 1.63" T2 + 0.012" T3
    p1 = 4.90822941839 + 0.0300025492114 * T + 7.85398163397e-06 * T2 + 5.3329504922e-08 * T3
    # e1 (e-one) Eccentricity of the Earth's orbit
    e1 = 0.01675104 - 0.00004180 * T - 0.000000126 * T2
    # χ1 (chi-one) right ascension of meridian of place of observations reckoned from the vernal equinox
    χ1 = t + h
    # l1 (lowercase-el(L) one) Longitude of sun in the ecliptic reckoned from the vernal equinox
    l1 = h + 2 * e1 * np.sin(h - p1)
    # cosθ (theta) θ represents the zenith angle of the moon
    cosθ = sinλ * np.sin(I) * np.sin(l) + cosλ * (np.cos(0.5 * I) ** 2 * np.cos(l - χ) + np.sin(0.5 * I) ** 2 *
                                                  np.cos(l + χ))
    # cosφ (phi) φ represents the zenith angle of the run
    cosφ = sinλ * np.sin(ω) * np.sin(l1) + cosλ * (np.cos(0.5 * ω) ** 2 * np.cos(l1 - χ1) + np.sin(0.5 * ω) ** 2 *
                                                np.cos(l1 + χ1))

    # Distance Calculations

    # (C) Distance parameter, equation 34
    # C**2 = 1/(1 + 0.006738 sinλ ** 2)
    C = np.sqrt(1. / (1 + 0.006738 * sinλ ** 2))
    # (r) Distance from point P to the center of the Earth
    r = C * a + H
    # a' (a prime) Distance parameter, equation 31
    aprime = 1. / (c * (1 - e * e))
    # a1' (a-one prime) Distance parameter, equation 31
    aprime1 = 1. / (c1 * (1 - e1 * e1))
    # (d) Distance between centers of the Earth and the moon
    d = 1. / ((1. / c) + aprime * e * np.cos(s - p) + aprime * e ** 2 * np.cos(2 * (s - p)) + (15. / 8) * aprime * m *
              e * np.cos(s - 2 * h + p) + aprime * m * m * np.cos(2 * (s - h)))
    # (D) Distance between centers of the Earth and the sun
    D = 1. / ((1. / c1) + aprime1 * e1 * np.cos(h - p1))

    # (gm) Vertical component of tidal acceleration due to the moon, equation (1):
    # gm = μMr/d^3 (3 * cos^2(θ) - 1) + 3/2 μMr/d^4 * (5 cos^3 θ - 3 cos(θ))
    gm = (μ * M * r / d ** 3) * (3 * cosθ ** 2 - 1) + (1.5 * (μ * M * r ** 2 / d ** 4) * (5 * cosθ ** 3 - 3 * cosθ))
    # (gs) Vertical component of tidal acceleration due to the sun
    gs = μ * S * r / D ** 3 * (3 * cosφ ** 2 - 1)

    # g0 Total vertical component due to Lunar and Solar forces
    g0 = (gm + gs) * 1e3 * love

    # Returns Lunar, Solar, Total corrections in mGals, as numpy 1d-arrays
    return gm * 1e3 * love, gs * 1e3 * love, g0

solve_longman_tide_scalar(lat, lon, alt, time)

Simple wrapper around solve_longman_tide that allows passing of singular scalar values instead of an array. This function simply wraps the scalar values within an ndarray and then extracts the result elements.

Parameters:

Name Type Description Default
lat float

Latitude as a scalar float value

required
lon float

Longitude as a scalar float value

required
alt float

Altitude as a scalar float value

required
time datetime

Time as a scalar datetime object

required

Returns:

Type Description
gm, gs, g0 : float, float, float

Where gm is the gravitational effect due to the moon gs is the gravitational effect due to the sun g0 is the total gravitational effect of the sun and moon (gm+gs)

Source code in gravtools/tides/longman1959.py
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
def solve_longman_tide_scalar(lat: float, lon: float, alt: float, time: datetime):
    """
    Simple wrapper around solve_longman_tide that allows passing of singular scalar values instead of an array.
    This function simply wraps the scalar values within an ndarray and then extracts the result elements.
    Parameters
    ----------
    lat : float
        Latitude as a scalar float value
    lon : float
        Longitude as a scalar float value
    alt : float
        Altitude as a scalar float value
    time : datetime
        Time as a scalar datetime object
    Returns
    -------
    gm, gs, g0 : float, float, float
        Where gm is the gravitational effect due to the moon
              gs is the gravitational effect due to the sun
              g0 is the total gravitational effect of the sun and moon (gm+gs)
    """
    lat_arr = np.array([lat])
    lon_arr = np.array([lon])
    alt_arr = np.array([alt])
    time_arr = np.array([time])

    gm, gs, g0 = solve_longman_tide(lat_arr, lon_arr, alt_arr, time_arr)
    return gm[0], gs[0], g0[0]

solve_point_corr(lat, lon, alt, t0=datetime.utcnow(), n=3600, increment='S')

Utility function to generate a tide correction DataFrame for a static lat/lon/alt given start time t0, an increment, and count (n) of datapoints to generate. Default parameters are supplied that will generate a correction series with one second increment over a one hour period, with start time being the time of execution.

Parameters:

Name Type Description Default
lat float

Latitude in decimal degrees

required
lon float

Longitude in decimal degrees

required
alt float

Altitude (height) above sea level in meters

required
t0 DateTime

Starting date/time

utcnow()
n int

Number of data points to generate

3600
increment String

Increment between data points, uses Pandas offset aliases: https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.Timedelta.html#pandas.Timedelta

Common options: * H : Hourly Frequency * T, min : Minutely frequency * S : Secondly frequency * L, ms : millisecond frequency * U, us : microsecond frequency * N : nanosecond frequency

'S'

Returns:

Name Type Description
df DataFrame

DataFrame (index=DateTime, lat, lon, alt, gm, gs, g0) shape: (n, 4) Pandas DataFrame indexed by DateTime containing the latitude, longitude, altitude, lunar, solar, and total corrections.

Notes

total_corr column renamed to g0 and gm/gs added in v0.3.0

Source code in gravtools/tides/longman1959.py
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
def solve_point_corr(lat: float, lon: float, alt: float, t0=datetime.utcnow(), n=3600, increment='S'):
    """
    Utility function to generate a tide correction DataFrame for a static lat/lon/alt given start time t0,
    an increment, and count (n) of datapoints to generate.
    Default parameters are supplied that will generate a correction series with one second increment over a one hour
    period, with start time being the time of execution.

    Parameters
    ----------
    lat : float
        Latitude in decimal degrees
    lon : float
        Longitude in decimal degrees
    alt : float
        Altitude (height) above sea level in meters
    t0 : DateTime
        Starting date/time
    n : int
        Number of data points to generate
    increment : String
        Increment between data points, uses Pandas offset aliases:
        https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.Timedelta.html#pandas.Timedelta

        Common options:
            * H : Hourly Frequency
            * T, min : Minutely frequency
            * S : Secondly frequency
            * L, ms : millisecond frequency
            * U, us : microsecond frequency
            * N : nanosecond frequency

    Returns
    -------
    df : pd.DataFrame
        DataFrame (index=DateTime, lat, lon, alt, gm, gs, g0) shape: (n, 4)
        Pandas DataFrame indexed by DateTime containing the latitude, longitude, altitude, lunar, solar,
        and total corrections.
    Notes
    -----
    total_corr column renamed to g0 and gm/gs added in v0.3.0
    """
    df = pd.DataFrame(data={'lat': np.repeat(lat, n), 'lon': np.repeat(lon, n), 'alt': np.repeat(alt, n)},
                      index=pd.date_range(start=t0, freq=increment, periods=n))

    gm, gs, g0 = solve_longman_tide(df.lat, df.lon, df.alt, df.index)
    df['gm'] = gm
    df['gs'] = gs
    df['g0'] = g0
    return df

solve_tide_df(df, lat='lat', lon='lon', alt='alt')

Solve tidal gravity corrections for a given Pandas DataFrame, returning the source DataFrame with an appended 'tide_corr' column. Source DataFrame should be indexed by datetime and must contain latitude, longitude, and altitude columns. Source column names can be specified by passing the applicable column name to the lat, lon, alt parameters as required. Time is assumed to be the index of the DataFrame

Returns:

Name Type Description
df DataFrame

Copy of source df with added columns 'gm' 'gs' and 'g0': lunar, solar and total corrections respectively

Source code in gravtools/tides/longman1959.py
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
def solve_tide_df(df: pd.DataFrame, lat='lat', lon='lon', alt='alt'):
    """
    Solve tidal gravity corrections for a given Pandas DataFrame, returning the source
    DataFrame with an appended 'tide_corr' column.
    Source DataFrame should be indexed by datetime and must contain latitude, longitude,
    and altitude columns.
    Source column names can be specified by passing the applicable column name to the
    lat, lon, alt parameters as required.
    Time is assumed to be the index of the DataFrame
    Returns
    -------
    df : pd.DataFrame
        Copy of source df with added columns 'gm' 'gs' and 'g0': lunar, solar and total corrections respectively
    """
    _lat = df[lat].values
    _lon = df[lon].values
    _alt = df[alt].values
    _time = df.index
    res_df = df.copy(deep=True)  # type: pd.DataFrame

    gm, gs, g0 = solve_longman_tide(_lat, _lon, _alt, _time)
    res_df['gm'] = gm
    res_df['gs'] = gs
    res_df['g0'] = g0
    return res_df