Skip to content

viewshed max_distance clips cells on anisotropic rasters #2691

@brendancol

Description

@brendancol

Description

viewshed(..., max_distance=...) drops cells that are within max_distance of the observer when the raster has anisotropic resolution (different cell spacing in x and y). Those cells come back INVISIBLE even though they should be evaluated.

In _viewshed_windowed() (xrspatial/viewshed.py) the analysis window is sized from a single cell radius:

cell_size = max(abs(ew_res), abs(ns_res))
radius_cells = int(np.ceil(max_distance / cell_size))
r_lo = max(0, obs_r - radius_cells)
r_hi = min(height, obs_r + radius_cells + 1)
c_lo = max(0, obs_c - radius_cells)
c_hi = min(width, obs_c + radius_cells + 1)

radius_cells is derived from the coarser of the two resolutions and then used to clip both rows and columns. When ew_res and ns_res differ, the window ends up too small along the finer axis, and any cell beyond that truncated window never gets computed.

There is a circular distance mask applied afterward, but it can only remove cells, not bring back ones already cut from the window.

Reproduce

Raster with ew_res=1 (x) and ns_res=10 (y). max_distance=8 reaches 8 cells east of the observer, but radius_cells = ceil(8/10) = 1, so only one column on each side survives.

import numpy as np, xarray as xa
from xrspatial import viewshed

ny, nx = 21, 21
t = np.zeros((ny, nx))
xs = np.arange(nx, dtype=float) * 1.0    # ew_res = 1
ys = np.arange(ny, dtype=float) * 10.0   # ns_res = 10
r = xa.DataArray(t, coords=dict(x=xs, y=ys), dims=['y', 'x'])

full = viewshed(r, x=xs[10], y=ys[10], observer_elev=50).values
win = viewshed(r.copy(), x=xs[10], y=ys[10], observer_elev=50, max_distance=8.0).values

# cell 5 columns east is 5 units away (< 8), should be visible
print(full[10, 15])  # 5.71...
print(win[10, 15])   # -1.0  (wrongly clipped)

Expected

Cells within max_distance should be evaluated no matter how anisotropic the resolution is. The windowed result should match the full viewshed inside the radius, the same way test_viewshed_max_distance_matches_full already checks for isotropic rasters.

Fix

Use a separate cell radius per axis: radius_rows from ns_res, radius_cols from ew_res.

Impact

_viewshed_windowed() handles max_distance for every backend (numpy, cupy, dask+numpy, dask+cupy), so all of them are affected.

Metadata

Metadata

Assignees

No one assigned

    Labels

    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