Fix hillshade gradient to use Horn's method (#1175)

hillshade was using np.gradient (2-point central difference) for
gradients while slope and aspect already used Horn's 8-neighbor Sobel
kernel. This caused measurable divergence from GDAL on noisy terrain
(up to 0.1 shading units on rough DEMs).

Replace np.gradient in _run_numpy with an @ngjit Horn kernel matching
slope.py. Update the CUDA kernel (_gpu_calc_numba) with the same
8-neighbor stencil. Both dask backends inherit the fix automatically
through their chunk functions.

Update test_hillshade_gdal_equivalence to validate against Horn's
method reference instead of np.gradient. Add
test_hillshade_horn_vs_central_diff to confirm the two methods
produce meaningfully different results on noisy data.

Author: brendancol

xrspatial/hillshade.py

@@ -16,10 +16,46 @@
 from .utils import (_boundary_to_dask, _pad_array, _validate_boundary,
                     _validate_raster, _validate_scalar,
                     calc_cuda_dims, get_dataarray_resolution,
-                    has_cuda_and_cupy, is_cupy_array, is_cupy_backed)
+                    has_cuda_and_cupy, is_cupy_array, is_cupy_backed,
+                    ngjit)
 from .dataset_support import supports_dataset
 
 
+@ngjit
+def _horn_hillshade(data, cellsize_x, cellsize_y, sin_alt, cos_alt,
+                    sin_az, cos_az):
+    """Hillshade using Horn's method (weighted 8-neighbor Sobel kernel).
+
+    Gradient kernel matches slope.py and GDAL gdaldem hillshade.
+    """
+    rows, cols = data.shape
+    out = np.empty(data.shape, dtype=np.float32)
+    out[:] = np.nan
+    for y in range(1, rows - 1):
+        for x in range(1, cols - 1):
+            a = data[y + 1, x - 1]
+            b = data[y + 1, x]
+            c = data[y + 1, x + 1]
+            d = data[y, x - 1]
+            f = data[y, x + 1]
+            g = data[y - 1, x - 1]
+            h = data[y - 1, x]
+            i = data[y - 1, x + 1]
+
+            dz_dx = ((c + 2 * f + i) - (a + 2 * d + g)) / (8 * cellsize_x)
+            dz_dy = ((g + 2 * h + i) - (a + 2 * b + c)) / (8 * cellsize_y)
+
+            xx_plus_yy = dz_dx * dz_dx + dz_dy * dz_dy
+            shaded = (sin_alt + cos_alt * (dz_dy * cos_az - dz_dx * sin_az)) \
+                / math.sqrt(1.0 + xx_plus_yy)
+            if shaded < 0.0:
+                shaded = 0.0
+            elif shaded > 1.0:
+                shaded = 1.0
+            out[y, x] = shaded
+    return out
+
+
 def _run_numpy(data, azimuth=225, angle_altitude=25,
                cellsize_x=1.0, cellsize_y=1.0, boundary='nan'):
     if boundary != 'nan':
