Repository: fatiando/harmonica
Pull request #61 (2019-05-08 13:06:23)

Refactor coordinates conversion
Geographic coordinates conversion functions now get the full coordinates
as arguments rather than only the `latitude` and `radius`/`height`, and returns
the full set of converted coordinates.
The previous implementation was inconvenient for real applications.

**Reminders**

- [ ] Run `make format` and `make check` to make sure the code follows the style guide.
- [ ] Add tests for new features or tests that would have caught the bug that you're fixing.
- [ ] Add new public functions/methods/classes to `doc/api/index.rst`.
- [ ] Write detailed docstrings for all functions/methods.
- [ ] If adding new functionality, add an example to the docstring, gallery, and/or tutorials.


leouieda: @santisoler just noticed that if we want these functions to work with Verde gridders, they should actually take `longitude, latitude, height` instead of `coordinates`: https://github.com/fatiando/verde/blob/master/verde/base/base_classes.py#L294

For example, if we make a spherical equivalent layer using Verde, we could use these functions to generate geodetic grids easily by passing them in as the `projection` argument to the `grid` method.

I don't think it's too inconvenient since the user code would change from `geodetic_to_spherical(coordinates)` to `geodetic_to_spherical(*coordinates)`.
santisoler: I'm a little bit confused about your proposal.

> just noticed that if we want these functions to work with Verde gridders

On what cases you want Verde gridders to convert from spherical to geodetic or from geodetic to spherical?
Maybe I'm missing an obvious example for that...

The `coordinates` tuple on the line you linked is expected to have only `longitude`, `latitude` coordinates, not `longitude`, `latitude`, `height`. Remeber that the `height`/`radius` array is crucial for the spherical to geodetic conversion. I feel I'm missing something around this...
leouieda: The use case I had in mind what something like this:

```python
# Load observed gravity in geodetic coordinates
lon, lat, height, gravity = np.loadtxt(...)

# Convert to geocentric spherical so we can use a spherical gridder
coordinates = hm.geodetic_to_spherical((lon, lat, height))

# Grid using equivalent sources in spherical coordinates
eql = hm.EQLHarmonic(kernel="spherical", ...)
eql.fit(coordinates, gravity)

# Generate a geodetic grid using our spherical gridder
grid = eql.grid(region=(...), spacing=0.1, extra_coords=[1000],
                projection=hm.spherical_to_geodetic)
# extra_coords is the ellipsoidal height of the grid
```

In this case, we can use a spherical gridder to generate a geodetic grid without much effort. This is the same way we can make geodetic grids from Cartesian gridders in Verde.

The expectation of `coordinates=(latitude, longitude)` is only in the docstring. The code itself doesn't care as long as the number of coordinates `projection` accepts is the same as the number `grid_coordinates` returns. Using `extra_coords` will make it return 3 instead of 2.

santisoler: Thanks for that example @leouieda.
Now I understand what you meant.
I agree and just made those changes.
I think it's ready to be merged.

Patch:
diff --git a/harmonica/coordinates.py b/harmonica/coordinates.py
index 227f3cf..4eb6202 100644
--- a/harmonica/coordinates.py
+++ b/harmonica/coordinates.py
@@ -6,7 +6,7 @@ import numpy as np
 from . import get_ellipsoid
 
 
