bordado.rolling_window_spherical

bordado.rolling_window_spherical(coordinates, window_size, overlap, *, region=None)

Split points into overlapping equal area windows on the sphere.

A window of the given latitudinal size is moved across the region at a given step (specified by the amount of overlap between adjacent windows). Windows are not regularly distributed on the sphere and are not “square” in longitude, latitude. They are evenly spaced in latitude but their longitudinal dimension varies to compensate for the convergence of meridians at the polar regions. The overlap will also wrap around the 360-0 longitude discontinuity (see examples below).

Returns the indices of points falling inside each window step. You can use the indices to select points falling inside a given window.

Parameters:

coordinatestuple = (longitude, latitude)

Tuple of arrays with the longitude and latitude coordinates of each point. Arrays can be Python lists or any numpy-compatible array type. Arrays can be of any shape but must all have the same shape.

window_sizefloat

The size of the windows along latitude in decimal degrees. The longitudinal window size is adjusted to retain equal area between windows. Must be > 0.

overlapfloat

The amount of overlap between adjacent windows. Should be within the range 1 > overlap ≥ 0. For example, an overlap of 0.5 means 50% overlap. The overlap may have to be adjusted to make sure windows fit inside the given region exactly.

regiontuple = (W, E, S, N)

The boundaries of a given region in geographic coordinates. Should have a lower and an upper boundary for each dimension of the coordinate system. If region is not given, will use the bounding region of the given coordinates.

Returns:

window_coordinatestuple = (longitude, latitude)

1D coordinate arrays for the center of each window.

indicesarray

1D array with each element of the array being a tuple of 2 arrays corresponding to the indices of the points that fall inside that particular window. Use these indices to index the given coordinates and select points from a window.

Notes

Uses the method of [Malkin2016] to divide the region into overlapping windows of equal area. The windows will have the specified window size in latitude but their longitudinal dimensions will be adjusted to account for the convergence of meridians at the poles.

Examples

Generate a set of sample coordinates on a grid to make it easier to visualize the windows:

>>> import bordado as bd
>>> import numpy as np

>>> coordinates = bd.grid_coordinates((0, 40, 60, 90), spacing=5)
>>> print(coordinates[0])
[[ 0.  5. 10. 15. 20. 25. 30. 35. 40.]
 [ 0.  5. 10. 15. 20. 25. 30. 35. 40.]
 [ 0.  5. 10. 15. 20. 25. 30. 35. 40.]
 [ 0.  5. 10. 15. 20. 25. 30. 35. 40.]
 [ 0.  5. 10. 15. 20. 25. 30. 35. 40.]
 [ 0.  5. 10. 15. 20. 25. 30. 35. 40.]
 [ 0.  5. 10. 15. 20. 25. 30. 35. 40.]]
>>> print(coordinates[1])
[[60. 60. 60. 60. 60. 60. 60. 60. 60.]
 [65. 65. 65. 65. 65. 65. 65. 65. 65.]
 [70. 70. 70. 70. 70. 70. 70. 70. 70.]
 [75. 75. 75. 75. 75. 75. 75. 75. 75.]
 [80. 80. 80. 80. 80. 80. 80. 80. 80.]
 [85. 85. 85. 85. 85. 85. 85. 85. 85.]
 [90. 90. 90. 90. 90. 90. 90. 90. 90.]]

Get the coordinates of the centers of rolling windows with 50% overlap and an indexer that allows us to select points from each window:

>>> window_size = 10  # degrees
>>> window_coords, indices = rolling_window_spherical(
...     coordinates, window_size=window_size, overlap=0.5,
... )

Window coordinates will be 1D arrays since the windows aren’t regular. Their longitudinal size is calculated to preserve their area. Their shape is the number of windows generated:

>>> print(window_coords[0].shape, window_coords[1].shape)
(8,) (8,)

The values of these arrays are the coordinates for the center of each rolling window. The latitude coordinates will be all at regular intervals dictated by the window size and the overlap:

>>> print(window_coords[1])
[65. 65. 70. 70. 75. 75. 80. 85.]

But in longitude, the window sizes (and thus their centers) will spread out as latitude increases to balance the convergence of meridians. The window size in longitude will be:

>>> window_size_lon = window_size / np.cos(np.radians(window_coords[1]))
>>> print(np.array_str(window_size_lon, precision=1))
[ 23.7  23.7  29.2  29.2  38.6  38.6  57.6 114.7]

Notice that as we get closer to the pole the window size is larger than the region, so in these cases the windows cannot be guaranteed to have equal area. The center of each window will be:

>>> print(np.array_str(window_coords[0], precision=1))
[11.8 28.2 14.6 25.4 19.3 20.7 20.  20. ]

If you look closely, you’ll see that the amount of overlap between the windows isn’t exactly 50%. This is because the overlap and window size do not result in multiples of the region. So we have to adjust the overlap to make things fit. This effect is less evident for smaller windows.

The indices of points falling on each window will have the same shape as the window center coordinates:

