While validating 3dep.deriv.slope for terrain-based geolocation analysis over CONUS, I found that the returned slope decreases with slope_scale_length beyond what terrain smoothing can explain: at 40 m on 1 m pixels it reads about 2/3 of the least-squares plane-fit slope over the same window.
Field measurement (3DEP 1 m tile USGS_1M_11_x65y418_NV_EastCentral_2021_D21, n=152 points, slope_scale_length=40): server slope = 0.68 x plane-fit slope with correlation 0.997, i.e. a pure scale factor rather than scatter.
Mechanism (GdalRaster::computeSlopeAspect, unchanged on main at v5.5.2): the accumulated sum(w * val * c) is divided by wsum_dx * dx * kHalf with wsum_dx = sum(w * |c|). That treats every column offset as if it sat at distance kHalf, so the estimator is only exact when all non-zero offsets are +/-kHalf, i.e. the 3x3 case. Emulating the exact loop on a planar ramp (true gradient unambiguous and scale-independent) recovers a fraction (2k+1)/(3k) of the gradient: k=1 gives 1.000, k=5 gives 0.733, k=20 (the 41x41 window from slope_scale_length=40 at 1 m) gives 0.683, and the limit for large windows is 2/3. The k=20 value matches the field measurement, and it explains why small-window tests pass while large windows read low.
Fix: the weighted least-squares gradient over the same symmetric window normalises by sum(w * c^2) * dx, which recovers a plane exactly for every k and is identical to the current formula at k=1 (both reduce to Horn's divisor 8*dx). Concretely, accumulate wsum_dx += w * c * c and wsum_dy += w * r * r instead of w * |c| / w * |r|, and drop the * kHalf from both divisors. Downstream quantities such as dh/tan(slope) currently inherit the inverse factor (about 1.46x at 40 m) silently. Happy to open a PR with the two-line change plus a planar-ramp selftest; if the current kernel is intentional, the k-dependence at least deserves a note on the field.
While validating
3dep.deriv.slopefor terrain-based geolocation analysis over CONUS, I found that the returned slope decreases withslope_scale_lengthbeyond what terrain smoothing can explain: at 40 m on 1 m pixels it reads about 2/3 of the least-squares plane-fit slope over the same window.Field measurement (3DEP 1 m tile USGS_1M_11_x65y418_NV_EastCentral_2021_D21, n=152 points, slope_scale_length=40): server slope = 0.68 x plane-fit slope with correlation 0.997, i.e. a pure scale factor rather than scatter.
Mechanism (GdalRaster::computeSlopeAspect, unchanged on main at v5.5.2): the accumulated
sum(w * val * c)is divided bywsum_dx * dx * kHalfwithwsum_dx = sum(w * |c|). That treats every column offset as if it sat at distance kHalf, so the estimator is only exact when all non-zero offsets are +/-kHalf, i.e. the 3x3 case. Emulating the exact loop on a planar ramp (true gradient unambiguous and scale-independent) recovers a fraction (2k+1)/(3k) of the gradient: k=1 gives 1.000, k=5 gives 0.733, k=20 (the 41x41 window from slope_scale_length=40 at 1 m) gives 0.683, and the limit for large windows is 2/3. The k=20 value matches the field measurement, and it explains why small-window tests pass while large windows read low.Fix: the weighted least-squares gradient over the same symmetric window normalises by
sum(w * c^2) * dx, which recovers a plane exactly for every k and is identical to the current formula at k=1 (both reduce to Horn's divisor 8*dx). Concretely, accumulatewsum_dx += w * c * candwsum_dy += w * r * rinstead ofw * |c|/w * |r|, and drop the* kHalffrom both divisors. Downstream quantities such as dh/tan(slope) currently inherit the inverse factor (about 1.46x at 40 m) silently. Happy to open a PR with the two-line change plus a planar-ramp selftest; if the current kernel is intentional, the k-dependence at least deserves a note on the field.