diff --git a/.github/workflows/auto-update-files.yml b/.github/workflows/auto-update-files.yml deleted file mode 100644 index db9e2fe8..00000000 --- a/.github/workflows/auto-update-files.yml +++ /dev/null @@ -1,70 +0,0 @@ -# This workflow will install Python dependencies and update the time and EOP files - -name: Auto-Update Files - -on: - schedule: - # * is a special character in YAML so you have to quote this string - - cron: '0 0 1 * *' - workflow_dispatch: - -jobs: - build: - - runs-on: ${{ matrix.os }} - strategy: - matrix: - os: [ubuntu-latest] - python-version: [3.11] - env: - OS: ${{ matrix.os }} - PYTHON: ${{ matrix.python-version }} - defaults: - run: - shell: bash -l {0} - - steps: - - uses: actions/checkout@v4 - - name: Set up mamba ${{ matrix.python-version }} - uses: mamba-org/setup-micromamba@v1 - with: - micromamba-version: 'latest' - environment-file: .binder/environment.yml - init-shell: bash - environment-name: gravity_toolkit - cache-environment: true - post-cleanup: 'all' - create-args: >- - python=${{ matrix.python-version }} - flake8 - pytest - pytest-cov - cython - - name: Lint with flake8 - run: | - # stop the build if there are Python syntax errors or undefined names - flake8 . --count --select=E9,F63,F7,F82 --show-source --statistics - # exit-zero treats all errors as warnings. The GitHub editor is 127 chars wide - flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics - - name: Test with pytest - run: | - pip install --no-deps . - pytest test/test_leap_seconds.py test/test_eop.py \ - --username=${{ secrets.EARTHDATA_USERNAME }} \ - --password=${{ secrets.EARTHDATA_PASSWORD }} - - name: Check for changes - id: changes - run: | - if [ -n "$(git status --porcelain)" ] ; then - echo "DETECTED=true" >> $GITHUB_OUTPUT; - echo "::debug::Changes detected"; - else - echo "DETECTED=false" >> $GITHUB_OUTPUT; - echo "::debug::No changes detected"; - fi - - name: Create pull request - if: steps.changes.outputs.DETECTED == 'true' - uses: peter-evans/create-pull-request@v5 - with: - assignees: ${{ github.actor }} - title: "data: automatic time updates" diff --git a/doc/source/api_reference/time.rst b/doc/source/api_reference/time.rst index 9b077af6..339666e9 100644 --- a/doc/source/api_reference/time.rst +++ b/doc/source/api_reference/time.rst @@ -82,12 +82,3 @@ General Methods .. autofunction:: gravity_toolkit.time.convert_calendar_decimal .. autofunction:: gravity_toolkit.time.convert_julian - -.. autoclass:: gravity_toolkit.time.timescale - :members: - -.. autofunction:: gravity_toolkit.time.count_leap_seconds - -.. autofunction:: gravity_toolkit.time.get_leap_seconds - -.. autofunction:: gravity_toolkit.time.update_leap_seconds diff --git a/doc/source/api_reference/time_series/lomb_scargle.rst b/doc/source/api_reference/time_series/lomb_scargle.rst new file mode 100644 index 00000000..f88f5a39 --- /dev/null +++ b/doc/source/api_reference/time_series/lomb_scargle.rst @@ -0,0 +1,11 @@ +======================== +time_series.lomb_scargle +======================== + +- Wrapper function for computing Lomb-Scargle periodograms using ``scipy.signal.lombscargle`` + +`Source code`__ + +.. __: https://github.com/tsutterley/gravity-toolkit/blob/main/gravity_toolkit/time_series/lomb_scargle.py + +.. autofunction:: gravity_toolkit.time_series.lomb_scargle diff --git a/doc/source/api_reference/time_series/regress.rst b/doc/source/api_reference/time_series/regress.rst index 892574aa..75fc1612 100644 --- a/doc/source/api_reference/time_series/regress.rst +++ b/doc/source/api_reference/time_series/regress.rst @@ -10,7 +10,7 @@ Calling Sequence .. code-block:: python import gravity_toolkit.time_series - tsbeta = gravity_toolkit.time_series.regress(t_in, d_in, ORDER=1, CYCLES=[0.5,1.0], CONF=0.95) + tsbeta = gravity_toolkit.time_series.regress(t_in, d_in, ORDER=1, CYCLES=[0.5,1.0], CONF=0.95) `Source code`__ diff --git a/doc/source/index.rst b/doc/source/index.rst index 5a80ba12..d60e5b78 100644 --- a/doc/source/index.rst +++ b/doc/source/index.rst @@ -74,6 +74,7 @@ missions api_reference/time.rst api_reference/time_series/amplitude.rst api_reference/time_series/fit.rst + api_reference/time_series/lomb_scargle.rst api_reference/time_series/piecewise.rst api_reference/time_series/regress.rst api_reference/time_series/savitzky_golay.rst diff --git a/gravity_toolkit/geocenter.py b/gravity_toolkit/geocenter.py index 4c4994c0..4926bad0 100644 --- a/gravity_toolkit/geocenter.py +++ b/gravity_toolkit/geocenter.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" geocenter.py -Written by Tyler Sutterley (09/2023) +Written by Tyler Sutterley (06/2024) Data class for reading and processing geocenter data PYTHON DEPENDENCIES: @@ -15,6 +15,8 @@ https://github.com/yaml/pyyaml UPDATE HISTORY: + Updated 06/2024: use wrapper to importlib for optional dependencies + Updated 05/2024: make subscriptable and allow item assignment Updated 09/2023: add group option to netCDF read function add functions to return variables or class attributes prevent double printing of filenames when using debug @@ -47,17 +49,13 @@ import time import uuid import yaml -import logging import pathlib -import warnings import numpy as np import gravity_toolkit.time +from gravity_toolkit.utilities import import_dependency # attempt imports -try: - import netCDF4 -except (AttributeError, ImportError, ModuleNotFoundError) as exc: - warnings.warn("netCDF4 not available", ImportWarning) +netCDF4 = import_dependency('netCDF4') class geocenter(object): """ @@ -1300,3 +1298,9 @@ def __next__(self): # add to index self.__index__ += 1 return temp + + def __getitem__(self, key): + return getattr(self, key) + + def __setitem__(self, key, value): + setattr(self, key, value) diff --git a/gravity_toolkit/harmonics.py b/gravity_toolkit/harmonics.py index c10fd9c1..c60a0c5c 100644 --- a/gravity_toolkit/harmonics.py +++ b/gravity_toolkit/harmonics.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" harmonics.py -Written by Tyler Sutterley (10/2023) +Written by Tyler Sutterley (06/2024) Contributions by Hugo Lecomte Spherical harmonic data class for processing GRACE/GRACE-FO Level-2 data @@ -25,6 +25,8 @@ destripe_harmonics.py: filters spherical harmonics for correlated errors UPDATE HISTORY: + Updated 06/2024: use wrapper to importlib for optional dependencies + Updated 05/2024: make subscriptable and allow item assignment Updated 10/2023: place time and month variables in try/except block Updated 09/2023: prevent double printing of filenames when using debug Updated 08/2023: add string representation of the harmonics object @@ -100,28 +102,19 @@ import logging import pathlib import zipfile -import warnings import numpy as np import gravity_toolkit.version from gravity_toolkit.time import adjust_months,calendar_to_grace from gravity_toolkit.destripe_harmonics import destripe_harmonics from gravity_toolkit.read_gfc_harmonics import read_gfc_harmonics from gravity_toolkit.read_GRACE_harmonics import read_GRACE_harmonics -from gravity_toolkit.utilities import reify +from gravity_toolkit.utilities import import_dependency, reify # attempt imports -try: - from geoid_toolkit.read_ICGEM_harmonics import read_ICGEM_harmonics -except (AttributeError, ImportError, ModuleNotFoundError) as exc: - warnings.warn("geoid_toolkit not available", ImportWarning) -try: - import h5py -except (AttributeError, ImportError, ModuleNotFoundError) as exc: - warnings.warn("h5py not available", ImportWarning) -try: - import netCDF4 -except (AttributeError, ImportError, ModuleNotFoundError) as exc: - warnings.warn("netCDF4 not available", ImportWarning) +geoidtk = import_dependency('geoid_toolkit') +h5py = import_dependency('h5py') +netCDF4 = import_dependency('netCDF4') +sparse = import_dependency('sparse') class harmonics(object): """ @@ -520,7 +513,7 @@ def from_gfc(self, filename, **kwargs): Ylms = read_gfc_harmonics(self.filename, TIDE=kwargs['tide']) else: - Ylms = read_ICGEM_harmonics(self.filename, + Ylms = geoidtk.read_ICGEM_harmonics(self.filename, TIDE=kwargs['tide']) # Output file information logging.info(self.filename) @@ -1188,6 +1181,18 @@ def to_masked_array(self): self.squeeze() # return the triangular matrix return Ylms + + def to_coo_array(self): + """ + Convert data arrays to a COO sparse matrices + """ + # assign degree and order fields + self.update_dimensions() + # create COO sparse matrices + self.clm = sparse.COO(self.clm) + self.slm = sparse.COO(self.slm) + # return the harmonics object + return self def update_dimensions(self): """ @@ -1930,3 +1935,9 @@ def __next__(self): # add to index self.__index__ += 1 return temp + + def __getitem__(self, key): + return getattr(self, key) + + def __setitem__(self, key, value): + setattr(self, key, value) diff --git a/gravity_toolkit/read_gfc_harmonics.py b/gravity_toolkit/read_gfc_harmonics.py index 5f44d06d..455f68e8 100644 --- a/gravity_toolkit/read_gfc_harmonics.py +++ b/gravity_toolkit/read_gfc_harmonics.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" read_gfc_harmonics.py -Written by Tyler Sutterley (05/2023) +Written by Tyler Sutterley (06/2023) Contributions by Hugo Lecomte Reads gfc files and extracts spherical harmonics for Swarm and @@ -55,6 +55,7 @@ calculate_tidal_offset.py: calculates the C20 offset for a tidal system UPDATE HISTORY: + Updated 06/2024: use wrapper to importlib for optional dependencies Updated 05/2023: use pathlib to define and operate on paths Updated 03/2023: improve typing for variables in docstrings Updated 04/2022: updated docstrings to numpy documentation format @@ -69,15 +70,12 @@ """ import re import pathlib -import warnings import numpy as np import gravity_toolkit.time +from gravity_toolkit.utilities import import_dependency # attempt imports -try: - from geoid_toolkit.read_ICGEM_harmonics import read_ICGEM_harmonics -except (AttributeError, ImportError, ModuleNotFoundError) as exc: - warnings.warn("geoid_toolkit not available", ImportWarning) +geoidtk = import_dependency('geoid_toolkit') # PURPOSE: read spherical harmonic coefficients of a gravity model def read_gfc_harmonics(input_file, TIDE=None, FLAG='gfc'): @@ -199,7 +197,7 @@ def read_gfc_harmonics(input_file, TIDE=None, FLAG='gfc'): # python dictionary with model input and headers ZIP = bool(re.search('ZIP', SFX, re.IGNORECASE)) - model_input = read_ICGEM_harmonics(input_file, TIDE=TIDE, + model_input = geoidtk.read_ICGEM_harmonics(input_file, TIDE=TIDE, FLAG=FLAG, ZIP=ZIP) # start and end day of the year diff --git a/gravity_toolkit/read_love_numbers.py b/gravity_toolkit/read_love_numbers.py index 16cf158b..c8742fbe 100755 --- a/gravity_toolkit/read_love_numbers.py +++ b/gravity_toolkit/read_love_numbers.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" read_love_numbers.py -Written by Tyler Sutterley (08/2023) +Written by Tyler Sutterley (05/2024) Reads sets of load Love numbers from PREM and applies isomorphic parameters Linearly interpolates load Love/Shida numbers for missing degrees @@ -56,6 +56,7 @@ 103(B12), 30205-30229, (1998) UPDATE HISTORY: + Updated 05/2024: make subscriptable and allow item assignment Updated 08/2023: add string representation of the love_numbers class Updated 05/2023: use pathlib to define and operate on paths Updated 03/2023: improve typing for variables in docstrings @@ -621,3 +622,9 @@ def __iter__(self): yield self.hl yield self.kl yield self.ll + + def __getitem__(self, key): + return getattr(self, key) + + def __setitem__(self, key, value): + setattr(self, key, value) diff --git a/gravity_toolkit/spatial.py b/gravity_toolkit/spatial.py index c68e5af3..ad34a75e 100644 --- a/gravity_toolkit/spatial.py +++ b/gravity_toolkit/spatial.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" spatial.py -Written by Tyler Sutterley (10/2023) +Written by Tyler Sutterley (06/2024) Data class for reading, writing and processing spatial data @@ -20,6 +20,8 @@ time.py: utilities for calculating time operations UPDATE HISTORY: + Updated 06/2024: use wrapper to importlib for optional dependencies + Updated 05/2024: make subscriptable and allow item assignment Updated 10/2023: place time and month variables in try/except block Updated 09/2023: prevent double printing of filenames when using debug Updated 08/2023: add string representation of the spatial object @@ -80,20 +82,14 @@ import logging import pathlib import zipfile -import warnings import numpy as np import gravity_toolkit.version from gravity_toolkit.time import adjust_months, calendar_to_grace +from gravity_toolkit.utilities import import_dependency # attempt imports -try: - import h5py -except (AttributeError, ImportError, ModuleNotFoundError) as exc: - warnings.warn("h5py not available", ImportWarning) -try: - import netCDF4 -except (AttributeError, ImportError, ModuleNotFoundError) as exc: - warnings.warn("netCDF4 not available", ImportWarning) +h5py = import_dependency('h5py') +netCDF4 = import_dependency('netCDF4') class spatial(object): """ @@ -2019,3 +2015,9 @@ def update_mask(self): if getattr(self, 'magnitude') is not None: self.magnitude[self.mask] = self.fill_value return self + + def __getitem__(self, key): + return getattr(self, key) + + def __setitem__(self, key, value): + setattr(self, key, value) diff --git a/gravity_toolkit/time.py b/gravity_toolkit/time.py index 90491f5e..e28abba0 100644 --- a/gravity_toolkit/time.py +++ b/gravity_toolkit/time.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" time.py -Written by Tyler Sutterley (06/2023) +Written by Tyler Sutterley (06/2024) Contributions by Hugo Lecomte Utilities for calculating time operations @@ -13,6 +13,7 @@ https://dateutil.readthedocs.io/en/stable/ UPDATE HISTORY: + Updated 06/2024: remove timescale class and leap second calculations Updated 06/2023: add timescale class for converting between time scales add functions for retrieving leap seconds from NIST or IERF servers Updated 05/2023: convert a string with time zone information to datetime @@ -854,311 +855,3 @@ def convert_julian(JD, **kwargs): return (year, month, day, hour, minute, second) elif (kwargs['format'] == 'zip'): return zip(year, month, day, hour, minute, second) - -class timescale: - """ - Class for converting between time scales - - Attributes - ---------- - leaps: np.ndarray - Number of leap seconds - MJD: np.ndarray - Modified Julian Days - day: float - Seconds in a day - """ - def __init__(self, MJD=None): - # leap seconds - self.leaps = None - # modified Julian Days - self.MJD = MJD - # seconds per day - self.day = 86400.0 - # iterator - self.__index__ = 0 - - def from_deltatime(self, - delta_time: np.ndarray, - epoch: str | tuple | list | np.ndarray, - standard: str = 'UTC' - ): - """ - Converts a delta time array and into a ``timescale`` object - - Parameters - ---------- - delta_time: np.ndarray - seconds since ``epoch`` - epoch: str, uuple, list or np.ndarray - epoch for input ``delta_time`` - standard: str, default 'UTC' - time standard for input ``delta_time`` - """ - # assert delta time is an array - delta_time = np.atleast_1d(delta_time) - # calculate leap seconds if specified - if (standard.upper() == 'GPS'): - GPS_Epoch_Time = convert_delta_time(0, epoch1=epoch, - epoch2= _gps_epoch, scale=1.0) - GPS_Time = convert_delta_time(delta_time, epoch1=epoch, - epoch2=_gps_epoch, scale=1.0) - # calculate difference in leap seconds from start of epoch - self.leaps = count_leap_seconds(GPS_Time) - \ - count_leap_seconds(np.atleast_1d(GPS_Epoch_Time)) - elif (standard.upper() == 'LORAN'): - # LORAN time is ahead of GPS time by 9 seconds - GPS_Epoch_Time = convert_delta_time(-9.0, epoch1=epoch, - epoch2=_gps_epoch, scale=1.0) - GPS_Time = convert_delta_time(delta_time - 9.0, epoch1=epoch, - epoch2=_gps_epoch, scale=1.0) - # calculate difference in leap seconds from start of epoch - self.leaps = count_leap_seconds(GPS_Time) - \ - count_leap_seconds(np.atleast_1d(GPS_Epoch_Time)) - elif (standard.upper() == 'TAI'): - # TAI time is ahead of GPS time by 19 seconds - GPS_Epoch_Time = convert_delta_time(-19.0, epoch1=epoch, - epoch2=_gps_epoch, scale=1.0) - GPS_Time = convert_delta_time(delta_time-19.0, epoch1=epoch, - epoch2=_gps_epoch, scale=1.0) - # calculate difference in leap seconds from start of epoch - self.leaps = count_leap_seconds(GPS_Time) - \ - count_leap_seconds(np.atleast_1d(GPS_Epoch_Time)) - else: - self.leaps = 0.0 - # convert time to days relative to Modified Julian days in UTC - self.MJD = convert_delta_time(delta_time - self.leaps, - epoch1=epoch, epoch2=_mjd_epoch, scale=(1.0/self.day)) - return self - - def from_datetime(self, dtime: np.ndarray): - """ - Reads a ``datetime`` array and converts into a ``timescale`` object - - Parameters - ---------- - dtime: np.ndarray - ``numpy.datetime64`` array - """ - # convert delta time array from datetime object - # to days relative to 1992-01-01T00:00:00 - self.MJD = convert_datetime(dtime, epoch=_mjd_epoch)/self.day - return self - - def to_deltatime(self, - epoch: str | tuple | list | np.ndarray, - scale: float = 1.0 - ): - """ - Convert a ``timescale`` object to a delta time array - - Parameters - ---------- - epoch: str, tuple, list, or np.ndarray - epoch for output ``delta_time`` - scale: float, default 1.0 - scaling factor for converting time to output units - - Returns - ------- - delta_time: np.ndarray - time since epoch - """ - # convert epochs to numpy datetime variables - epoch1 = np.datetime64(datetime.datetime(*_mjd_epoch)) - if isinstance(epoch, (tuple, list)): - epoch = np.datetime64(datetime.datetime(*epoch)) - elif isinstance(epoch, str): - epoch = np.datetime64(parse(epoch)) - # calculate the difference in epochs in days - delta_time_epochs = (epoch - epoch1) / np.timedelta64(1, 'D') - # return the date in time (default days) since epoch - return scale*np.array(self.MJD - delta_time_epochs, dtype=np.float64) - - def to_datetime(self): - """ - Convert a ``timescale`` object to a ``datetime`` array - - Returns - ------- - dtime: np.ndarray - ``numpy.datetime64`` array - """ - # convert Modified Julian Day epoch to datetime variable - epoch = np.datetime64(datetime.datetime(*_mjd_epoch)) - # use nanoseconds to keep as much precision as possible - delta_time = np.atleast_1d(self.MJD*self.day*1e9).astype(np.int64) - # return the datetime array - return np.array(epoch + delta_time.astype('timedelta64[ns]')) - - @gravity_toolkit.utilities.reify - def ut1(self): - """Universal Time (UT) as Julian Days - """ - return self.MJD + 2400000.5 - - @gravity_toolkit.utilities.reify - def year(self): - """Universal Time (UT) as calendar year - """ - Y, M, D, h, m, s = convert_julian(self.ut1, format='tuple') - return convert_calendar_decimal(Y, M, D, hour=h, minute=m, second=s) - - @property - def dtype(self): - """Main data type of ``timescale`` object""" - return self.MJD.dtype - - @property - def shape(self): - """Dimensions of ``timescale`` object - """ - return np.shape(self.MJD) - - @property - def ndim(self): - """Number of dimensions in ``timescale`` object - """ - return np.ndim(self.MJD) - - def __len__(self): - """Number of time values - """ - return len(np.atleast_1d(self.MJD)) - - def __iter__(self): - """Iterate over time values - """ - self.__index__ = 0 - return self - - def __next__(self): - """Get the next time step - """ - temp = timescale() - try: - temp.MJD = np.atleast_1d(self.MJD)[self.__index__].copy() - except IndexError as exc: - raise StopIteration from exc - # add to index - self.__index__ += 1 - return temp - -# PURPOSE: Count number of leap seconds that have passed for each GPS time -def count_leap_seconds(GPS_Time, truncate=True): - """ - Counts the number of leap seconds between a given GPS time and UTC - - Parameters - ---------- - GPS_Time: np.ndarray or float - seconds since January 6, 1980 at 00:00:00 - truncate: bool, default True - Reduce list of leap seconds to positive GPS times - - Returns - ------- - n_leaps: float - number of elapsed leap seconds - """ - # get the valid leap seconds - leaps = get_leap_seconds(truncate=truncate) - # number of leap seconds prior to GPS_Time - n_leaps = np.zeros_like(GPS_Time,dtype=np.float64) - for i,leap in enumerate(leaps): - count = np.count_nonzero(GPS_Time >= leap) - if (count > 0): - indices = np.nonzero(GPS_Time >= leap) - n_leaps[indices] += 1.0 - # return the number of leap seconds for converting to UTC - return n_leaps - -# PURPOSE: Define GPS leap seconds -def get_leap_seconds(truncate=True): - """ - Gets a list of GPS times for when leap seconds occurred - - Parameters - ---------- - truncate: bool, default True - Reduce list of leap seconds to positive GPS times - - Returns - ------- - GPS time: float - GPS seconds when leap seconds occurred - """ - leap_secs = gravity_toolkit.utilities.get_data_path( - ['data','leap-seconds.list']) - # find line with file expiration as delta time - with leap_secs.open(mode='r', encoding='utf8') as fid: - secs, = [re.findall(r'\d+',i).pop() for i in fid.read().splitlines() - if re.match(r'^(?=#@)',i)] - # check that leap seconds file is still valid - expiry = datetime.datetime(*_ntp_epoch) + datetime.timedelta(seconds=int(secs)) - today = datetime.datetime.utcnow() - update_leap_seconds() if (expiry < today) else None - # get leap seconds - leap_UTC,TAI_UTC = np.loadtxt(leap_secs).T - # TAI time is ahead of GPS by 19 seconds - TAI_GPS = 19.0 - # convert leap second epochs from NTP to GPS - # convert from time of 2nd leap second to time of 1st leap second - leap_GPS = convert_delta_time(leap_UTC + TAI_UTC - TAI_GPS - 1, - epoch1=_ntp_epoch, epoch2=_gps_epoch) - # return the GPS times of leap second occurrence - if truncate: - return leap_GPS[leap_GPS >= 0].astype(np.float64) - else: - return leap_GPS.astype(np.float64) - -# PURPOSE: connects to servers and downloads leap second files -def update_leap_seconds(timeout=20, verbose=False, mode=0o775): - """ - Connects to servers to download leap-seconds.list files from NIST servers - - - https://www.nist.gov/pml/time-and-frequency-division/leap-seconds-faqs - - Servers and Mirrors - - - ftp://ftp.nist.gov/pub/time/leap-seconds.list - - https://www.ietf.org/timezones/data/leap-seconds.list - - Parameters - ---------- - timeout: int or None, default 20 - timeout in seconds for blocking operations - verbose: bool, default False - print file information about output file - mode: oct, default 0o775 - permissions mode of output file - """ - # local version of file - FILE = 'leap-seconds.list' - LOCAL = gravity_toolkit.utilities.get_data_path(['data',FILE]) - HASH = gravity_toolkit.utilities.get_hash(LOCAL) - - # try downloading from NIST ftp servers - HOST = ['ftp.nist.gov','pub','time',FILE] - try: - gravity_toolkit.utilities.check_ftp_connection(HOST[0]) - gravity_toolkit.utilities.from_ftp(HOST, - timeout=timeout, local=LOCAL, hash=HASH, - verbose=verbose, mode=mode) - except Exception as exc: - logging.debug(traceback.format_exc()) - pass - else: - return - - # try downloading from Internet Engineering Task Force (IETF) mirror - REMOTE = ['https://www.ietf.org','timezones','data',FILE] - try: - gravity_toolkit.utilities.from_http(REMOTE, - timeout=timeout, local=LOCAL, hash=HASH, - verbose=verbose, mode=mode) - except Exception as exc: - logging.debug(traceback.format_exc()) - pass - else: - return diff --git a/gravity_toolkit/time_series/__init__.py b/gravity_toolkit/time_series/__init__.py index 9702dad9..56df5574 100644 --- a/gravity_toolkit/time_series/__init__.py +++ b/gravity_toolkit/time_series/__init__.py @@ -1,5 +1,6 @@ from .amplitude import * from .fit import * +from .lomb_scargle import * from .piecewise import * from .regress import * from .savitzky_golay import * diff --git a/gravity_toolkit/time_series/lomb_scargle.py b/gravity_toolkit/time_series/lomb_scargle.py new file mode 100755 index 00000000..9a00e6cc --- /dev/null +++ b/gravity_toolkit/time_series/lomb_scargle.py @@ -0,0 +1,181 @@ +#!/usr/bin/env python +u""" +lomb_scargle.py +Written by Tyler Sutterley (06/2024) + +Wrapper function for computing Lomb-Scargle periodograms + +Probability is calculated using a calculation of the independent + frequencies from Horne and Baliunas (1986) + +INPUTS: + t_in: input time array + d_in: input data array + +OUTPUTS: + PowerDensity: normalized spectral power density + probability: probability of each frequency + frequency: considered frequencies array (omega/2pi) + period: considered periods array (1/f) + peak: period at peak power density + centroid: centroid of power density and period + +OPTIONS: + NORMALIZE: compute normalized periodogram + OMEGA: angular frequency range [min, max] + FREQUENCY: temporal frequency range [min, max] + PERIOD: range of periods [min, max] + N: number of frequencies + p: probability levels for contours + +PYTHON DEPENDENCIES: + numpy: Scientific Computing Tools For Python + https://numpy.org + https://numpy.org/doc/stable/user/numpy-for-matlab-users.html + scipy: Scientific Tools for Python + https://docs.scipy.org/doc/ + +REFERENCES: + N. R. Lomb, "Least-squares frequency analysis of unequally spaced data", + Astrophysics and Space Science, 39, 447--462, 1976. + https://doi.org/10.1007/BF00648343 + J. D. Scargle, "Studies in astronomical time series analysis. II - + Statistical aspects of spectral analysis of unevenly spaced data", + The Astrophysical Journal, 263, 835--853, 1982. + https://doi.org/10.1086/160554 + J. H. Horne and S. L. Baliunas, "A Prescription for Period + Analysis of Unevenly Sampled Time Series", + The Astrophysical Journal, 302, 757--763, 1986. + https://doi.org/10.1086/164037 + +UPDATE HISTORY: + Updated 06/2024: add function docstrings + Updated 01/2015: added centroid output + Written 08/2013 +""" +import numpy as np +import scipy.signal + +def lomb_scargle(t_in, d_in, **kwargs): + """ + Computes periodograms for least-squares spectral analysis following + [Lomb1976]_ [Scargle1982]_ and computes the frequency probabilities + following [Horne1986]_ + + Parameters + ---------- + t_in: float + input time array + d_in: floatto consider + input data array + NORMALIZE: bool, default False + Compute normalized periodogram + OMEGA: list, default [] + Angular frequency range + FREQUENCY: list, default [] + Temporal frequency range + PERIOD: list, default [] + Temporal period range + N: int, default 1000 + Number of frequencies + p: list, default [0.05, 0.01, 0.001, 1e-4, 1e-5, 1e-6] + Probability levels for contours + + Returns + ------- + PowerDensity: float + spectral power density + probability: float + probability of each frequency + frequency: float + considered frequencies array + period: float + periods array + peak: float + period at peak power density + centroid: float + centroid of power density and period + + References + ---------- + .. [Lomb1976] N. R. Lomb, "Least-squares frequency analysis of + unequally spaced data", *Astrophysics and Space Science*, + 39, 447--462, 1976. `doi: 10.1007/BF00648343 + `_ + .. [Scargle1982] J. D. Scargle, "Studies in astronomical time series + analysis. II - Statistical aspects of spectral analysis of + unevenly spaced data", *The Astrophysical Journal*, + 263, 835--853, 1982. `doi: 10.1086/160554 + `_ + .. [Horne1986] J. H. Horne and S. L. Baliunas, "A Prescription for + Period Analysis of Unevenly Sampled Time Series", + *The Astrophysical Journal*, 302, 757--763, 1986. + `doi: 10.1086/164037 `_ + """ + # default keyword arguments + kwargs.setdefault('NORMALIZE', False) + kwargs.setdefault('OMEGA', []) + kwargs.setdefault('FREQUENCY', []) + kwargs.setdefault('PERIOD', []) + kwargs.setdefault('N', 1000) + kwargs.setdefault('p', [0.05, 0.01, 0.001, 1e-4, 1e-5, 1e-6]) + + # remove singleton dimensions + t_in = np.squeeze(t_in) + d_in = np.squeeze(d_in) + + # number of independent measurements + nmax = np.count_nonzero(np.isfinite(d_in)) + nyquist = 1.0/(2.0*np.mean(t_in[1:] - t_in[0:-1])) + + # angular frequency range + if kwargs['OMEGA']: + OMEGA = np.atleast_1d(kwargs['OMEGA']) + elif kwargs['FREQUENCY']: + OMEGA = np.atleast_1d(kwargs['FREQUENCY'])/(2.0*np.pi) + elif kwargs['PERIOD']: + OMEGA = (2.0*np.pi)/np.atleast_1d(kwargs['PERIOD']) + else: + raise ValueError('Frequency range must be defined') + + # number of angular frequencies + N = int(kwargs['N']) + assert len(OMEGA) == 2, 'Angular frequency range must have 2 values' + # array of angular frequencies + angular_freq = np.linspace(OMEGA[0], OMEGA[1], N) + + # Estimate of number of independent frequencies in Lomb-Scargle + # analysis based on sample size. From Horne and Baliunas, + # "A Prescription for Period Analysis of Unevenly Sampled Time Series", + # The Astrophysical Journal, 392: 757-763, 1986. + independent_freq = np.round(-6.362 + 1.193*nmax + 0.00098*nmax**2) + # if less than 1 independent frequency: set equal to 1 + independent_freq = np.maximum(independent_freq, 1) + + # scaling the date (t[0] = 0) + t = t_in - t_in[0] + # periods and frequencies considered + frequency = angular_freq/(2.0*np.pi) + period = 2.0*np.pi/angular_freq + # scaling the data to be mean 0 with variance 1 + data_norm = (d_in - d_in.mean())/d_in.std() + # computing the lomb-scargle periodogram + # "normalized" spectral density refers to variance term in denominator + # PowerDensity has exponential probability distribution with unit mean + # can calculate normalized as described in Scipy reference + PowerDensity = scipy.signal.lombscargle(t, data_norm, angular_freq, + normalize=kwargs['NORMALIZE']) + # probability of frequencies (NULL test, significance of peak) + probability = 1.0 - (1.0-np.exp(-PowerDensity))**independent_freq + # probability contours + p = np.atleast_1d(kwargs['p']) + contour = -np.log(1.0 - (1 - p)**(1.0/independent_freq)) + # period at peak (maximum probability) + ipeak = np.argmax(PowerDensity) + peak = period[ipeak] + # period at signal centroid + centroid = np.sum(period*PowerDensity)/np.sum(PowerDensity) + + return {'PowerDensity':PowerDensity, 'Probability':probability, + 'frequency':frequency, 'period':period, 'contour':contour, + 'Nyquist':nyquist, 'peak':peak, 'centroid':centroid} diff --git a/gravity_toolkit/time_series/regress.py b/gravity_toolkit/time_series/regress.py index 0022683f..3795dba0 100755 --- a/gravity_toolkit/time_series/regress.py +++ b/gravity_toolkit/time_series/regress.py @@ -31,6 +31,7 @@ NRMSE: normalized root mean square error AIC: Akaike information criterion (Second-Order, AICc) BIC: Bayesian information criterion (Schwarz criterion) + LOGLIK: log likelihood model: modeled timeseries simple: modeled timeseries without oscillating components residual: model residual @@ -200,7 +201,7 @@ def regress(t_in, d_in, ORDER=1, CYCLES=[0.5,1.0], TERMS=[], # calculate epoch for calculating relative times if isinstance(RELATIVE, (list, np.ndarray)): t_rel = t_in[RELATIVE].mean() - elif isinstance(RELATIVE, (float, int, np.float_, np.int_)): + elif isinstance(RELATIVE, (float, int, np.float64, np.int_)): t_rel = np.copy(RELATIVE) elif (RELATIVE == Ellipsis): t_rel = t_in[RELATIVE].mean() diff --git a/gravity_toolkit/time_series/savitzky_golay.py b/gravity_toolkit/time_series/savitzky_golay.py index 4c509128..10fa3f44 100644 --- a/gravity_toolkit/time_series/savitzky_golay.py +++ b/gravity_toolkit/time_series/savitzky_golay.py @@ -40,6 +40,7 @@ polynomial of high order over an odd-sized window centered at the point. +REFERENCES: A. Savitzky, M. J. E. Golay, Smoothing and Differentiation of Data by Simplified Least Squares Procedures. Analytical Chemistry, 1964, 36 (8), pp 1627-1639. diff --git a/gravity_toolkit/tools.py b/gravity_toolkit/tools.py index 8a89160f..5794e586 100644 --- a/gravity_toolkit/tools.py +++ b/gravity_toolkit/tools.py @@ -11,7 +11,7 @@ scipy: Scientific Tools for Python https://docs.scipy.org/doc/ netCDF4: Python interface to the netCDF C library - https://unidata.github.io/netcdf4-python/netCDF4/index.html + https://unidata.github.io/netcdf4-python/netCDF4/index.html tkinter: Python interface to the Tcl/Tk GUI toolkit https://docs.python.org/3/library/tkinter.html ipywidgets: interactive HTML widgets for Jupyter notebooks and IPython @@ -28,6 +28,7 @@ UPDATE HISTORY: Updated 04/2024: add widget for setting endpoint for accessing PODAAC data + place colormap registration within try/except to check for existing Updated 05/2023: use pathlib to define and operate on paths Updated 03/2023: add wrap longitudes function to change convention improve typing for variables in docstrings @@ -1020,7 +1021,10 @@ def from_cpt(filename, use_extremes=True, **kwargs): if use_extremes: cmap = cmap.with_extremes(**extremes) # register colormap to be recognizable by cm.get_cmap() - cm.register_cmap(name=name, cmap=cmap) + try: + cm.register_cmap(name=name, cmap=cmap) + except: + pass # return the colormap return cmap @@ -1101,7 +1105,10 @@ def custom_colormap(N, map_name, **kwargs): # create colormap for use in matplotlib cmap = colors.LinearSegmentedColormap(map_name, cdict, **kwargs) # register colormap to be recognizable by cm.get_cmap() - cm.register_cmap(name=map_name, cmap=cmap) + try: + cm.register_cmap(name=map_name, cmap=cmap) + except: + pass # return the colormap return cmap diff --git a/gravity_toolkit/units.py b/gravity_toolkit/units.py index 9807f250..af1ac59a 100644 --- a/gravity_toolkit/units.py +++ b/gravity_toolkit/units.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" units.py -Written by Tyler Sutterley (09/2023) +Written by Tyler Sutterley (05/2024) Contributions by Hugo Lecomte Class for converting GRACE/GRACE-FO Level-2 data to specific units @@ -10,6 +10,7 @@ numpy: Scientific Computing Tools For Python (https://numpy.org) UPDATE HISTORY: + Updated 05/2024: make subscriptable and allow item assignment Updated 09/2023: added property for the approximate mass of the Earth Updated 03/2023: include option to not compensate for elastic deformation include option to include effects for Earth's oblateness @@ -308,4 +309,10 @@ def get_attributes(var: str) -> str: try: return named_attributes[var] except Exception as exc: - raise ValueError(f'Unknown units {var}') from exc \ No newline at end of file + raise ValueError(f'Unknown units {var}') from exc + + def __getitem__(self, key): + return getattr(self, key) + + def __setitem__(self, key, value): + setattr(self, key, value) diff --git a/gravity_toolkit/utilities.py b/gravity_toolkit/utilities.py index fb903450..a9e68010 100644 --- a/gravity_toolkit/utilities.py +++ b/gravity_toolkit/utilities.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" utilities.py -Written by Tyler Sutterley (04/2024) +Written by Tyler Sutterley (06/2024) Download and management utilities for syncing time and auxiliary files PYTHON DEPENDENCIES: @@ -9,6 +9,8 @@ https://pypi.python.org/pypi/lxml UPDATE HISTORY: + Updated 06/2024: added wrapper to importlib for optional dependencies + make default case for an import exception be a class Updated 04/2024: added argument for products in CMR shortname query Updated 11/2023: updated ssl context to fix deprecation error Updated 10/2023: add capability to download CSR LRI solutions @@ -71,6 +73,7 @@ import builtins import dateutil import warnings +import importlib import posixpath import lxml.etree import subprocess @@ -84,12 +87,6 @@ from urllib.parse import urlencode import urllib.request as urllib2 -# attempt imports -try: - import boto3 -except (AttributeError, ImportError, ModuleNotFoundError) as exc: - warnings.warn("boto3 not available", ImportWarning) - # PURPOSE: get absolute path within a package from a relative path def get_data_path(relpath: list | str | pathlib.Path): """ @@ -109,6 +106,48 @@ def get_data_path(relpath: list | str | pathlib.Path): elif isinstance(relpath, str): return filepath.joinpath(relpath) + +def import_dependency( + name: str, + extra: str = "", + raise_exception: bool = False + ): + """ + Import an optional dependency + + Adapted from ``pandas.compat._optional::import_optional_dependency`` + + Parameters + ---------- + name: str + Module name + extra: str, default "" + Additional text to include in the ``ImportError`` message + raise_exception: bool, default False + Raise an ``ImportError`` if the module is not found + + Returns + ------- + module: obj + Imported module + """ + # check if the module name is a string + msg = f"Invalid module name: '{name}'; must be a string" + assert isinstance(name, str), msg + # default error if module cannot be imported + err = f"Missing optional dependency '{name}'. {extra}" + module = type('module', (), {}) + # try to import the module + try: + module = importlib.import_module(name) + except (ImportError, ModuleNotFoundError) as exc: + if raise_exception: + raise ImportError(err) from exc + else: + logging.debug(err) + # return the module + return module + class reify(object): """Class decorator that puts the result of the method it decorates into the instance""" @@ -1154,6 +1193,7 @@ def s3_region(): region_name: str AWS region name """ + boto3 = import_dependency('boto3') region_name = boto3.session.Session().region_name return region_name @@ -1184,6 +1224,7 @@ def s3_client( response = urllib2.urlopen(request, timeout=timeout) cumulus = json.loads(response.read()) # get AWS client object + boto3 = import_dependency('boto3') client = boto3.client('s3', aws_access_key_id=cumulus['accessKeyId'], aws_secret_access_key=cumulus['secretAccessKey'], diff --git a/scripts/calc_degree_one.py b/scripts/calc_degree_one.py index 4301c13d..8d9861d8 100755 --- a/scripts/calc_degree_one.py +++ b/scripts/calc_degree_one.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" calc_degree_one.py -Written by Tyler Sutterley (10/2023) +Written by Tyler Sutterley (06/2024) Calculates degree 1 variations using GRACE coefficients of degree 2 and greater, and ocean bottom pressure variations from ECCO and OMCT/MPIOM @@ -167,6 +167,7 @@ https://doi.org/10.1029/2007JB005338 UPDATE HISTORY: + Updated 06/2024: use wrapper to importlib for optional dependencies Updated 10/2023: generalize mission variable to be GRACE/GRACE-FO Updated 09/2023: output comprehensive netCDF4 files with all components simplify I-matrix and G-matrix calculations @@ -257,7 +258,6 @@ import logging import pathlib import argparse -import warnings import traceback import collections import numpy as np @@ -265,17 +265,11 @@ import gravity_toolkit as gravtk # attempt imports -try: - import matplotlib.pyplot as plt - import matplotlib.cm as cm - import matplotlib.offsetbox - import matplotlib.ticker as ticker -except (AttributeError, ImportError, ModuleNotFoundError) as exc: - warnings.warn("matplotlib not available", ImportWarning) -try: - import netCDF4 -except (AttributeError, ImportError, ModuleNotFoundError) as exc: - warnings.warn("netCDF4 not available", ImportWarning) +plt = gravtk.utilities.import_dependency('matplotlib.pyplot') +cm = gravtk.utilities.import_dependency('matplotlib.cm') +offsetbox = gravtk.utilities.import_dependency('matplotlib.offsetbox') +ticker = gravtk.utilities.import_dependency('matplotlib.ticker') +netCDF4 = gravtk.utilities.import_dependency('netCDF4') # PURPOSE: keep track of threads def info(args): @@ -1064,7 +1058,7 @@ def calc_degree_one(base_dir, PROC, DREL, MODEL, LMAX, RAD, ax[i].set_ylabel('[mm]', fontsize=14) # add axis labels and adjust font sizes for axis ticks # axis label - artist = matplotlib.offsetbox.AnchoredText(key, pad=0., + artist = offsetbox.AnchoredText(key, pad=0., prop=dict(size=16,weight='bold'), frameon=False, loc=2) ax[i].add_artist(artist) # axes tick adjustments @@ -1115,7 +1109,7 @@ def calc_degree_one(base_dir, PROC, DREL, MODEL, LMAX, RAD, for i,key in enumerate(iteration_mmwe.fields): ax[i].set_ylabel('mm', fontsize=14) # axis label - artist = matplotlib.offsetbox.AnchoredText(key, pad=0., + artist = offsetbox.AnchoredText(key, pad=0., prop=dict(size=16,weight='bold'), frameon=False, loc=2) ax[i].add_artist(artist) # axes tick adjustments diff --git a/scripts/grace_raster_grids.py b/scripts/grace_raster_grids.py index b50b5772..781a6f70 100644 --- a/scripts/grace_raster_grids.py +++ b/scripts/grace_raster_grids.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" grace_raster_grids.py -Written by Tyler Sutterley (03/2024) +Written by Tyler Sutterley (06/2024) Reads in GRACE/GRACE-FO spherical harmonic coefficients and exports projected spatial fields @@ -147,6 +147,7 @@ utilities.py: download and management utilities for files UPDATE HISTORY: + Updated 06/2024: use wrapper to importlib for optional dependencies Updated 03/2024: increase mask buffer to twice the smoothing radius Written 08/2023 """ @@ -157,7 +158,6 @@ import time import logging import pathlib -import warnings import numpy as np import argparse import traceback @@ -165,14 +165,8 @@ import gravity_toolkit as gravtk # attempt imports -try: - import geoid_toolkit as geoidtk -except (AttributeError, ImportError, ModuleNotFoundError) as exc: - warnings.warn("geoid_toolkit not available", ImportWarning) -try: - import pyproj -except (AttributeError, ImportError, ModuleNotFoundError) as exc: - warnings.warn("pyproj not available", ImportWarning) +geoidtk = gravtk.utilities.import_dependency('geoid_toolkit') +pyproj = gravtk.utilities.import_dependency('pyproj') # PURPOSE: keep track of threads def info(args): @@ -430,7 +424,8 @@ def grace_raster_grids(base_dir, PROC, DREL, DSET, LMAX, RAD, # projection attributes attributes['crs'] = {} # add projection attributes - attributes['crs']['standard_name'] = crs1.to_cf()['grid_mapping_name'].title() + attributes['crs']['standard_name'] = \ + crs1.to_cf()['grid_mapping_name'].title() attributes['crs']['spatial_epsg'] = crs1.to_epsg() attributes['crs']['spatial_ref'] = crs1.to_wkt() attributes['crs']['proj4_params'] = crs1.to_proj4() diff --git a/scripts/monte_carlo_degree_one.py b/scripts/monte_carlo_degree_one.py index 4d4d8d1f..e801ecbc 100644 --- a/scripts/monte_carlo_degree_one.py +++ b/scripts/monte_carlo_degree_one.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" monte_carlo_degree_one.py -Written by Tyler Sutterley (10/2023) +Written by Tyler Sutterley (06/2024) Calculates degree 1 errors using GRACE coefficients of degree 2 and greater, and ocean bottom pressure variations from OMCT/MPIOM in a Monte Carlo scheme @@ -157,6 +157,7 @@ https://doi.org/10.1029/2005GL025305 UPDATE HISTORY: + Updated 06/2024: use wrapper to importlib for optional dependencies Updated 10/2023: generalize mission variable to be GRACE/GRACE-FO Updated 09/2023: add more root level attributes to output netCDF4 files simplify I-matrix and G-matrix calculations @@ -215,7 +216,6 @@ import logging import pathlib import argparse -import warnings import traceback import collections import numpy as np @@ -223,17 +223,11 @@ import gravity_toolkit as gravtk # attempt imports -try: - import matplotlib.pyplot as plt - import matplotlib.cm as cm - import matplotlib.offsetbox - from matplotlib.ticker import MultipleLocator -except (AttributeError, ImportError, ModuleNotFoundError) as exc: - warnings.warn("matplotlib not available", ImportWarning) -try: - import netCDF4 -except (AttributeError, ImportError, ModuleNotFoundError) as exc: - warnings.warn("netCDF4 not available", ImportWarning) +plt = gravtk.utilities.import_dependency('matplotlib.pyplot') +cm = gravtk.utilities.import_dependency('matplotlib.cm') +offsetbox = gravtk.utilities.import_dependency('matplotlib.offsetbox') +ticker = gravtk.utilities.import_dependency('matplotlib.ticker') +netCDF4 = gravtk.utilities.import_dependency('netCDF4') # PURPOSE: keep track of threads def info(args): @@ -901,12 +895,12 @@ def monte_carlo_degree_one(base_dir, PROC, DREL, LMAX, RAD, ax[2].set_ylabel('mm', fontsize=14) ax[2].set_xlabel('Grace Month', fontsize=14) ax[2].set_xlim(np.floor(months[0]/10.)*10.,np.ceil(months[-1]/10.)*10.) - ax[2].xaxis.set_minor_locator(MultipleLocator(5)) + ax[2].xaxis.set_minor_locator(ticker.MultipleLocator(5)) ax[2].xaxis.get_major_formatter().set_useOffset(False) # add axis labels and adjust font sizes for axis ticks for i,lbl in enumerate(['C10','C11','S11']): # axis label - artist = matplotlib.offsetbox.AnchoredText(lbl, pad=0.0, + artist = offsetbox.AnchoredText(lbl, pad=0.0, frameon=False, loc=2, prop=dict(size=16,weight='bold')) ax[i].add_artist(artist) # axes tick adjustments diff --git a/scripts/podaac_cumulus.py b/scripts/podaac_cumulus.py index 7dad5282..df050a8e 100644 --- a/scripts/podaac_cumulus.py +++ b/scripts/podaac_cumulus.py @@ -184,7 +184,7 @@ def podaac_cumulus(client, DIRECTORY, PROC=[], DREL=[], VERSION=[], if (ENDPOINT == 's3'): # get shortname for CMR query cmr_shortname, = gravtk.utilities.cmr_product_shortname( - mission=mi, center='GFZ', release=rl, level='L1B') + mission='grace', center='GFZ', release=rl, level='L1B') # attempt to list objects in s3 bucket try: objects = client.list_objects(Bucket=bucket, diff --git a/scripts/quick_mascon_plot.py b/scripts/quick_mascon_plot.py index 3763eaec..a2fc9336 100644 --- a/scripts/quick_mascon_plot.py +++ b/scripts/quick_mascon_plot.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" quick_mascon_plot.py -Written by Tyler Sutterley (08/2023) +Written by Tyler Sutterley (06/2024) Plots a mascon time series file for a particular format COMMAND LINE OPTIONS: @@ -35,6 +35,7 @@ time_series.regress.py: calculates trend coefficients using least-squares UPDATE HISTORY: + Updated 06/2024: use wrapper to importlib for optional dependencies Updated 08/2023: add option for changing drawing order of time series Updated 05/2023: use pathlib to define and operate on paths use fit module for getting tidal aliasing terms @@ -56,15 +57,11 @@ import sys import pathlib import argparse -import warnings import numpy as np import gravity_toolkit as gravtk # attempt imports -try: - import matplotlib.pyplot as plt -except ModuleNotFoundError: - warnings.warn("matplotlib not available", ImportWarning) +plt = gravtk.utilities.import_dependency('matplotlib.pyplot') # PURPOSE: read mascon time series file and create plot def run_plot(i, input_file, diff --git a/scripts/run_sea_level_equation.py b/scripts/run_sea_level_equation.py index 72cb8c41..61928e01 100644 --- a/scripts/run_sea_level_equation.py +++ b/scripts/run_sea_level_equation.py @@ -43,7 +43,7 @@ numpy: Scientific Computing Tools For Python https://numpy.org netCDF4: Python interface to the netCDF C library - https://unidata.github.io/netcdf4-python/netCDF4/index.html + https://unidata.github.io/netcdf4-python/netCDF4/index.html h5py: Pythonic interface to the HDF5 binary data format. http://www.h5py.org/ diff --git a/test/test_time.py b/test/test_time.py index ea78019b..13de3316 100644 --- a/test/test_time.py +++ b/test/test_time.py @@ -128,58 +128,6 @@ def test_parse_date_string(): assert np.all(epoch == [2000,1,1,18,0,0]) assert (to_secs == 0.0) -# PURPOSE: verify forward and backwards delta time conversions -@pytest.mark.parametrize("delta_time", np.random.randint(1,31536000,size=4)) -def test_delta_time(delta_time, gps_epoch=1198800018.0): - # convert to array if single value - delta_time = np.atleast_1d(delta_time) - # calculate gps time from delta_time - gps_seconds = gps_epoch + delta_time - time_leaps = gravtk.time.count_leap_seconds(gps_seconds) - # compare output delta times with original values - output_time = gravtk.time.convert_delta_time(gps_seconds - time_leaps, - epoch1=gravtk.time._gps_epoch, epoch2=(2018,1,1,0,0,0), - scale=1.0) - assert (delta_time == output_time) - # compare with original values using string epochs - output_time = gravtk.time.convert_delta_time(gps_seconds - time_leaps, - epoch1='1980-01-06T00:00:00', epoch2='2018-01-01T00:00:00', - scale=1.0) - assert (delta_time == output_time) - # calculate and compare with timescale values - GPS = gravtk.time.timescale().from_deltatime(gps_seconds, - epoch=gravtk.time._gps_epoch, standard='GPS') - output_time = GPS.to_deltatime((2018,1,1,0,0,0), scale=GPS.day) - assert np.isclose(delta_time, output_time) - -# PURPOSE: test timescale class conversions and constants -def test_timescale(): - # ATLAS Standard Data Epoch - atlas_sdp_epoch = np.datetime64('2018-01-01T00:00:00') - # from datetime - ATLAS = gravtk.time.timescale().from_datetime(atlas_sdp_epoch) - assert np.all(ATLAS.MJD == 58119) - delta_time_epochs = (ATLAS.to_datetime() - atlas_sdp_epoch) - assert np.all(delta_time_epochs/np.timedelta64(1, 'ns') == 0) - # from deltatime - ATLAS = gravtk.time.timescale().from_deltatime(0, epoch=(2018,1,1)) - assert np.all(ATLAS.MJD == 58119) - delta_time_epochs = (ATLAS.to_datetime() - atlas_sdp_epoch) - assert np.all(delta_time_epochs/np.timedelta64(1, 'ns') == 0) - # from MJD - ATLAS = gravtk.time.timescale(MJD=58119) - assert np.all(ATLAS.ut1 == 2458119.5) - assert np.all((ATLAS.MJD - 51544.5) == (ATLAS.ut1 - 2451545.0)) - # from MJD hourly array - delta_time = np.arange(0, 365, 1.0/24.0) - ATLAS = gravtk.time.timescale(MJD=58119 + delta_time) - delta_time_epochs = (ATLAS.to_datetime() - atlas_sdp_epoch) - assert np.allclose(delta_time_epochs/np.timedelta64(1, 'D'), delta_time) - assert np.allclose(ATLAS.to_deltatime(epoch=atlas_sdp_epoch), delta_time) - assert np.allclose(ATLAS.to_deltatime(epoch='2018-01-01'), delta_time) - # check constants - assert (ATLAS.day == 86400.0) - # PURPOSE: test months adjustment for special cases # parameterize calendar dates @pytest.mark.parametrize("PROC", ['CSR','GFZ','GSFC','JPL'])