Use pymap3d for coordinate conversions

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

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)

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()
+ fig.show()

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()

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

env/requirements-tests.txt

@@ -1,5 +1,6 @@
 # Requirements for running the tests
 boule
+pymap3d
 pytest
 pytest-cov
 coverage

environment.yml

@@ -28,6 +28,7 @@ dependencies:
  - pytest-cov
  - coverage
  - boule
+ - pymap3d
  # Documentation requirements
  - sphinx==4.5.*
  - sphinx-book-theme==0.2.*

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)

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)