refactor: modified custom units in spherical caps and disc loads (#129)

doc/source/getting_started/Spatial-Maps.rst

@@ -129,7 +129,7 @@ which can be spectrally filtered following [Swenson2006]_.
 References
 ##########
 
-.. [Bandikova2019] T. Bandikova, C. McCullough, G. L. Kruizinga, H. Save, and B. Christophe, "GRACE accelerometer data transplant", *Advances in Space Research*, 64(3), 623--644, (2019). `doi: 10.1016/j.asr.2019.05.021 <10.1016/j.asr.2019.05.021>`_
+.. [Bandikova2019] T. Bandikova, C. McCullough, G. L. Kruizinga, H. Save, and B. Christophe, "GRACE accelerometer data transplant", *Advances in Space Research*, 64(3), 623--644, (2019). `doi: 10.1016/j.asr.2019.05.021 <https://doi.org/10.1016/j.asr.2019.05.021>`_
 
 .. [Blewett2003] G. Blewitt, "Self-consistency in reference frames, geocenter definition, and surface loading of the solid Earth", *Journal of Geophysical Research: Solid Earth*, 108(B2), 2103, (2003). `doi: 10.1029/2002JB002082 <https://doi.org/10.1029/2002JB002082>`_
 

gravity_toolkit/gen_disc_load.py

@@ -1,7 +1,7 @@
 #!/usr/bin/env python
 u"""
 gen_disc_load.py
-Written by Tyler Sutterley (03/2023)
+Written by Tyler Sutterley (06/2023)
 Calculates gravitational spherical harmonic coefficients for a uniform disc load
 
 CALLING SEQUENCE:
@@ -55,6 +55,7 @@
         https://doi.org/10.1007/s00190-011-0522-7
 
 UPDATE HISTORY:
+    Updated 06/2023: modified custom units case to not convert to cmwe
     Updated 03/2023: simplified unit degree factors using units class
         improve typing for variables in docstrings
     Updated 02/2023: set custom units as top option in if/else statements
@@ -153,7 +154,7 @@ def gen_disc_load(data, lon, lat, area, LMAX=60, MMAX=None, UNITS=2,
     th = (90.0 - lat)*np.pi/180.0# Colatitude in radians
 
     # Earth Parameters
-    factors = gravity_toolkit.units(lmax=LMAX).spatial(*LOVE)
+    factors = gravity_toolkit.units(lmax=LMAX)
 
     # convert input area into cm^2 and then divide by area of a half sphere
     # alpha will be 1 - the ratio of the input area with the half sphere
@@ -162,26 +163,32 @@ def gen_disc_load(data, lon, lat, area, LMAX=60, MMAX=None, UNITS=2,
     # Calculate factor to convert from input units into g/cm^2
     if isinstance(UNITS, (list, np.ndarray)):
         # custom units
-        unit_conv = np.copy(UNITS)
+        unit_conv = 1.0
+        dfactor = np.copy(UNITS)
     elif (UNITS == 1):
         # Input data is in cm water equivalent (cmwe)
         unit_conv = 1.0
+        # degree dependent factors to convert from coefficients
+        # of mass into normalized geoid coefficients
+        dfactor = 4.0*np.pi*factors.spatial(*LOVE).cmwe/(1.0 + 2.0*factors.l)
     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)
+        # degree dependent factors to convert from coefficients
+        # of mass into normalized geoid coefficients
+        dfactor = 4.0*np.pi*factors.spatial(*LOVE).cmwe/(1.0 + 2.0*factors.l)
     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
+        # degree dependent factors to convert from coefficients
+        # of mass into normalized geoid coefficients
+        dfactor = 4.0*np.pi*factors.spatial(*LOVE).cmwe/(1.0 + 2.0*factors.l)
     else:
         raise ValueError(f'Unknown units {UNITS}')
 
-    # calculate SH degree dependent factors to convert from coefficients
-    # of mass into normalized geoid coefficients
-    dfactor = 4.0*np.pi*factors.cmwe/(1.0 + 2.0*factors.l)
-
     # Calculating plms of the disc
     # allocating for constructed array
     pl_alpha = np.zeros((LMAX+1))

gravity_toolkit/gen_spherical_cap.py

@@ -1,12 +1,9 @@
 #!/usr/bin/env python
 u"""
 gen_spherical_cap.py
-Written by Tyler Sutterley (03/2023)
+Written by Tyler Sutterley (06/2023)
 Calculates gravitational spherical harmonic coefficients for a spherical cap
 
-Spherical cap derivation from Longman (1962), Farrell (1972), Pollack (1973)
-    and Jacob (2012)
-
 Creating a spherical cap with generating angle alpha is a 2 step process:
     1) obtain harmonics when cap is located at the north pole
     2) rotate the cap to an arbitrary latitude and longitude
@@ -64,6 +61,7 @@
         https://doi.org/10.1007/s00190-011-0522-7
 
 UPDATE HISTORY:
+    Updated 06/2023: modified custom units case to not convert to cmwe
     Updated 03/2023: simplified unit degree factors using units class
         improve typing for variables in docstrings
     Updated 02/2023: set custom units as top option in if/else statements