-def geodetic_to_spherical(latitude, height):
+def geodetic_to_spherical(longitude, latitude, height):
     """
     Convert from geodetic to geocentric spherical coordinates.
 
@@ -16,23 +16,50 @@ def geodetic_to_spherical(latitude, height):
 
     Parameters
     ----------
-    latitude : float or array
-        The geodetic latitude (in degrees).
-    height : float or array
-        The ellipsoidal (geometric or geodetic) height (in meters).
+    longitude : array
+        Longitude coordinates on geodetic coordinate system in degrees.
+    latitude : array
+        Latitude coordinates on geodetic coordinate system in degrees.
+    height : array
+        Ellipsoidal heights in meters.
 
     Returns
     -------
-    geocentric_latitude : float or array
-        The latitude coordinate in the geocentric spherical reference system
-        (in degrees).
-    radius : float or array
-        The radial coordinate in the geocentric spherical reference system (in meters).
+    longitude : array
+        Longitude coordinates on geocentric spherical coordinate system in degrees.
+        The longitude coordinates are not modified during this conversion.
+    spherical_latitude : array
+        Converted latitude coordinates on geocentric spherical coordinate system in
+        degrees.
+    radius : array
+        Converted spherical radius coordinates in meters.
 
     See also
     --------
     spherical_to_geodetic : Convert from geocentric spherical to geodetic coordinates.
+
+    Examples
+    --------
+
+    In the poles, the radius should be the reference ellipsoid's semi-minor axis:
+
+    >>> import harmonica as hm
+    >>> spherical = hm.geodetic_to_spherical(longitude=0, latitude=90, height=0)
+    >>> print(", ".join("{:.4f}".format(i) for i in spherical))
+    0.0000, 90.0000, 6356752.3142
+    >>> print("{:.4f}".format(hm.get_ellipsoid().semiminor_axis))
+    6356752.3142
+
+    In  the equator, it should be the semi-major axis:
+
+    >>> spherical = hm.geodetic_to_spherical(longitude=0, latitude=0, height=0)
+    >>> print(", ".join("{:.4f}".format(i) for i in spherical))
+    0.0000, 0.0000, 6378137.0000
+    >>> print("{:.4f}".format(hm.get_ellipsoid().semimajor_axis))
+    6378137.0000
+
     """
+    # Get ellipsoid
     ellipsoid = get_ellipsoid()
     # Convert latitude to radians
     latitude_rad = np.radians(latitude)
@@ -46,11 +73,11 @@ def geodetic_to_spherical(latitude, height):
         height + (1 - ellipsoid.first_eccentricity ** 2) * prime_vertical_radius
     ) * np.sin(latitude_rad)
     radius = np.sqrt(xy_projection ** 2 + z_cartesian ** 2)
-    geocentric_latitude = np.degrees(np.arcsin(z_cartesian / radius))
-    return geocentric_latitude, radius
+    spherical_latitude = np.degrees(np.arcsin(z_cartesian / radius))
+    return longitude, spherical_latitude, radius
 
 
-def spherical_to_geodetic(geocentric_latitude, radius):
+def spherical_to_geodetic(longitude, spherical_latitude, radius):
     """
     Convert from geocentric spherical to geodetic coordinates.
 
@@ -60,25 +87,56 @@ def spherical_to_geodetic(geocentric_latitude, radius):
 
     Parameters
     ----------
-    geocentric_latitude : float or array
-        The latitude coordinate in the geocentric spherical reference system
-        (in degrees).
-    radius : float or array
-        The radial coordinate in the geocentric spherical reference system (in meters).
+    longitude : array
+        Longitude coordinates on geocentric spherical coordinate system in degrees.
+    spherical_latitude : array
+        Latitude coordinates on geocentric spherical coordinate system in degrees.
+    radius : array
+        Spherical radius coordinates in meters.
 
     Returns
     -------
-    latitude : float or array
-        The geodetic latitude (in degrees).
-    height : float or array
-        The ellipsoidal (geometric or geodetic) height (in meters).
+    longitude : array
+        Longitude coordinates on geodetic coordinate system in degrees.
+        The longitude coordinates are not modified during this conversion.
+    latitude : array
+        Converted latitude coordinates on geodetic coordinate system in degrees.
+    height : array
+        Converted ellipsoidal height coordinates in meters.
 
     See also
     --------
     geodetic_to_spherical : Convert from geodetic to geocentric spherical coordinates.
+
+    Examples
+    --------
+
+    In the poles and equator, using the semi-minor or semi-major axis of the ellipsoid
+    as the radius should yield 0 height:
+
+    >>> import harmonica as hm
+    >>> geodetic = hm.spherical_to_geodetic(
+    ...     longitude=0, spherical_latitude=90, radius=hm.get_ellipsoid().semiminor_axis
+    ... )
+    >>> print(", ".join("{:.1f}".format(i) for i in geodetic))
+    0.0, 90.0, 0.0
+    >>> geodetic = hm.spherical_to_geodetic(
+    ...     longitude=0, spherical_latitude=0, radius=hm.get_ellipsoid().semimajor_axis
+    ... )
+    >>> print(", ".join("{:.1f}".format(i) for i in geodetic))
+    0.0, 0.0, 0.0
+    >>> geodetic = hm.spherical_to_geodetic(
+    ...     longitude=0,
+    ...     spherical_latitude=-90,
+    ...     radius=hm.get_ellipsoid().semiminor_axis + 2
+    ... )
+    >>> print(", ".join("{:.1f}".format(i) for i in geodetic))
+    0.0, -90.0, 2.0
+
     """
