hillshade computes gradients with np.gradient (2-point central difference) while GDAL's gdaldem hillshade uses Horn's method (weighted 8-neighbor Sobel kernel). The shade formula itself is correct; it's only the gradient kernel that differs.
This is inconsistent within xrspatial: slope and aspect already use Horn's method (the explicit Sobel kernel in their _run_numpy functions), but hillshade does not.
Measured impact
Tested on synthetic DEMs with realistic noise (σ=2 m on 30 m cells):
| Metric |
64×64 |
256×256 |
| Max abs diff vs GDAL |
0.076 |
0.100 |
| Mean abs diff |
0.013 |
0.016 |
| P95 abs diff |
0.034 |
0.040 |
On smooth (noise-free) data the two methods converge (max diff < 0.001). The gap grows with terrain roughness because np.gradient only uses 2 neighbors per axis while Horn averages over 6, making it more noise-sensitive.
The comment at hillshade.py:40 says "matches GDAL Horn method." It doesn't.
Suggested fix
Replace np.gradient(data, cellsize_y, cellsize_x) with the same Horn kernel used in slope.py. The GPU kernel (_gpu_calc_numba) already uses simple central differences and needs the same update.
All other terrain functions (slope, aspect, roughness, TPI, TRI) match GDAL within float precision.
hillshadecomputes gradients withnp.gradient(2-point central difference) while GDAL'sgdaldem hillshadeuses Horn's method (weighted 8-neighbor Sobel kernel). The shade formula itself is correct; it's only the gradient kernel that differs.This is inconsistent within xrspatial:
slopeandaspectalready use Horn's method (the explicit Sobel kernel in their_run_numpyfunctions), buthillshadedoes not.Measured impact
Tested on synthetic DEMs with realistic noise (σ=2 m on 30 m cells):
On smooth (noise-free) data the two methods converge (max diff < 0.001). The gap grows with terrain roughness because
np.gradientonly uses 2 neighbors per axis while Horn averages over 6, making it more noise-sensitive.The comment at
hillshade.py:40says "matches GDAL Horn method." It doesn't.Suggested fix
Replace
np.gradient(data, cellsize_y, cellsize_x)with the same Horn kernel used inslope.py. The GPU kernel (_gpu_calc_numba) already uses simple central differences and needs the same update.All other terrain functions (slope, aspect, roughness, TPI, TRI) match GDAL within float precision.