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… .