Skip to content

kde returns all-zero or wrong output on dask backends with descending-coordinate templates #3627

Description

@brendancol

Describe the bug

kde() with a dask-backed template whose y or x coordinates descend (row 0 = north, the usual layout for geospatial rasters and the default in our own create_test_raster helper) returns all zeros for the compact kernels (epanechnikov, quartic). The gaussian kernel returns nonzero but wrong values, which is worse: it looks plausible. dask+cupy fails the same way. Eager numpy and cupy are fine.

This is #1198 again, one layer up. #1199 fixed the sign handling in the CPU kernels, but the dask paths pre-filter points per tile in _filter_points_to_tile (xrspatial/kde.py:384) and that filter still assumes positive pixel spacing:

tile_x1 = tile_x0 + tile_cols * dx
tile_y1 = tile_y0 + tile_rows * dy
mask = ((xs >= tile_x0 - cutoff) & (xs <= tile_x1 + cutoff) &
        (ys >= tile_y0 - cutoff) & (ys <= tile_y1 + cutoff))

With negative dx/dy, tile_x1 < tile_x0, so the interval is inverted and most or all points are dropped before the kernel runs. The gaussian kernel's 4*bw cutoff is wide enough for a few points to survive, so its output is partially populated rather than empty.

Expected behavior

Dask results match numpy for any coordinate orientation, same as the eager path already guarantees via the #1199 regression tests.

Reproduction (this host, worktree at b9ba8a8)

import numpy as np
import xarray as xr
import dask.array as dask_array
from xrspatial.kde import kde

rng = np.random.default_rng(42)
x = rng.uniform(-3, 3, 50)
y = rng.uniform(-3, 3, 50)

def template(desc_y, desc_x, dask=False):
    ys = np.linspace(4, -4, 16) if desc_y else np.linspace(-4, 4, 16)
    xs = np.linspace(4, -4, 16) if desc_x else np.linspace(-4, 4, 16)
    t = xr.DataArray(np.zeros((16, 16)), dims=['y', 'x'], coords={'y': ys, 'x': xs})
    if dask:
        t = t.copy(data=dask_array.from_array(t.values, chunks=(8, 8)))
    return t

for kernel in ['quartic', 'gaussian']:
    for desc_y, desc_x, label in [(False, False, 'ascending'), (True, False, 'desc-y'),
                                  (False, True, 'desc-x'), (True, True, 'desc-both')]:
        np_res = kde(x, y, bandwidth=1.0, kernel=kernel, template=template(desc_y, desc_x)).values
        dk_res = kde(x, y, bandwidth=1.0, kernel=kernel, template=template(desc_y, desc_x, dask=True)).values
        print(kernel, label, 'numpy.sum=%.4f' % np_res.sum(), 'dask.sum=%.4f' % dk_res.sum(),
              'max|diff|=%.3e' % np.abs(np_res - dk_res).max())

Output:

quartic ascending numpy.sum=175.9665 dask.sum=175.9665 max|diff|=0.000e+00
quartic desc-y    numpy.sum=175.9665 dask.sum=0.0000   max|diff|=3.475e+00
quartic desc-x    numpy.sum=175.9665 dask.sum=0.0000   max|diff|=3.475e+00
quartic desc-both numpy.sum=175.9665 dask.sum=0.0000   max|diff|=3.475e+00
gaussian ascending numpy.sum=171.5630 dask.sum=171.5630 max|diff|=0.000e+00
gaussian desc-y    numpy.sum=171.5640 dask.sum=134.8469 max|diff|=1.409e+00
gaussian desc-x    numpy.sum=171.5631 dask.sum=135.9322 max|diff|=9.542e-01
gaussian desc-both numpy.sum=171.5640 dask.sum=107.2278 max|diff|=1.588e+00

dask+cupy gives the same zeros (verified on this host with CUDA).

Additional context

The existing TestDaskParity and TestCuPyParity tests only use an ascending template, which is how this slipped past the #1199 regression tests. Found by the 2026-07-03 accuracy sweep of the kde module.

Metadata

Metadata

Assignees

No one assigned

    Labels

    backend-coverageAdding missing dask/cupy/dask+cupy backend supportbugSomething isn't workingdaskDask backend / chunked arraysseverity:highSweep finding: HIGHsweep-accuracyFound by /sweep-accuracy

    Type

    No type

    Fields

    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