feat: update geocenter plot programs for AGU (#60)

feat: added netCDF4 reader for UCI iteration files
feat: added custom colormap function for some common scales
feat: added UNITS list option for converting from custom units

doc/source/user_guide/clenshaw_summation.rst

@@ -31,12 +31,13 @@ Keyword arguments
 - ``RAD``: Gaussian smoothing radius (km)
 - ``UNITS``: output data units
 
-   * ``1`` cm of water thickness
-   * ``2`` mm of geoid height
-   * ``3`` mm of elastic crustal deformation [Davis2004]_
-   * ``4`` microGal gravitational perturbation
-   * ``5`` millibar equivalent surface pressure
-   * ``6`` cm of viscoelastic crustal uplift (GIA) [Wahr2000]_
+   * ``1``: cm of water thickness
+   * ``2``: mm of geoid height
+   * ``3``: mm of elastic crustal deformation [Davis2004]_
+   * ``4``: microGal gravitational perturbation
+   * ``5``: millibar equivalent surface pressure
+   * ``6``: cm of viscoelastic crustal uplift (GIA) [Wahr2000]_
+   * ``list``: custom degree-dependent unit conversion factor
 - ``LMAX``: Upper bound of Spherical Harmonic Degrees
 - ``LOVE``: input load Love numbers up to degree of truncation (``hl``, ``kl``, ``ll``)
 - ``ASTYPE``: floating point precision for calculating Clenshaw summation

doc/source/user_guide/gen_disc_load.md

@@ -21,6 +21,11 @@ Ylms = gen_disc_load(data, lon, lat, area, LMAX=LMAX, PLM=PLM, LOVE=(hl,kl,ll))
 #### Keyword arguments
 - `LMAX`:  maximum spherical harmonic degree of the output harmonics
 - `MMAX`: maximum spherical harmonic order of the output harmonics
+- `UNITS`: input data units
+   * `1` cm water equivalent thickness (cm w.e., g/cm<sup>2</sup>)
+   * `2` gigatonnes of mass (Gt)
+   * `3` mm water equivalent thickness (mm w.e., kg/m<sup>2</sup>)
+   * `list`: custom unit conversion factor
 - `PLM`: input Legendre polynomials for `cos(theta)` (disc center)
 - `LOVE`: input load Love numbers up to degree of truncation
 

doc/source/user_guide/gen_harmonics.md

@@ -24,7 +24,7 @@ Ylms = gen_harmonics(data, lon, lat, LMAX=LMAX, PLM=PLM)
 - `PLM`: input fully normalized associated Legendre polynomials or Fourier coefficients of Legendre polynomials
 - `METHOD`: conversion method for calculating harmonics
     * `'integration'`: for global grids
-    *  `'fourier'`: for regional or global grids
+    * `'fourier'`: for regional or global grids
 
 #### Returns
 - `clm`: Cosine spherical harmonic coefficients (4-pi normalized)

doc/source/user_guide/gen_point_load.md

@@ -19,6 +19,7 @@ Ylms = gen_point_load(data, lon, lat, UNITS=1, LMAX=LMAX, LOVE=(hl,kl,ll))
 - `UNITS`: input data units
    * `1`: grams of mass (default)
    * `2`: gigatonnes of mass (Gt)
+   * `list`: custom degree-dependent unit conversion factor
 - `LMAX`:  maximum spherical harmonic degree of the output harmonics
 - `MMAX`: maximum spherical harmonic order of the output harmonics
 - `LOVE`: input load Love numbers up to degree of truncation

doc/source/user_guide/gen_spherical_cap.md

@@ -25,8 +25,9 @@ Ylms = gen_spherical_cap(data, lon, lat, UNITS=1, LMAX=LMAX, PLM=PLM, LOVE=(hl,k
 - `AREA`: spherical cap area (cm<sup>2</sup>)
 - `UNITS`: input data units
    * `1`: cm water equivalent thickness (cm w.e., g/cm<sup>2</sup>)
-   * `2` gigatonnes of mass (Gt)
-   * `3` mm water equivalent thickness (mm w.e., kg/m<sup>2</sup>)
+   * `2`: gigatonnes of mass (Gt)
+   * `3`: mm water equivalent thickness (mm w.e., kg/m<sup>2</sup>)
+   * `list`: custom unit conversion factor
 - `PLM`: input Legendre polynomials for spherical cap centers
 - `LOVE`: input load Love numbers up to degree of truncation
 

doc/source/user_guide/gen_stokes.md

@@ -19,13 +19,14 @@ Ylms = gen_stokes(data, lon, lat, UNITS=1, LMAX=LMAX, PLM=PLM, LOVE=(hl,kl,ll))
 
 #### Keyword arguments
 - `UNITS`: input data units
-   * `1` cm water equivalent thickness (cm w.e., g/cm<sup>2</sup>)
-   * `2` gigatonnes of mass (Gt)
-   * `3` mm water equivalent thickness (mm w.e., kg/m<sup>2</sup>)
+   * `1`: cm water equivalent thickness (cm w.e., g/cm<sup>2</sup>)
+   * `2`: gigatonnes of mass (Gt)
+   * `3`: mm water equivalent thickness (mm w.e., kg/m<sup>2</sup>)
+   * `list`: custom degree-dependent unit conversion factor
 - `LMIN`: minimum spherical harmonic degree of the output harmonics
 - `LMAX`:  maximum spherical harmonic degree of the output harmonics
 - `MMAX`: maximum spherical harmonic order of the output harmonics
-- `PLM`: input Legendre polynomials (for improving computational time)
+- `PLM`: input Legendre polynomials
 - `LOVE`: input load Love numbers up to degree of truncation
 
 #### Returns

doc/source/user_guide/geocenter.rst

@@ -64,7 +64,7 @@ General Attributes and Methods
             ``release``: GRACE/GRACE-FO/Swarm data release for dealiasing product
 
             ``calendar_year``: calendar year of data
-            
+
             ``calendar_month``: calendar month of data
 
 
@@ -133,17 +133,17 @@ General Attributes and Methods
     .. method:: object.from_tellus(geocenter_file, **kwargs)
 
         - Reads monthly geocenter spherical harmonic data files from `GRACE Tellus Technical Notes <https://podaac-tools.jpl.nasa.gov/drive/files/allData/tellus/L2/degree_1>`_
-  
+
         - Estimates calculated using GRACE measurements and Ocean Models of Degree 1 [Swenson2008]_ [Sun2016]_
 
         Arguments:
 
             ``geocenter_file``: degree 1 file
-            
+
                 * CSR: ``TN-13_GEOC_CSR_RL06.txt``
 
                 * GFZ: ``TN-13_GEOC_GFZ_RL06.txt``
-                
+
                 * JPL: ``TN-13_GEOC_JPL_RL06.txt``
 
         Keyword arguments:
@@ -153,6 +153,19 @@ General Attributes and Methods
             ``JPL``: use JPL TN-13 geocenter files calculated following [Sun2016]_
 
 
+    .. method:: object.from_netCDF4(geocenter_file, **kwargs)
+
+        - Reads geocenter file and extracts dates and spherical harmonic data from a netCDF4 file [Sutterley2019]_
+
+        Arguments:
+
+            ``geocenter_file``: degree 1 netCDF4 file
+
+        Keyword arguments:
+
+            ``compression``: netCDF4 file is compressed or streaming as bytes
+
+
     .. method:: object.copy(**kwargs)
 
         Copy a geocenter object to a new geocenter object
@@ -228,7 +241,7 @@ General Attributes and Methods
 
         Keyword arguments:
 
-            ``kl``: gravitational load love number of degree 1     
+            ``kl``: gravitational load love number of degree 1
 
 
     .. method:: object.to_mmwe(kl=0.0):
@@ -256,8 +269,8 @@ General Attributes and Methods
         Keyword arguments:
 
             ``kl``: gravitational load love number of degree 1
-       
-        
+
+
     .. method:: object.from_mmwe(kl=0.0):
 
         Normalizes spherical harmonics from millimeters water equivalent (mmwe)
@@ -291,15 +304,15 @@ General Attributes and Methods
         Subtract one geocenter object from another
 
         Arguments:
-        
+
             ``temp``: geocenter object to be subtracted
 
     .. method:: object.multiply(self, temp):
 
         Multiply two geocenter objects
 
         Arguments:
-        
+
             ``temp``: geocenter object to be multiplied
 
 
@@ -308,7 +321,7 @@ General Attributes and Methods
         Divide one geocenter object from another
 
         Arguments:
-        
+
             ``temp``: geocenter object to be divided
 
 
@@ -317,7 +330,7 @@ General Attributes and Methods
         Multiply a geocenter object by a constant
 
         Arguments:
-        
+
             ``var``: scalar value to which the geocenter object will be multiplied
 
 
@@ -326,7 +339,7 @@ General Attributes and Methods
         Raise a geocenter object to a power
 
         Arguments:
-        
+
             ``power``: power to which the geocenter object will be raised
 
 

doc/source/user_guide/tools.rst

@@ -164,12 +164,34 @@ General Attributes and Methods
 
     Reads GMT color palette table files and registers the colormap to be recognizable by ``plt.cm.get_cmap()``
 
+.. function:: custom_colormap(N, map_name)
+
+    Calculates a custom colormap and registers it to be recognizable by ``plt.cm.get_cmap()``
+
+    Arguments:
+
+        - ``N``: number of slices in initial HSV color map
+
+        - ``map_name``: name of color map
+
+            * ``'Joughin'``: [Joughin2018]_ standard velocity colormap
+
+            * ``'Rignot'``: [Rignot2011]_ standard velocity colormap
+
+            * ``'Seroussi'``:  [Seroussi2011]_ velocity divergence colormap
+
 
 References
 ##########
 
 .. [Fagiolini2015] E. Fagiolini, F. Flechtner, M. Horwath, and H. Dobslaw, "Correction of inconsistencies in ECMWF's operational analysis data during de-aliasing of GRACE gravity models", *Geophysical Journal International*, 202(3), 2150--2158, (2015). `doi: 10.1093/gji/ggv276 <https://doi.org/10.1093/gji/ggv276>`_
 
+.. [Joughin2018] I. Joughin, B. E. Smith, and I. Howat, "Greenland Ice Mapping Project: ice flow velocity variation at sub-monthly to decadal timescales", *The Cryosphere*, 12, 2211--2227, (2018). `doi: 10.5194/tc-12-2211-2018 <https://doi.org/10.5194/tc-12-2211-2018>`_
+
+.. [Rignot2011] E. Rignot J. Mouginot, and B. Scheuchl, "Ice Flow of the Antarctic Ice Sheet", *Science*, 333(6048), 1427--1430, (2011). `doi: 10.1126/science.1208336 <https://doi.org/10.1126/science.1208336>`_
+
+.. [Seroussi2011] H. Seroussi, M. Morlighem, E. Rignot, E. Larour, D. Aubry, H. Ben Dhia, and S. S. Kristensen, "Ice flux divergence anomalies on 79north Glacier, Greenland", *Geophysical Research Letters*, 38(L09501), (2011). `doi: 10.1029/2011GL047338 <https://doi.org/10.1029/2011GL047338>`_
+
 .. [Swenson2006] S. Swenson and J. Wahr, "Post‐processing removal of correlated errors in GRACE data", *Geophysical Research Letters*, 33(L08402), (2006). `doi: 10.1029/2005GL025285 <https://doi.org/10.1029/2005GL025285>`_
 
 .. [Wahr2015] J. Wahr, R. S. Nerem, and S. V. Bettadpur, "The pole tide and its effect on GRACE time‐variable gravity measurements: Implications for estimates of surface mass variations". *Journal of Geophysical Research: Solid Earth*, 120(6), 4597--4615, (2015). `doi: 10.1002/2015JB011986 <https://doi.org/10.1002/2015JB011986>`_

gravity_toolkit/clenshaw_summation.py

@@ -1,7 +1,7 @@
 #!/usr/bin/env python
 u"""
 clenshaw_summation.py
-Written by Tyler Sutterley (06/2021)
+Written by Tyler Sutterley (11/2021)
 Calculates the spatial field for a series of spherical harmonics for a
     sequence of ungridded points
 
@@ -24,6 +24,7 @@
         4: microGal gravitational perturbation
         5: mbar equivalent surface pressure
         6: cm of viscoelastic crustal uplift (GIA) [See Wahr 1995 or Wahr 2000]
+        list: custom degree-dependent unit conversion factor
     LMAX: Upper bound of Spherical Harmonic Degrees
     LOVE: input load Love numbers up to degree LMAX (hl,kl,ll)
     ASTYPE: floating point precision for calculating Clenshaw summation
@@ -48,6 +49,7 @@
         Bollettino di Geodesia e Scienze (1982)
 
 UPDATE HISTORY:
+    Updated 11/2021: added UNITS list option for converting to custom units
     Updated 09/2021: fix passing SCALE keyword argument to clenshaw_s_m
     Updated 06/2021: output equivalent pressure in pascals
     Updated 08/2020: parameterize float precision to improve computational time
@@ -85,6 +87,7 @@ def clenshaw_summation(clm, slm, lon, lat, RAD=0, UNITS=0, LMAX=0, LOVE=None,
         4: microGal gravitational perturbation
         5: mbar equivalent surface pressure
         6: cm of viscoelastic crustal uplift (GIA)
+        list: custom degree-dependent unit conversion factor
     LMAX: Upper bound of Spherical Harmonic Degrees
     LOVE: input load Love numbers up to degree LMAX (hl,kl,ll)
 
@@ -139,8 +142,11 @@ def clenshaw_summation(clm, slm, lon, lat, RAD=0, UNITS=0, LMAX=0, LOVE=None,
     elif (UNITS == 6):
         #-- 6: cmVCU, cm viscoelastic  crustal uplift (GIA ONLY)
         dfactor = factors.cmVCU
+    elif isinstance(UNITS,(list,np.ndarray)):
+        #-- custom units 
+        dfactor = np.copy(UNITS)
     else:
-        raise ValueError('Invalid units code {0:d}'.format(UNITS))
+        raise ValueError('Unknown units {0}'.format(UNITS))
 
     #-- calculate arrays for clenshaw summations over colatitudes
     s_m_c = np.zeros((npts,LMAX*2+2))

gravity_toolkit/gen_disc_load.py

@@ -1,6 +1,6 @@
 #!/usr/bin/env python
 u"""
-gen_disc_load.py (01/2021)
+gen_disc_load.py (11/2021)
 Calculates gravitational spherical harmonic coefficients for a uniform disc load
 
 CALLING SEQUENCE:
@@ -21,7 +21,12 @@
 OPTIONS:
     LMAX: Upper bound of Spherical Harmonic Degrees
     MMAX: Upper bound of Spherical Harmonic Orders
-    PLM: input Legendre polynomials (for improving computational time)
+    UNITS: input data units
+        1: cm of water thickness
+        2: Gigatonnes of mass
+        3: kg/m^2
+        list: custom unit conversion factor
+    PLM: input Legendre polynomials
     LOVE: input load Love numbers up to degree LMAX (hl,kl,ll)
 
 PYTHON DEPENDENCIES:
@@ -52,6 +57,7 @@
         https://doi.org/10.1007/s00190-011-0522-7
 
 UPDATE HISTORY:
+    Updated 11/2021: added UNITS option for converting from different inputs
     Updated 01/2021: use harmonics class for spherical harmonic operations
     Updated 07/2020: added function docstrings
     Updated 05/2020: vectorize calculation over degrees to improve compute time
@@ -69,7 +75,8 @@
 from gravity_toolkit.plm_holmes import plm_holmes
 from gravity_toolkit.legendre_polynomials import legendre_polynomials
 
-def gen_disc_load(data,lon,lat,area,LMAX=60,MMAX=None,PLM=None,LOVE=None):
+def gen_disc_load(data, lon, lat, area, LMAX=60, MMAX=None, UNITS=2,
+    PLM=None, LOVE=None):
     """
     Calculates spherical harmonic coefficients for a uniform disc load
 
@@ -84,6 +91,11 @@ def gen_disc_load(data,lon,lat,area,LMAX=60,MMAX=None,PLM=None,LOVE=None):
     -----------------
     LMAX: Upper bound of Spherical Harmonic Degrees
     MMAX: Upper bound of Spherical Harmonic Orders
+    UNITS: input data units
+        1: cm of water thickness
+        2: Gigatonnes of mass
+        3: kg/m^2
+        list: custom unit conversion factor
     PLM: input Legendre polynomials
     LOVE: input load Love numbers up to degree LMAX (hl,kl,ll)
 
@@ -112,9 +124,24 @@ def gen_disc_load(data,lon,lat,area,LMAX=60,MMAX=None,PLM=None,LOVE=None):
     #-- alpha will be 1 - the ratio of the input area with the half sphere
     alpha = (1.0 - 1e10*area/(2.0*np.pi*rad_e**2))
 
-    #-- Input data is in gigatonnes (Gt)
-    #-- 1e15 converts from Gt to grams, 1e10 converts from km^2 to cm^2
-    unit_conv = 1e15/(1e10*area)
+    #-- Calculate factor to convert from input units into g/cm^2
+    if (UNITS == 1):
+        #-- Input data is in cm water equivalent (cmH2O)
+        unit_conv = 1.0
+    elif (UNITS == 2):
+        #-- Input data is in gigatonnes (Gt)
+        #-- 1e15 converts from Gt to grams, 1e10 converts from km^2 to cm^2
+        unit_conv = 1e15/(1e10*area)
+    elif (UNITS == 3):
+        #-- Input data is in kg/m^2
+        #-- 1 kg = 1000 g
+        #-- 1 m^2 = 100*100 cm^2 = 1e4 cm^2
+        unit_conv = 0.1
+    elif isinstance(UNITS,(list,np.ndarray)):
+        #-- custom units 
+        unit_conv = np.copy(UNITS)
+    else:
+        raise ValueError('Unknown units {0}'.format(UNITS))
 
     #-- Coefficient for calculating Stokes coefficients for a disc load
     #-- From Jacob et al (2012), Farrell (1972) and Longman (1962)