>>> print(indices.shape)
(8,)

Each element of the indices array is a tuple of arrays, one for each element in the coordinates:

>>> print(len(indices[0]))
2

They are indices of the points that fall inside the selected window (window 0). The first element indexes the axis 0 of the coordinate arrays and so forth. So this:

>>> print(indices[0][0])
[0 0 0 0 0 1 1 1 1 1 2 2 2 2 2]

corresponds to the rows of the coordinate arrays that belong to the first window, and this

>>> print(indices[0][1])
[0 1 2 3 4 0 1 2 3 4 0 1 2 3 4]

corresponds to the columns of the coordinate arrays that belong to the first window.

The same can be showed for the second window:

>>> print(indices[1][0])
[0 0 0 0 0 1 1 1 1 1 2 2 2 2 2]
>>> print(indices[1][1])
[4 5 6 7 8 4 5 6 7 8 4 5 6 7 8]

and the third window:

>>> print(indices[-1][0])
[4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 6 6 6 6 6 6 6 6 6]
>>> print(indices[-1][1])
[0 1 2 3 4 5 6 7 8 0 1 2 3 4 5 6 7 8 0 1 2 3 4 5 6 7 8]

Use these indices to select the coordinates the points that fall inside a window:

>>> points_window_0 = [c[indices[0]] for c in coordinates]
>>> print(points_window_0[0])
[ 0.  5. 10. 15. 20.  0.  5. 10. 15. 20.  0.  5. 10. 15. 20.]
>>> print(points_window_0[1])
[60. 60. 60. 60. 60. 65. 65. 65. 65. 65. 70. 70. 70. 70. 70.]

>>> points_window_1 = [c[indices[1]] for c in coordinates]
>>> print(points_window_1[0])
[20. 25. 30. 35. 40. 20. 25. 30. 35. 40. 20. 25. 30. 35. 40.]
>>> print(points_window_1[1])
[60. 60. 60. 60. 60. 65. 65. 65. 65. 65. 70. 70. 70. 70. 70.]

If the coordinates are 1D, the indices will also be 1D:

>>> coordinates1d = [c.ravel() for c in coordinates]
>>> window_coords, indices = rolling_window_spherical(
...     coordinates1d, window_size=window_size, overlap=0.5,
... )
>>> print(len(indices[0]))
1

In this case, the indices will refer to the raveled coordinate array. The indexer for the first window will be:

>>> print(indices[0][0])
[ 0  1  2  3  4  9 10 11 12 13 18 19 20 21 22]

And for the second window:

>>> print(indices[1][0])
[ 4  5  6  7  8 13 14 15 16 17 22 23 24 25 26]

The returned indices can be used in the same way as before to get the same coordinates:

>>> print(coordinates1d[0][indices[0]])
[ 0.  5. 10. 15. 20.  0.  5. 10. 15. 20.  0.  5. 10. 15. 20.]
>>> print(coordinates1d[1][indices[0]])
[60. 60. 60. 60. 60. 65. 65. 65. 65. 65. 70. 70. 70. 70. 70.]

By default, the windows will span the entire data region. You can also control the specific region you’d like the windows to cover:

>>> window_coords, indices = rolling_window_spherical(
...     coordinates,
...     window_size=window_size,
...     overlap=0.5,
...     region=(0, 20, 60, 75),
... )

The windows will now try to fit the smaller region instead of the full extent of the coordinates:

>>> print(coordinates[0][indices[0]])
[ 0.  5. 10. 15. 20.  0.  5. 10. 15. 20.  0.  5. 10. 15. 20.]
>>> print(coordinates[1][indices[0]])
[60. 60. 60. 60. 60. 65. 65. 65. 65. 65. 70. 70. 70. 70. 70.]

>>> print(coordinates[0][indices[1]])
[ 0.  5. 10. 15. 20.  0.  5. 10. 15. 20.  0.  5. 10. 15. 20.]
>>> print(coordinates[1][indices[1]])
[65. 65. 65. 65. 65. 70. 70. 70. 70. 70. 75. 75. 75. 75. 75.]

If the longitude range is a full 360 degrees, the windows will wrap around the 360-0 discontinuity:

>>> coordinates = bd.grid_coordinates((0, 360, 0, 0), spacing=30)
>>> print(coordinates[0])
[[  0.  30.  60.  90. 120. 150. 180. 210. 240. 270. 300. 330. 360.]]
>>> print(coordinates[1])
[[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]]

>>> window_coords, indices = rolling_window_spherical(
...     coordinates, window_size=90, overlap=0.5,
... )
>>> print(indices.shape)
(8,)
>>> print(coordinates[0][indices[0]])
[  0.  30.  60.  90. 360.]

The 360 point is there because it’s the same as the 0 point and so they are both inside the window.

The last window will be centered at 0 and will wrap around the 360-0 divide:

>>> print(coordinates[0][indices[-1]])
[  0.  30. 330. 360.]

This way, the windows wrap around the globe and overlap all the way around.