Skip to content

GPU viewshed ignores x/y resolution: anisotropic rasters diverge from CPU #2861

@brendancol

Description

@brendancol

Describe the bug

The GPU ray-tracing viewshed builds its terrain mesh on integer grid coordinates and ignores the real cell resolution. In xrspatial/gpu_rtx/mesh_utils.py, the mesh vertices sit at (w, h) integer indices (_triangulate_terrain_kernel, around line 77/84). The real ew_res/ns_res only come into play later, in xrspatial/gpu_rtx/viewshed.py around line 243, when the output vertical angle is computed.

So on an anisotropic raster (ew_res != ns_res), the GPU traces a different terrain than the CPU sweep does. The CPU path (xrspatial/viewshed.py) scales horizontal distances by the real ew_res/ns_res wherever it computes gradients and distances, so occlusion gets decided in real-world units. The GPU decides occlusion in integer grid units. A ridge that blocks line of sight in real units may not block it in grid units, or the other way around, and the two backends end up reporting different sets of visible cells.

Expected behavior

GPU and CPU viewshed should agree on an anisotropic raster, within the small float32 ray-tracing tolerance. The GPU mesh geometry should reflect the real cell resolution so the ray tracer sees the same terrain the CPU sweep does.

Reproduction

import numpy as np
import xarray as xa
import cupy as cp
from xrspatial import viewshed

ny, nx = 15, 15
terrain = np.full((ny, nx), 1.0)
terrain[:, 7] = 8.0  # ridge that occludes

xs = np.arange(nx, dtype=float) * 1.0    # ew_res = 1
ys = np.arange(ny, dtype=float) * 10.0   # ns_res = 10  (anisotropic)

obs_x, obs_y, obs_elev = xs[3], ys[7], 2

v_np = viewshed(xa.DataArray(terrain, coords=dict(x=xs, y=ys), dims=["y","x"]),
                x=obs_x, y=obs_y, observer_elev=obs_elev)
v_gpu = viewshed(xa.DataArray(cp.asarray(terrain), coords=dict(x=xs, y=ys), dims=["y","x"]),
                 x=obs_x, y=obs_y, observer_elev=obs_elev)

n = v_np.values
g = v_gpu.data.get()
print("visibility disagreements:", ((n > -1) != (g > -1)).sum(), "of", n.size)
# -> 105 of 225

Additional context

The fix is to build the GPU mesh (and the matching camera/viewpoint ray origins) in real-resolution coordinates instead of integer grid indices, and keep the downstream angle calculation consistent so resolution is not applied twice. create_triangulation is shared with the RTX hillshade path, so that path has to stay consistent too. GPU-only code path.

Metadata

Metadata

Assignees

No one assigned

    Labels

    bugSomething isn't workinggpuCuPy / 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