Repository: fatiando/verde (python)

Issue:
Improve support for geographic coordinates
**Description of the desired feature**

Some functions from `verde` are intended to work with both projected and geographic coordinates.
Nevertheless, working with the latter imposes some requirements.
For example, one must ensure continuity around the globe: setting a grid between longitudes 170° and -170° or between 350° and 10° is not possible at the time because numerically the `east` coordinate is lower than the `west` one.
One way to fix this is to add a `latlon=True` argument to those functions in order to make clear that we are passing geographic coordinates.

I'll try to make a complete list of the functions that needs to be enhanced:
- `inside()`
- `grid_coordinates()`
- `scatter_coordinates()`
- `profile_coordinates()`
- `get_region()`
- `pad_region()`

**Are you willing to help implement and maintain this feature?** Yes


Pull request:
Add enhanced support for geographic coordinates
Add `latlon_continuity()` function that moves a set of longitudinal coordinates either
to the `[0, 360)` or `[-180, 180)` degrees intervals, depending which one is more suited for
that particular set of coordinates.
It also allow to modify an extra set of coordinates based on the decision made on the
first set.
Modify the `inside()` function to enhance how the inside points are selected, even if
their longitudinal coordinates are not defined on the same degree interval. 

Fixes #171 

**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` and `verde/__init__.py`.
- [ ] Write detailed docstrings for all functions/classes/methods.
- [ ] If adding new functionality, add an example to the docstring, gallery, and/or tutorials.


Patch:
diff --git a/verde/__init__.py b/verde/__init__.py
index fc96f69..b11cac4 100644
--- a/verde/__init__.py
+++ b/verde/__init__.py
@@ -11,6 +11,7 @@ from .coordinates import (
     get_region,
     pad_region,
     project_region,
+    longitude_continuity,
 )
 from .mask import distance_mask
 from .utils import variance_to_weights, maxabs, grid_to_table
diff --git a/verde/coordinates.py b/verde/coordinates.py
index f64724f..945b79e 100644
--- a/verde/coordinates.py
+++ b/verde/coordinates.py
@@ -32,7 +32,9 @@ def check_region(region):
     w, e, s, n = region
     if w > e:
         raise ValueError(
-            "Invalid region '{}' (W, E, S, N). Must have W =< E.".format(region)
+            "Invalid region '{}' (W, E, S, N). Must have W =< E. ".format(region)
+            + "If working with geographic coordinates, don't forget to match geographic"
+            + " region with coordinates using 'verde.longitude_continuity'."
         )
     if s > n:
         raise ValueError(
@@ -635,6 +637,17 @@ def inside(coordinates, region):
      [False  True  True]
      [False False False]]
 
+    Geographic coordinates are also supported using :func:`verde.longitude_continuity`:
+
+    >>> from verde import longitude_continuity
+    >>> east, north = grid_coordinates([0, 350, -20, 20], spacing=10)
+    >>> region = [-10, 10, -10, 10]
+    >>> are_inside = inside(*longitude_continuity([east, north], region))
+    >>> print(east[are_inside])
+    [  0.  10. 350.   0.  10. 350.   0.  10. 350.]
+    >>> print(north[are_inside])
+    [-10. -10. -10.   0.   0.   0.  10.  10.  10.]
+
     """
     check_region(region)
     w, e, s, n = region
@@ -751,3 +764,123 @@ def block_split(coordinates, spacing=None, adjust="spacing", region=None, shape=
     tree = kdtree(block_coords)
     labels = tree.query(np.transpose(n_1d_arrays(coordinates, 2)))[1]
     return block_coords, labels
+
+
+def longitude_continuity(coordinates, region):
+    """
+    Modify coordinates and region boundaries to ensure longitude continuity.
+
+    Longitudinal boundaries of the region are moved to the ``[0, 360)`` or ``[-180, 180)``
+    degrees interval depending which one is better suited for that specific region.
+
+    Parameters
+    ----------
+    coordinates : list or array
+        Set of geographic coordinates that will be moved to the same degrees
+        interval as the one of the modified region.
+    region : list or array
+        List or array containing the boundary coordinates `w`, `e`, `s`, `n` of the
+        region in degrees.
+
+    Returns
+    -------
+    modified_coordinates : array
+        Modified set of extra geographic coordinates.
+    modified_region : array
+        List containing the modified boundary coordinates `w, `e`, `s`, `n` of the
+        region.
+
+    Examples
+    --------
+
+    >>> # Modify region with west > east
+    >>> w, e, s, n = 350, 10, -10, 10
+    >>> print(longitude_continuity(coordinates=None, region=[w, e, s, n]))
+    [-10  10 -10  10]
+    >>> # Modify region and extra coordinates
+    >>> from verde import grid_coordinates
+    >>> region = [-70, -60, -40, -30]
+    >>> coordinates = grid_coordinates([270, 320, -50, -20], spacing=5)
+    >>> [longitude, latitude], region = longitude_continuity(coordinates, region)
+    >>> print(region)
+    [290 300 -40 -30]
+    >>> print(longitude.min(), longitude.max())
+    270.0 320.0
+    >>> # Another example
+    >>> region = [-20, 20, -20, 20]
+    >>> coordinates = grid_coordinates([0, 350, -90, 90], spacing=10)
+    >>> [longitude, latitude], region = longitude_continuity(coordinates, region)
+    >>> print(region)
+    [-20  20 -20  20]
+    >>> print(longitude.min(), longitude.max())
+    -180.0 170.0
+    """
+    # Get longitudinal boundaries and check region
+    w, e, s, n = region[:4]
+    # Run sanity checks for region
+    _check_geographic_region([w, e, s, n])
+    # Check if region is defined all around the globe
+    all_globe = np.allclose(abs(e - w), 360)
+    # Move coordinates to [0, 360)
+    interval_360 = True
+    w = w % 360
+    e = e % 360
+    # Move west=0 and east=360 if region longitudes goes all around the globe
+    if all_globe:
+        w, e = 0, 360
+    # Check if the [-180, 180) interval is better suited
+    if w > e:
+        interval_360 = False
+        e = ((e + 180) % 360) - 180
+        w = ((w + 180) % 360) - 180
+    region = np.array(region)
+    region[:2] = w, e
+    # Modify extra coordinates if passed
+    if coordinates:
+        # Run sanity checks for coordinates
+        _check_geographic_coordinates(coordinates)
+        longitude = coordinates[0]
+        if interval_360:
+            longitude = longitude % 360
+        else:
+            longitude = ((longitude + 180) % 360) - 180
+        coordinates = np.array(coordinates)
+        coordinates[0] = longitude
+        return coordinates, region
+    return region
+
+
+def _check_geographic_coordinates(coordinates):
+    "Check if geographic coordinates are within accepted degrees intervals"
+    longitude, latitude = coordinates[:2]
+    if np.any(longitude > 360) or np.any(longitude < -180):
+        raise ValueError(
+            "Invalid longitude coordinates. They should be < 360 and > -180 degrees."
+        )
+    if np.any(latitude > 90) or np.any(latitude < -90):
+        raise ValueError(
+            "Invalid latitude coordinates. They should be < 90 and > -90 degrees."
+        )
+
+
+def _check_geographic_region(region):
+    "Check if region in geographic coordinates are within accepted degree intervals"
+    w, e, s, n = region[:4]
+    # Check if coordinates are within accepted degrees intervals
+    if np.any(np.array([w, e]) > 360) or np.any(np.array([w, e]) < -180):
+        raise ValueError(
+            "Invalid region '{}' (W, E, S, N). ".format(region)
+            + "Longitudinal coordinates should be < 360 and > -180 degrees."
+        )
+    if np.any(np.array([s, n]) > 90) or np.any(np.array([s, n]) < -90):
+        raise ValueError(
+            "Invalid region '{}' (W, E, S, N). ".format(region)
+            + "Latitudinal coordinates should be < 90 and > -90 degrees."
+        )
+    # Check if longitude boundaries do not involve more than one spin around the globe
+    if abs(e - w) > 360:
+        raise ValueError(
+            "Invalid region '{}' (W, E, S, N). ".format(region)
+            + "East and West boundaries must not be separated by an angle greater "
+            + "than 360 degrees."
+        )


Test patch:
diff --git a/verde/tests/test_coordinates.py b/verde/tests/test_coordinates.py
index d62b4e8..33cffda 100644
--- a/verde/tests/test_coordinates.py
+++ b/verde/tests/test_coordinates.py
@@ -1,6 +1,7 @@
 """
 Test the coordinate generation functions
 """
+import numpy as np
 import numpy.testing as npt
 import pytest
 
@@ -9,6 +10,7 @@ from ..coordinates import (
     spacing_to_shape,
     profile_coordinates,
     grid_coordinates,
+    longitude_continuity,
 )
 
 
@@ -92,3 +94,108 @@ def test_profile_coordiantes_fails():
         profile_coordinates((0, 1), (1, 2), size=0)
     with pytest.raises(ValueError):
         profile_coordinates((0, 1), (1, 2), size=-10)
+
+
+def test_longitude_continuity():
+    "Test continuous boundary conditions in geographic coordinates."
+    # Define longitude coordinates around the globe for [0, 360) and [-180, 180)
+    longitude_360 = np.linspace(0, 350, 36)
+    longitude_180 = np.hstack((longitude_360[:18], longitude_360[18:] - 360))
+    latitude = np.linspace(-90, 90, 36)
+    s, n = -90, 90
+    # Check w, e in [0, 360)
+    w, e = 10.5, 20.3
+    for longitude in [longitude_360, longitude_180]:
+        coordinates = [longitude, latitude]
+        coordinates_new, region_new = longitude_continuity(coordinates, (w, e, s, n))
+        w_new, e_new = region_new[:2]
+        assert w_new == w
+        assert e_new == e
+        npt.assert_allclose(coordinates_new[0], longitude_360)
+    # Check w, e in [-180, 180)
+    w, e = -20, 20
+    for longitude in [longitude_360, longitude_180]:
+        coordinates = [longitude, latitude]
+        coordinates_new, region_new = longitude_continuity(coordinates, (w, e, s, n))
+        w_new, e_new = region_new[:2]
+        assert w_new == -20
+        assert e_new == 20
+        npt.assert_allclose(coordinates_new[0], longitude_180)
+    # Check region around the globe
+    for w, e in [[0, 360], [-180, 180], [-20, 340]]:
+        for longitude in [longitude_360, longitude_180]:
+            coordinates = [longitude, latitude]
+            coordinates_new, region_new = longitude_continuity(
+                coordinates, (w, e, s, n)
+            )
+            w_new, e_new = region_new[:2]
+            assert w_new == 0
+            assert e_new == 360
+            npt.assert_allclose(coordinates_new[0], longitude_360)
+    # Check w == e
+    w, e = 20, 20
+    for longitude in [longitude_360, longitude_180]:
+        coordinates = [longitude, latitude]
+        coordinates_new, region_new = longitude_continuity(coordinates, (w, e, s, n))
+        w_new, e_new = region_new[:2]
+        assert w_new == 20
+        assert e_new == 20
+        npt.assert_allclose(coordinates_new[0], longitude_360)
+    # Check angle greater than 180
+    w, e = 0, 200
+    for longitude in [longitude_360, longitude_180]:
+        coordinates = [longitude, latitude]
+        coordinates_new, region_new = longitude_continuity(coordinates, (w, e, s, n))
+        w_new, e_new = region_new[:2]
+        assert w_new == 0
+        assert e_new == 200
+        npt.assert_allclose(coordinates_new[0], longitude_360)
+    w, e = -160, 160
+    for longitude in [longitude_360, longitude_180]:
+        coordinates = [longitude, latitude]
+        coordinates_new, region_new = longitude_continuity(coordinates, (w, e, s, n))
+        w_new, e_new = region_new[:2]
+        assert w_new == -160
+        assert e_new == 160
+        npt.assert_allclose(coordinates_new[0], longitude_180)
+
+
+def test_invalid_geographic_region():
+    "Check if passing invalid region to longitude_continuity raises a ValueError"
+    # Region with latitude over boundaries
+    w, e = -10, 10
+    for s, n in [[-200, 90], [-90, 200]]:
+        with pytest.raises(ValueError):
+            longitude_continuity(None, [w, e, s, n])
+    # Region with longitude over boundaries
+    s, n = -10, 10
+    for w, e in [[-200, 0], [0, 380]]:
+        with pytest.raises(ValueError):
+            longitude_continuity(None, [w, e, s, n])
+    # Region with longitudinal difference greater than 360
+    w, e, s, n = -180, 200, -10, 10
+    with pytest.raises(ValueError):
+        longitude_continuity(None, [w, e, s, n])
+
+
+def test_invalid_geographic_coordinates():
+    "Check if passing invalid coordinates to longitude_continuity raises a ValueError"
+    boundaries = [0, 360, -90, 90]
+    spacing = 10
+    region = [-20, 20, -20, 20]
+    # Region with longitude point over boundaries
+    longitude, latitude = grid_coordinates(boundaries, spacing=spacing)
+    longitude[0] = -200
+    with pytest.raises(ValueError):
+        longitude_continuity([longitude, latitude], region)
+    longitude[0] = 400
+    with pytest.raises(ValueError):
+        longitude_continuity([longitude, latitude], region)
+    # Region with latitude point over boundaries
+    longitude, latitude = grid_coordinates(boundaries, spacing=spacing)
+    latitude[0] = -100
+    with pytest.raises(ValueError):
+        longitude_continuity([longitude, latitude], region)
+    latitude[0] = 100
+    with pytest.raises(ValueError):
+        longitude_continuity([longitude, latitude], region)