# Copyright (c) 2017 The Verde Developers.
# Distributed under the terms of the BSD 3-Clause License.
# SPDX-License-Identifier: BSD-3-Clause
#
# This code is part of the Fatiando a Terra project (https://www.fatiando.org)
#
"""
Functions for generating and manipulating coordinates.
"""
import warnings
import numpy as np
from sklearn.utils import check_random_state
from .base.utils import check_coordinates, n_1d_arrays
from .utils import kdtree
def check_region(region):
"""
Check that the given region dimensions are valid.
For example, the west limit should not be greater than the east and there
must be exactly 4 values given.
Parameters
----------
region : list = [W, E, S, N]
The boundaries of a given region in Cartesian or geographic
coordinates.
Raises
------
ValueError
If the region doesn't have exactly 4 entries, W > E, or S > N.
"""
if len(region) != 4:
raise ValueError("Invalid region '{}'. Only 4 values allowed.".format(region))
w, e, s, n = region
if w > e:
raise ValueError(
"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(
"Invalid region '{}' (W, E, S, N). Must have S =< N.".format(region)
)
[docs]
def get_region(coordinates):
"""
Get the bounding region of the given coordinates.
Parameters
----------
coordinates : tuple of arrays
Arrays with the coordinates of each data point. Should be in the
following order: (easting, northing, vertical, ...). Only easting and
northing will be used, all subsequent coordinates will be ignored.
Returns
-------
region : tuple = (W, E, S, N)
The boundaries of a given region in Cartesian or geographic
coordinates.
Examples
--------
>>> coords = grid_coordinates((0, 1, -10, -6), shape=(10, 10))
>>> print(get_region(coords))
(0.0, 1.0, -10.0, -6.0)
"""
easting, northing = coordinates[:2]
region = (np.min(easting), np.max(easting), np.min(northing), np.max(northing))
return region
[docs]
def pad_region(region, pad):
"""
Extend the borders of a region by the given amount.
Parameters
----------
region : list = [W, E, S, N]
The boundaries of a given region in Cartesian or geographic
coordinates.
pad : float or tuple = (pad_north, pad_east)
The amount of padding to add to the region. If it's a single number,
add this to all boundaries of region equally. If it's a tuple of
numbers, then will add different padding to the North-South and
East-West dimensions.
Returns
-------
padded_region : list = [W, E, S, N]
The padded region.
Examples
--------
>>> pad_region((0, 1, -5, -3), 1)
(-1, 2, -6, -2)
>>> pad_region((0, 1, -5, -3), (3, 2))
(-2, 3, -8, 0)
"""
if np.isscalar(pad):
pad = (pad, pad)
w, e, s, n = region
padded = (w - pad[1], e + pad[1], s - pad[0], n + pad[0])
return padded
[docs]
def scatter_points(region, size, random_state=None, extra_coords=None):
"""
Generate the coordinates for a random scatter of points.
The points are drawn from a uniform distribution.
Parameters
----------
region : list = [W, E, S, N]
The boundaries of a given region in Cartesian or geographic
coordinates.
size : int
The number of points to generate.
random_state : numpy.random.RandomState or an int seed
A random number generator used to define the state of the random
permutations. Use a fixed seed to make sure computations are
reproducible. Use ``None`` to choose a seed automatically (resulting in
different numbers with each run).
extra_coords : None, scalar, or list
If not None, then value(s) of extra coordinate arrays to be generated.
These extra arrays will have the same *size* as the others but will
contain a constant value. Will generate an extra array per value given
in *extra_coords*. Use this to generate arrays of constant heights or
times, for example, that might be needed to evaluate a gridder.
Returns
-------
coordinates : tuple of arrays
Arrays with coordinates of each point in the grid. Each array contains
values for a dimension in the order: easting, northing, vertical, and
any extra dimensions given in *extra_coords*. All arrays will have the
specified *size*.
Examples
--------
>>> # We'll use a seed value will ensure that the same will be generated
>>> # every time.
>>> easting, northing = scatter_points((0, 10, -2, -1), 4, random_state=0)
>>> print(', '.join(['{:.4f}'.format(i) for i in easting]))
5.4881, 7.1519, 6.0276, 5.4488
>>> print(', '.join(['{:.4f}'.format(i) for i in northing]))
-1.5763, -1.3541, -1.5624, -1.1082
>>> easting, northing, height = scatter_points(
... (0, 10, -2, -1), 4, random_state=0, extra_coords=12
... )
>>> print(height)
[12. 12. 12. 12.]
>>> easting, northing, height, time = scatter_points(
... (0, 10, -2, -1), 4, random_state=0, extra_coords=[12, 1986])
>>> print(height)
[12. 12. 12. 12.]
>>> print(time)
[1986. 1986. 1986. 1986.]
See also
--------
grid_coordinates : Generate coordinates for each point on a regular grid
profile_coordinates : Coordinates for a profile between two points
"""
check_region(region)
random = check_random_state(random_state)
coordinates = []
for lower, upper in np.array(region).reshape((len(region) // 2, 2)):
coordinates.append(random.uniform(lower, upper, size))
if extra_coords is not None:
for value in np.atleast_1d(extra_coords):
coordinates.append(np.ones_like(coordinates[0]) * value)
return tuple(coordinates)
[docs]
def line_coordinates(
start, stop, size=None, spacing=None, adjust="spacing", pixel_register=False
):
"""
Generate evenly spaced points between two values.
Able to handle either specifying the number of points required (*size*) or
the size of the interval between points (*spacing*). If using *size*, the
output will be similar to using :func:`numpy.linspace`. When using
*spacing*, if the interval is not divisible by the desired spacing, either
the interval or the spacing will have to be adjusted. By default, the
spacing will be rounded to the nearest multiple. Optionally, the *stop*
value can be adjusted to fit the exact spacing given.
Parameters
----------
start : float
The starting value of the sequence.
stop : float
The end value of the sequence.
num : int or None
The number of points in the sequence. If None, *spacing* must be
provided.
spacing : float or None
The step size (interval) between points in the sequence. If None,
*size* must be provided.
adjust : {'spacing', 'region'}
Whether to adjust the spacing or the interval/region if required.
Ignored if *size* is given instead of *spacing*. Defaults to adjusting
the spacing.
pixel_register : bool
If True, the points will refer to the center of each interval (pixel)
instead of the boundaries. In practice, this means that there will be
one less element in the sequence if *spacing* is provided. If *size* is
provided, the requested number of elements is respected. Default is
False.
Returns
-------
sequence : array
The generated sequence of values.
Examples
--------
>>> # Lower printing precision to shorten this example
>>> import numpy as np; np.set_printoptions(precision=2, suppress=True)
>>> values = line_coordinates(0, 5, spacing=2.5)
>>> print(values.shape)
(3,)
>>> print(values)
[0. 2.5 5. ]
>>> print(line_coordinates(0, 10, size=5))
[ 0. 2.5 5. 7.5 10. ]
>>> print(line_coordinates(0, 10, spacing=2.5))
[ 0. 2.5 5. 7.5 10. ]
The spacing is adjusted to fit the interval by default but this can be
changed to adjusting the interval/region instead:
>>> print(line_coordinates(0, 10, spacing=2.4))
[ 0. 2.5 5. 7.5 10. ]
>>> print(line_coordinates(0, 10, spacing=2.4, adjust="region"))
[0. 2.4 4.8 7.2 9.6]
>>> print(line_coordinates(0, 10, spacing=2.6))
[ 0. 2.5 5. 7.5 10. ]
>>> print(line_coordinates(0, 10, spacing=2.6, adjust="region"))
[ 0. 2.6 5.2 7.8 10.4]
Optionally, return values at the center of the intervals instead of their
boundaries:
>>> print(line_coordinates(0, 10, spacing=2.5, pixel_register=True))
[1.25 3.75 6.25 8.75]
Notice that this produces one value less than the non-pixel registered
version. If using *size* instead of *spacing*, the number of values will be
*size* regardless and the spacing will therefore be different from the
non-pixel registered version:
>>> print(line_coordinates(0, 10, size=5, pixel_register=True))
[1. 3. 5. 7. 9.]
"""
if size is not None and spacing is not None:
raise ValueError("Both size and spacing provided. Only one is allowed.")
if size is None and spacing is None:
raise ValueError("Either a size or a spacing must be provided.")
if spacing is not None:
size, stop = spacing_to_size(start, stop, spacing, adjust)
elif pixel_register:
# Starts by generating grid-line registered coordinates and shifting
# them to the center of the pixel. Need 1 more point if given a size
# instead of spacing so that we can do that because we discard the last
# point when shifting the coordinates.
size = size + 1
values = np.linspace(start, stop, size)
if pixel_register:
values = values[:-1] + (values[1] - values[0]) / 2
return values
[docs]
def grid_coordinates(
region,
shape=None,
spacing=None,
adjust="spacing",
pixel_register=False,
extra_coords=None,
meshgrid=True,
):
"""
Generate the coordinates for each point on a regular grid.
The grid can be specified by either the number of points in each dimension
(the *shape*) or by the grid node spacing.
If the given region is not divisible by the desired spacing, either the
region or the spacing will have to be adjusted. By default, the spacing
will be rounded to the nearest multiple. Optionally, the East and North
boundaries of the region can be adjusted to fit the exact spacing given.
See the examples below.
Parameters
----------
region : list = [W, E, S, N]
The boundaries of a given region in Cartesian or geographic
coordinates.
shape : tuple = (n_north, n_east) or None
The number of points in the South-North and West-East directions,
respectively.
spacing : float, tuple = (s_north, s_east), or None
The grid spacing in the South-North and West-East directions,
respectively. A single value means that the spacing is equal in both
directions.
adjust : {'spacing', 'region'}
Whether to adjust the spacing or the region if required. Ignored if
*shape* is given instead of *spacing*. Defaults to adjusting the
spacing.
pixel_register : bool
If True, the coordinates will refer to the center of each grid pixel
instead of the grid lines. In practice, this means that there will be
one less element per dimension of the grid when compared to grid line
registered (only if given *spacing* and not *shape*). Default is False.
extra_coords : None, scalar, or list
If not None, then value(s) of extra coordinate arrays to be generated.
These extra arrays will have the same *shape* as the others but will
contain a constant value. Will generate an extra array per value given
in *extra_coords*. Use this to generate arrays of constant heights or
times, for example, that might be needed to evaluate a gridder.
meshgrid : bool
If True, will call :func:`numpy.meshgrid` on the generated coordinates
and return 2D arrays (useful if you need the coordinate values for
every single point on a grid). Otherwise, will return 1D coordinate
arrays (useful if you're making a :class:`xarray.DataArray` or looping
over grid coordinates with two ``for`` loops). Passing False to
*meshgrid* is **incompatible with *extra_coords* and an exception will
be raised** if used together (:class:`ValueError`).
Returns
-------
coordinates : tuple of arrays
Arrays with coordinates of each point in the grid. Each array contains
values for a dimension in the order: easting, northing, vertical, and
any extra dimensions given in *extra_coords*. All arrays will have the
specified *shape*.
Examples
--------
>>> east, north = grid_coordinates(region=(0, 5, 0, 10), shape=(5, 3))
>>> print(east.shape, north.shape)
(5, 3) (5, 3)
>>> # Lower printing precision to shorten this example
>>> import numpy as np; np.set_printoptions(precision=1, suppress=True)
>>> print(east)
[[0. 2.5 5. ]
[0. 2.5 5. ]
[0. 2.5 5. ]
[0. 2.5 5. ]
[0. 2.5 5. ]]
>>> print(north)
[[ 0. 0. 0. ]
[ 2.5 2.5 2.5]
[ 5. 5. 5. ]
[ 7.5 7.5 7.5]
[10. 10. 10. ]]
The grid can also be specified using the spacing between points instead of
the shape:
>>> east, north = grid_coordinates(region=(0, 5, 0, 10), spacing=2.5)
>>> print(east.shape, north.shape)
(5, 3) (5, 3)
>>> print(east)
[[0. 2.5 5. ]
[0. 2.5 5. ]
[0. 2.5 5. ]
[0. 2.5 5. ]
[0. 2.5 5. ]]
>>> print(north)
[[ 0. 0. 0. ]
[ 2.5 2.5 2.5]
[ 5. 5. 5. ]
[ 7.5 7.5 7.5]
[10. 10. 10. ]]
If you don't need the 2D arrays, then use ``meshgrid=False``:
>>> east, north = grid_coordinates(
... region=(0, 5, 0, 10), spacing=2.5, meshgrid=False,
... )
>>> print(east.shape, north.shape)
(3,) (5,)
>>> print(east)
[0. 2.5 5. ]
>>> print(north)
[ 0. 2.5 5. 7.5 10. ]
The spacing can be different for northing and easting, respectively:
>>> east, north = grid_coordinates(region=(-5, 1, 0, 10), spacing=(2.5, 1))
>>> print(east.shape, north.shape)
(5, 7) (5, 7)
>>> print(east)
[[-5. -4. -3. -2. -1. 0. 1.]
[-5. -4. -3. -2. -1. 0. 1.]
[-5. -4. -3. -2. -1. 0. 1.]
[-5. -4. -3. -2. -1. 0. 1.]
[-5. -4. -3. -2. -1. 0. 1.]]
>>> print(north)
[[ 0. 0. 0. 0. 0. 0. 0. ]
[ 2.5 2.5 2.5 2.5 2.5 2.5 2.5]
[ 5. 5. 5. 5. 5. 5. 5. ]
[ 7.5 7.5 7.5 7.5 7.5 7.5 7.5]
[10. 10. 10. 10. 10. 10. 10. ]]
If the region can't be divided into the desired spacing, the spacing will
be adjusted to conform to the region:
>>> east, north = grid_coordinates(region=(-5, 0, 0, 5), spacing=2.6)
>>> print(east.shape, north.shape)
(3, 3) (3, 3)
>>> print(east)
[[-5. -2.5 0. ]
[-5. -2.5 0. ]
[-5. -2.5 0. ]]
>>> print(north)
[[0. 0. 0. ]
[2.5 2.5 2.5]
[5. 5. 5. ]]
>>> east, north = grid_coordinates(region=(-5, 0, 0, 5), spacing=2.4)
>>> print(east.shape, north.shape)
(3, 3) (3, 3)
>>> print(east)
[[-5. -2.5 0. ]
[-5. -2.5 0. ]
[-5. -2.5 0. ]]
>>> print(north)
[[0. 0. 0. ]
[2.5 2.5 2.5]
[5. 5. 5. ]]
You can choose to adjust the East and North boundaries of the region
instead:
>>> east, north = grid_coordinates(region=(-5, 0, 0, 5), spacing=2.6,
... adjust='region')
>>> print(east.shape, north.shape)
(3, 3) (3, 3)
>>> print(east)
[[-5. -2.4 0.2]
[-5. -2.4 0.2]
[-5. -2.4 0.2]]
>>> print(north)
[[0. 0. 0. ]
[2.6 2.6 2.6]
[5.2 5.2 5.2]]
>>> east, north = grid_coordinates(region=(-5, 0, 0, 5), spacing=2.4,
... adjust='region')
>>> print(east.shape, north.shape)
(3, 3) (3, 3)
>>> print(east)
[[-5. -2.6 -0.2]
[-5. -2.6 -0.2]
[-5. -2.6 -0.2]]
>>> print(north)
[[0. 0. 0. ]
[2.4 2.4 2.4]
[4.8 4.8 4.8]]
We can optionally generate coordinates for the center of each grid pixel
instead of the corner (default):
>>> east, north = grid_coordinates(region=(0, 5, 0, 10), spacing=2.5,
... pixel_register=True)
>>> # Raise the printing precision for this example
>>> np.set_printoptions(precision=2, suppress=True)
>>> # Notice that the shape is 1 less than when pixel_register=False
>>> print(east.shape, north.shape)
(4, 2) (4, 2)
>>> print(east)
[[1.25 3.75]
[1.25 3.75]
[1.25 3.75]
[1.25 3.75]]
>>> print(north)
[[1.25 1.25]
[3.75 3.75]
[6.25 6.25]
[8.75 8.75]]
>>> east, north = grid_coordinates(region=(0, 5, 0, 10), shape=(4, 2),
... pixel_register=True)
>>> print(east)
[[1.25 3.75]
[1.25 3.75]
[1.25 3.75]
[1.25 3.75]]
>>> print(north)
[[1.25 1.25]
[3.75 3.75]
[6.25 6.25]
[8.75 8.75]]
Generate arrays for other coordinates that have a constant value:
>>> east, north, height = grid_coordinates(
... region=(0, 5, 0, 10), spacing=2.5, extra_coords=57
... )
>>> print(east.shape, north.shape, height.shape)
(5, 3) (5, 3) (5, 3)
>>> print(height)
[[57. 57. 57.]
[57. 57. 57.]
[57. 57. 57.]
[57. 57. 57.]
[57. 57. 57.]]
>>> east, north, height, time = grid_coordinates(
... region=(0, 5, 0, 10), spacing=2.5, extra_coords=[57, 0.1]
... )
>>> print(east.shape, north.shape, height.shape, time.shape)
(5, 3) (5, 3) (5, 3) (5, 3)
>>> print(height)
[[57. 57. 57.]
[57. 57. 57.]
[57. 57. 57.]
[57. 57. 57.]
[57. 57. 57.]]
>>> print(time)
[[0.1 0.1 0.1]
[0.1 0.1 0.1]
[0.1 0.1 0.1]
[0.1 0.1 0.1]
[0.1 0.1 0.1]]
See also
--------
scatter_points : Generate the coordinates for a random scatter of points
profile_coordinates : Coordinates for a profile between two points
line_coordinates: Generate evenly spaced points between two values
"""
check_region(region)
if shape is not None and spacing is not None:
raise ValueError("Both grid shape and spacing provided. Only one is allowed.")
if shape is None and spacing is None:
raise ValueError("Either a grid shape or a spacing must be provided.")
if shape is None:
shape = (None, None)
# Make sure the spacing is a tuple of 2 numbers
spacing = np.atleast_1d(spacing)
if len(spacing) == 1:
spacing = (spacing[0], spacing[0])
elif len(spacing) > 2:
raise ValueError(f"Only two values allowed for grid spacing: {spacing}")
else:
spacing = (None, None)
east = line_coordinates(
region[0],
region[1],
size=shape[1],
spacing=spacing[1],
adjust=adjust,
pixel_register=pixel_register,
)
north = line_coordinates(
region[2],
region[3],
size=shape[0],
spacing=spacing[0],
adjust=adjust,
pixel_register=pixel_register,
)
coordinates = [east, north]
if meshgrid:
coordinates = list(np.meshgrid(east, north))
if extra_coords is not None:
if not meshgrid:
raise ValueError(
"Using 'meshgrid=False' and 'extra_coords' is 'verde.grid_coordinates' is not allowed."
)
for value in np.atleast_1d(extra_coords):
coordinates.append(np.ones_like(coordinates[0]) * value)
return tuple(coordinates)
def spacing_to_size(start, stop, spacing, adjust):
"""
Convert a spacing to the number of points between start and stop.
Takes into account if the spacing or the interval needs to be adjusted.
Parameters
----------
start : float
The starting value of the sequence.
stop : float
The end value of the sequence.
spacing : float
The step size (interval) between points in the sequence.
adjust : {'spacing', 'region'}
Whether to adjust the spacing or the interval/region if required.
Defaults to adjusting the spacing.
Returns
-------
size : int
The number of points between start and stop.
stop : float
The end of the interval, which may or may not have been adjusted.
"""
if adjust not in ["spacing", "region"]:
raise ValueError(
"Invalid value for *adjust* '{}'. Should be 'spacing' or 'region'".format(
adjust
)
)
# Add 1 to get the number of nodes, not segments
size = int(round((stop - start) / spacing)) + 1
# If the spacing >= 2 * (stop - start), it rounds to zero so we'd be
# generating a single point, which isn't equivalent to adjusting the
# spacing or the region. To get the appropriate behaviour of decreasing the
# spacing until it fits the region or increasing the region until it fits
# at least 1 spacing, we need to always round to at least 1 in the code
# above.
if size == 1:
size += 1
if adjust == "region":
# The size is the same but we adjust the interval so that the spacing
# isn't altered when we do the linspace.
stop = start + (size - 1) * spacing
return size, stop
def shape_to_spacing(region, shape, pixel_register=False):
"""
Calculate the spacing of a grid given region and shape.
Parameters
----------
region : list = [W, E, S, N]
The boundaries of a given region in Cartesian or geographic
coordinates.
shape : tuple = (n_north, n_east) or None
The number of points in the South-North and West-East directions,
respectively.
pixel_register : bool
If True, the coordinates will refer to the center of each grid pixel
instead of the grid lines. In practice, this means that there will be
one less element per dimension of the grid when compared to grid line
registered (only if given *spacing* and not *shape*). Default is False.
Returns
-------
spacing : tuple = (s_north, s_east)
The grid spacing in the South-North and West-East directions,
respectively.
Examples
--------
>>> spacing = shape_to_spacing([0, 10, -5, 1], (7, 11))
>>> print("{:.1f}, {:.1f}".format(*spacing))
1.0, 1.0
>>> spacing = shape_to_spacing([0, 10, -5, 1], (14, 11))
>>> print("{:.1f}, {:.1f}".format(*spacing))
0.5, 1.0
>>> spacing = shape_to_spacing([0, 10, -5, 1], (7, 21))
>>> print("{:.1f}, {:.1f}".format(*spacing))
1.0, 0.5
>>> spacing = shape_to_spacing(
... [-0.5, 10.5, -5.5, 1.5], (7, 11), pixel_register=True,
... )
>>> print("{:.1f}, {:.1f}".format(*spacing))
1.0, 1.0
>>> spacing = shape_to_spacing(
... [-0.25, 10.25, -5.5, 1.5], (7, 21), pixel_register=True,
... )
>>> print("{:.1f}, {:.1f}".format(*spacing))
1.0, 0.5
"""
spacing = []
for i, n_points in enumerate(reversed(shape)):
if not pixel_register:
n_points -= 1
spacing.append((region[2 * i + 1] - region[2 * i]) / n_points)
return tuple(reversed(spacing))
[docs]
def profile_coordinates(point1, point2, size, extra_coords=None):
"""
Coordinates for a profile along a straight line between two points.
Parameters
----------
point1 : tuple or list
``(easting, northing)`` West-East and South-North coordinates of the
first point, respectively.
point2 : tuple or list
``(easting, northing)`` West-East and South-North coordinates of the
second point, respectively.
size : int
Number of points to sample along the line.
extra_coords : None, scalar, or list
If not None, then value(s) of extra coordinate arrays to be generated.
These extra arrays will have the same *size* as the others but will
contain a constant value. Will generate an extra array per value given
in *extra_coords*. Use this to generate arrays of constant heights or
times, for example, that might be needed to evaluate a gridder.
Returns
-------
coordinates, distances : tuple and 1d array
The coordinates of points along the straight line and the distances
from the first point.
Examples
--------
>>> (east, north), dist = profile_coordinates((1, 10), (1, 20), size=11)
>>> print('easting:', ', '.join('{:.1f}'.format(i) for i in east))
easting: 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0
>>> print('northing:', ', '.join('{:.1f}'.format(i) for i in north))
northing: 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0, 20.0
>>> print('distance:', ', '.join('{:.1f}'.format(i) for i in dist))
distance: 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0
>>> (east, north, height), dist = profile_coordinates(
... (1, 10), (1, 20), size=11, extra_coords=35)
>>> print(height)
[35. 35. 35. 35. 35. 35. 35. 35. 35. 35. 35.]
>>> (east, north, height, time), dist = profile_coordinates(
... (1, 10), (1, 20), size=11, extra_coords=[35, 0.1])
>>> print(height)
[35. 35. 35. 35. 35. 35. 35. 35. 35. 35. 35.]
>>> print(time)
[0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1]
See also
--------
scatter_points : Generate the coordinates for a random scatter of points
grid_coordinates : Generate coordinates for each point on a regular grid
"""
if size <= 0:
raise ValueError("Invalid profile size '{}'. Must be > 0.".format(size))
diffs = [i - j for i, j in zip(point2, point1)]
separation = np.hypot(*diffs)
distances = np.linspace(0, separation, size)
angle = np.arctan2(*reversed(diffs))
coordinates = [
point1[0] + distances * np.cos(angle),
point1[1] + distances * np.sin(angle),
]
if extra_coords is not None:
for value in np.atleast_1d(extra_coords):
coordinates.append(np.ones_like(coordinates[0]) * value)
return tuple(coordinates), distances
[docs]
def inside(coordinates, region):
"""
Determine which points fall inside a given region.
Points at the boundary are counted as being outsize.
Parameters
----------
coordinates : tuple of arrays
Arrays with the coordinates of each data point. Should be in the
following order: (easting, northing, vertical, ...). Only easting and
northing will be used, all subsequent coordinates will be ignored.
region : list = [W, E, S, N]
The boundaries of a given region in Cartesian or geographic
coordinates.
Returns
-------
are_inside : array of booleans
An array of booleans with the same shape as the input coordinate
arrays. Will be ``True`` if the respective coordinates fall inside the
area, ``False`` otherwise.
Examples
--------
>>> import numpy as np
>>> east = np.array([1, 2, 3, 4, 5, 6])
>>> north = np.array([10, 11, 12, 13, 14, 15])
>>> region = [2.5, 5.5, 12, 15]
>>> print(inside((east, north), region))
[False False True True True False]
>>> # This also works for 2D-arrays
>>> east = np.array([[1, 1, 1],
... [2, 2, 2],
... [3, 3, 3]])
>>> north = np.array([[5, 7, 9],
... [5, 7, 9],
... [5, 7, 9]])
>>> region = [0.5, 2.5, 6, 9]
>>> print(inside((east, north), region))
[[False True True]
[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
easting, northing = coordinates[:2]
# Allocate temporary arrays to minimize memory allocation overhead
out = np.empty_like(easting, dtype=bool)
tmp = tuple(np.empty_like(easting, dtype=bool) for i in range(4))
# Using the logical functions is a lot faster than & > < for some reason
# Plus, this way avoids repeated allocation of intermediate arrays
in_we = np.logical_and(
np.greater_equal(easting, w, out=tmp[0]),
np.less_equal(easting, e, out=tmp[1]),
out=tmp[2],
)
in_ns = np.logical_and(
np.greater_equal(northing, s, out=tmp[0]),
np.less_equal(northing, n, out=tmp[1]),
out=tmp[3],
)
are_inside = np.logical_and(in_we, in_ns, out=out)
return are_inside
[docs]
def block_split(coordinates, spacing=None, adjust="spacing", region=None, shape=None):
"""
Split a region into blocks and label points according to where they fall.
The labels are integers corresponding to the index of the block. Also
returns the coordinates of the center of each block (following the same
index as the labels).
The size of the blocks can be specified by the *spacing* parameter.
Alternatively, the number of blocks in the South-North and West-East
directions can be specified using the *shape* parameter.
.. note::
If installed, package ``pykdtree`` will be used instead of
:class:`scipy.spatial.cKDTree` for better performance.
Parameters
----------
coordinates : tuple of arrays
Arrays with the coordinates of each data point. Should be in the
following order: (easting, northing, vertical, ...). Only easting and
northing will be used, all subsequent coordinates will be ignored.
shape : tuple = (n_north, n_east) or None
The number of blocks in the South-North and West-East directions,
respectively.
spacing : float, tuple = (s_north, s_east), or None
The block size in the South-North and West-East directions,
respectively. A single value means that the size is equal in both
directions.
adjust : {'spacing', 'region'}
Whether to adjust the spacing or the region if required. Ignored if
*shape* is given instead of *spacing*. Defaults to adjusting the
spacing.
region : list = [W, E, S, N]
The boundaries of a given region in Cartesian or geographic
coordinates. If not region is given, will use the bounding region of
the given points.
Returns
-------
block_coordinates : tuple of arrays
(easting, northing) arrays with the coordinates of the center of each
block.
labels : array
integer label for each data point. The label is the index of the block
to which that point belongs.
See also
--------
BlockReduce : Apply a reduction operation to the data in blocks (windows).
rolling_window : Select points on a rolling (moving) window.
expanding_window : Select points on windows of changing size.
Examples
--------
>>> from verde import grid_coordinates
>>> coords = grid_coordinates((-5, 0, 5, 10), spacing=1)
>>> block_coords, labels = block_split(coords, spacing=2.5)
>>> for coord in block_coords:
... print(', '.join(['{:.2f}'.format(i) for i in coord]))
-3.75, -1.25, -3.75, -1.25
6.25, 6.25, 8.75, 8.75
>>> print(labels.reshape(coords[0].shape))
[[0 0 0 1 1 1]
[0 0 0 1 1 1]
[0 0 0 1 1 1]
[2 2 2 3 3 3]
[2 2 2 3 3 3]
[2 2 2 3 3 3]]
>>> # Use the shape instead of the block size
>>> block_coords, labels = block_split(coords, shape=(4, 2))
>>> for coord in block_coords:
... print(', '.join(['{:.3f}'.format(i) for i in coord]))
-3.750, -1.250, -3.750, -1.250, -3.750, -1.250, -3.750, -1.250
5.625, 5.625, 6.875, 6.875, 8.125, 8.125, 9.375, 9.375
>>> print(labels.reshape(coords[0].shape))
[[0 0 0 1 1 1]
[0 0 0 1 1 1]
[2 2 2 3 3 3]
[4 4 4 5 5 5]
[6 6 6 7 7 7]
[6 6 6 7 7 7]]
"""
# Select the coordinates after checking to make sure indexing will still
# work on the ignored coordinates.
coordinates = check_coordinates(coordinates)[:2]
if region is None:
region = get_region(coordinates)
block_coords = grid_coordinates(
region, spacing=spacing, shape=shape, adjust=adjust, pixel_register=True
)
tree = kdtree(block_coords)
labels = tree.query(np.transpose(n_1d_arrays(coordinates, 2)))[1]
return n_1d_arrays(block_coords, len(block_coords)), labels
[docs]
def rolling_window(
coordinates, size, spacing=None, shape=None, region=None, adjust="spacing"
):
"""
Select points on a rolling (moving) window.
A window of the given size is moved across the region at a given step
(specified by *spacing* or *shape*). Returns the indices of points falling
inside each window step. You can use the indices to select points falling
inside a given window.
The size of the step when moving the windows can be specified by the
*spacing* parameter. Alternatively, the number of windows in the
South-North and West-East directions can be specified using the *shape*
parameter. **One of the two must be given.**
Parameters
----------
coordinates : tuple of arrays
Arrays with the coordinates of each data point. Should be in the
following order: (easting, northing, vertical, ...). Only easting and
northing will be used, all subsequent coordinates will be ignored.
size : float
The size of the windows. Units should match the units of *coordinates*.
spacing : float, tuple = (s_north, s_east), or None
The window size in the South-North and West-East directions,
respectively. A single value means that the size is equal in both
directions.
shape : tuple = (n_north, n_east) or None
The number of blocks in the South-North and West-East directions,
respectively.
region : list = [W, E, S, N]
The boundaries of a given region in Cartesian or geographic
coordinates. If not region is given, will use the bounding region of
the given points.
adjust : {'spacing', 'region'}
Whether to adjust the spacing or the region if required. Ignored if
*shape* is given instead of *spacing*. Defaults to adjusting the
spacing.
Returns
-------
window_coordinates : tuple of arrays
Coordinate arrays for the center of each window.
indices : array
Each element of the array corresponds the indices of points falling
inside a window. The array will have the same shape as the
*window_coordinates*. Use the array elements to index the coordinates
for each window. The indices will depend on the number of dimensions in
the input coordinates. For example, if the coordinates are 2D arrays,
each window will contain indices for 2 dimensions (row, column).
See also
--------
block_split : Split a region into blocks and label points accordingly.
expanding_window : Select points on windows of changing size.
Examples
--------
Generate a set of sample coordinates on a grid and determine the indices
of points for each rolling window:
>>> from verde import grid_coordinates
>>> coords = grid_coordinates((-5, -1, 6, 10), spacing=1)
>>> print(coords[0])
[[-5. -4. -3. -2. -1.]
[-5. -4. -3. -2. -1.]
[-5. -4. -3. -2. -1.]
[-5. -4. -3. -2. -1.]
[-5. -4. -3. -2. -1.]]
>>> print(coords[1])
[[ 6. 6. 6. 6. 6.]
[ 7. 7. 7. 7. 7.]
[ 8. 8. 8. 8. 8.]
[ 9. 9. 9. 9. 9.]
[10. 10. 10. 10. 10.]]
>>> # Get the rolling window indices
>>> window_coords, indices = rolling_window(coords, size=2, spacing=2)
>>> # Window coordinates will be 2D arrays. Their shape is the number of
>>> # windows in each dimension
>>> print(window_coords[0].shape, window_coords[1].shape)
(2, 2) (2, 2)
>>> # The there are the easting and northing coordinates for the center of
>>> # each rolling window
>>> for coord in window_coords:
... print(coord)
[[-4. -2.]
[-4. -2.]]
[[7. 7.]
[9. 9.]]
>>> # The indices of points falling on each window will have the same shape
>>> # as the window center coordinates
>>> print(indices.shape)
(2, 2)
>>> # The points in the first window. Indices are 2D positions because the
>>> # coordinate arrays are 2D.
>>> print(len(indices[0, 0]))
2
>>> for dimension in indices[0, 0]:
... print(dimension)
[0 0 0 1 1 1 2 2 2]
[0 1 2 0 1 2 0 1 2]
>>> for dimension in indices[0, 1]:
... print(dimension)
[0 0 0 1 1 1 2 2 2]
[2 3 4 2 3 4 2 3 4]
>>> for dimension in indices[1, 0]:
... print(dimension)
[2 2 2 3 3 3 4 4 4]
[0 1 2 0 1 2 0 1 2]
>>> for dimension in indices[1, 1]:
... print(dimension)
[2 2 2 3 3 3 4 4 4]
[2 3 4 2 3 4 2 3 4]
>>> # To get the coordinates for each window, use indexing
>>> print(coords[0][indices[0, 0]])
[-5. -4. -3. -5. -4. -3. -5. -4. -3.]
>>> print(coords[1][indices[0, 0]])
[6. 6. 6. 7. 7. 7. 8. 8. 8.]
If the coordinates are 1D, the indices will also be 1D:
>>> coords1d = [coord.ravel() for coord in coords]
>>> window_coords, indices = rolling_window(coords1d, size=2, spacing=2)
>>> print(len(indices[0, 0]))
1
>>> print(indices[0, 0][0])
[ 0 1 2 5 6 7 10 11 12]
>>> print(indices[0, 1][0])
[ 2 3 4 7 8 9 12 13 14]
>>> print(indices[1, 0][0])
[10 11 12 15 16 17 20 21 22]
>>> print(indices[1, 1][0])
[12 13 14 17 18 19 22 23 24]
>>> # The returned indices can be used in the same way as before
>>> print(coords1d[0][indices[0, 0]])
[-5. -4. -3. -5. -4. -3. -5. -4. -3.]
>>> print(coords1d[1][indices[0, 0]])
[6. 6. 6. 7. 7. 7. 8. 8. 8.]
By default, the windows will span the entire data region. You can also
control the specific region you'd like the windows to cover:
>>> # Coordinates on a larger region but with the same spacing as before
>>> coords = grid_coordinates((-10, 5, 0, 20), spacing=1)
>>> # Get the rolling window indices but limited to the region from before
>>> window_coords, indices = rolling_window(
... coords, size=2, spacing=2, region=(-5, -1, 6, 10),
... )
>>> # The windows should still be in the same place as before
>>> for coord in window_coords:
... print(coord)
[[-4. -2.]
[-4. -2.]]
[[7. 7.]
[9. 9.]]
>>> # And indexing the coordinates should also provide the same result
>>> print(coords[0][indices[0, 0]])
[-5. -4. -3. -5. -4. -3. -5. -4. -3.]
>>> print(coords[1][indices[0, 0]])
[6. 6. 6. 7. 7. 7. 8. 8. 8.]
Only the first 2 coordinates are considered (assumed to be the horizontal
ones). All others will be ignored by the function.
>>> coords = grid_coordinates((-5, -1, 6, 10), spacing=1, extra_coords=20)
>>> print(coords[2])
[[20. 20. 20. 20. 20.]
[20. 20. 20. 20. 20.]
[20. 20. 20. 20. 20.]
[20. 20. 20. 20. 20.]
[20. 20. 20. 20. 20.]]
>>> window_coords, indices = rolling_window(coords, size=2, spacing=2)
>>> # The windows would be the same in this case since coords[2] is ignored
>>> for coord in window_coords:
... print(coord)
[[-4. -2.]
[-4. -2.]]
[[7. 7.]
[9. 9.]]
>>> print(indices.shape)
(2, 2)
>>> for dimension in indices[0, 0]:
... print(dimension)
[0 0 0 1 1 1 2 2 2]
[0 1 2 0 1 2 0 1 2]
>>> for dimension in indices[0, 1]:
... print(dimension)
[0 0 0 1 1 1 2 2 2]
[2 3 4 2 3 4 2 3 4]
>>> for dimension in indices[1, 0]:
... print(dimension)
[2 2 2 3 3 3 4 4 4]
[0 1 2 0 1 2 0 1 2]
>>> for dimension in indices[1, 1]:
... print(dimension)
[2 2 2 3 3 3 4 4 4]
[2 3 4 2 3 4 2 3 4]
>>> # The indices can still be used with the third coordinate
>>> print(coords[0][indices[0, 0]])
[-5. -4. -3. -5. -4. -3. -5. -4. -3.]
>>> print(coords[1][indices[0, 0]])
[6. 6. 6. 7. 7. 7. 8. 8. 8.]
>>> print(coords[2][indices[0, 0]])
[20. 20. 20. 20. 20. 20. 20. 20. 20.]
"""
# Check if shape or spacing were passed
if shape is None and spacing is None:
raise ValueError("Either a shape or a spacing must be provided.")
# Select the coordinates after checking to make sure indexing will still
# work on the ignored coordinates.
coordinates = check_coordinates(coordinates)[:2]
if region is None:
region = get_region(coordinates)
# Check if window size is bigger than the minimum dimension of the region
region_min_width = min(region[1] - region[0], region[3] - region[2])
if region_min_width < size:
raise ValueError(
"Window size '{}' is larger ".format(size)
+ "than dimensions of the region '{}'.".format(region)
)
# Calculate the region spanning the centers of the rolling windows
window_region = [
dimension + (-1) ** (i % 2) * size / 2 for i, dimension in enumerate(region)
]
_check_rolling_window_overlap(window_region, size, shape, spacing)
centers = grid_coordinates(
window_region, spacing=spacing, shape=shape, adjust=adjust
)
# pykdtree doesn't support query_ball_point yet and we need that
tree = kdtree(coordinates, use_pykdtree=False)
# Coordinates must be transposed because the kd-tree wants them as columns
# of a matrix
# Use p=inf (infinity norm) to get square windows instead of circular ones
indices1d = tree.query_ball_point(
np.transpose(n_1d_arrays(centers, 2)), r=size / 2, p=np.inf
)
# Make the indices array the same shape as the center coordinates array.
# That preserves the information of the number of windows in each
# dimension. Need to first create an empty array of object type because
# otherwise numpy tries to use the index tuples as dimensions (even if
# given ndim=1 explicitly). Can't make it 1D and then reshape because the
# reshape is ignored for some reason. The workaround is to create the array
# with the correct shape and assign the values to a raveled view of the
# array.
indices = np.empty(centers[0].shape, dtype="object")
# Need to convert the indices to int arrays because unravel_index doesn't
# like empty lists but can handle empty integer arrays in case a window has
# no points inside it.
indices.ravel()[:] = [
np.unravel_index(np.array(i, dtype="int"), shape=coordinates[0].shape)
for i in indices1d
]
return centers, indices
def _check_rolling_window_overlap(region, size, shape, spacing):
"""
Warn the user if there is no overlap between neighboring windows.
"""
if shape is not None:
ndims = len(shape)
dimensions = [region[i * ndims + 1] - region[i * ndims] for i in range(ndims)]
# The - 1 is because we need to divide by the number of intervals, not
# the number of nodes.
spacing = tuple(dim / (n - 1) for dim, n in zip(dimensions, shape))
spacing = np.atleast_1d(spacing)
if np.any(spacing > size):
warnings.warn(
f"Rolling windows do not overlap (size '{size}' and spacing '{spacing}'). "
"Some data points may not be included in any window. "
"Increase size or decrease spacing to avoid this.",
stacklevel=2,
)
[docs]
def expanding_window(coordinates, center, sizes):
"""
Select points on windows of changing size around a center point.
Returns the indices of points falling inside each window.
Parameters
----------
coordinates : tuple of arrays
Arrays with the coordinates of each data point. Should be in the
following order: (easting, northing, vertical, ...). Only easting and
northing will be used, all subsequent coordinates will be ignored.
center : tuple
The coordinates of the center of the window. Should be in the
following order: (easting, northing, vertical, ...).
sizes : array
The sizes of the windows. Does not have to be in any particular order.
The order of indices returned will match the order of window sizes
given. Units should match the units of *coordinates* and *center*.
Returns
-------
indices : list
Each element of the list corresponds to the indices of points falling
inside a window. Use them to index the coordinates for each window. The
indices will depend on the number of dimensions in the input
coordinates. For example, if the coordinates are 2D arrays, each window
will contain indices for 2 dimensions (row, column).
See also
--------
block_split : Split a region into blocks and label points accordingly.
rolling_window : Select points on a rolling (moving) window.
Examples
--------
Generate a set of sample coordinates on a grid and determine the indices
of points for each expanding window:
>>> from verde import grid_coordinates
>>> coords = grid_coordinates((-5, -1, 6, 10), spacing=1)
>>> print(coords[0])
[[-5. -4. -3. -2. -1.]
[-5. -4. -3. -2. -1.]
[-5. -4. -3. -2. -1.]
[-5. -4. -3. -2. -1.]
[-5. -4. -3. -2. -1.]]
>>> print(coords[1])
[[ 6. 6. 6. 6. 6.]
[ 7. 7. 7. 7. 7.]
[ 8. 8. 8. 8. 8.]
[ 9. 9. 9. 9. 9.]
[10. 10. 10. 10. 10.]]
>>> # Get the expanding window indices
>>> indices = expanding_window(coords, center=(-3, 8), sizes=[1, 2, 4])
>>> # There is one index per window
>>> print(len(indices))
3
>>> # The points in the first window. Indices are 2D positions because the
>>> # coordinate arrays are 2D.
>>> print(len(indices[0]))
2
>>> for dimension in indices[0]:
... print(dimension)
[2]
[2]
>>> for dimension in indices[1]:
... print(dimension)
[1 1 1 2 2 2 3 3 3]
[1 2 3 1 2 3 1 2 3]
>>> for dimension in indices[2]:
... print(dimension)
[0 0 0 0 0 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3 4 4 4 4 4]
[0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4]
>>> # To get the coordinates for each window, use indexing
>>> print(coords[0][indices[0]])
[-3.]
>>> print(coords[1][indices[0]])
[8.]
>>> print(coords[0][indices[1]])
[-4. -3. -2. -4. -3. -2. -4. -3. -2.]
>>> print(coords[1][indices[1]])
[7. 7. 7. 8. 8. 8. 9. 9. 9.]
If the coordinates are 1D, the indices will also be 1D:
>>> coords1d = [coord.ravel() for coord in coords]
>>> indices = expanding_window(coords1d, center=(-3, 8), sizes=[1, 2, 4])
>>> print(len(indices))
3
>>> # Since coordinates are 1D, there is only one index
>>> print(len(indices[0]))
1
>>> print(indices[0][0])
[12]
>>> print(indices[1][0])
[ 6 7 8 11 12 13 16 17 18]
>>> # The returned indices can be used in the same way as before
>>> print(coords1d[0][indices[0]])
[-3.]
>>> print(coords1d[1][indices[0]])
[8.]
Only the first 2 coordinates are considered (assumed to be the horizontal
ones). All others will be ignored by the function.
>>> coords = grid_coordinates((-5, -1, 6, 10), spacing=1, extra_coords=15)
>>> print(coords[2])
[[15. 15. 15. 15. 15.]
[15. 15. 15. 15. 15.]
[15. 15. 15. 15. 15.]
[15. 15. 15. 15. 15.]
[15. 15. 15. 15. 15.]]
>>> indices = expanding_window(coords, center=(-3, 8), sizes=[1, 2, 4])
>>> # The returned indices should be the same as before, ignoring coords[2]
>>> print(len(indices[0]))
2
>>> for dimension in indices[0]:
... print(dimension)
[2]
[2]
>>> for dimension in indices[1]:
... print(dimension)
[1 1 1 2 2 2 3 3 3]
[1 2 3 1 2 3 1 2 3]
>>> for dimension in indices[2]:
... print(dimension)
[0 0 0 0 0 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3 4 4 4 4 4]
[0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4]
>>> # The indices can be used to index all 3 coordinates
>>> print(coords[0][indices[0]])
[-3.]
>>> print(coords[1][indices[0]])
[8.]
>>> print(coords[2][indices[0]])
[15.]
"""
# Select the coordinates after checking to make sure indexing will still
# work on the ignored coordinates.
coordinates = check_coordinates(coordinates)[:2]
shape = coordinates[0].shape
center = np.atleast_2d(center)
# pykdtree doesn't support query_ball_point yet and we need that
tree = kdtree(coordinates, use_pykdtree=False)
indices = []
for size in sizes:
# Use p=inf (infinity norm) to get square windows instead of circular
index1d = tree.query_ball_point(center, r=size / 2, p=np.inf)[0]
# Convert indices to an array to avoid errors when the index is empty
# (no points in the window). unravel_index doesn't like empty lists.
indices.append(np.unravel_index(np.array(index1d, dtype="int"), shape=shape))
return indices
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 is in geographic coordinates are within accepted 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."
)