Skip to content

Commit

Permalink
Add longitude_continuity func to support slicing longitudes (#181)
Browse files Browse the repository at this point in the history
The function moves a set of longitudinal coordinates either to the `[0, 360)`
or `[-180, 180)` degrees intervals to ensure that longitude is continuous
throughout the array (no 360 -> 0 jump). Optionally change a given region
to match the longitude values.
This way, longitude arrays can be sliced to a sub region with simple `>= <=`
comparisons.

Fixes #171
  • Loading branch information
santisoler authored and leouieda committed Jul 16, 2019
1 parent 96bfb1c commit 50877d8
Show file tree
Hide file tree
Showing 3 changed files with 242 additions and 1 deletion.
1 change: 1 addition & 0 deletions verde/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@
get_region,
pad_region,
project_region,
longitude_continuity,
)
from .mask import distance_mask
from .utils import variance_to_weights, maxabs, grid_to_table
Expand Down
135 changes: 134 additions & 1 deletion verde/coordinates.py
Original file line number Diff line number Diff line change
Expand Up @@ -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(
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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."
)
107 changes: 107 additions & 0 deletions verde/tests/test_coordinates.py
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
"""
Test the coordinate generation functions
"""
import numpy as np
import numpy.testing as npt
import pytest

Expand All @@ -9,6 +10,7 @@
spacing_to_shape,
profile_coordinates,
grid_coordinates,
longitude_continuity,
)


Expand Down Expand Up @@ -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)

0 comments on commit 50877d8

Please sign in to comment.