From 06141b9a008458735bc8528850b969f80ab93f57 Mon Sep 17 00:00:00 2001
From: Brendan Collins <brendancol@gmail.com>
Date: Wed, 4 Mar 2026 16:56:01 -0800
Subject: [PATCH 1/4] Add sky-view factor function (#962)

New xrspatial.sky_view_factor module with all four backends:
NumPy (@ngjit), CuPy (@cuda.jit), Dask+NumPy, and Dask+CuPy.
---
 xrspatial/__init__.py        |   1 +
 xrspatial/sky_view_factor.py | 258 +++++++++++++++++++++++++++++++++++
 2 files changed, 259 insertions(+)
 create mode 100644 xrspatial/sky_view_factor.py

diff --git a/xrspatial/__init__.py b/xrspatial/__init__.py
index bb30a13a..56dc235b 100644
--- a/xrspatial/__init__.py
+++ b/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
diff --git a/xrspatial/sky_view_factor.py b/xrspatial/sky_view_factor.py
new file mode 100644
index 00000000..df8b0bf8
--- /dev/null
+++ b/xrspatial/sky_view_factor.py
@@ -0,0 +1,258 @@
+"""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(max_radius, rows - max_radius):
+        for x in range(max_radius, cols - max_radius):
+            center = data[y, x]
+            if center != center:  # NaN check
+                continue
+
+            svf_sum = 0.0
+            valid_dirs = 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)
+                valid_dirs += 1
+
+            if valid_dirs > 0:
+                out[y, x] = 1.0 - svf_sum / valid_dirs
+    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 < max_radius or y >= rows - max_radius:
+        return
+    if x < max_radius or x >= cols - max_radius:
+        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 within *max_radius*
+        of the boundary are NaN.
+
+    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,
+    )

From 0e1cbdb268dbf17258178c2d845b4bd4f7a5931d Mon Sep 17 00:00:00 2001
From: Brendan Collins <brendancol@gmail.com>
Date: Wed, 4 Mar 2026 16:58:55 -0800
Subject: [PATCH 2/4] Add tests for sky_view_factor, fix edge cell handling
 (#962)

25 tests covering correctness, NaN handling, range validation,
cross-backend equivalence (NumPy/Dask/CuPy/Dask+CuPy), and
Dataset support. Edge cells now compute with truncated rays
instead of returning NaN, keeping numpy and dask paths consistent.
---
 xrspatial/sky_view_factor.py            |  17 +-
 xrspatial/tests/test_sky_view_factor.py | 238 ++++++++++++++++++++++++
 2 files changed, 244 insertions(+), 11 deletions(-)
 create mode 100644 xrspatial/tests/test_sky_view_factor.py

diff --git a/xrspatial/sky_view_factor.py b/xrspatial/sky_view_factor.py
index df8b0bf8..d83a12de 100644
--- a/xrspatial/sky_view_factor.py
+++ b/xrspatial/sky_view_factor.py
@@ -60,14 +60,13 @@ def _svf_cpu(data, max_radius, n_directions):
     out = np.empty((rows, cols), dtype=np.float64)
     out[:] = np.nan
 
-    for y in range(max_radius, rows - max_radius):
-        for x in range(max_radius, cols - max_radius):
+    for y in range(rows):
+        for x in range(cols):
             center = data[y, x]
             if center != center:  # NaN check
                 continue
 
             svf_sum = 0.0
-            valid_dirs = 0
             for d in range(n_directions):
                 angle = 2.0 * _pi * d / n_directions
                 dx = _cos(angle)
@@ -89,10 +88,8 @@ def _svf_cpu(data, max_radius, n_directions):
                         max_elev_angle = elev_angle
 
                 svf_sum += _sin(max_elev_angle)
-                valid_dirs += 1
 
-            if valid_dirs > 0:
-                out[y, x] = 1.0 - svf_sum / valid_dirs
+            out[y, x] = 1.0 - svf_sum / n_directions
     return out
 
 
@@ -106,9 +103,7 @@ def _svf_gpu(data, out, max_radius, n_directions):
     y, x = cuda.grid(2)
     rows, cols = data.shape
 
-    if y < max_radius or y >= rows - max_radius:
-        return
-    if x < max_radius or x >= cols - max_radius:
+    if y >= rows or x >= cols:
         return
 
     center = data[y, x]
@@ -223,8 +218,8 @@ def sky_view_factor(
     -------
     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 within *max_radius*
-        of the boundary are NaN.
+        dims, and attrs as the input.  Edge cells use truncated rays
+        that stop at the raster boundary.
 
     References
     ----------
diff --git a/xrspatial/tests/test_sky_view_factor.py b/xrspatial/tests/test_sky_view_factor.py
new file mode 100644
index 00000000..00c771b5
--- /dev/null
+++ b/xrspatial/tests/test_sky_view_factor.py
@@ -0,0 +1,238 @@
+import numpy as np
+import pytest
+import xarray as xr
+
+from xrspatial import sky_view_factor
+from xrspatial.tests.general_checks import (
+    assert_numpy_equals_cupy,
+    assert_numpy_equals_dask_cupy,
+    assert_numpy_equals_dask_numpy,
+    create_test_raster,
+    cuda_and_cupy_available,
+    dask_array_available,
+    general_output_checks,
+)
+
+
+# ---------------------------------------------------------------------------
+# Fixtures
+# ---------------------------------------------------------------------------
+
+@pytest.fixture
+def flat_surface():
+    """Constant elevation -- SVF should be 1.0 everywhere in the interior."""
+    return np.full((30, 30), 100.0, dtype=np.float64)
+
+
+@pytest.fixture
+def valley_surface():
+    """V-shaped valley along the x axis.
+
+    Cross-section:  z = abs(x - center) * slope_factor
+    Cells at the valley floor have obstructed sky on both sides,
+    so SVF < 1.
+    """
+    rows, cols = 40, 40
+    x = np.arange(cols, dtype=np.float64)
+    z = np.abs(x - cols / 2.0) * 5.0
+    return np.broadcast_to(z, (rows, cols)).copy()
+
+
+@pytest.fixture
+def ridge_surface():
+    """Cone-shaped peak in the center.
+
+    Cells at the summit have unobstructed sky (SVF close to 1),
+    cells at the base have partial obstruction.
+    """
+    rows, cols = 50, 50
+    y, x = np.mgrid[:rows, :cols].astype(np.float64)
+    cy, cx = rows / 2.0, cols / 2.0
+    dist = np.sqrt((y - cy)**2 + (x - cx)**2)
+    return np.maximum(0.0, 200.0 - dist * 10.0)
+
+
+# ---------------------------------------------------------------------------
+# Correctness tests
+# ---------------------------------------------------------------------------
+
+def test_flat_surface_svf(flat_surface):
+    """Flat terrain: no horizon obstruction, SVF should be 1.0."""
+    agg = create_test_raster(flat_surface)
+    result = sky_view_factor(agg, max_radius=5, n_directions=16)
+    general_output_checks(agg, result)
+    interior = result.data[5:-5, 5:-5]
+    np.testing.assert_allclose(interior, 1.0, atol=1e-10)
+
+
+def test_valley_svf_less_than_one(valley_surface):
+    """Valley floor cells should have SVF < 1 due to obstruction."""
+    agg = create_test_raster(valley_surface)
+    result = sky_view_factor(agg, max_radius=8, n_directions=16)
+    # Valley floor at col=20, interior rows
+    floor_svf = result.data[20, 20]
+    assert np.isfinite(floor_svf)
+    assert floor_svf < 1.0, f"Expected SVF < 1 at valley floor, got {floor_svf}"
+
+
+def test_ridge_summit_high_svf(ridge_surface):
+    """Summit of a cone has nearly unobstructed sky."""
+    agg = create_test_raster(ridge_surface)
+    result = sky_view_factor(agg, max_radius=10, n_directions=16)
+    summit_svf = result.data[25, 25]
+    assert np.isfinite(summit_svf)
+    assert summit_svf > 0.9, f"Expected SVF > 0.9 at summit, got {summit_svf}"
+
+
+def test_svf_range():
+    """All valid SVF values should be in [0, 1]."""
+    rng = np.random.default_rng(962)
+    data = rng.random((30, 30)).astype(np.float64) * 500
+    agg = create_test_raster(data)
+    result = sky_view_factor(agg, max_radius=5, n_directions=8)
+    valid = result.data[~np.isnan(result.data)]
+    assert np.all(valid >= 0.0), "SVF values below 0"
+    assert np.all(valid <= 1.0), "SVF values above 1"
+
+
+def test_more_directions_same_flat(flat_surface):
+    """Different n_directions on a flat surface should still give SVF=1."""
+    agg = create_test_raster(flat_surface)
+    for n_dir in [4, 8, 16, 32]:
+        result = sky_view_factor(agg, max_radius=5, n_directions=n_dir)
+        interior = result.data[5:-5, 5:-5]
+        np.testing.assert_allclose(interior, 1.0, atol=1e-10,
+                                   err_msg=f"Failed for n_directions={n_dir}")
+
+
+# ---------------------------------------------------------------------------
+# Edge / NaN handling
+# ---------------------------------------------------------------------------
+
+def test_edge_cells_computed():
+    """Edge cells should be computed (with truncated rays), not NaN."""
+    data = np.random.default_rng(962).random((20, 20)).astype(np.float64) * 100
+    agg = create_test_raster(data)
+    radius = 3
+    result = sky_view_factor(agg, max_radius=radius, n_directions=8)
+    # All cells should have valid values (no NaN except where input is NaN)
+    assert np.all(np.isfinite(result.data))
+
+
+def test_nan_in_input():
+    """NaN in the DEM: center cell NaN should produce NaN output."""
+    data = np.random.default_rng(962).random((20, 20)).astype(np.float64) * 100
+    data[10, 10] = np.nan
+    agg = create_test_raster(data)
+    result = sky_view_factor(agg, max_radius=3, n_directions=8)
+    assert np.isnan(result.data[10, 10])
+
+
+def test_single_cell_raster():
+    """A 1x1 raster should return SVF=1 (no neighbors to obstruct)."""
+    data = np.array([[42.0]])
+    agg = create_test_raster(data)
+    result = sky_view_factor(agg, max_radius=1, n_directions=4)
+    np.testing.assert_allclose(result.data[0, 0], 1.0, atol=1e-10)
+
+
+# ---------------------------------------------------------------------------
+# Validation tests
+# ---------------------------------------------------------------------------
+
+def test_invalid_max_radius():
+    data = np.ones((10, 10), dtype=np.float64)
+    agg = create_test_raster(data)
+    with pytest.raises(ValueError, match="max_radius"):
+        sky_view_factor(agg, max_radius=0)
+
+
+def test_invalid_n_directions():
+    data = np.ones((10, 10), dtype=np.float64)
+    agg = create_test_raster(data)
+    with pytest.raises(ValueError, match="n_directions"):
+        sky_view_factor(agg, n_directions=0)
+
+
+def test_invalid_input_type():
+    with pytest.raises(TypeError, match="DataArray"):
+        sky_view_factor(np.ones((5, 5)))
+
+
+# ---------------------------------------------------------------------------
+# Cross-backend tests
+# ---------------------------------------------------------------------------
+
+@dask_array_available
+@pytest.mark.parametrize("size", [(20, 20), (30, 40)])
+@pytest.mark.parametrize("dtype", [np.float32, np.float64])
+def test_numpy_equals_dask(random_data, size, dtype):
+    numpy_agg = create_test_raster(random_data, backend='numpy')
+    dask_agg = create_test_raster(random_data, backend='dask', chunks=(10, 10))
+    func = lambda agg: sky_view_factor(agg, max_radius=3, n_directions=8)
+    assert_numpy_equals_dask_numpy(numpy_agg, dask_agg, func, nan_edges=False)
+
+
+@cuda_and_cupy_available
+@pytest.mark.parametrize("size", [(20, 20), (30, 40)])
+@pytest.mark.parametrize("dtype", [np.float32, np.float64])
+def test_numpy_equals_cupy(random_data, size, dtype):
+    numpy_agg = create_test_raster(random_data, backend='numpy')
+    cupy_agg = create_test_raster(random_data, backend='cupy')
+    func = lambda agg: sky_view_factor(agg, max_radius=3, n_directions=8)
+    assert_numpy_equals_cupy(numpy_agg, cupy_agg, func, nan_edges=False)
+
+
+@dask_array_available
+@cuda_and_cupy_available
+@pytest.mark.parametrize("size", [(20, 20), (30, 40)])
+@pytest.mark.parametrize("dtype", [np.float32, np.float64])
+def test_numpy_equals_dask_cupy(random_data, size, dtype):
+    numpy_agg = create_test_raster(random_data, backend='numpy')
+    dask_cupy_agg = create_test_raster(random_data, backend='dask+cupy',
+                                       chunks=(10, 10))
+    func = lambda agg: sky_view_factor(agg, max_radius=3, n_directions=8)
+    assert_numpy_equals_dask_cupy(numpy_agg, dask_cupy_agg, func,
+                                  nan_edges=False, atol=1e-6, rtol=1e-6)
+
+
+# ---------------------------------------------------------------------------
+# Dataset support
+# ---------------------------------------------------------------------------
+
+def test_dataset_support():
+    data = np.random.default_rng(962).random((20, 20)).astype(np.float64) * 100
+    da1 = xr.DataArray(data, dims=['y', 'x'], attrs={'res': (0.5, 0.5)})
+    da1['y'] = np.linspace(9.5, 0, 20)
+    da1['x'] = np.linspace(0, 9.5, 20)
+    ds = xr.Dataset({'elev1': da1, 'elev2': da1 * 2})
+    result = sky_view_factor(ds, max_radius=3, n_directions=8)
+    assert isinstance(result, xr.Dataset)
+    assert set(result.data_vars) == {'elev1', 'elev2'}
+    for var in result.data_vars:
+        expected = sky_view_factor(ds[var], max_radius=3, n_directions=8,
+                                   name=var)
+        np.testing.assert_allclose(
+            result[var].data, expected.data, equal_nan=True)
+
+
+# ---------------------------------------------------------------------------
+# Known reference value test
+# ---------------------------------------------------------------------------
+
+def test_known_svf_simple_ramp():
+    """A tilted plane: SVF should be < 1 for cells looking uphill.
+
+    For a plane tilted at 45 degrees (z = x), the horizon angle
+    along the uphill direction is 45 degrees, so SVF should be
+    strictly less than 1 at interior cells.
+    """
+    rows, cols = 30, 30
+    x = np.arange(cols, dtype=np.float64)
+    data = np.broadcast_to(x * 1.0, (rows, cols)).copy()
+    agg = create_test_raster(data)
+    result = sky_view_factor(agg, max_radius=5, n_directions=16)
+    center = result.data[15, 15]
+    assert np.isfinite(center)
+    assert center < 1.0, f"Expected SVF < 1 on tilted plane, got {center}"
+    assert center > 0.5, f"Expected SVF > 0.5 on gentle tilt, got {center}"

From 49874670804e4880a4e3f13ae4057879917e8449 Mon Sep 17 00:00:00 2001
From: Brendan Collins <brendancol@gmail.com>
Date: Wed, 4 Mar 2026 17:54:38 -0800
Subject: [PATCH 3/4] Add docs and user guide notebook for sky_view_factor
 (#962)

- Add Sky-View Factor entry to docs/source/reference/surface.rst
- Add examples/user_guide/22_Sky_View_Factor.ipynb with parameter
  exploration and hillshade comparison
---
 docs/source/reference/surface.rst            |   7 ++
 examples/user_guide/22_Sky_View_Factor.ipynb | 101 +++++++++++++++++++
 2 files changed, 108 insertions(+)
 create mode 100644 examples/user_guide/22_Sky_View_Factor.ipynb

diff --git a/docs/source/reference/surface.rst b/docs/source/reference/surface.rst
index 5985b309..5632db2e 100644
--- a/docs/source/reference/surface.rst
+++ b/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::
diff --git a/examples/user_guide/22_Sky_View_Factor.ipynb b/examples/user_guide/22_Sky_View_Factor.ipynb
new file mode 100644
index 00000000..b69004aa
--- /dev/null
+++ b/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
+}
\ No newline at end of file

From 9b75726bd9d7456ee57cc00365e5adf862466e44 Mon Sep 17 00:00:00 2001
From: Brendan Collins <brendancol@gmail.com>
Date: Wed, 4 Mar 2026 17:54:52 -0800
Subject: [PATCH 4/4] Add sky_view_factor to README feature matrix (#962)

---
 README.md | 1 +
 1 file changed, 1 insertion(+)

diff --git a/README.md b/README.md
index 18ac4b8a..f9a7e784 100644
--- a/README.md
+++ b/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) | ✅️ | ✅️ | ✅️ | ✅️ |