xrspatial.viewshed.viewshed#

xrspatial.viewshed.viewshed(raster: DataArray, x: int | float, y: int | float, observer_elev: float = 0, target_elev: float = 0, max_distance: float = None, name: str | None = 'viewshed') → DataArray[source]#

Calculate viewshed of a raster (the visible cells in the raster) for the given viewpoint (observer) location.

Parameters:

  • raster (xr.DataArray) – Input raster image.

  • x (int, float) – x-coordinate in data space of observer location.

  • y (int, float) – y-coordinate in data space of observer location.

  • observer_elev (float) – Observer elevation above the terrain.

  • target_elev (float) – Target elevation offset above the terrain, which is the height in surface units to be added to the z-value of each pixel when it is being analyzed for visibility.

  • max_distance (float, optional) – Maximum analysis distance from the observer in surface units. Must be a finite number >= 0; a negative or non-finite value raises ValueError. Cells beyond this distance are marked INVISIBLE without being evaluated. When set and the raster is dask-backed, only the chunks within the distance window are loaded — this is the most efficient way to run viewshed on very large dask rasters.

  • name (str, default='viewshed') – Name of the output DataArray. Set on every backend so the result name does not depend on which backend ran.

Returns:

viewshed – A cell x in the visibility grid is recorded as follows: If it is invisible, then x is set to INVISIBLE. If it is visible, then x is set to the vertical angle w.r.t the viewpoint. The value returned is in [0, 180]. A value of 0 is directly below the specified viewing position, 90 is due horizontal, and 180 is directly above the observer.

Return type:

xarray.DataArray

Notes

The CPU (numpy), GPU (cupy with RTX), and dask backends use different algorithms and may produce slightly different results for the same input.

  • CPU: Angular sweep algorithm ported from GRASS GIS r.viewshed. Operates directly on the grid and produces exact visibility angles.

  • GPU: Ray tracing via NVIDIA OptiX RTX against a triangulated mesh of the terrain. The mesh discretisation can introduce small angular errors (typically < 0.03 degrees for visible cells).

  • Dask: When max_distance is set or the grid fits in memory, the exact CPU algorithm is used on the relevant window, so results match the numpy backend. For very large grids that exceed memory and have no max_distance, an out-of-core horizon-profile distance-sweep algorithm is used instead. This is a different, approximate visibility model from the exact GRASS sweep, and the two do not agree cell-for-cell. On rough terrain the visibility mask can differ for a substantial fraction of cells (measured at up to ~20% on small random rasters), with both false positives and false negatives relative to the exact sweep. The error is geometric and does not shrink with finer angle discretisation. When this path runs, viewshed() emits a UserWarning so the approximation is not silent. If you need results that match the exact sweep on a large dask raster, set max_distance to restrict the analysis to a window that fits in memory.

The CPU and GPU backends, and the dask exact-window path, agree on which cells are visible vs invisible in the vast majority of cases; reported vertical angles may differ by a small amount near cell boundaries. The dask out-of-core distance sweep is the exception described above.

Examples

>>> import numpy as np
>>> import xarray as xr
>>> from xrspatial import viewshed
>>> data = np.array([
...     [0, 0, 1, 0, 0],
...     [1, 3, 0, 0, 0],
...     [10, 2, 5, 2, -1],
...     [11, 1, 2, 9, 0]])
>>> terrain = xr.DataArray(data, dims=['y', 'x'])
>>> h, w = data.shape
>>> terrain['y'] = np.linspace(1, h, h)
>>> terrain['x'] = np.linspace(1, w, w)
>>> terrain
<xarray.DataArray (y: 4, x: 5)>
array([[ 0,  0,  1,  0,  0],
       [ 1,  3,  0,  0,  0],
       [10,  2,  5,  2, -1],
       [11,  1,  2,  9,  0]])
Coordinates:
  * y        (y) float64 1.0 2.0 3.0 4.0
  * x        (x) float64 1.0 2.0 3.0 4.0 5.0
>>> viewshed(terrain, x=3, y=2)
<xarray.DataArray (y: 4, x: 5)>
array([[ -1.        ,  90.        , 135.        ,  90.        , -1.],
       [ -1.        , 161.56505118, 180.        ,  90.        , 90.],
       [167.39561735, 144.73561032, 168.69006753, 144.73561032, -1.],
       [165.57993189,  -1.        ,  -1.        , 166.0472636 , -1.]])
Coordinates:
  * x        (x) float64 1.0 2.0 3.0 4.0 5.0
  * y        (y) float64 1.0 2.0 3.0 4.0