From 50877d84776bfecb656e27d325def7f83a16554c Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Mon, 15 Jul 2019 23:58:27 -0300 Subject: [PATCH] Add longitude_continuity func to support slicing longitudes (#181) 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 --- verde/__init__.py | 1 + verde/coordinates.py | 135 +++++++++++++++++++++++++++++++- verde/tests/test_coordinates.py | 107 +++++++++++++++++++++++++ 3 files changed, 242 insertions(+), 1 deletion(-) diff --git a/verde/__init__.py b/verde/__init__.py index fc96f69c3..b11cac4c8 100644 --- a/verde/__init__.py +++ b/verde/__init__.py @@ -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 diff --git a/verde/coordinates.py b/verde/coordinates.py index f64724f3f..945b79ea8 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." + ) diff --git a/verde/tests/test_coordinates.py b/verde/tests/test_coordinates.py index d62b4e8ca..33cffdad6 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 @@ 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)