Skip to content

Bounded dask GREAT_CIRCLE proximity misses targets across the antimeridian #3108

@brendancol

Description

@brendancol

Describe the bug

proximity/allocation/direction with distance_metric='GREAT_CIRCLE' and a finite max_distance disagree between dask and numpy when the raster spans the antimeridian.

The bounded dask path sizes its map_overlap halo in array space: _halo_depth divides max_distance by the densest per-column step distance. That works for planar metrics, but great-circle distance is periodic in longitude. Haversine takes the short way around, so a pixel at lon 179.5 is about 111 km from a target at lon -179.5. In array space that target is at the opposite end of the axis, far outside any per-chunk halo, and the chunk never sees it.

Repro:

import numpy as np, xarray as xr, dask.array as da
from xrspatial import proximity

lon = np.arange(-179.5, 180.0, 1.0)
lat = np.arange(-5.0, 6.0, 1.0)
data = np.zeros((lat.size, lon.size))
data[5, 0] = 1.0  # target at (lat 0, lon -179.5)

r = xr.DataArray(data, dims=['lat', 'lon'], coords={'lat': lat, 'lon': lon})
kw = dict(x='lon', y='lat', distance_metric='GREAT_CIRCLE', max_distance=250_000.0)

proximity(r, **kw).values[5, -1]                                                   # 111319.49
proximity(r.copy(data=da.from_array(data, chunks=(11, 60))), **kw).values[5, -1]   # nan

numpy and cupy run a whole-raster brute-force search and find the target. dask+numpy and dask+cupy both return NaN.

The linear column halo has a second problem at high latitude: a great-circle chord between two points on a parallel is shorter than the sum of the per-column steps, so the halo can be undersized even without a wrap. The worst-latitude sizing usually rescues this by blowing up the pad until the axis folds, but that is not guaranteed.

Expected behavior

All four backends return the same result. The chord bound dist >= 2R asin(cos(lat_max) |sin(dlon/2)|) holds for any pair of grid points, so the column halo can be derived from it. When the seam gap (or the 180-degree chord at the worst-case latitude) is still within max_distance, no array-space halo works and the x axis has to be folded into a single chunk, which _fit_halo_to_chunks already does for oversized pads.

Additional context

_halo_depth is shared by the dask+numpy and dask+cupy bounded paths, so one fix covers both.

Metadata

Metadata

Assignees

No one assigned

    Labels

    bugSomething isn't working

    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