GDAL使用DEM数据计算地形指数

零、 前言 本文是接上文《GDAL使用DEM数据计算坡度坡向》，还是一篇关于DEM计算地形指数的一篇文章。这里所要计算的地形指数主要包括以下三个指数：地形耐用指数（Terrain Ruggedness Index，TPI）、地形位置指数（Topographic Position Index，TRI）和粗糙度（roughness）。 上面三个地形指数都是在一个3×3的窗口中计算的，3×3的窗口定义如图1所示，下文中的公式中的符号与图1中的符号一致。下文中的算法及其含义具体可以参考Wilson等人2007年的著作《Marine Geodesy》。

图1 3×3的地形窗口及其编号一、 地形指数简介与计算公式 1、地形耐用指数（Terrain Ruggedness Index，TPI） TRI是指中心点高程与周围高程的差的平均值。如公式（1）所示，式（1）中的e表示中心点高程值，ei如图1所示，也就是周围的高程值。


（1） 其中关于TRI的pdf可以在这里（http://download.osgeo.org/qgis/doc/reference-docs/Terrain_Ruggedness_InInd.pdf）下载到。一个简单的示例如图2所示：

图2 TRI计算示例 2、地形位置指数（Topographic Position Index，TRI） TPI是中心点高程值减去周围高程的平均值。计算公式如式（2）所示，其中式中的符号与式（1）中的意义相同，此处不再赘述。


（2）3、粗糙度（roughness）粗糙度的定义是指中心点高程与周围高程中的最大值和最小值的高程差。计算公式如式（3）所示，符号同上。


（3）二、 算法编写 有了上面的计算公式，那么我们就可以利用之前的3×3模板算子来编写地形指数的函数代码了。从上面的公式看，这三种地形指数的计算非常简单，代码基本都不超过10行，分别如下： 1、地形耐用指数（Terrain Ruggedness Index，TPI）

float GDALTRIAlg (float* afWin, float fDstNoDataValue,void* pData)
{
         return(fabs(afWin[0]-afWin[4]) +
                   fabs(afWin[1]-afWin[4])+
                   fabs(afWin[2]-afWin[4])+
                   fabs(afWin[3]-afWin[4])+
                   fabs(afWin[5]-afWin[4])+
                   fabs(afWin[6]-afWin[4])+
                   fabs(afWin[7]-afWin[4])+
                   fabs(afWin[8]-afWin[4]))/ 8;
}

2、地形位置指数（Topographic Position Index，TRI）

float GDALTPIAlg (float* afWin, float fDstNoDataValue,void* pData)
{
         returnafWin[4] -
                   ((afWin[0]+ afWin[1] + afWin[2] + afWin[3] +
                   afWin[5]+ afWin[6] + afWin[7] + afWin[8]) / 8);
}

3、粗糙度（roughness）

float GDALRoughnessAlg (float* afWin, floatfDstNoDataValue, void* pData)
{
         floatpafRoughnessMin = afWin[0];
         floatpafRoughnessMax = afWin[0];
 
         for ( int k =1; k < 9; k++)
         {
                   if(afWin[k] > pafRoughnessMax)
                            pafRoughnessMax= afWin[k];
 
                   if(afWin[k] < pafRoughnessMin)
                            pafRoughnessMin= afWin[k];
         }
 
         returnpafRoughnessMax - pafRoughnessMin;
}

三、 处理截图 使用的数据依旧是SRTM的DEM数据，我裁切了一块用于快速计算。图3是原始的DEM显示结果，图4、图5、图6分别是TRI、TPI和粗糙度的显示截图。

图3 DEM数据


图4 TRI

图5 TPI


图6 粗糙度 四、 参考文献[1] http://download.osgeo.org/qgis/doc/reference-docs/Terrain_Ruggedness_Index.pdf
[2] http://www.gdal.org/gdaldem.html