+    # Get ellipsoid
     ellipsoid = get_ellipsoid()
-    k, big_d, big_z = _spherical_to_geodetic_parameters(geocentric_latitude, radius)
+    k, big_d, big_z = _spherical_to_geodetic_parameters(spherical_latitude, radius)
     latitude = np.degrees(
         2 * np.arctan(big_z / (big_d + np.sqrt(big_d ** 2 + big_z ** 2)))
     )
@@ -87,18 +145,19 @@ def spherical_to_geodetic(geocentric_latitude, radius):
         / k
         * np.sqrt(big_d ** 2 + big_z ** 2)
     )
-    return latitude, height
+    return longitude, latitude, height
 
 
-def _spherical_to_geodetic_parameters(geocentric_latitude, radius):
+def _spherical_to_geodetic_parameters(spherical_latitude, radius):
     "Compute parameters for spherical to geodetic coordinates conversion"
+    # Get ellipsoid
     ellipsoid = get_ellipsoid()
     # Convert latitude to radians
-    geocentric_latitude_rad = np.radians(geocentric_latitude)
-    big_z = radius * np.sin(geocentric_latitude_rad)
+    spherical_latitude_rad = np.radians(spherical_latitude)
+    big_z = radius * np.sin(spherical_latitude_rad)
     p_0 = (
         radius ** 2
-        * np.cos(geocentric_latitude_rad) ** 2
+        * np.cos(spherical_latitude_rad) ** 2
         / ellipsoid.semimajor_axis ** 2
     )
     q_0 = (
@@ -116,7 +175,7 @@ def _spherical_to_geodetic_parameters(geocentric_latitude, radius):
     big_d = (
         k
         * radius
-        * np.cos(geocentric_latitude_rad)
+        * np.cos(spherical_latitude_rad)
         / (k + ellipsoid.first_eccentricity ** 2)
     )
     return k, big_d, big_z


Test patch:
diff --git a/harmonica/tests/test_coordinates.py b/harmonica/tests/test_coordinates.py
index d660b0b..8b53d07 100644
--- a/harmonica/tests/test_coordinates.py
+++ b/harmonica/tests/test_coordinates.py
@@ -22,37 +22,52 @@ def test_geodetic_to_spherical_with_spherical_ellipsoid():
     spherical_ellipsoid = ReferenceEllipsoid(
         "unit_sphere", sphere_radius, np.infty, 0, 0
     )
+    size = 5
+    longitude = np.linspace(0, 180, size)
     with set_ellipsoid(spherical_ellipsoid):
-        latitude = np.linspace(-90, 90, 5)
-        height = np.linspace(-0.2, 0.2, 5)
-        geocentric_latitude, radius = geodetic_to_spherical(latitude, height)
-        npt.assert_allclose(geocentric_latitude, latitude, rtol=rtol)
+        latitude = np.linspace(-90, 90, size)
+        height = np.linspace(-0.2, 0.2, size)
+        spherical_longitude, spherical_latitude, radius = geodetic_to_spherical(
+            longitude, latitude, height
+        )
+        npt.assert_allclose(spherical_longitude, longitude, rtol=rtol)
+        npt.assert_allclose(spherical_latitude, latitude, rtol=rtol)
         npt.assert_allclose(radius, sphere_radius + height, rtol=rtol)
 
 
 def test_geodetic_to_spherical_on_equator():
     "Test geodetic to geocentric coordinates conversion on equator."
     rtol = 1e-10
-    height = np.linspace(-1e4, 1e4, 5)
-    latitude = np.full(height.shape, 0)
+    size = 5
+    longitude = np.linspace(0, 180, size)
+    height = np.linspace(-1e4, 1e4, size)
+    latitude = np.zeros_like(size)
     for ellipsoid_name in KNOWN_ELLIPSOIDS:
         with set_ellipsoid(ellipsoid_name):
             ellipsoid = get_ellipsoid()
-            geocentric_latitude, radius = geodetic_to_spherical(latitude, height)
-            npt.assert_allclose(geocentric_latitude, latitude, rtol=rtol)
+            spherical_longitude, spherical_latitude, radius = geodetic_to_spherical(
+                longitude, latitude, height
+            )
+            npt.assert_allclose(spherical_longitude, longitude, rtol=rtol)
+            npt.assert_allclose(spherical_latitude, latitude, rtol=rtol)
             npt.assert_allclose(radius, height + ellipsoid.semimajor_axis, rtol=rtol)
 
 
 def test_geodetic_to_spherical_on_poles():
     "Test geodetic to geocentric coordinates conversion on poles."
     rtol = 1e-10
-    height = np.hstack([np.linspace(-1e4, 1e4, 5)] * 2)
-    latitude = np.array([90.0] * 5 + [-90.0] * 5)
+    size = 5
+    longitude = np.hstack([np.linspace(0, 180, size)] * 2)
+    height = np.hstack([np.linspace(-1e4, 1e4, size)] * 2)
+    latitude = np.array([90.0] * size + [-90.0] * size)
     for ellipsoid_name in KNOWN_ELLIPSOIDS:
         with set_ellipsoid(ellipsoid_name):
             ellipsoid = get_ellipsoid()
-            geocentric_latitude, radius = geodetic_to_spherical(latitude, height)
-            npt.assert_allclose(geocentric_latitude, latitude, rtol=rtol)
+            spherical_longitude, spherical_latitude, radius = geodetic_to_spherical(
+                longitude, latitude, height
+            )
+            npt.assert_allclose(spherical_longitude, longitude, rtol=rtol)
+            npt.assert_allclose(spherical_latitude, latitude, rtol=rtol)
             npt.assert_allclose(radius, height + ellipsoid.semiminor_axis, rtol=rtol)
 
 
@@ -65,11 +80,16 @@ def test_spherical_to_geodetic_with_spherical_ellipsoid():
     spherical_ellipsoid = ReferenceEllipsoid(
         "unit_sphere", sphere_radius, np.infty, 0, 0
     )
+    size = 5
+    spherical_longitude = np.linspace(0, 180, size)
     with set_ellipsoid(spherical_ellipsoid):
-        geocentric_latitude = np.linspace(-90, 90, 5)
-        radius = np.linspace(0.8, 1.2, 5)
-        latitude, height = spherical_to_geodetic(geocentric_latitude, radius)
-        npt.assert_allclose(geocentric_latitude, latitude, rtol=rtol)
+        spherical_latitude = np.linspace(-90, 90, size)
+        radius = np.linspace(0.8, 1.2, size)
+        longitude, latitude, height = spherical_to_geodetic(
+            spherical_longitude, spherical_latitude, radius
+        )
+        npt.assert_allclose(spherical_longitude, longitude, rtol=rtol)
+        npt.assert_allclose(spherical_latitude, latitude, rtol=rtol)
         npt.assert_allclose(radius, sphere_radius + height, rtol=rtol)
 
 
@@ -77,13 +97,17 @@ def test_spherical_to_geodetic_on_equator():
     "Test spherical to geodetic coordinates conversion on equator."
     rtol = 1e-10
     size = 5
-    geocentric_latitude = np.zeros(size)
+    spherical_latitude = np.zeros(size)
     for ellipsoid_name in KNOWN_ELLIPSOIDS:
         with set_ellipsoid(ellipsoid_name):
             ellipsoid = get_ellipsoid()
+            spherical_longitude = np.linspace(0, 180, size)
             radius = np.linspace(-1e4, 1e4, size) + ellipsoid.semimajor_axis
-            latitude, height = spherical_to_geodetic(geocentric_latitude, radius)
-            npt.assert_allclose(geocentric_latitude, latitude, rtol=rtol)
+            longitude, latitude, height = spherical_to_geodetic(
+                spherical_longitude, spherical_latitude, radius
+            )
+            npt.assert_allclose(spherical_longitude, longitude, rtol=rtol)
+            npt.assert_allclose(spherical_latitude, latitude, rtol=rtol)
             npt.assert_allclose(radius, height + ellipsoid.semimajor_axis, rtol=rtol)
 
 
@@ -91,29 +115,31 @@ def test_spherical_to_geodetic_on_poles():
     "Test spherical to geodetic coordinates conversion on poles."
     rtol = 1e-10
     size = 5
-    geocentric_latitude = np.array([90.0] * size + [-90.0] * size)
+    spherical_longitude = np.hstack([np.linspace(0, 180, size)] * 2)
+    spherical_latitude = np.array([90.0] * size + [-90.0] * size)
     for ellipsoid_name in KNOWN_ELLIPSOIDS:
         with set_ellipsoid(ellipsoid_name):
             ellipsoid = get_ellipsoid()
             radius = np.hstack(
                 [np.linspace(-1e4, 1e4, size) + ellipsoid.semiminor_axis] * 2
             )
-            latitude, height = spherical_to_geodetic(geocentric_latitude, radius)
-            npt.assert_allclose(geocentric_latitude, latitude, rtol=rtol)
+            longitude, latitude, height = spherical_to_geodetic(
+                spherical_longitude, spherical_latitude, radius
+            )
+            npt.assert_allclose(spherical_longitude, longitude, rtol=rtol)
+            npt.assert_allclose(spherical_latitude, latitude, rtol=rtol)
             npt.assert_allclose(radius, height + ellipsoid.semiminor_axis, rtol=rtol)
 
 
-def test_identity():
+def test_spherical_to_geodetic_identity():
     "Test if geodetic_to_spherical and spherical_to_geodetic is the identity operator"
     rtol = 1e-10
+    longitude = np.linspace(0, 350, 36)
     latitude = np.linspace(-90, 90, 19)
     height = np.linspace(-1e4, 1e4, 8)
-    latitude, height = np.meshgrid(latitude, height)
+    coordinates = np.meshgrid(longitude, latitude, height)
     for ellipsoid_name in KNOWN_ELLIPSOIDS:
         with set_ellipsoid(ellipsoid_name):
-            geocentric_latitude, radius = geodetic_to_spherical(latitude, height)
-            converted_latitude, converted_height = spherical_to_geodetic(
-                geocentric_latitude, radius
-            )
-            npt.assert_allclose(latitude, converted_latitude, rtol=rtol)
-            npt.assert_allclose(height, converted_height, rtol=rtol)
+            spherical_coordinates = geodetic_to_spherical(*coordinates)
+            reconverted_coordinates = spherical_to_geodetic(*spherical_coordinates)
+            npt.assert_allclose(coordinates, reconverted_coordinates, rtol=rtol)