Skip to content

Planar aspect() cupy backend snaps near-360 aspects to 0, diverging from numpy #2827

@brendancol

Description

@brendancol

Description

The planar aspect() GPU device function clamps aspect values just below 360 down to 0, but the numpy/numba CPU kernel does not. On the same input the cupy and dask+cupy backends return a value about 360 degrees away from the numpy and dask+numpy backends.

In xrspatial/aspect.py, the cupy device function _gpu (lines 126-129) ends with:

if _aspect > 359.999:  # lame float equality check...
    return 0
else:
    return _aspect

The numpy/numba _cpu kernel has no equivalent clamp. The CPU compass conversion produces output in the half-open range [0, 360) and never reaches 360, so there is nothing on the CPU side for this clamp to canonicalize. The clamp lives only on the GPU side, so instead of bringing the backends into agreement it pulls them apart.

Reproduce

import numpy as np, xarray as xr, cupy
from xrspatial import aspect

# Near-degenerate gradient plane: tiny east gradient, strong north gradient.
n = 5
xs = np.arange(n); ys = np.arange(n)
z = (1e-6 * xs[None, :] + 1.0 * ys[:, None]).astype(np.float64)

def mk(backend):
    a = xr.DataArray(z.copy(), dims=['y', 'x'],
                     coords={'y': np.arange(n), 'x': np.arange(n)})
    if backend == 'cupy':
        a.data = cupy.asarray(a.data)
    return a

np_r = aspect(mk('numpy')).data
cu_r = aspect(mk('cupy')).data.get()
print(np_r[2, 2])  # 359.99994
print(cu_r[2, 2])  # 0.0

numpy returns 359.99994; cupy returns 0.0. Both point north, so the physical direction is the same, but the stored numeric values differ by about 360. That breaks any downstream code that compares backends or does arithmetic on the raw degrees.

Why the existing tests miss it

The cross-backend random-data tests never land in the (359.999, 360) band. Across 400 random rasters at the test sizes, zero hit it. The band is narrower than a thousandth of a degree and only reachable with a near-degenerate gradient (east gradient near zero, north gradient positive), so random elevation data almost never triggers it.

Scope

This is planar-only. The geodesic path canonicalizes the same way on CPU and GPU (both do if aspect_deg >= 360.0: aspect_deg -= 360.0), so it is fine. The planar cellsize handling fixed in #2780 is unrelated and should be left alone.

Proposed fix

Remove the asymmetric clamp from _gpu so the cupy output matches numpy and stays in [0, 360). Add a cross-backend test that builds the near-360 gradient plane and checks that numpy, cupy, dask+numpy, and dask+cupy agree.

Metadata

Metadata

Assignees

No one assigned

    Labels

    bugSomething isn't workingdaskDask backend / chunked arraysgpuCuPy / CUDA GPU support

    Type

    No type
    No fields configured for issues without a type.

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions