diff --git a/doc/conf.py b/doc/conf.py index 5a2f6dd0b..906e3f731 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -61,6 +61,7 @@ None, ), "pygmt": ("https://www.pygmt.org/latest/", None), + "pymap3d": ("https://geospace-code.github.io/pymap3d/", None), } # Autosummary pages will be generated by sphinx-autogen instead of sphinx-build diff --git a/doc/user_guide/coordinate_systems.rst b/doc/user_guide/coordinate_systems.rst index eafe37e5b..30f6730ec 100644 --- a/doc/user_guide/coordinate_systems.rst +++ b/doc/user_guide/coordinate_systems.rst @@ -210,23 +210,30 @@ the same radius equal to the *mean radius of the Earth*. print("radius:", radius) We can convert spherical coordinates to geodetic ones through -:meth:`boule.Ellipsoid.spherical_to_geodetic`: +:func:`pymap3d.spherical2geodetic`: .. jupyter-execute:: - coordinates_geodetic = ellipsoid.spherical_to_geodetic(*coordinates) - longitude, latitude, height = coordinates_geodetic[:] + import pymap3d + + latitude, longitude, height = pymap3d.spherical2geodetic( + sph_latitude, longitude, radius, ell=ellipsoid + ) print("longitude:", longitude) print("latitude:", latitude) print("height:", height) While the conversion of spherical coordinates into geodetic ones can be -carried out through :meth:`boule.Ellipsoid.geodetic_to_spherical`: +carried out through :func:`pymap3d.geodetic2spherical`: .. jupyter-execute:: - coordinates_spherical = ellipsoid.geodetic_to_spherical(*coordinates_geodetic) - longitude, sph_latitude, radius = coordinates_spherical[:] + sph_latitude, longitude, radius = pymap3d.geodetic2spherical( + latitude, + longitude, + height, + ell=ellipsoid, + ) print("longitude:", longitude) print("spherical latitude:", sph_latitude) print("radius:", radius) diff --git a/doc/user_guide/equivalent_sources/eq-sources-spherical.rst b/doc/user_guide/equivalent_sources/eq-sources-spherical.rst index 0bebeccb1..0b96d8bd0 100644 --- a/doc/user_guide/equivalent_sources/eq-sources-spherical.rst +++ b/doc/user_guide/equivalent_sources/eq-sources-spherical.rst @@ -6,7 +6,7 @@ Equivalent sources in spherical coordinates When interpolating gravity or magnetic data over a very large region that covers full continents we need to take into account the curvature of the Earth. In these cases projecting the data to plain Cartesian coordinates may introduce -errors due to the distorsions caused by it. +errors due to the distortions caused by it. Therefore using the :class:`harmonica.EquivalentSources` class is not well suited for it. @@ -76,14 +76,16 @@ some first guess values for the ``damping`` and ``depth`` parameters. Before we can fit the sources' coefficients we need to convert the data given in geographical coordinates to spherical ones. We can do it through the -:meth:`boule.Ellipsoid.geodetic_to_spherical` method of the WGS84 ellipsoid -defined in :mod:`boule`. +:func:`pymap3d.geodetic2spherical` function. .. jupyter-execute:: - coordinates = ellipsoid.geodetic_to_spherical( - data.longitude, data.latitude, data.height_sea_level_m + import pymap3d + + latitude, longitude, height = pymap3d.geodetic2spherical( + data.latitude, data.longitude, data.height_sea_level_m, ell=ellipsoid ) + coordinates = (longitude, latitude, height) And then use them to fit the sources: @@ -114,7 +116,10 @@ field we need to convert the grid coordinates to spherical. .. jupyter-execute:: - grid_coords_sph = ellipsoid.geodetic_to_spherical(*grid_coords) + grid_lat_sph, grid_lon_sph, grid_upward_sph = pymap3d.geodetic2spherical( + grid_coords[1], grid_coords[0], grid_coords[2], ell=ellipsoid + ) + grid_coords_sph = (grid_lon_sph, grid_lat_sph, grid_upward_sph) And then predict the gravity disturbance on the grid points: @@ -183,4 +188,4 @@ Lets plot it: fig.coast(shorelines="0.5p,black", area_thresh=1e4) fig.colorbar(cmap=True, frame=["a50f25", "x+lmGal"]) - fig.show() \ No newline at end of file + fig.show() diff --git a/doc/user_guide/forward_modelling/point.rst b/doc/user_guide/forward_modelling/point.rst index d4aaefdb4..385595d43 100644 --- a/doc/user_guide/forward_modelling/point.rst +++ b/doc/user_guide/forward_modelling/point.rst @@ -93,7 +93,7 @@ Lets plot this gravitational field: .. jupyter-execute:: - import pygmt + import pygmt grid = vd.make_xarray_grid( coordinates, g_z, data_names="g_z", extra_coords_names="extra") @@ -167,8 +167,8 @@ Working with sources defined in geodetic coordinates ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ If our point sources and computation points are defined in geodetic -coordinates, we can use the :meth:`boule.Ellipsoid.geodetic_to_spherical` -method to convert them to spherical coordinates and use them to compute their +coordinates, we can use the :func:`pymap3d.geodetic2spherical` function to +convert them to spherical coordinates and use them to compute their gravitational field. Lets define a set point sources in geodetic coordinates and their masses: @@ -196,12 +196,21 @@ coordinates and located at 1km above the ellipsoid: Before we can start forward modelling these sources, we need to convert them to spherical coordinates. To do so, we can use the -:meth:`boule.Ellipsoid.geodetic_to_spherical` method: +:func:`pymap3d.geodetic2spherical` function: .. jupyter-execute:: - points_spherical = ellipsoid.geodetic_to_spherical(*points) - coordinates_spherical = ellipsoid.geodetic_to_spherical(*coordinates) + import pymap3d + + lat_sph, lon, radius = pymap3d.geodetic2spherical( + points[1], points[0], points[2], ell=ellipsoid + ) + points_spherical = (lon, lat_sph, radius) + + lat_sph, lon, radius = pymap3d.geodetic2spherical( + coordinates[1], coordinates[0], coordinates[2], ell=ellipsoid + ) + coordinates_spherical = (lon, lat_sph, radius) We can finally use these converted coordinates to compute the gravitational field the source generate on every computation point: @@ -220,7 +229,7 @@ Lets plot these results using :mod:`pygmt`: .. jupyter-execute:: - import pygmt + import pygmt grid = vd.make_xarray_grid( coordinates_spherical, g_z, data_names="g_z", extra_coords_names="extra") @@ -236,6 +245,6 @@ Lets plot these results using :mod:`pygmt`: grid=grid.g_z, frame=[f"WSne+t{title}", "x", "y"], cmap=True,) - + fig.colorbar(cmap=True, position="JMR", frame=["a0.000000005", "x+lmGal"]) fig.show() diff --git a/env/requirements-docs.txt b/env/requirements-docs.txt index a20b7d151..fd42b7f5e 100644 --- a/env/requirements-docs.txt +++ b/env/requirements-docs.txt @@ -6,6 +6,7 @@ sphinx-design==0.2.* sphinx-copybutton==0.5.* jupyter-sphinx==0.3.* boule +pymap3d pyproj ensaio netcdf4 diff --git a/env/requirements-tests.txt b/env/requirements-tests.txt index 32f8f1b22..7c166f96e 100644 --- a/env/requirements-tests.txt +++ b/env/requirements-tests.txt @@ -1,5 +1,6 @@ # Requirements for running the tests boule +pymap3d pytest pytest-cov coverage diff --git a/environment.yml b/environment.yml index cae320ea2..bda9bb773 100644 --- a/environment.yml +++ b/environment.yml @@ -28,6 +28,7 @@ dependencies: - pytest-cov - coverage - boule + - pymap3d # Documentation requirements - sphinx==4.5.* - sphinx-book-theme==0.2.* diff --git a/examples/equivalent_sources/spherical.py b/examples/equivalent_sources/spherical.py index d265d0347..7fc7fb38a 100644 --- a/examples/equivalent_sources/spherical.py +++ b/examples/equivalent_sources/spherical.py @@ -21,6 +21,7 @@ import boule as bl import numpy as np import pygmt +import pymap3d import verde as vd import harmonica as hm @@ -44,7 +45,10 @@ # Convert data coordinates from geodetic (longitude, latitude, height) to # spherical (longitude, spherical_latitude, radius). -coordinates = ellipsoid.geodetic_to_spherical(longitude, latitude, elevation) +lat_sph, lon, radius = pymap3d.geodetic2spherical( + latitude, longitude, elevation, ell=ellipsoid +) +coordinates = (lon, lat_sph, radius) # Create the equivalent sources eqs = hm.EquivalentSourcesSph(damping=1e-3, relative_depth=10000) diff --git a/harmonica/tests/test_forward_utils.py b/harmonica/tests/test_forward_utils.py index ef77cd480..877c7835d 100644 --- a/harmonica/tests/test_forward_utils.py +++ b/harmonica/tests/test_forward_utils.py @@ -11,6 +11,7 @@ import numpy as np import numpy.testing as npt import pytest +from pymap3d import geodetic2spherical from ..forward.utils import check_coordinate_system, distance @@ -67,8 +68,10 @@ def test_geodetic_distance_vs_spherical(): # Compute distance using closed-form formula dist = distance(point_a, point_b, coordinate_system="geodetic", ellipsoid=ellipsoid) # Convert points to spherical coordinates - point_a_sph = ellipsoid.geodetic_to_spherical(*point_a) - point_b_sph = ellipsoid.geodetic_to_spherical(*point_b) + lat, lon, up = geodetic2spherical(point_a[1], point_a[0], point_a[2], ell=ellipsoid) + point_a_sph = (lon, lat, up) + lat, lon, up = geodetic2spherical(point_b[1], point_b[0], point_b[2], ell=ellipsoid) + point_b_sph = (lon, lat, up) # Compute distance using these converted points dist_sph = distance(point_a_sph, point_b_sph, coordinate_system="spherical") npt.assert_allclose(dist, dist_sph)