diff --git a/README.md b/README.md index 2c25a11..78fd941 100644 --- a/README.md +++ b/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 diff --git a/examples/scripts/TeMaari_ocean_loading_pyhardisp.py b/examples/scripts/TeMaari_ocean_loading_pyhardisp.py index fb8c5e2..f8a33e8 100644 --- a/examples/scripts/TeMaari_ocean_loading_pyhardisp.py +++ b/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) - diff --git a/pyproject.toml b/pyproject.toml index f6cf800..c263736 100644 --- a/pyproject.toml +++ b/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", diff --git a/src/gsolve/reductions/anomalies.py b/src/gsolve/reductions/anomalies.py index fd505fb..f0346ad 100644 --- a/src/gsolve/reductions/anomalies.py +++ b/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 diff --git a/src/gsolve/reductions/corrections.py b/src/gsolve/reductions/corrections.py index 6ecf2ca..bd10a7c 100644 --- a/src/gsolve/reductions/corrections.py +++ b/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) diff --git a/tests/test_corrections.py b/tests/test_corrections.py index 07a112e..85b312e 100644 --- a/tests/test_corrections.py +++ b/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 diff --git a/uv.lock b/uv.lock index 0de95b7..bad9896 100644 --- a/uv.lock +++ b/uv.lock @@ -92,16 +92,16 @@ wheels = [ [[package]] name = "boule" -version = "0.5.0" +version = "0.6.0" source = { registry = "https://pypi.org/simple" } dependencies = [ { name = "attrs" }, { name = "numpy" }, { name = "scipy" }, ] -sdist = { url = "https://files.pythonhosted.org/packages/6f/cb/ba840a6f38c045310ac223edc799a69842a959f549aeea3fee0644472e83/boule-0.5.0.tar.gz", hash = "sha256:fde5908b8a88f0f05ccbf23f2aa422b9d81a2a3e6ed6c1ce65afda0bb6254717", size = 33514, upload-time = "2024-10-22T21:35:43.393Z" } +sdist = { url = "https://files.pythonhosted.org/packages/05/3e/36215c05560b9daadaccbd05950cbe30484c91605dbc15c4ea1b690e9ffd/boule-0.6.0.tar.gz", hash = "sha256:ee74bec713b9694154aaffafc8ba99f448e6cd1bd59b73b9986804570c7fcad3", size = 32133, upload-time = "2026-03-30T13:53:37.336Z" } wheels = [ - { url = "https://files.pythonhosted.org/packages/3e/a1/fff65b40476cbdd2626aca3f744347921aca9656638277302679c2a608d2/boule-0.5.0-py3-none-any.whl", hash = "sha256:374129df477a4551db9cdfae7fb1578f9d3f1978f826f5d60bd0eab98ba31f15", size = 37180, upload-time = "2024-10-22T21:35:42.026Z" }, + { url = "https://files.pythonhosted.org/packages/3f/c2/680c115640f145466db2d9988484141d0822856a13ca9eade47f6301d846/boule-0.6.0-py3-none-any.whl", hash = "sha256:8b4812d088c12d64bc2ab9e7cdfa6a33d624cf7227de9d281fc4ff976dc73f39", size = 32516, upload-time = "2026-03-30T13:53:36.023Z" }, ] [[package]] @@ -658,7 +658,7 @@ docs = [ [package.metadata] requires-dist = [ - { name = "boule", specifier = "==0.5.*" }, + { name = "boule", specifier = "==0.6.*" }, { name = "harmonica", specifier = "<0.8" }, { name = "pandas", extras = ["excel", "plot"], specifier = "<3" }, { name = "pygtide", specifier = ">=0.8.2" },