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) → 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. 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.

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. For very large grids that exceed memory, a horizon-profile distance-sweep algorithm is used instead. This algorithm discretises angles and may produce minor visibility differences near the boundary of occluded regions.

Both backends agree on which cells are visible vs invisible in the vast majority of cases, but the reported vertical angles may differ by a small amount near cell boundaries.

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