@@ -37,25 +73,8 @@ def _run_numpy(data, azimuth=225, angle_altitude=25,
     sin_az = np.sin(az_rad)
     cos_az = np.cos(az_rad)
 
-    # Gradient with actual cell spacing (matches GDAL Horn method)
-    dy, dx = np.gradient(data, cellsize_y, cellsize_x)
-    xx_plus_yy = dx * dx + dy * dy
-
-    # GDAL-equivalent hillshade formula (simplified from the original
-    # trig-heavy version; see issue #748 and GDAL gdaldem_lib.cpp):
-    #   shaded = (sin(alt) + cos(alt) * sqrt(xx+yy) * sin(aspect - az))
-    #            / sqrt(1 + xx+yy)
-    # where aspect = atan2(dy, dx), expanded inline:
-    #   sin(aspect - az) = (dy*cos(az) - dx*sin(az)) / sqrt(xx+yy)
-    # so sqrt(xx+yy) cancels, giving:
-    shaded = (sin_alt + cos_alt * (dy * cos_az - dx * sin_az)) \
-        / np.sqrt(1.0 + xx_plus_yy)
-
-    # Clamp negatives (shadow) then scale to [0, 1]
-    result = np.clip(shaded, 0.0, 1.0)
-    result[(0, -1), :] = np.nan
-    result[:, (0, -1)] = np.nan
-    return result
+    return _horn_hillshade(data, cellsize_x, cellsize_y,
+                           sin_alt, cos_alt, sin_az, cos_az)
 
 
 def _run_dask_numpy(data, azimuth, angle_altitude,
@@ -86,15 +105,27 @@ def _gpu_calc_numba(
 
     i, j = cuda.grid(2)
     if i > 0 and i < data.shape[0]-1 and j > 0 and j < data.shape[1] - 1:
-        dx = (data[i, j+1] - data[i, j-1]) / (2.0 * cellsize_x)
-        dy = (data[i+1, j] - data[i-1, j]) / (2.0 * cellsize_y)
-
-        xx_plus_yy = dx * dx + dy * dy
-        shaded = (sin_alt + cos_alt * (dy * cos_az - dx * sin_az)) \
+        # Horn's method (weighted 8-neighbor Sobel kernel)
+        a = data[i + 1, j - 1]
+        b = data[i + 1, j]
+        c = data[i + 1, j + 1]
+        d = data[i, j - 1]
+        f = data[i, j + 1]
+        g = data[i - 1, j - 1]
+        h = data[i - 1, j]
+        ii = data[i - 1, j + 1]
+
+        dz_dx = ((c + 2.0 * f + ii) - (a + 2.0 * d + g)) / (8.0 * cellsize_x)
+        dz_dy = ((g + 2.0 * h + ii) - (a + 2.0 * b + c)) / (8.0 * cellsize_y)
+
+        xx_plus_yy = dz_dx * dz_dx + dz_dy * dz_dy
+        shaded = (sin_alt + cos_alt * (dz_dy * cos_az - dz_dx * sin_az)) \
             / math.sqrt(1.0 + xx_plus_yy)
 
         if shaded < 0.0:
             shaded = 0.0
+        elif shaded > 1.0:
+            shaded = 1.0
         output[i, j] = shaded
 
 

xrspatial/tests/test_hillshade.py

@@ -96,40 +96,93 @@ def test_hillshade_resolution_sensitivity():
         "hillshade should be sensitive to cell resolution"
 
 
-def test_hillshade_gdal_equivalence():
+def _horn_reference(data, cellsize_x, cellsize_y, azimuth, altitude):
+    """Pure-numpy Horn's method reference for testing.
+
+    Computes gradients with the 8-neighbor Sobel kernel and applies the
+    GDAL hillshade formula.
     """
-    Verify that _run_numpy matches the GDAL hillshade formula directly.
+    az_rad = azimuth * np.pi / 180.
+    alt_rad = altitude * np.pi / 180.
 
-    GDAL formula (gdaldem_lib.cpp):
-        aspect = atan2(dy, dx)
-        shaded = (sin(alt) + cos(alt)*sqrt(xx+yy)*sin(aspect-az))
-                 / sqrt(1 + xx+yy)
-    where dx, dy are gradients divided by cell spacing.
+    # Horn / Sobel kernel gradients (same as GDAL gdaldem)
+    a = data[2:, :-2]
+    b = data[2:, 1:-1]
+    c = data[2:, 2:]
+    d = data[1:-1, :-2]
+    f = data[1:-1, 2:]
+    g = data[:-2, :-2]
+    h = data[:-2, 1:-1]
+    i = data[:-2, 2:]
+
+    dz_dx = ((c + 2 * f + i) - (a + 2 * d + g)) / (8 * cellsize_x)
+    dz_dy = ((g + 2 * h + i) - (a + 2 * b + c)) / (8 * cellsize_y)
+
+    xx_plus_yy = dz_dx**2 + dz_dy**2
+    shaded = (
+        np.sin(alt_rad) + np.cos(alt_rad) * (dz_dy * np.cos(az_rad) - dz_dx * np.sin(az_rad))
+    ) / np.sqrt(1 + xx_plus_yy)
+    return np.clip(shaded, 0.0, 1.0)
+
+
+def test_hillshade_gdal_equivalence():
+    """
+    Verify that _run_numpy matches the GDAL hillshade formula with
+    Horn's method (weighted 8-neighbor Sobel kernel) for gradients.
     """
     rng = np.random.default_rng(99)
     data = rng.random((20, 20)).astype(np.float32) * 500
     cellsize_x, cellsize_y = 30.0, 30.0
     azimuth, altitude = 315.0, 45.0
 
-    # Reference: direct GDAL formula
+    ref = _horn_reference(data, cellsize_x, cellsize_y, azimuth, altitude)
+
+    result = _run_numpy(data, azimuth=azimuth, angle_altitude=altitude,
+                        cellsize_x=cellsize_x, cellsize_y=cellsize_y)
+
+    interior = slice(1, -1)
+    assert_allclose(result[interior, interior], ref, atol=1e-6)
+
+
+def test_hillshade_horn_vs_central_diff():
+    """Horn's method should differ from simple central differences on noisy data.
+
+    This validates that the implementation actually uses the 8-neighbor Sobel
+    kernel rather than 2-point np.gradient, which was the bug in #1175.
+    """
+    rng = np.random.default_rng(1175)
+    # Smooth ramp with noise to amplify the difference between methods
+    rows, cols = 64, 64
+    base = np.tile(np.linspace(0, 500, cols), (rows, 1))
+    noise = rng.normal(0, 2, (rows, cols))
+    data = (base + noise).astype(np.float32)
+    cellsize_x, cellsize_y = 30.0, 30.0
+    azimuth, altitude = 315.0, 45.0
+
+    # Our result (should be Horn)
+    result = _run_numpy(data, azimuth=azimuth, angle_altitude=altitude,
+                        cellsize_x=cellsize_x, cellsize_y=cellsize_y)
+    interior = slice(1, -1)
+
+    # Horn reference
+    horn_ref = _horn_reference(data, cellsize_x, cellsize_y, azimuth, altitude)
+    assert_allclose(result[interior, interior], horn_ref, atol=1e-6)
+
+    # np.gradient reference (old method) -- should NOT match
     az_rad = azimuth * np.pi / 180.
     alt_rad = altitude * np.pi / 180.
     dy, dx = np.gradient(data, cellsize_y, cellsize_x)
     xx_plus_yy = dx**2 + dy**2
-    aspect = np.arctan2(dy, dx)
-    gdal_shaded = (
-        np.sin(alt_rad)
-        + np.cos(alt_rad) * np.sqrt(xx_plus_yy) * np.sin(aspect - az_rad)
+    old_shaded = (
+        np.sin(alt_rad) + np.cos(alt_rad) * (dy * np.cos(az_rad) - dx * np.sin(az_rad))
     ) / np.sqrt(1 + xx_plus_yy)
-    gdal_ref = np.clip(gdal_shaded, 0.0, 1.0)
+    old_ref = np.clip(old_shaded, 0.0, 1.0)
 
-    # Our implementation
-    result = _run_numpy(data, azimuth=azimuth, angle_altitude=altitude,
-                        cellsize_x=cellsize_x, cellsize_y=cellsize_y)
-
-    interior = slice(1, -1)
-    assert_allclose(result[interior, interior],
-                    gdal_ref[interior, interior], atol=1e-6)
+    max_diff = np.max(np.abs(result[interior, interior] - old_ref[interior, interior]))
+    assert max_diff > 0.01, (
+        f"Expected meaningful difference between Horn and central-diff, "
+        f"got max diff {max_diff}"
+    )
 
 
 @dask_array_available