updates for boule v0.6.0

README.md

@@ -1,8 +1,10 @@
 ![gSolve logo](docs/source/_static/gsolve_logo.png)
 
 [![codecov](https://codecov.io/gh/GNS-Science/gsolve/branch/main/graph/badge.svg)](https://codecov.io/gh/GNS-Science/gsolve)
-![GitHub License](https://img.shields.io/github/license/GNS-Science/gsolve)
-[![Publish to PyPI](https://github.com/GNS-Science/gsolve/actions/workflows/publish.yml/badge.svg)](https://github.com/GNS-Science/gsolve/actions/workflows/publish.yml)
+[![GitHub License](https://img.shields.io/github/license/GNS-Science/gsolve)](https://github.com/GNS-Science/gsolve/blob/main/LICENSE)
+![Publish to PyPI](https://github.com/GNS-Science/gsolve/actions/workflows/publish.yml/badge.svg)
+[![Pypi version](https://img.shields.io/pypi/v/gsolve)](https://pypi.org/project/gsolve/)
+
 
 
 # gSolve
@@ -24,10 +26,11 @@ Process gravity data for time change microgravity and Bouguer surveys.
 ## Algorithm
 
 * calculate the meter calibration correction, beta.
-* corrects for earth tide using Longman and ETERNA formula.
-* option to correct for ocean loading using pyhardisp
+* corrects for earth tide using Longman and [ETERNA](https://github.com/hydrogeoscience/pygtide/) formula.
+* option to correct for ocean loading using [pyhardisp](https://github.com/craigmillernz/pyhardisp).
 * correct for drift across loops or whole survey.
 * network adjustment with three different network adjustment algorithms depending on user requirements.
+# calibrate meters to absolute values.
 * residuals can be filtered using a percentile cut filter.  
 
 ## Corrections
@@ -58,7 +61,7 @@ Example scripts are provided to demonstrate the software capability.
 
 # Support
 
-Full documentation is available here. [gSolve](https://silver-spork-5vqw3le.pages.github.io/fundamentals.html)
+Full documentation is available here. [gSolve](https://gns-science.github.io/gsolve/)
 
 # Authors and acknowledgment
 

examples/scripts/TeMaari_ocean_loading_pyhardisp.py

@@ -24,8 +24,6 @@
 import pathlib
 
 from gsolve import (
-    GravityAnomalies,
-    GravityCorrectionParameters,
     GravityObservations,
     GravitySites,
     GravitySurvey,
@@ -35,6 +33,7 @@
 from gsolve.reports import GSolveReport
 from gsolve.tide.earth_tide import LongmanTidalCorrection
 from gsolve.tide.ocean_load import HardispOceanLoadCorrector
+
 # %%
 data_path = pathlib.Path(__file__).parent.parent
 
@@ -52,13 +51,13 @@
 # %%
 # Read in observations
 obs = GravityObservations.from_excel(
-    survey_file, sheet_name='obs', parse_split_datetime=False
-    )
+    survey_file, sheet_name="obs", parse_split_datetime=False
+)
 
 # Read in site location information
 sites = GravitySites.from_excel(survey_file, sheet_name="Locations")
 
-# set which reference stations are used in this survey (must be in sites). In 
+# set which reference stations are used in this survey (must be in sites). In
 # this case i just set TGKB to be 0.0 mGal
 ref_sites = ReferenceGravity(site_id="TGKB", gravity=0.0)
 _ = sites.set_reference_gravity(ref_sites)
@@ -111,4 +110,3 @@
 
 
 report = GSolveReport(observations=obs, sites=sites, results=results)
-

pyproject.toml

@@ -25,7 +25,7 @@ classifiers = [
     "Topic :: Scientific/Engineering",
 ]
 dependencies = [
-    "boule==0.5.*",
+    "boule==0.6.*",
     "harmonica<0.8",
     "pandas[excel,plot]<3",
     "pygtide>=0.8.2",

src/gsolve/reductions/anomalies.py

@@ -345,8 +345,8 @@ class GravityAnomalies(GSolveTable):
         "free_air_correction": DataFieldSpecification(
             "free_air_correction", float, required=False, default=_np.nan
         ),
-        "bouguer_slab": DataFieldSpecification(
-            "bouguer_slab", float, required=False, default=_np.nan
+        "bouguer_slab_correction": DataFieldSpecification(
+            "bouguer_slab_correction", float, required=False, default=_np.nan
         ),
         "bouguer_slab_curvature_corrected": DataFieldSpecification(
             "bouguer_slab_curvature_corrected", float, required=False, default=_np.nan

src/gsolve/reductions/corrections.py

@@ -47,6 +47,7 @@
 
 
 def normal_gravity_at_stn_elevation(
+    longitude: ArrayLike,
     latitude: ArrayLike,
     height_ellipsoidal: ArrayLike,
     ellipsoid: Literal["WGS84", "GRS80"] | boule.Ellipsoid = "GRS80",
@@ -69,9 +70,11 @@ def normal_gravity_at_stn_elevation(
 
     Parameters
     ----------
+    longitude : array-like
+        The geodetic longitude in decimal degrees.
     latitude : array-like
         The geodetic latitude in decimal degrees.
-    height_ellipsoidal : array_likefloat or array-like
+    height_ellipsoidal : float or array-like
         The ellipsoidal height in meters.
     ellipsoid: 'WGS84', 'GRS80' or boule.Ellipsoid, default 'GRS80'
         The ellipsoid to use for normal gravity calculation.
@@ -106,7 +109,7 @@ def normal_gravity_at_stn_elevation(
         raise ValueError(f"Unknown ellipsoid '{ellipsoid}': must be 'WGS84' or 'GRS80'")
 
     return e.normal_gravity(
-        latitude=latitude, height=height_ellipsoidal, si_units=si_units
+        coordinates=(longitude, latitude, height_ellipsoidal), si_units=si_units
     )
 
 
@@ -577,6 +580,7 @@ class GravityCorrections(GSolveTable):
 
     _known_fields = {
         "site_id": DataFieldSpecification("site_id", str, required=True),
+        "longitude": DataFieldSpecification("longitude", float, required=False),
         "latitude": DataFieldSpecification("latitude", float, required=False),
         "height_ellipsoidal": DataFieldSpecification(
             "height_ellipsoidal", float, required=False, legacy_name="height"
@@ -697,17 +701,17 @@ def compute(
         Parameters
         ----------
         sites : GravitySites | DataFrame
-            An object providing site latitude and ellipsoidal height, and indexed
+            An object providing site longitude, latitude and ellipsoidal height, and indexed
             by ``'site_id'``. If `sites` is a DataFrame, it is expected to have columns named
-            ``'latitude'`` and ``'height_ellipsoidal'``, unless alternative columns are
+            ``'longitude'``, ``'latitude'`` and ``'height_ellipsoidal'``, unless alternative columns are
             specified using the `column_names` argument.
         corrections : str | Sequence[str], optional
             An array or string of corrections to compute.  By default compute all
             corrections required for generating a Bouguer anomaly as specified in
             ``self.params``.
         column_names : dict[str, str] | None, optional
-            A dictionary mapping dexpected columns ``latitude`` and ``height_ellipsoidal``
-            to alternative column names. E.g. ``{'latitude': 'lat', 'height_ellipsoidal': 'height'}``.
+            A dictionary mapping expected columns ``longitude``, ``latitude`` and ``height_ellipsoidal``
+            to alternative column names. E.g. ``{'longitude': 'lon', 'latitude': 'lat', 'height_ellipsoidal': 'height'}``.
         include_coords : bool, default False
             If True, include site latitude and height in output.
 
@@ -717,7 +721,11 @@ def compute(
             Object containing computed gravity corrections and the correction parameters.
 
         """
-        _cols = {"latitude": "latitude", "height_ellipsoidal": "height_ellipsoidal"}
+        _cols = {
+            "longitude": "longitude",
+            "latitude": "latitude",
+            "height_ellipsoidal": "height_ellipsoidal",
+        }
         if column_names is not None:
             _cols.update(column_names)
 
@@ -731,6 +739,7 @@ def compute(
                 f"'{type(sites)}'"
             )
 
+        lon = sites_df[_cols["longitude"]].to_numpy()
         lat = sites_df[_cols["latitude"]].to_numpy()
         ht = sites_df[_cols["height_ellipsoidal"]].to_numpy()
         idx = sites_df.index.copy()
@@ -749,13 +758,17 @@ def compute(
         else:
             _corrs = self.params.bouguer_correction_fields()
 
-        df = pd.DataFrame(index=idx, data={"latitude": lat, "height_ellipsoidal": ht})
+        df = pd.DataFrame(
+            index=idx,
+            data={"longitude": lon, "latitude": lat, "height_ellipsoidal": ht},
+        )
         for c in _corrs:
             df[c] = np.nan
 
         k = "normal_gravity_at_stn_elevation"
         if k in _corrs:
             df[k] = normal_gravity_at_stn_elevation(
+                longitude=lon,
                 latitude=lat,
                 height_ellipsoidal=ht,
                 ellipsoid=self.params.ellipsoid,  # type: ignore[arg-type]
@@ -799,7 +812,13 @@ def compute(
             )
 
         if not include_coords:
-            df = df.drop(columns=[_cols["latitude"], _cols["height_ellipsoidal"]])
+            df = df.drop(
+                columns=[
+                    _cols["longitude"],
+                    _cols["latitude"],
+                    _cols["height_ellipsoidal"],
+                ]
+            )
         _df_dict = {str(k): v for k, v in df.to_dict("list").items()}
         return GravityCorrections(params=self.params, site_id=idx, **_df_dict)
 

tests/test_corrections.py

@@ -65,42 +65,45 @@ def _make_sites(
 
 class TestNormalGravityAtStnElevation:
     def test_scalar_at_equator_sea_level(self):
-        g = normal_gravity_at_stn_elevation(latitude=0.0, height_ellipsoidal=0.0)
+        g = normal_gravity_at_stn_elevation(
+            longitude=0.0, latitude=0.0, height_ellipsoidal=0.0
+        )
         # GRS80 equatorial normal gravity is ~978032.7 mGal
         assert pytest.approx(g, rel=1e-4) == 978032.67715
 
     def test_si_units(self):
-        g_mgal = normal_gravity_at_stn_elevation(0.0, 0.0, si_units=False)
-        g_si = normal_gravity_at_stn_elevation(0.0, 0.0, si_units=True)
+        g_mgal = normal_gravity_at_stn_elevation(0.0, 0.0, 0.0, si_units=False)
+        g_si = normal_gravity_at_stn_elevation(0.0, 0.0, 0.0, si_units=True)
         assert pytest.approx(g_si, rel=1e-6) == g_mgal * 1e-5
 
     def test_wgs84_ellipsoid(self):
-        g = normal_gravity_at_stn_elevation(0.0, 0.0, ellipsoid="WGS84")
+        g = normal_gravity_at_stn_elevation(0.0, 0.0, 0.0, ellipsoid="WGS84")
         assert g > 0
 
     def test_boule_ellipsoid_object(self):
-        g = normal_gravity_at_stn_elevation(0.0, 0.0, ellipsoid=boule.GRS80)
+        g = normal_gravity_at_stn_elevation(0.0, 0.0, 0.0, ellipsoid=boule.GRS80)
         assert g > 0
 
     def test_negative_height_raises(self):
         with pytest.raises(ValueError, match="heights must be >= 0"):
-            normal_gravity_at_stn_elevation(0.0, -1.0)
+            normal_gravity_at_stn_elevation(0.0, 0.0, -1.0)
 
     def test_bad_ellipsoid_raises(self):
         with pytest.raises(ValueError, match="Unknown ellipsoid"):
-            normal_gravity_at_stn_elevation(0.0, 0.0, ellipsoid="GRS67")
+            normal_gravity_at_stn_elevation(0.0, 0.0, 0.0, ellipsoid="GRS67")
 
     def test_array_input(self):
+        lons = np.array([0.0, 0.0, 0.0])
         lats = np.array([0.0, -45.0, -90.0])
         hts = np.array([0.0, 0.0, 0.0])
-        g = normal_gravity_at_stn_elevation(lats, hts)
+        g = normal_gravity_at_stn_elevation(lons, lats, hts)
         assert g.shape == (3,)
         # gravity at pole > gravity at equator
         assert g[2] > g[0]
 
     def test_gravity_decreases_with_height(self):
-        g0 = normal_gravity_at_stn_elevation(0.0, 0.0)
-        g100 = normal_gravity_at_stn_elevation(0.0, 100.0)
+        g0 = normal_gravity_at_stn_elevation(0.0, 0.0, 0.0)
+        g100 = normal_gravity_at_stn_elevation(0.0, 0.0, 100.0)
         assert g100 < g0
 
 
@@ -416,7 +419,11 @@ def test_compute_with_gravity_sites(self):
 
     def test_compute_with_dataframe(self):
         df = pd.DataFrame(
-            {"latitude": [-45.0, 0.0], "height_ellipsoidal": [100.0, 0.0]},
+            {
+                "longitude": [0.0, 0.0],
+                "latitude": [-45.0, 0.0],
+                "height_ellipsoidal": [100.0, 0.0],
+            },
             index=pd.Index(["A", "B"], name="site_id"),
         )
         provider = GravityCorrectionProvider()
@@ -450,6 +457,7 @@ def test_compute_include_coords(self):
         gc = provider.compute(
             sites, corrections=["free_air_correction"], include_coords=True
         )
+        assert "longitude" in gc.data.columns
         assert "latitude" in gc.data.columns
         assert "height_ellipsoidal" in gc.data.columns