Add sky-view factor function (#962)

README.md

@@ -266,6 +266,7 @@ In the GIS world, rasters are used for representing continuous phenomena (e.g. e
 | [Curvature](xrspatial/curvature.py) | Measures rate of slope change (concavity/convexity) at each cell | ✅️ |✅️ |✅️ | ✅️  |
 | [Hillshade](xrspatial/hillshade.py) | Simulates terrain illumination from a given sun angle and azimuth | ✅️ | ✅️ | ✅️ | ✅️ |
 | [Roughness](xrspatial/terrain_metrics.py) | Computes local relief as max minus min elevation in a 3×3 window | ✅️ | ✅️ | ✅️ | ✅️ |
+| [Sky-View Factor](xrspatial/sky_view_factor.py) | Measures the fraction of visible sky hemisphere at each cell | ✅️ | ✅️ | ✅️ | ✅️ |
 | [Slope](xrspatial/slope.py) | Computes terrain gradient steepness at each cell in degrees | ✅️  | ✅️  | ✅️ | ✅️  |
 | [Terrain Generation](xrspatial/terrain.py) | Generates synthetic terrain from fBm or ridged fractal noise with optional domain warping, Worley blending, and hydraulic erosion | ✅️ | ✅️ | ✅️ | ✅️ |
 | [TPI](xrspatial/terrain_metrics.py) | Computes Topographic Position Index (center minus mean of neighbors) | ✅️ | ✅️ | ✅️ | ✅️ |

docs/source/reference/surface.rst

@@ -39,6 +39,13 @@ Terrain Generation
 
     xrspatial.terrain.generate_terrain
 
+Sky-View Factor
+===============
+.. autosummary::
+    :toctree: _autosummary
+
+    xrspatial.sky_view_factor.sky_view_factor
+
 Viewshed
 ========
 .. autosummary::

examples/user_guide/22_Sky_View_Factor.ipynb

@@ -0,0 +1,101 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "0fnki2gcn3md",
+   "source": "# Sky-View Factor\n\nSky-view factor (SVF) measures the fraction of the sky hemisphere visible from each cell in a digital elevation model. Values range from 0 (fully obstructed, e.g. a deep canyon) to 1 (flat open terrain with no horizon obstruction).\n\nSVF is used in:\n- **LiDAR archaeology** to reveal subtle terrain features hidden under canopy\n- **Urban heat island studies** to quantify street canyon geometry\n- **Solar energy modeling** as a proxy for diffuse sky irradiance\n- **Terrain visualization** as an illumination-independent alternative to hillshade\n\nThe algorithm casts rays at evenly spaced azimuths from each cell and records the maximum elevation angle to the horizon along each ray. SVF is then:\n\n$$\\text{SVF} = 1 - \\frac{1}{N}\\sum_{d=1}^{N} \\sin(\\theta_d)$$\n\nwhere $\\theta_d$ is the maximum horizon angle in direction $d$.",
+   "metadata": {}
+  },
+  {
+   "cell_type": "code",
+   "id": "050wobw30odo",
+   "source": "import numpy as np\nimport xarray as xr\nimport matplotlib.pyplot as plt\n\nfrom xrspatial import sky_view_factor",
+   "metadata": {},
+   "execution_count": null,
+   "outputs": []
+  },
+  {
+   "cell_type": "markdown",
+   "id": "9w2idcivtrw",
+   "source": "## Generate synthetic terrain\n\nWe'll create a terrain surface with valleys, ridges, and a flat plateau to show how SVF responds to different landforms.",
+   "metadata": {}
+  },
+  {
+   "cell_type": "code",
+   "id": "yubj4xlrxd",
+   "source": "# Create a 200x200 terrain with ridges and valleys\nrows, cols = 200, 200\ny = np.arange(rows, dtype=np.float64)\nx = np.arange(cols, dtype=np.float64)\nY, X = np.meshgrid(y, x, indexing='ij')\n\n# Sinusoidal ridges + a central peak\nterrain = (\n    40 * np.sin(X / 15) * np.cos(Y / 20)\n    + 80 * np.exp(-((Y - 100)**2 + (X - 100)**2) / (2 * 30**2))\n    + 20 * np.sin(Y / 10)\n    + 200\n)\n\ndem = xr.DataArray(\n    terrain,\n    dims=['y', 'x'],\n    coords={'y': np.arange(rows, dtype=float), 'x': np.arange(cols, dtype=float)},\n    attrs={'res': (1.0, 1.0)},\n)\n\nfig, ax = plt.subplots(figsize=(8, 7))\ndem.plot(ax=ax, cmap='terrain', cbar_kwargs={'label': 'Elevation'})\nax.set_title('Synthetic DEM')\nax.set_aspect('equal')\nplt.tight_layout()\nplt.show()",
+   "metadata": {},
+   "execution_count": null,
+   "outputs": []
+  },
+  {
+   "cell_type": "markdown",
+   "id": "35xi1fqq831",
+   "source": "## Compute SVF with default parameters\n\nThe default uses `max_radius=10` cells and `n_directions=16` azimuth rays.",
+   "metadata": {}
+  },
+  {
+   "cell_type": "code",
+   "id": "mrt8dfxrni",
+   "source": "svf = sky_view_factor(dem, max_radius=10, n_directions=16)\n\nfig, ax = plt.subplots(figsize=(8, 7))\nsvf.plot(ax=ax, cmap='gray', vmin=0, vmax=1,\n         cbar_kwargs={'label': 'Sky-View Factor'})\nax.set_title('SVF (max_radius=10, n_directions=16)')\nax.set_aspect('equal')\nplt.tight_layout()\nplt.show()",
+   "metadata": {},
+   "execution_count": null,
+   "outputs": []
+  },
+  {
+   "cell_type": "markdown",
+   "id": "ot8twzbdqgq",
+   "source": "## Effect of `max_radius`\n\nA larger search radius captures more distant obstructions. This matters in terrain with broad features like wide valleys or distant ridgelines.",
+   "metadata": {}
+  },
+  {
+   "cell_type": "code",
+   "id": "pq6q4dvd2dd",
+   "source": "radii = [5, 15, 30]\nfig, axes = plt.subplots(1, 3, figsize=(18, 5))\n\nfor ax, r in zip(axes, radii):\n    result = sky_view_factor(dem, max_radius=r, n_directions=16)\n    result.plot(ax=ax, cmap='gray', vmin=0, vmax=1, add_colorbar=False)\n    ax.set_title(f'max_radius={r}')\n    ax.set_aspect('equal')\n\nfig.suptitle('SVF at different search radii', fontsize=14, y=1.02)\nplt.tight_layout()\nplt.show()",
+   "metadata": {},
+   "execution_count": null,
+   "outputs": []
+  },
+  {
+   "cell_type": "markdown",
+   "id": "l1rcmgtk8",
+   "source": "## Effect of `n_directions`\n\nMore azimuth directions give a smoother result at the cost of longer computation. For most uses, 16 directions is a good balance.",
+   "metadata": {}
+  },
+  {
+   "cell_type": "code",
+   "id": "b33ppidyce",
+   "source": "directions = [4, 16, 64]\nfig, axes = plt.subplots(1, 3, figsize=(18, 5))\n\nfor ax, n in zip(axes, directions):\n    result = sky_view_factor(dem, max_radius=15, n_directions=n)\n    result.plot(ax=ax, cmap='gray', vmin=0, vmax=1, add_colorbar=False)\n    ax.set_title(f'n_directions={n}')\n    ax.set_aspect('equal')\n\nfig.suptitle('SVF at different direction counts', fontsize=14, y=1.02)\nplt.tight_layout()\nplt.show()",
+   "metadata": {},
+   "execution_count": null,
+   "outputs": []
+  },
+  {
+   "cell_type": "markdown",
+   "id": "n7hl0m8uhmq",
+   "source": "## Comparison: SVF vs Hillshade\n\nHillshade depends on a specific sun angle, which creates directional bias. SVF provides omnidirectional illumination information, making features visible regardless of orientation.",
+   "metadata": {}
+  },
+  {
+   "cell_type": "code",
+   "id": "v8qltgseoba",
+   "source": "from xrspatial import hillshade\n\nhs = hillshade(dem, azimuth=315, angle_altitude=45)\nsvf_result = sky_view_factor(dem, max_radius=15, n_directions=16)\n\nfig, axes = plt.subplots(1, 2, figsize=(14, 6))\n\nhs.plot(ax=axes[0], cmap='gray', add_colorbar=False)\naxes[0].set_title('Hillshade (azimuth=315)')\naxes[0].set_aspect('equal')\n\nsvf_result.plot(ax=axes[1], cmap='gray', vmin=0, vmax=1, add_colorbar=False)\naxes[1].set_title('Sky-View Factor')\naxes[1].set_aspect('equal')\n\nplt.tight_layout()\nplt.show()",
+   "metadata": {},
+   "execution_count": null,
+   "outputs": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "name": "python",
+   "version": "3.10.0"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

xrspatial/__init__.py

@@ -72,6 +72,7 @@
 from xrspatial.snap_pour_point import snap_pour_point  # noqa
 from xrspatial.stream_link import stream_link  # noqa
 from xrspatial.stream_order import stream_order  # noqa
+from xrspatial.sky_view_factor import sky_view_factor  # noqa
 from xrspatial.slope import slope  # noqa
 from xrspatial.surface_distance import surface_allocation  # noqa
 from xrspatial.surface_distance import surface_direction  # noqa

xrspatial/sky_view_factor.py

@@ -0,0 +1,253 @@
+"""Sky-view factor (SVF) for digital elevation models.
+
+SVF measures the fraction of the sky hemisphere visible from each cell,
+ranging from 0 (fully obstructed) to 1 (flat open terrain).  For each
+cell, rays are cast at *n_directions* evenly spaced azimuths out to
+*max_radius* cells.  Along each ray the maximum elevation angle to the
+horizon is recorded and SVF is computed as:
+
+    SVF(i,j) = 1 - mean(sin(max_horizon_angle(d)) for d in directions)
+
+References
+----------
+Zakek, K., Ostir, K., Kokalj, Z. (2011).  Sky-View Factor as a Relief
+Visualization Technique.  Remote Sensing, 3(2), 398-415.
+
+Yokoyama, R., Sidle, R.C., Noguchi, S. (2002).  Application of
+Topographic Shading to Terrain Visualization.  International Journal
+of Geographical Information Science, 16(5), 489-502.
+"""
+from __future__ import annotations
+
+from functools import partial
+from math import atan2 as _atan2, cos as _cos, pi as _pi, sin as _sin, sqrt as _sqrt
+from typing import Optional
+
+import numpy as np
+import xarray as xr
+from numba import cuda
+
+try:
+    import cupy
+except ImportError:
+    class cupy:
+        ndarray = False
+
+try:
+    import dask.array as da
+except ImportError:
+    da = None
+
+from xrspatial.dataset_support import supports_dataset
+from xrspatial.utils import (
+    ArrayTypeFunctionMapping,
+    _boundary_to_dask,
+    _validate_raster,
+    _validate_scalar,
+    cuda_args,
+    ngjit,
+)
+
+
+# ---------------------------------------------------------------------------
+# CPU kernel
+# ---------------------------------------------------------------------------
+
+@ngjit
+def _svf_cpu(data, max_radius, n_directions):
+    """Compute SVF over an entire 2-D array on the CPU."""
+    rows, cols = data.shape
+    out = np.empty((rows, cols), dtype=np.float64)
+    out[:] = np.nan
+
+    for y in range(rows):
+        for x in range(cols):
+            center = data[y, x]
+            if center != center:  # NaN check
+                continue
+
+            svf_sum = 0.0
+            for d in range(n_directions):
+                angle = 2.0 * _pi * d / n_directions
+                dx = _cos(angle)
+                dy = _sin(angle)
+
+                max_elev_angle = 0.0
+                for r in range(1, max_radius + 1):
+                    sx = x + int(round(r * dx))
+                    sy = y + int(round(r * dy))
+                    if sx < 0 or sx >= cols or sy < 0 or sy >= rows:
+                        break
+                    elev = data[sy, sx]
+                    if elev != elev:  # NaN
+                        break
+                    dz = elev - center
+                    dist = float(r)
+                    elev_angle = _atan2(dz, dist)
+                    if elev_angle > max_elev_angle:
+                        max_elev_angle = elev_angle
+
+                svf_sum += _sin(max_elev_angle)
+
+            out[y, x] = 1.0 - svf_sum / n_directions
+    return out
+
+
+# ---------------------------------------------------------------------------
+# GPU kernels
+# ---------------------------------------------------------------------------
+
+@cuda.jit
+def _svf_gpu(data, out, max_radius, n_directions):
+    """CUDA global kernel: one thread per cell."""
+    y, x = cuda.grid(2)
+    rows, cols = data.shape
+
+    if y >= rows or x >= cols:
+        return
+
+    center = data[y, x]
+    if center != center:  # NaN
+        return
+
+    svf_sum = 0.0
+    pi2 = 2.0 * 3.141592653589793
+
+    for d in range(n_directions):
+        angle = pi2 * d / n_directions
+        dx = _cos(angle)
+        dy = _sin(angle)
+
+        max_elev_angle = 0.0
+        for r in range(1, max_radius + 1):
+            sx = x + int(round(r * dx))
+            sy = y + int(round(r * dy))
+            if sx < 0 or sx >= cols or sy < 0 or sy >= rows:
+                break
+            elev = data[sy, sx]
+            if elev != elev:  # NaN
+                break
+            dz = elev - center
+            dist = float(r)
+            elev_angle = _atan2(dz, dist)
+            if elev_angle > max_elev_angle:
+                max_elev_angle = elev_angle
+
+        svf_sum += _sin(max_elev_angle)
+
+    out[y, x] = 1.0 - svf_sum / n_directions
+
+
+# ---------------------------------------------------------------------------
+# Backend wrappers
+# ---------------------------------------------------------------------------
+
+def _run_numpy(data, max_radius, n_directions):
+    data = data.astype(np.float64)
+    return _svf_cpu(data, max_radius, n_directions)
+
+
+def _run_cupy(data, max_radius, n_directions):
+    data = data.astype(cupy.float64)
+    out = cupy.full(data.shape, cupy.nan, dtype=cupy.float64)
+    griddim, blockdim = cuda_args(data.shape)
+    _svf_gpu[griddim, blockdim](data, out, max_radius, n_directions)
+    return out
+
+
+def _run_dask_numpy(data, max_radius, n_directions):
+    data = data.astype(np.float64)
+    _func = partial(_svf_cpu, max_radius=max_radius, n_directions=n_directions)
+    out = data.map_overlap(
+        _func,
+        depth=(max_radius, max_radius),
+        boundary=np.nan,
+        meta=np.array(()),
+    )
+    return out
+
+
+def _run_dask_cupy(data, max_radius, n_directions):
+    data = data.astype(cupy.float64)
+    _func = partial(_run_cupy, max_radius=max_radius, n_directions=n_directions)
+    out = data.map_overlap(
+        _func,
+        depth=(max_radius, max_radius),
+        boundary=cupy.nan,
+        meta=cupy.array(()),
+    )
+    return out
+
+
+# ---------------------------------------------------------------------------
+# Public API
+# ---------------------------------------------------------------------------
+
+@supports_dataset
+def sky_view_factor(
+    agg: xr.DataArray,
+    max_radius: int = 10,
+    n_directions: int = 16,
+    name: Optional[str] = 'sky_view_factor',
+) -> xr.DataArray:
+    """Compute the sky-view factor for each cell of a DEM.
+
+    SVF measures the fraction of the sky hemisphere visible from each
+    cell on a scale from 0 (fully obstructed) to 1 (flat open terrain).
+    Rays are cast at *n_directions* evenly spaced azimuths out to
+    *max_radius* cells, and the maximum elevation angle along each
+    ray determines the horizon obstruction.
+
+    Parameters
+    ----------
+    agg : xarray.DataArray or xr.Dataset
+        2D NumPy, CuPy, NumPy-backed Dask, or CuPy-backed Dask
+        xarray DataArray of elevation values.
+        If a Dataset is passed, the operation is applied to each
+        data variable independently.
+    max_radius : int, default=10
+        Maximum search distance in cells along each ray direction.
+        Cells within *max_radius* of the raster edge will be NaN.
+    n_directions : int, default=16
+        Number of azimuth directions to sample, evenly spaced
+        around 360 degrees.
+    name : str, default='sky_view_factor'
+        Name of the output DataArray.
+
+    Returns
+    -------
+    xarray.DataArray or xr.Dataset
+        2D array of SVF values in [0, 1] with the same shape, coords,
+        dims, and attrs as the input.  Edge cells use truncated rays
+        that stop at the raster boundary.
+
+    References
+    ----------
+    Zakek, K., Ostir, K., Kokalj, Z. (2011).  Sky-View Factor as a
+    Relief Visualization Technique.  Remote Sensing, 3(2), 398-415.
+
+    Examples
+    --------
+    >>> from xrspatial import sky_view_factor
+    >>> svf = sky_view_factor(dem, max_radius=100, n_directions=16)
+    """
+    _validate_raster(agg, func_name='sky_view_factor', name='agg')
+    _validate_scalar(max_radius, func_name='sky_view_factor',
+                     name='max_radius', dtype=int, min_val=1)
+    _validate_scalar(n_directions, func_name='sky_view_factor',
+                     name='n_directions', dtype=int, min_val=1)
+
+    mapper = ArrayTypeFunctionMapping(
+        numpy_func=_run_numpy,
+        cupy_func=_run_cupy,
+        dask_func=_run_dask_numpy,
+        dask_cupy_func=_run_dask_cupy,
+    )
+    out = mapper(agg)(agg.data, max_radius, n_directions)
+    return xr.DataArray(
+        out,
+        name=name,
+        coords=agg.coords,
+        dims=agg.dims,
+        attrs=agg.attrs,
+    )