diff --git a/doc/source/getting_started/Citations.rst b/doc/source/getting_started/Citations.rst index 50c8c8b..7a0b4cc 100644 --- a/doc/source/getting_started/Citations.rst +++ b/doc/source/getting_started/Citations.rst @@ -49,7 +49,7 @@ Problem Reports ############### If you have found a problem in these programs, or you would like to suggest an improvement or modification, please submit a -`GitHub issue `_ +`GitHub issue `_ and we will get back to you. Dependencies diff --git a/doc/source/getting_started/Install.rst b/doc/source/getting_started/Install.rst index 8bda705..4fa5bc9 100644 --- a/doc/source/getting_started/Install.rst +++ b/doc/source/getting_started/Install.rst @@ -3,28 +3,35 @@ Installation ============ ``gravity-toolkit`` is available for download from the `GitHub repository `_, -and the `Python Package Index (pypi) `_, -The contents of the repository can be downloaded as a `zipped file `_ or cloned. +the `Python Package Index (pypi) `_, +and from `conda-forge `_. -To use this repository, please fork into your own account and then clone onto your system: + +The simplest installation for most users will likely be using ``conda`` or ``mamba``: .. code-block:: bash - git clone https://github.com/tsutterley/gravity-toolkit.git + conda install -c conda-forge gravity-toolkit -Can then install using ``setuptools``: +``conda`` installed versions of ``gravity-toolkit`` can be upgraded to the latest stable release: .. code-block:: bash - python3 setup.py install + conda update gravity-toolkit + +To use the development repository, please fork ``gravity-toolkit`` into your own account and then clone onto your system: + +.. code-block:: bash + + git clone https://github.com/tsutterley/gravity-toolkit.git -or ``pip`` +``gravity-toolkit`` can then be installed within the package directory using ``pip``: .. code-block:: bash python3 -m pip install --user . -Alternatively can install the ``gravity_toolkit`` utilities directly from GitHub with ``pip``: +The development version of ``gravity-toolkit`` can also be installed directly from GitHub using ``pip``: .. code-block:: bash diff --git a/gravity_toolkit/gen_stokes.py b/gravity_toolkit/gen_stokes.py index 5856ee3..6b43f3f 100755 --- a/gravity_toolkit/gen_stokes.py +++ b/gravity_toolkit/gen_stokes.py @@ -43,6 +43,7 @@ and filters the GRACE/GRACE-FO coefficients for striping errors UPDATE HISTORY: + Updated 06/2025: copy latitude and longitude as float64 for numpy 2.0 stability Updated 04/2023: allow love numbers to be None for custom units case Updated 03/2023: improve typing for variables in docstrings Updated 02/2023: set custom units as top option in if/else statements @@ -136,11 +137,13 @@ def gen_stokes(data, lon, lat, LMIN=0, LMAX=60, MMAX=None, UNITS=1, # colatitude degree spacing in radians dth = dlat*np.pi/180.0 + # convert latitude and longitude to float if integers + lon = lon.astype(np.float64) + lat = lat.astype(np.float64) # reformatting longitudes to range 0:360 (if previously -180:180) lon = np.squeeze(lon.copy()) if np.any(lon < 0): - lon_ind, = np.nonzero(lon < 0) - lon[lon_ind] += 360.0 + lon[lon < 0] += 360.0 # Longitude in radians phi = lon[np.newaxis,:]*np.pi/180.0 # Colatitude in radians diff --git a/gravity_toolkit/sea_level_equation.py b/gravity_toolkit/sea_level_equation.py index 002ed8c..b0ec7b3 100644 --- a/gravity_toolkit/sea_level_equation.py +++ b/gravity_toolkit/sea_level_equation.py @@ -1,6 +1,6 @@ #!/usr/bin/env python u""" -sea_level_equation.py (03/2023) +sea_level_equation.py (06/2025) Solves the sea level equation with the option of including polar motion feedback Uses a Clenshaw summation to calculate the spherical harmonic summation @@ -32,6 +32,7 @@ 2: Munk and MacDonald (1960) fluid love number 3: Lambeck (1980) fluid love number list or tuple: custom value (klf) + DENSITY: Density of water [g/cm^3] POLAR: Include polar feedback ITERATIONS: maximum number of iterations for the solver PLM: input Legendre polynomials @@ -90,6 +91,7 @@ https://doi.org/10.1029/JB090iB11p09363 UPDATE HISTORY: + Updated 06/2025: added option to set the density of sea water (g/cm^3) Updated 03/2023: improve typing for variables in docstrings Updated 01/2023: refactored associated legendre polynomials Updated 11/2022: use f-strings for formatting verbose or ascii output @@ -128,8 +130,9 @@ # PURPOSE: Computes Sea Level Fingerprints including polar motion feedback def sea_level_equation(loadClm, loadSlm, glon, glat, land_function, LMAX=0, - LOVE=None, BODY_TIDE_LOVE=0, FLUID_LOVE=0, POLAR=True, ITERATIONS=6, - PLM=None, FILL_VALUE=0, ASTYPE=np.longdouble, SCALE=1e-280, **kwargs): + LOVE=None, BODY_TIDE_LOVE=0, FLUID_LOVE=0, DENSITY=1.0, POLAR=True, + ITERATIONS=6, PLM=None, FILL_VALUE=0, ASTYPE=np.longdouble, SCALE=1e-280, + **kwargs): """ Solves the sea level equation with the option of including polar motion feedback :cite:p:`Farrell:1976hm,Kendall:2005ds,Mitrovica:2003cq` @@ -167,6 +170,8 @@ def sea_level_equation(loadClm, loadSlm, glon, glat, land_function, LMAX=0, - ``2``: :cite:p:`Munk:1960uk` fluid love number - ``3``: :cite:p:`Lambeck:1980um` fluid love number - list or tuple: custom value ``(klf)`` + DENSITY: float, default 1.0 + Density of water [g/cm\ :sup:`3`] POLAR: bool, default True Include polar feedback ITERATIONS: int, default 6 @@ -199,7 +204,7 @@ def sea_level_equation(loadClm, loadSlm, glon, glat, land_function, LMAX=0, # extract arrays of kl, hl, and ll Love Numbers hl,kl,ll = LOVE # density of water [g/cm^3] - rho_water = 1.0 + rho_water = np.float64(DENSITY) # Earth Parameters factors = units(lmax=LMAX) # Average Density of the Earth [g/cm^3] diff --git a/requirements.txt b/requirements.txt index 86591e4..dab10e5 100644 --- a/requirements.txt +++ b/requirements.txt @@ -2,6 +2,7 @@ boto3 future lxml matplotlib +netCDF4 numpy python-dateutil pyyaml diff --git a/scripts/convert_harmonics.py b/scripts/convert_harmonics.py index bc9b97a..482ade8 100644 --- a/scripts/convert_harmonics.py +++ b/scripts/convert_harmonics.py @@ -251,14 +251,17 @@ def arguments(): parser.add_argument('--reference', type=str.upper, default='CF', choices=['CF','CM','CE'], help='Reference frame for load Love numbers') - # output units + # input units + # 1: cm of water thickness (cmwe) + # 2: Gigatonnes (Gt) + # 3: mm of water thickness kg/m^2 parser.add_argument('--units','-U', type=int, default=1, choices=[1,2,3], - help='Output units') + help='Input units of spatial fields') # output grid parameters parser.add_argument('--spacing','-S', type=float, nargs='+', default=[0.5,0.5], metavar=('dlon','dlat'), - help='Spatial resolution of output data') + help='Spatial resolution of input data') parser.add_argument('--interval','-I', type=int, default=2, choices=[1,2], help='Input grid interval (1: global, 2: centered global)') diff --git a/scripts/run_sea_level_equation.py b/scripts/run_sea_level_equation.py index 61928e0..e4a0e27 100644 --- a/scripts/run_sea_level_equation.py +++ b/scripts/run_sea_level_equation.py @@ -1,6 +1,6 @@ #!/usr/bin/env python u""" -run_sea_level_equation.py (05/2023) +run_sea_level_equation.py (06/2025) Solves the sea level equation with the option of including polar motion feedback Uses a Clenshaw summation to calculate the spherical harmonic summation @@ -10,7 +10,7 @@ --verbose --mode 0o775 input_file output_file INPUTS: - input_file: input load file + input_file: input load file (harmonics or spatial field) output_file: output sea level fingerprints file COMMAND LINE OPTIONS: @@ -28,6 +28,7 @@ 1: Munk and MacDonald (1960) secular love number 2: Munk and MacDonald (1960) fluid love number 3: Lambeck (1980) fluid love number + -d X, --density X: Density of water in g/cm^3 --polar-feedback: Include polar feedback --reference X: Reference frame for load love numbers -I X, --iterations X: maximum number of iterations for the solver @@ -35,7 +36,14 @@ ascii netCDF4 HDF5 + -T X, --input-type X: Input data type for load files + harmonics: spherical harmonic coefficients + spatial: spatial fields -D, --date: input and output files have date information + -U X, --units X: input units for spatial files + 1: cm of water thickness (cmwe) + 2: Gigatonnes (Gt) + 3: mm of water thickness kg/m^2 -V, --verbose: verbose output of processing run -M X, --mode X: permissions mode of the output files @@ -68,6 +76,9 @@ Bollettino di Geodesia e Scienze (1982) UPDATE HISTORY: + Updated 06/2025: added options to run from input spatial fields + added attributes for lineage to track input files + added option to set the density of water in g/cm^3 Updated 05/2023: use pathlib to define and operate on paths Updated 03/2023: add root attributes to output netCDF4 and HDF5 files Updated 02/2023: use love numbers class with additional attributes @@ -111,6 +122,7 @@ import sys import os import re +import copy import logging import pathlib import argparse @@ -135,11 +147,14 @@ def run_sea_level_equation(INPUT_FILE, OUTPUT_FILE, LOVE_NUMBERS=0, BODY_TIDE_LOVE=0, FLUID_LOVE=0, + DENSITY=1.0, REFERENCE=None, ITERATIONS=0, POLAR=False, DATAFORM=None, + INPUT_TYPE=None, DATE=False, + UNITS=None, MODE=0o775): # set default paths @@ -192,14 +207,49 @@ def run_sea_level_equation(INPUT_FILE, OUTPUT_FILE, if POLAR: attributes['polar_motion_feedback'] = 'Kendall et al. (2005)' - # read spherical harmonic coefficients from input file format - load_Ylms = gravtk.harmonics().from_file(INPUT_FILE, - format=DATAFORM, date=DATE) - # truncate harmonics to degree and order LMAX - load_Ylms.truncate(lmax=LMAX, mmax=LMAX) - # expand dimensions to iterate over slices - load_Ylms.expand_dims() - l1,m1,nt = load_Ylms.shape + # read input spherical harmonic coefficients from file + single_file_formats = ('ascii', 'netCDF4', 'HDF5') + index_file_formats = ('index-ascii', 'index-netCDF4', 'index-HDF5') + if DATAFORM in single_file_formats and (INPUT_TYPE == 'spatial'): + # read spatial data from input file format + dataform = copy.copy(DATAFORM) + load_spatial = gravtk.spatial().from_file(INPUT_FILE, + format=DATAFORM, date=DATE) + attributes['lineage'] = load_spatial.filename.name + elif DATAFORM in index_file_formats and (INPUT_TYPE == 'spatial'): + # read spatial data from index file + _,dataform = DATAFORM.split('-') + load_spatial = gravtk.spatial().from_index(INPUT_FILE, + format=dataform, date=DATE) + attributes['lineage'] = [f.name for f in load_spatial.filename] + elif DATAFORM in single_file_formats: + dataform = copy.copy(DATAFORM) + # read spherical harmonic coefficients from input file format + load_Ylms = gravtk.harmonics().from_file(INPUT_FILE, + format=DATAFORM, date=DATE) + attributes['lineage'] = load_Ylms.filename.name + elif DATAFORM in index_file_formats: + # read spherical harmonic coefficients from index file + _,dataform = DATAFORM.split('-') + load_Ylms = gravtk.harmonics().from_index(INPUT_FILE, + format=dataform, date=DATE) + attributes['lineage'] = [f.name for f in load_Ylms.filename] + else: + raise ValueError(f'Unknown input data format {DATAFORM:s} for {INPUT_TYPE:s}') + + # convert input data to be iterable over time slices + if (INPUT_TYPE == 'spatial'): + # expand dimensions to iterate over slices + load_spatial.expand_dims() + # number of time slices + nt = load_spatial.shape[2] + else: + # truncate harmonics to degree and order LMAX + load_Ylms.truncate(lmax=LMAX, mmax=LMAX) + # expand dimensions to iterate over slices + load_Ylms.expand_dims() + # number of time slices + nt = load_Ylms.shape[2] # calculate the legendre functions using Holmes and Featherstone relation PLM, dPLM = gravtk.plm_holmes(LMAX, np.cos(th)) @@ -211,20 +261,35 @@ def run_sea_level_equation(INPUT_FILE, OUTPUT_FILE, for i in range(nt): # print iteration if running a series if (nt > 1): - logging.info('Index {0:d} of {1:d}'.format(i+1,nt)) - # subset harmonics to indice - Ylms = load_Ylms.index(i, date=DATE) + logging.info(f'Index {i+1:d} of {nt:d}') + # subset harmonics/spatial fields to indice + if (INPUT_TYPE == 'spatial'): + spatial_data = load_spatial.index(i, date=DATE) + # convert missing values to zero + spatial_data.replace_invalid(0.0) + # convert spatial field to spherical harmonics + Ylms = gravtk.gen_stokes(spatial_data.data.T, + spatial_data.lon, spatial_data.lat, UNITS=UNITS, + LMIN=0, LMAX=LMAX, LOVE=LOVE) + else: + Ylms = load_Ylms.index(i, date=DATE) # run pseudo-spectral sea level equation solver sea_level.data[:,:,i] = gravtk.sea_level_equation(Ylms.clm, Ylms.slm, landsea.lon, landsea.lat, land_function.T, LMAX=LMAX, LOVE=LOVE, BODY_TIDE_LOVE=BODY_TIDE_LOVE, - FLUID_LOVE=FLUID_LOVE, POLAR=POLAR, PLM=PLM, - ITERATIONS=ITERATIONS, FILL_VALUE=0).T + FLUID_LOVE=FLUID_LOVE, DENSITY=DENSITY, POLAR=POLAR, + PLM=PLM, ITERATIONS=ITERATIONS, FILL_VALUE=0).T sea_level.mask[:,:,i] = (sea_level.data[:,:,i] == 0) # copy dimensions sea_level.lon = np.copy(landsea.lon) sea_level.lat = np.copy(landsea.lat) - sea_level.time = np.copy(load_Ylms.time) if DATE else None + # copy date variables + if (INPUT_TYPE == 'spatial') and DATE: + # copy time from spatial data + sea_level.time = np.copy(load_spatial.time) + elif DATE: + # copy time from load Ylms + sea_level.time = np.copy(load_Ylms.time) # remove singleton dimensions if necessary sea_level.squeeze() # add attributes to output spatial field @@ -236,16 +301,16 @@ def run_sea_level_equation(INPUT_FILE, OUTPUT_FILE, kwargs['units'] = 'centimeters' kwargs['longname'] = 'Equivalent_Water_Thickness' # save as output DATAFORM - if (DATAFORM == 'ascii'): + if (dataform == 'ascii'): # ascii (.txt) # only print ocean points sea_level.fill_value = 0 sea_level.update_mask() sea_level.to_ascii(OUTPUT_FILE, date=DATE) - elif (DATAFORM == 'netCDF4'): + elif (dataform == 'netCDF4'): # netCDF4 (.nc) sea_level.to_netCDF4(OUTPUT_FILE, date=DATE, **kwargs) - elif (DATAFORM == 'HDF5'): + elif (dataform == 'HDF5'): # HDF5 (.H5) sea_level.to_HDF5(OUTPUT_FILE, date=DATE, **kwargs) # set the permissions mode of the output file @@ -292,6 +357,10 @@ def arguments(): parser.add_argument('--body','-b', type=int, default=0, choices=[0,1], help='Treatment of the body tide Love number') + # density of water in g/cm^3 + parser.add_argument('--density','-d', + type=float, default=1.0, + help='Density of water in g/cm^3') # different treatments of the fluid Love number of gravitational potential # 0: Han and Wahr (1989) fluid love number # 1: Munk and MacDonald (1960) secular love number @@ -315,13 +384,27 @@ def arguments(): type=str.upper, default='CF', choices=['CF','CM','CE'], help='Reference frame for load Love numbers') # input and output data format (ascii, netCDF4, HDF5) + choices = [] + choices.extend(['ascii','netCDF4','HDF5']) + choices.extend(['index-ascii','index-netCDF4','index-HDF5']) parser.add_argument('--format','-F', - type=str, default='netCDF4', choices=['ascii','netCDF4','HDF5'], + type=str, default='netCDF4', choices=choices, help='Input and output data format') + # define the input data type for the load files + parser.add_argument('--input-type','-T', + type=str, default='harmonics', choices=['harmonics','spatial'], + help='Input data type for load fields') # Input and output files have date information parser.add_argument('--date','-D', default=False, action='store_true', help='Input and output files have date information') + # input units + # 1: cm of water thickness (cmwe) + # 2: Gigatonnes (Gt) + # 3: mm of water thickness kg/m^2 + parser.add_argument('--units','-U', + type=int, default=1, choices=[1,2,3], + help='Input units of spatial fields') # print information about processing run parser.add_argument('--verbose','-V', action='count', default=0, @@ -353,11 +436,14 @@ def main(): LOVE_NUMBERS=args.love, BODY_TIDE_LOVE=args.body, FLUID_LOVE=args.fluid, + DENSITY=args.density, REFERENCE=args.reference, ITERATIONS=args.iterations, POLAR=args.polar_feedback, DATAFORM=args.format, + INPUT_TYPE=args.input_type, DATE=args.date, + UNITS=args.units, MODE=args.mode) except Exception as exc: # if there has been an error exception