D8算法的matlab实现 身的编程实践介绍D8算法的matlab实现过程，与大家分享。 D8算法是当今非常成熟的提取水系（或分水岭）的计算机程序算法。这里笔者结合自 D8算法分为两部分实现：每个栅格点水流方向的计算；每个栅格点汇水面积的计算。 1、水流方向的计算 在水流方向的计算中，本人结合matlab内嵌的滤波工具，利用一个简单的非线性滤波 来实现，具体内容包括一个自己编辑的确定水流方向的函数flowdirection和若干脚本命 令，如下： 以下是自己编辑的flowdirection函数代码： 1 2 3 4 5 6 7 8 9 function output=flowdirection(a) n=size(a); for i=1:n(2) k=0; b=-inf; for j=[1:4,6:9] if rem(j,2)==0 r(j)=a(5,i)-a(j,i); r(j)=(a(5,i)-a(j,i))/sqrt(2); else end if r(j)>b b=r(j); k=j; end end output(i)=k; end 执行的脚本代码如下： >> uiopen('F:\课件\2009田淑芳\正射校正\DEM.tif',1)%打开DEM高程数据，生成的数据矩阵 名可自行定义，这里默认为DEM >>DEM=mat2gray(DEM);%高程数据归一化处理 >>DEM=padarray(DEM,[1 1],'replicate');%对数据边缘进行自动填充 >>direc8=colfilt(DEM,[3 3],'sliding',@flowdirection);%执行非线性滤波，得到方向矩阵，得到 的值分别为{1 ,2,3,4,6,7,8,9}.其方向意义如上图矩阵所示，如1代表北西向，8代表正东 1、汇水面积的计算 这里笔者根据自己的编程实践经验，给出一个函数来计算不同窗口大小下汇水面积的大 小。函数代码如下： function f=flowareax(a,x) %参数a为水流方向矩阵，x为表征窗口大小的参数， b=a; [m n]=size(b); 

f=ones([m+2 n+2 x]);%默认每点有一初始水量 for i=1:m for j=1:n if a(i,j)==1 f(i,j,1)=f(i,j,1)+1; elseif a(i,j)==2 f(i+1,j,1)=f(i+1,j,1)+1; elseif a(i,j)==3 f(i+2,j,1)=f(i+2,j,1)+1; elseif a(i,j)==4 f(i,j+1,1)=f(i,j+1,1)+1; elseif a(i,j)==6 f(i+2,j+1,1)=f(i+2,j+1,1)+1; elseif a(i,j)==7 f(i,j+2,1)=f(i,j+2,1)+1; elseif a(i,j)==8 f(i+1,j+2,1)=f(i+1,j+2,1)+1; f(i+2,j+2,1)=f(i+2,j+2,1)+1; else end end end%通过该循环确定经过一步流动时某一点水流量 for s=2:x for i=1:m for j=1:n if a(i,j)==1 f(i,j,s)=f(i,j,s)+f(i+1,j+1,s-1); f(i+1,j,s)=f(i+1,j,s)+f(i+1,j+1,s-1); f(i+2,j,s)=f(i+2,j,s)+f(i+1,j+1,s-1); f(i,j+1,s)=f(i,j+1,s)+f(i+1,j+1,s-1); elseif a(i,j)==2 elseif a(i,j)==3 elseif a(i,j)==4 elseif a(i,j)==6 elseif a(i,j)==7 elseif a(i,j)==8 f(i+2,j+1,s)=f(i+2,j+1,s)+f(i+1,j+1,s-1); f(i,j+2,s)=f(i,j+2,s)+f(i+1,j+1,s-1); f(i+1,j+2,s)=f(i+1,j+2,s)+f(i+1,j+1,s-1); f(i+2,j+2,s)=f(i+2,j+2,s)+f(i+1,j+1,s-1); else end end end 

end f=f(2:m+1,2:n+1,:); 该函数的执行结果为一个三维数据体，其中f（:,:,1）、f(:,:,2)、。。。。f(:,:,x)依次代表了不同窗 口下每个栅格点的汇水面积，可以根据对比他们的效果并且根据需要选择不同的数据层。 需要帮助 1045105061@qq.com