Open Source software is not, of course just about warm and fuzzies, great support, rapid development cycles, shared costs, etc., it’s also about getting your hands dirty with someone else’s code and implementing stuff more quickly and more intelligently because of it, and hopefully learning something in the process.  You don’t have to poke under the hood to drive the car, but sometimes it’s nice to anyway.

The two normal approaches to calculating slope in GIS raster tools (used in applications from ESRI’s Spatial Analyst and GDAL’s gdaldem) are Horn’s formula and Zevenbergen & Thorne’s formula, both integrated calculations of slope that do a bit of generalization in the process of calculating slope.  This is particularly useful if you are using slope for more than a localized calculations, for example in doing estimates of local light conditions.  That said, if you want to match local hydrologic slope (something which will more closely approximate a field measurement of slope), you are better off using a different formula altogether.  An article a colleague sent me: “When GIS Slope is Not What You Think” in this month’s The Forestry Source tipped me off to this, and I thought I’d do some quick and dirty digging in the GDAL code to see if I could implement something a little different.

Here is the standard Horn approach as currently implemented in GDAL:

float GDALSlopeHornAlg (float* afWin, float fDstNoDataValue, void* pData)
{
    const double radiansToDegrees = 180.0 / M_PI;
    GDALSlopeAlgData* psData = (GDALSlopeAlgData*)pData;
    double dx, dy, key;

    using std::max;
    dx = (max((afWin[0], afWin[3], afWin[3], afWin[6])) -
          max((afWin[2], afWin[5], afWin[5], afWin[8])))/psData->ewres;

    dy = ((afWin[6], afWin[7], afWin[7], afWin[8]) -
          max((afWin[0], afWin[1], afWin[1], afWin[2])))/psData->nsres;

    key = (dx * dx + dy * dy);

    if (psData->slopeFormat == 1)
        return (float) (atan(sqrt(key) / (8*psData->scale)) * radiansToDegrees);
    else
        return (float) (100*(sqrt(key) / (8*psData->scale)));
}

Borrowing from gdaldem roughness we implement this:

float GDALSlopeHydroAlg (float* afWin, float fDstNoDataValue, void* pData)

{
    // Hydrologic Slope is the max
    //  local slope btw center cell and adjacent cells

    const double radiansToDegrees = 180.0 / M_PI;
    GDALSlopeAlgData* psData = (GDALSlopeAlgData*)pData;
                double key;

    float pafLocalMin = afWin[0];
     float pafLocalMax = afWin[0];

    for ( int k = 1; k < 9; k++)
    {
        if (afWin[k] < pafLocalMin)
        {
            pafLocalMin=afWin[k];
        }

        if (afWin[k] > pafLocalMax)
        {
            pafLocalMax=afWin[k];
        }
    }

    key = afWin[4] - pafLocalMin;

    if (psData->slopeFormat == 1)
        return (float) (atan(sqrt(key) / (2*psData->scale)) * radiansToDegrees);
    else
        return (float) (100*(sqrt(key) / (2*psData->scale)));
}

Output is rougher than the original, as expected:

And looking at the stats we see a larger standard deviation for local hydrologic slope:

Z:\>gdalinfo Hydrologic_slope.tif
Driver: GTiff/GeoTIFF
Files: Hydrologic_slope.tif
...

Metadata:
STATISTICS_MINIMUM=0
STATISTICS_MAXIMUM=46.168109893799
STATISTICS_MEAN=10.206419514979
STATISTICS_MEDIAN=2.8041767686283e-266
STATISTICS_MODE=1.4463930749435e-307
STATISTICS_STDDEV=5.3794018830967
LAYER_TYPE=athematic

Z:\>gdalinfo Zevenbergen_Thorne.tif
Driver: GTiff/GeoTIFF
Files: Zevenbergen_Thorne.tif
...
Metadata:
STATISTICS_MINIMUM=0
STATISTICS_MAXIMUM=52.424499511719
STATISTICS_MEAN=3.0666069784504
STATISTICS_MEDIAN=4.331070843283e-184
STATISTICS_MODE=-4
STATISTICS_STDDEV=3.567882216098
LAYER_TYPE=athematic

A couple of additional notes.  First, this is 2.5ft resolution data.  At this resolution, the differences in algorithms may be somewhat academic.  Also, of note– the range for the Zevenbergen-Thorn slope algorithm is greater than that for the local hydrologic slope.  This is an outcome I did not anticipate.

Update:

Hmmm, looks like I have more revisions and questions from the GDAL dev board. May have to revisit this… .  Forgot about the corners of my 3×3 window… .