@@ -174,7 +172,7 @@ def gen_spherical_cap(data, lon, lat, LMAX=60, MMAX=None,
     th = (90.0 - lat)*np.pi/180.0# Colatitude in radians
 
     # Earth Parameters
-    factors = gravity_toolkit.units(lmax=LMAX).spatial(*LOVE)
+    factors = gravity_toolkit.units(lmax=LMAX)
 
     # Converting input area into an equivalent spherical cap radius
     # Following Jacob et al. (2012) Equation 4 and 5
@@ -197,13 +195,17 @@ def gen_spherical_cap(data, lon, lat, LMAX=60, MMAX=None,
     else:
         raise ValueError('Input RAD_CAP, AREA or RAD_KM of spherical cap')
 
-    # Calculate factor to convert from input units into cmwe
+    # Calculate factor to convert from input units
     if isinstance(UNITS, (list, np.ndarray)):
         # custom units
-        unit_conv = np.copy(UNITS)
+        unit_conv = 1.0
+        dfactor = np.copy(UNITS)
     elif (UNITS == 1):
         # Input data is in cm water equivalent (cmwe)
         unit_conv = 1.0
+        # degree dependent factors to convert from coefficients
+        # of mass into normalized geoid coefficients
+        dfactor = 4.0*np.pi*factors.spatial(*LOVE).cmwe/(1.0 + 2.0*factors.l)
     elif (UNITS == 2):
         # Input data is in gigatonnes (Gt)
         # calculate spherical cap area from angular radius
@@ -212,18 +214,20 @@ def gen_spherical_cap(data, lon, lat, LMAX=60, MMAX=None,
         # 1 g/cm^3 = 1000 kg/m^3 = density water
         # 1 Gt = 1 Pg = 1.e15 g
         unit_conv = 1.e15/area
+        # degree dependent factors to convert from coefficients
+        # of mass into normalized geoid coefficients
+        dfactor = 4.0*np.pi*factors.spatial(*LOVE).cmwe/(1.0 + 2.0*factors.l)
     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
+        # degree dependent factors to convert from coefficients
+        # of mass into normalized geoid coefficients
+        dfactor = 4.0*np.pi*factors.spatial(*LOVE).cmwe/(1.0 + 2.0*factors.l)
     else:
         raise ValueError(f'Unknown units {UNITS}')
 
-    # calculate SH degree dependent factors to convert from coefficients
-    # of mass into normalized geoid coefficients
-    dfactor = 4.0*np.pi*factors.cmwe/(1.0 + 2.0*factors.l)
-
     # Calculating plms of the spherical caps
     # From Longman et al. (1962)
     # pl_alpha = F(alpha) from Jacob 2011

gravity_toolkit/utilities.py

@@ -1255,7 +1255,7 @@ def drive_list(
     try:
         # Create and submit request.
         request = urllib2.Request(posixpath.join(*HOST))
-        tree = lxml.etree.parse(urllib2.urlopen(request,timeout=timeout),parser)
+        tree = lxml.etree.parse(urllib2.urlopen(request, timeout=timeout),parser)
     except (urllib2.HTTPError, urllib2.URLError) as exc:
         raise Exception('List error from {0}'.format(posixpath.join(*HOST)))
     else:
@@ -1349,7 +1349,7 @@ def from_drive(
     try:
         # Create and submit request.
         request = urllib2.Request(posixpath.join(*HOST))
-        response = urllib2.urlopen(request,timeout=timeout)
+        response = urllib2.urlopen(request, timeout=timeout)
     except (urllib2.HTTPError, urllib2.URLError) as exc:
         raise Exception('Download error from {0}'.format(posixpath.join(*HOST)))
     else:
@@ -2018,7 +2018,7 @@ def from_figshare(
     local_dir.mkdir(mode=mode, parents=True, exist_ok=True)
     # Create and submit request.
     request = urllib2.Request(posixpath.join(*HOST))
-    response = urllib2.urlopen(request,timeout=timeout,context=context)
+    response = urllib2.urlopen(request, timeout=timeout,context=context)
     resp = json.loads(response.read())
     # reduce list of geocenter files
     geocenter_files = [f for f in resp['files'] if re.match(pattern,f['name'])]
@@ -2352,7 +2352,7 @@ def icgem_list(
     try:
         # Create and submit request.
         request = urllib2.Request(host)
-        tree = lxml.etree.parse(urllib2.urlopen(request,timeout=timeout),parser)
+        tree = lxml.etree.parse(urllib2.urlopen(request, timeout=timeout),parser)
     except:
         raise Exception(f'List error from {host}')
     else:

notebooks/GRACE-Spatial-Maps.ipynb

@@ -11,7 +11,7 @@
     "\n",
     "The primary and secondary instrumentation onboard the GRACE/GRACE-FO satellites are the ranging instrument (GRACE has a microwave ranging instrument, GRACE-FO has both a microwave ranging instrument and a laser interferometer), the global positioning system (GPS), the accelerometers and the star cameras.\n",
     "Data from these instruments are combined to estimate the distance between the two satellites, the positions of the satellites in space, the pointing vector of the satellites and any non-gravitational accelerations the satellites experience.\n",
-    "These measurements combined with background gravity field estimates, atmospheric and oceanic variability estimates and tidal estimates are used to create the [Level-2 spherical harmonic product of GRACE and GRACE-FO](https://podaac-tools.jpl.nasa.gov/drive/files/GeodeticsGravity/gracefo/docs/GRACE-FO_L2-UserHandbook_v1.1.pdf).  \n",
+    "These measurements combined with background gravity field estimates, atmospheric and oceanic variability estimates and tidal estimates are used to create the [Level-2 spherical harmonic product of GRACE and GRACE-FO](https://archive.podaac.earthdata.nasa.gov/podaac-ops-cumulus-docs/gracefo/open/docs/GRACE-FO_L2_UserHandbook.pdf).  \n",
     "\n",
     "There are three main processing centers that create the Level-2 spherical harmonic data as part of the GRACE/GRACE-FO Science Data System (SDS): the [University of Texas Center for Space Research (CSR)](http://www2.csr.utexas.edu/grace/), the [German Research Centre for Geosciences (GeoForschungsZentrum, GFZ)](https://www.gfz-potsdam.de/en/grace/) and the [Jet Propulsion Laboratory (JPL)](https://grace.jpl.nasa.gov/).  \n",
     "\n",