三维场景下坡度、坡向图生成
在得到区域较大比例尺DEM数据后，通过各格网单元的领域高程差等运算，可生成此区域的三维场景下的坡度与坡向图。其用途较广，如坡度图一般用于退耕还林、土地适宜性评价等方面，坡向图为辐照度等进一步的计算提供基础。另外，坡度、坡向图也应用于自然条件下的水流分布、滑坡等地质灾害监测等方面。以下给出基本公式与程序代码参考。
此图引自：Tarboton, D.G. (1997): A new method for the determination of flow directions and upslope areas in grid digital elevation models, Water Resources Research, Vol.33, No.2, p.309-319
参考上图，将这个简单的三角形看成一坡度，则将坡度用矢量（S1,S2）的形式表示，有
S1 = (e0– e1)/d1；
S2 = (e1 – e2)/d2；
其中，e0，e1，e2代表格网单元上的高程数值，d1，d2代表单元之间的距离（有时就是1），
则此坡度的方向r与大小（模）s是：
r = tan-1(S2 / S1)；
s = sqrt(S1*S1+ S2*S2)；
参考代码如下：
作为第一步的入口函数：
bool SlopeCreate::SlopeStart(void) //以DEM为基本数据计算地形坡度、坡向
{
long gridNum = m_iRows*m_iColumns;
double* pSurfaceSlope = new double[gridNum];
long lIndex = 0;
double dSlopeMin = 3.14; //90角度等于.5707963弧度45角度等于.78弧度
double dSlopeMax = -3.14;
intx, y;
for(y=0; y<m_iRows; y++) //两重循环，取DEM上的每个元素 {
for(x=0; x<m_iColumns; x++)
{
double slope = ComputeSlope(x, y); //对格网上每个元素进行计算
lIndex = y*m_iColumns+x;
pSurfaceSlope[lIndex] = slope;
if(dSlopeMin>slope) dSlopeMin = slope;//得到较小的坡度数值, 记录下来
if(dSlopeMax<slope) dSlopeMax = slope;//得到较大的坡度数值，记录下来
}
}
//设置颜色表
m_byRsAnalysis = new BYTE[m_iRows * m_iColumns]; //
m_byGsAnalysis = new BYTE[m_iRows * m_iColumns];
m_byBsAnalysis = new BYTE[m_iRows * m_iColumns];
for(lIndex=0;lIndex<gridNum;lIndex++)
{
OutputSlopeColor(lIndex,pSurfaceSlope[lIndex]); //对格网上的每个单元格赋颜色
}
return true;
}
计算坡度、坡向的核心函数：
double SlopeCreate::ComputeSlope(int x, int y)
{
int i, ix, iy, j;
double z, zm[8], iSlope, iAspect, Slope, Aspect, G, H;
if( m_p3dDEMGrid->IsNodata(x,y))
{
long tIndex = y * m_iColumns + x;
m_byRsAnalysis[tIndex] = 200; //没有数值的地方填充灰白
m_byGsAnalysis[tIndex] = 200;
m_byBsAnalysis[tIndex] = 200;
return 0;
}
//-----------------------------------------------------
z = m_p3dDEMPoints[y*m_iColumns+x][2]; //得到本位置处的高程数值
for(i=0; i<8; i++) //对本位置的相邻个元素进行循环
{
ix = m_p3dDEMGrid->GetxTo(i, x); //通过方向判断得到邻近XY数值
iy = m_p3dDEMGrid->GetyTo(i, y);
if( m_p3dDEMGrid->IsInGrid(int(ix),int(iy))) //判断是否落在格网范围内
{
zm[i] = m_p3dDEMPoints[iy*m_iColumns+ix][2]; //记录下这个邻域单元的高程数值
}
else //如果没有落在格网范围之内
{
ix = m_p3dDEMGrid->GetxTo(i + 4, x);//得到邻近单元的索引
iy = m_p3dDEMGrid->GetyTo(i + 4, y);//得到邻近单元的索引
if( m_p3dDEMGrid->IsInGrid(int(ix),int(iy)) )
{
zm[i] = z - (m_p3dDEMPoints[iy*m_iColumns+ix][2] - z);
}
else
{
zm[i] = z;
}
}
}
//---------------------------------------------
Slope = 0.0; //初始数值
Aspect = -1.0;
for(i=0, j=1; i<8; i++, j=(j+1)%8) //相邻元素的循环
{
if( i % 2 ) {
G = (z - zm[j]) / m_iXCellSize;
H = (zm[j]- zm[i]) / m_iXCellSize;
}
else {
G = (z - zm[i]) / m_iXCellSize;
H = (zm[i]- zm[j]) / m_iXCellSize;
}
//参考上面的公式分析
if( H < 0.0 )
{
iAspect = 0.0;
iSlope = G;
}
else if( H > G )
{
iAspect = M_PI_045;
iSlope = (z - zm[i % 2 ? i : j]) / (sqrt(2.0) * m_iXCellSize);
}
else
{
iAspect = atan(H / G); //求反函数
iSlope = sqrt(G*G + H*H);
}
if( iSlope > Slope )
{
Aspect = i * M_PI_045 + (i % 2 ? M_PI_045 - iAspect : iAspect);
Slope = iSlope;
}
}
//---------------------------------------------
if( Aspect < 0.0 )
{
return 0;
}
else
{
return atan(Slope); //返回坡度数值
}
}
寻找邻近单元函数
int CGV3dDEMGrid::GetxTo(int Direction, int x)
{
static int ix[8] = { 0, 1, 1, 1, 0,-1,-1,-1 };
return( x + ix[Direction % 8] ); //得到邻近单元的索引
}
intCGV3dDEMGrid::GetyTo(int Direction, int y)
{
static int iy[8] = { 1, 1, 0,-1,-1,-1, 0, 1 };
return( y + iy[Direction % 8] ); //得到邻近单元的索引
}
原始图象(配图例)
三维场景下坡度图(大)
三维场景下坡向图(大)