MapX应用讲义第06章 绘等值线

1、第6章 绘等值线1576.1 分析数据是否为网格化数据1576.2 离散数据三角化1586.3 数据网格化1666.3.1 网格化方法1666.3.2 网格化程序1696.4 网格化数据的坡度、坡向分析1706.4.1 方法简介1706.4.2 坡度、坡向分析程序1716.5 用三角形法绘平面等值线1756.5.1 方法简介1756.5.2 三角法绘平面等值线子程序1766.6 用网格法绘平面等值线1826.6.1 方法简介1826.6.2 网格法绘平面等值线子程序1826.7 平面等值线填色1836.7.1 方法说明1836.7.2 绘平面等值线程序1836.8 用网格法绘立体等值线1926

2、.8.1 方法简介1926.8.2 网格法绘立体等值线程序1926.9 用网格法绘立体表面图1956.9.1 方法说明1956.9.2 用网格法绘立体表面图程序1956.10 绘等值线剖面图1986.10.1 方法简介1986.10.2 绘等值线剖面图程序1986.10.3 运行程序2026.11 用浮动水平线算法绘网格表面消隐图形2046.11.1 方法简介2046.11.2 浮动水平线算法绘网格表面程序2046.12 一个完整的绘等值线及三维表面图程序205第6章 绘等值线本章主要讲述如何用三角法绘平面等值线，如何用网格法绘平面等值线、立体等值线、立体表面图，以及如何进行等值线填色等，并给

3、出全部源程序。MapInfo系统无绘等值线及立体表面图的程序，本章给出全部源程序，用户将其稍作修改后便可用于自己的MapInfo应用程序系统。6.1 分析数据是否为网格化数据用户给出的数据有可能是已经网格化好了的数据，也有可能是任意散乱点数据，对这些数据必须加以区分，以便进行下一步处理。已经网格化好了的数据就可以直接用于绘平面图及立体图，而对于任意散乱点数据，除可直接用于绘平面图外，必须经过网格化后才能用于绘立体图。下面给出的就是用于分析数据是否为网格化数据的程序。源程序见光盘61CheckGridCheckGrid.vbp。1.子程序语句Private Sub CheckGrid(X() A

4、s Single, Y() As Single, N As Integer, NX As Integer, NY As Integer)2.参数说明l X()：存放用户数据的X坐标，实型一维数组，长度为N。l Y()：，存放用户数据的Y坐标，实型一维数组，长度为N。l N：数组X、Y的长度, 整型。l NX： X方向的网格数, 整型。l NY： Y方向的网格数，整型。3.子程序'判断是否为网格数据Private Sub CheckGrid(X() As Single, Y() As Single, N As Long, NX As Integer, NY As Integer)Dim

5、IX As Long, IY As LongDim RR() As Single, VV() As Single, V As Single, C As SingleDim Xtemp As Single, Ytemp As Single, Error1 As IntegerDim I As LongDataType = 0'判断是否是网格数据Xtemp = X(0)Ytemp = Y(0)NX = 0NY = 0DY = Y(2) - Y(1)For I = 0 To N If (X(I) <> Xtemp) Then Exit For NY = NY + 1Next II

6、f (NY < 3) Then 'X方向数据小于3个,肯定是非网格数据 DataType = 1 GoTo Error1End IfNX = Fix(N + 1) / NY)If (NY * NX = N + 1) Then '有可能是网格数据 DataType = 0 '判断X坐标是否等间距 DX = X(NY) - X(0) For IX = 2 To NX - 1 Xtemp = X(IX * NY) - X(IX - 1) * NY) If (Abs(Xtemp - DX) > 0.000001) Then DataType = 1 GoTo Err

7、or1 End If Next IX '判断Y是否等间距 DY = Y(1) - Y(0) For IY = 2 To NY - 1 Ytemp = Y(IY) - Y(IY - 1) If (Abs(Ytemp - DY) > 0.000001) Then DataType = 1 GoTo Error1 End If Next IYElse '肯定不是网格数据 DataType = 1End IfError1:End Sub6.2 离散数据三角化1.方法说明平面散乱点集的三角化是指将平面上的N个点用互不相交的直线段连接起来，构成一三角形网格。一般而言，平面上N个点集连

8、成三角形网格的方式有多种，不同的三角化准则对应有不同的三角化方式。1934年，俄国数学家Delaunay提出了三角形最小内角最大的三角化准则，并证明了在没有四点或四点以上共圆条件下的平面散乱点存在的一种三角化方式，使连成的三角形网中的三角形满足这一条件，最接近等边三角形，通常称这类三角形为Delaunay三角形。平面散乱点集的三角化的基本思想是：先构造一外接圆不包含第四点的初始三角形，然后将此形成的三角形的每条边作为基边，以基边所对应角度最大为准则，在其邻边区域搜索第三点，构成新的三角形，如此循环进行至所有点集都被三角化。经证明该算法所产生的三角形仍具有Delaunay性质。平面点集三角化的主

9、要步骤如下：(1)将平面散乱点按X、Y坐标从小到大分类排序，它将有助于消除重合点、快速搜索候选点，提高运行效率。(2)从点集中任取一点赋值给Pi作为初始边界的起点，并在Pi(Xi,Yi)的邻近区域寻找Pi最近的点Pj，即在Xi-j,Xi+jYi-j,Yi+j矩形区域搜索(j=1,2,)。(3)将PiPj赋给边界环表构成初始边界环。(4)从边界环表中弹出一有向边界PiPj作为基边。(5)在基边PiPj正侧的邻边区域搜索构成三角形的第三点Pk的候选点，即在Min(Xi-e,Xj+1),Max(Xi-e,Xj+1)、Min(Yi-e,Yj+1),Max(Yi-e,Yj+1)区域搜索（e=1,2,）。

10、(6)从搜索到的Pk候选点中选择使PiPkPj最大则者为Pk。(7)检查Pk点是否为边界环上的点，若Pk为边界环上的点转8，否则转9。(8)检查Pk是否与基边PiPj相邻，若不相邻则置当前基边的下一条边界边为新的基边，转5，若相邻则转9。(9)记录三角形PiPjPk，并修改边界环表。l 若Pk为基边PiPj的右邻，则边界边PiPj用PiPk置换。l 若Pk为基边PiPj的左邻，则边界边PiPj用PkPj置换。l 若Pk不为边界环上的点，则边界PiPj用PiPk.。同时在其后插入一新的边界PkPj。(10)判断点集中的点是否为空，若为空，则三角化结束，否则转5。2.子程序语句Private Su

11、b SJX()3.参数说明无。4.子程序源程序见光盘62SJXSJX.vbp，用VB运行出现如图6-1所示的界面。图6-1 三角化主要的源程序如下：Private Sub SJX()Dim R As Double, RT As DoubleDim J As Integer, K As IntegerDim A1 As Double, B1 As Double, C1 As Double, Tt As DoubleDim AA As Double, BB As Double, CC As Double, CosC As DoubleDim bCheck As Boolean, bEQBD As

12、BooleanDim Ni As Integer, Nj As Integer, Pi As Integer, Pj As Integer, Pk As IntegerDim M1 As Integer, M2 As IntegerDim X1 As Double, Y1 As Double, X2 As Double, Y2 As Double, X3 As Double, Y3 As DoubleDim bSJX As Boolean'边界段数Dim nBD As Integer'用于判断边界环是否搜索过Dim bBD() As Boolean'用于判断数据点是否在

13、边界环内Dim bPoint() As Byte'第i边界起点BD,第i边对应顶点Dim BD() As Integer, BDij() As Integer'第i边界的上一节点、下一节点Dim Nlast As Integer, Nnext As IntegerReDim ID1(0 To 2 * NContou), ID2(0 To 2 * NContou), ID3(0 To 2 * NContou)'Begin生成第一个三角形Pi = 0'找出距第一点最近的点2RT = 1E+20For J = 1 To NContou R = (Xcontou(J)

14、- Xcontou(Pi) 2 + (Ycontou(J) - Ycontou(Pi) 2 If (R < RT) Then RT = R Pj = J End IfNext J'找出第三点Tt = 0For J = 1 To NContou If (J <> Pj) Then AA = (Xcontou(Pj) - Xcontou(J) 2 + (Ycontou(Pj) - Ycontou(J) 2 BB = (Xcontou(Pi) - Xcontou(J) 2 + (Ycontou(Pi) - Ycontou(J) 2 CC = (Xcontou(Pi) - X

15、contou(Pj) 2 + (Ycontou(Pi) - Ycontou(Pj) 2 CosC = 1# - (AA + BB - CC) / (2# * Sqr(AA * BB) If (CosC > Tt + 0.00001) Then Tt = CosC Pk = J End If End IfNext JnSJX = 0ID1(nSJX) = PiID2(nSJX) = PjID3(nSJX) = Pk'End生成第一个三角形X1 = Xcontou(Pi)X2 = Xcontou(Pj)X3 = Xcontou(Pk) Y1 = Ycontou(Pi)Y2 = Yco

16、ntou(Pj)Y3 = Ycontou(Pk)PictureGrid.DrawMode = 13PictureGrid.ForeColor = QBColor(12)PictureGrid.Line (X1, Y1)-(X2, Y2)PictureGrid.Line (X2, Y2)-(X3, Y3)PictureGrid.Line (X3, Y3)-(X1, Y1)DoEvents'定义边界环数组ReDim BD(0 To NContou), BDij(0 To NContou), bBD(0 To NContou), bPoint(0 To NContou)For J = 0 T

17、o NContou bBD(J) = False bPoint(J) = 0Next J'生成三条边nBD = 3BD(1) = PiBDij(1) = PkbBD(1) = TrueBD(2) = PjBDij(2) = PibBD(2) = TrueBD(3) = PkBDij(3) = PjbBD(3) = TruebPoint(Pi) = 1bPoint(Pj) = 1bPoint(Pk) = 1Do Ni = Ni + 1 bEQBD = False For J = Ni To nBD If (bBD(J) = True) Then bEQBD = True Exit For

18、 End If Next J If (bEQBD = False) Then For J = 1 To Ni - 1 If (bBD(J) = True) Then bEQBD = True Exit For End If Next J End If Ni = J If (bEQBD = False) Then Exit Do '第Ni边的终点节点 If (Ni = nBD) Then Nj = 1 Else Nj = Ni + 1 End If '第Ni边的后邻 If (Ni = 1) Then Nlast = nBD Else Nlast = Ni - 1 End If &

19、#39;第Ni边的前邻 If (Nj = nBD) Then Nnext = 1 Else Nnext = Nj + 1 End If '第Ni边的三个顶点 Pi = BD(Ni) Pj = BD(Nj) Pk = BDij(Ni) 'Begin找下一点 bCheck = False If (Xcontou(Pj) = Xcontou(Pi) Then C1 = Xcontou(Pk) - Xcontou(Pi) Else A1 = (Ycontou(Pj) - Ycontou(Pi) / (Xcontou(Pj) - Xcontou(Pi) B1 = Ycontou(Pj)

20、- A1 * Xcontou(Pj) C1 = Ycontou(Pk) - Xcontou(Pk) * A1 - B1 End If If (C1 <= 0#) Then M1 = -1 Else M1 = 1 End If Tt = 0# For K = 1 To NContou If (K = Pi Or K = Pj Or bPoint(K) = 2) Then Else If (Xcontou(Pj) = Xcontou(Pi) Then C1 = Xcontou(K) - Xcontou(Pi) Else C1 = Ycontou(K) - Xcontou(K) * A1 -

21、B1 End If If (C1 <= 0#) Then M2 = -1 Else M2 = 1 End If If (M1 <> M2) Then AA = (Xcontou(Pj) - Xcontou(K) 2 + (Ycontou(Pj) - Ycontou(K) 2 BB = (Xcontou(Pi) - Xcontou(K) 2 + (Ycontou(Pi) - Ycontou(K) 2 CC = (Xcontou(Pi) - Xcontou(Pj) 2 + (Ycontou(Pi) - Ycontou(Pj) 2 CosC = 1# - (AA + BB - CC

22、) / (2# * Sqr(AA * BB) If (CosC > Tt + 0.00001) Then Tt = CosC Pk = K bCheck = True End If End If End If Next K 'End找下一点 If (bCheck = True) Then '找到下一点 bSJX = True '找前邻、后邻 If (Pk = BD(Nlast) Then '后邻,删除Pi点 '修改后邻点参数 BDij(Nlast) = Pi bBD(Nlast) = True bPoint(Pk) = 1 For J = Ni +

23、 1 To nBD BD(J - 1) = BD(J) BDij(J - 1) = BDij(J) bBD(J - 1) = bBD(J) Next J nBD = nBD - 1 '减少一个边 bPoint(Pi) = 2 'Pi点位于圈内 ElseIf (Pk = BD(Nnext) Then '前邻,删除Pj点 '修改Ni点参数 BDij(Ni) = Pj bBD(Ni) = True bPoint(Pk) = 1 For J = Nj + 1 To nBD BD(J - 1) = BD(J) BDij(J - 1) = BDij(J) bBD(J - 1

24、) = bBD(J) Next J nBD = nBD - 1 '减少一个边 bPoint(Pj) = 2 'Pj点位于圈内 ElseIf (bPoint(Pk) < 1) Then '下一点不在边界环上,插入Pk点 '修改Ni点参数 BDij(Ni) = Pj bBD(Ni) = True bPoint(Pk) = 1 '数据后推 For J = nBD To Ni + 1 Step -1 BD(J + 1) = BD(J) BDij(J + 1) = BDij(J) bBD(J + 1) = bBD(J) Next J nBD = nBD +

25、1 'Pk点参数 BD(Ni + 1) = Pk BDij(Ni + 1) = Pi bBD(Ni + 1) = True Ni = Ni + 1 Else bSJX = False End If If (bSJX = True) Then '生成新三角形 nSJX = nSJX + 1 ID1(nSJX) = Pi ID2(nSJX) = Pj ID3(nSJX) = Pk X1 = Xcontou(Pi) X2 = Xcontou(Pj) X3 = Xcontou(Pk) Y1 = Ycontou(Pi) Y2 = Ycontou(Pj) Y3 = Ycontou(Pk)

26、PictureGrid.Line (X1, Y1)-(X2, Y2) PictureGrid.Line (X2, Y2)-(X3, Y3) PictureGrid.Line (X3, Y3)-(X1, Y1) DoEvents End If Else '没找到下一点,该边界不再搜索 bBD(Ni) = False End IfLoop'整理边界nBorder = nBD - 1ReDim BorderX(0 To nBorder + 1), BorderY(0 To nBorder + 1)For J = 0 To nBorder BorderX(J) = Xcontou(BD

27、(J + 1) BorderY(J) = Ycontou(BD(J + 1)Next JBorderX(nBorder + 1) = BorderX(0)BorderY(nBorder + 1) = BorderY(0)End Sub6.3 数据网格化离散数据网格化的方法有很多，本节提供三角面平面插值、双二次三角曲面光滑插值、双三次三角曲面光滑插值、按距离平方反比加权插值、按方位取点加权插值、趋势面拟合及加权最小二乘曲面拟合等七种方法。用户可根据自己的需要选用不同的方法。6.3.1 网格化方法下面简述各网格化方法的基本原理。1三角平面插值该方法是基于对三维散列数据的定义域做三角剖分，然后在网格

28、点所在的三角形上做线性插值。2三角曲面光滑插值该方法是基于对三维散列数据的定义域做三角剖分，然后在每个三角形上做双三次分片插值。三维散列数据的光滑插值的基本思想是：已知一组3D数据(xi,yi,zi)，i=1,2,N，构造一个具有所需光滑度的曲面z=S(x,y),使之插值于这些数据，即满足S(xi,yi)=zi，i=1,2, ,N，对数据点(xi,yi)要求互不重合且不共线。下面我们导出分片三角曲面的Lagrange插值公式。设插值曲面S定义在xy平面上区域D中，(xi,yi)D,i=1,2,N。将Vi=(xi,yi)作为顶点，对区域D作三角剖分，在每个三角形上求插值函数S(x,y)。我们希望

29、插值曲面S具有C1连续性，即它的切平面连续变化。对xy平面的任意三角形T，顶点为V1、V2、V3，V是平面上任一点，可以定义V相对于T的重心坐标(u1,u2,u3),其中分量u1=A1/A，A1和A分别是VV2V3和T的面积，u1、u2类推。于是又因为一个双三次多项式有10个未知系数，所以需要10个插值条件才能确定。若在T上已知10个标记点处的数据，即z=S(vi)，i=1,10，则可以确定插值曲面是S(x,y)，为其中属于双三次多项式函数类，且满足，据此可以导出用重心坐标表示的双三次Lagrange插值函数：可以证明上式具有仿射不变性，因此适合用于三维图形的表示，且具有C1连续性。类似地，对

30、双二次插值曲面（有6个未知系数），若已知T上数据点zi=s(pi)，i=1,6，则有插值曲面的中心坐标式：其中，Q1=u1(2u1-1)，Q2=u2(2u2-1)，Q3=u3(2u3-1)，Q4=4u2u3，Q5=4u1u3，Q6=4u1u2，双二次Lagrange插值在分片三角曲面上是C1连续的，但在分片曲面连接处可能有切向跳跃，达不到全局C1连续。3按距离平方反比加权插值该方法是基于对三维散列数据的定义域作三角剖分，然后在网格点所在的三角形上，把距离平方反比作为权系数加权平均求网格点的值。三角形顶点坐标为（Xi,Yi）,i=1,2,3，网格点坐标为（X0,Y0），且网格点在该三角形内，则该

31、网格点的值为：4按方位取点加权插值该方法的基本原理是：欲求某个网格点(i,j)的函数时，则以(i,j)为原点将平面分成四个基本象限，再把每个象限分成n0份。这样就把全平面分成4n0等份。然后在每个等分角内寻找一个距(i,j)最近的数据点，其值为Zi1，它到(i,j)的距离为ri1。则网格(i,j)上的值为：式中参数由上式可以看出，4n0和Ci1之和即：因此，Ci1是符合权系数定义的。再看当ri1=0时，即Zi1就在网格(i,j)上时，由于中时，因此有且其他的Cj1(ji1)都为零，因此网格(i,j)上的值就是数据点的值Zi1。5趋势面拟合趋势面拟合可以把长周期的趋势变化和短周期的局部变化分开。

32、趋势面分析是用简单的幂级多项式来拟合复杂的曲面，有削平、填平实际曲面的作用。为了准确地反映实际曲面，不能单独采用这种方法，趋势面分析与其他方法一起使用可以得到较准确的拟合。令趋势面次数为n0,则趋势分析的公式为:n0趋势面的系数ak,i的个数为：若有Nl个测量值，现在要做n0次趋势面(应有NlS),需要对Nl个观测值Z1、Z2、ZNl用最小二乘法确定系数aK,I。先作观测值与理论差的平方和：再求对系数aK,i的偏导数，并令其为零：这样线性代数方程组有S个方程和S个未知数ak,I，求出aK,i后就确定了趋势面表达式。6加权最小二乘曲面拟合假定我们要计算网格点(a,b)上的曲面高度f(a,b)，为

33、此需要求出一个多项式p(x,y)，一般二次多项式对一般趋势分析应有为最小，式中n是数据点个数，(xi,yi)是数据点坐标，Zi是（xi,yi）上的观测值，拟合值 p(xi,yi)和观测值 zi的误差平方和为Q。考虑按距离加权，应对上式做些改动。这里仅以二次多项式p(x,y)为例：式中就是权，它是距离的函数。为了求出方程中的系数Crs,按最小二乘原理应有：(r,s=0,1,2)上式是有6个未知数和6个方程的方程组，可以用高斯消去法求解。程序设计时可用如下的简单方法。设E为6×6的对称矩阵，V是长度为6的列向量。对Crs行和Cra列对应的E的元素，我们加上对应于Crs的V的元素加上对所有

34、的数据点都做过后，再解线性代数方程组EC=V，就可求出系数C。加权方式可以是以下三种：l 型式中d为网点(a,b)到离散数据点（xi,yi）的距离，即d2=(xi-a)2+(yi-b)2。是一个很小的数，为了防止运算时算术溢出。l 型式中的d与的意义同型。l 型式中的是常数，其数量级相当于相邻数据点距离平方的倒数。d与的意义同型。6.3.2 网格化程序源程序见光盘63GridGrid.vbp，用VB运行出现如图6-2所示的界面。图6-2 网格化程序l 数据网格化主程序是：Private Sub Grid(iGrid As Integer, PictureGrid As Object)，该主程序

35、中含“三角平面插值网格化”、“距离平方反比插值网格化”及“加权最小二乘曲面拟合网格化”程序。l 双二次Lagrange、双三次Lagrange插值网格化程序是：Private Sub Lagrange(X1 As Double, Y1 As Double, Z1 As Double, X2 As Double, Y2 As Double, Z2 As Double, X3 As Double, Y3 As Double, Z3 As Double, XX As Double, YY As Double, ZZ As Double, b23 As Integer)l 按方位取点加权网格化程序是：

36、Private Sub GridAZimuth(NY As Integer, NX As Integer, N0 As Integer, Xmin As Double, Ymin As Double, DX As Double, DY As Double, Xcontou() As Double, Ycontou() As Double, Zcontou() As Double, NContou As Long, Zgrid() As Double, bZgrid() As Boolean, Dr0 As Double)l 趋势面拟合网格化程序是：Private Sub GridTrendSu

37、rface(NY As Integer, NX As Integer, N0 As Integer, Xmin As Double, Ymin As Double, DX As Double, DY As Double, Xcontou() As Double, Ycontou() As Double, Zcontou() As Double, NContou As Long, Zgrid() As Double, bZgrid() As Boolean, Dr0 As Double)6.4 网格化数据的坡度、坡向分析6.4.1 方法简介根据空间解析几何的原理，基本单位长度为X、Y，原点坐标为

38、(x0,y0)的格网模型，任一格网交点的标准向量为，即对于由四个相邻格网点确定的地表基本单元，其基本向量的计算公式如下通过基本向量就可以确定基本单元的空间特性。向量的向量积就是基本单元的法向量。利用，就可以进行地表单元坡度、坡向等的分析和计算。格网模型中地表基本单元的坡度等于其法向量与Z轴之夹角，而向量的夹角的余弦等于两向量的数量积与模的乘积之商，即地表基本单元的坡向为其法向量在xoy平面上的投影与x轴的夹角。计算公式为其中6.4.2 坡度、坡向分析程序1.功能计算坡度、坡向。2.语句Private Sub Slope(ZGrid() As Double, zGridSlope() As Do

39、uble, nRows As Integer, nCols As Integer, dx As Double, dy As Double, Vmin As Double, Vmax As Double)Private Sub Aspect(ZGrid() As Double, ZGridDir() As Double, nRows As Integer, nCols As Integer, dx As Double, dy As Double, Vmin As Double, Vmax As Double)3.参数说明l ZGrid()：输入网格点值，双精度，。l zGridSlope()：输

40、出网格点坡度，双精度。l ZGridDir()：输出网格点坡向，双精度。l nRows：输入参数，纵向网格数。l nCols：输入参数，横向网格数。l dx、dy：输入参数，横向、纵向网格步长。l Vmin、Vmax：输出参数，坡度、坡向最小值和最大值。.子程序源程序参见书附光盘64SlopeSlope.vbp，用VB运行出现如图6-3所示界面。图6- 坡度、坡向分析在图6-3中，选择网格化等值线数据，单击【开始计算】，在输入数据目录下生成坡度、坡向结果文件。下面为计算坡度、坡向的源程序。Private Sub Slope(ZGrid() As Double, zGridSlope() As

41、Double, nRows As Integer, nCols As Integer, dx As Double, dy As Double, Vmin As Double, Vmax As Double) Dim I As Integer, J As Integer Dim DYZx As Double, DXZy As Double Dim Zx As Double, Zy As Double Dim X As Double ReDim zGridSlope(1 To nCols, 1 To nRows) '计算非边界点 For J = 2 To nRows - 1 For I =

42、 2 To nCols - 1 Zx = ZGrid(I + 1, J - 1) - ZGrid(I - 1, J - 1) + ZGrid(I + 1, J) - ZGrid(I - 1, J) + ZGrid(I + 1, J + 1) - ZGrid(I - 1, J + 1) Zy = ZGrid(I - 1, J + 1) - ZGrid(I - 1, J - 1) + ZGrid(I, J + 1) - ZGrid(I, J - 1) + ZGrid(I + 1, J + 1) - ZGrid(I + 1, J - 1) DYZx = dy * Zx DXZy = dx * Zy

43、X = 2# * dx * dy / Sqr(DYZx 2 + DXZy 2 + (2# * dx * dy) 2) zGridSlope(I, J) = ArcCos(X) Next I Next J '边界点赋值 For J = 1 To nRows zGridSlope(1, J) = zGridSlope(2, J) zGridSlope(nCols, J) = zGridSlope(nCols - 1, J) Next J For I = 1 To nCols zGridSlope(I, 1) = zGridSlope(I, 2) zGridSlope(I, nRows) =

44、 zGridSlope(I, nRows - 1) Next I Vmin = zGridSlope(1, 1) Vmax = zGridSlope(1, 1) For J = 1 To nRows For I = 1 To nCols If (zGridSlope(I, J) > Vmax) Then Vmax = zGridSlope(I, J) If (zGridSlope(I, J) < Vmin) Then Vmin = zGridSlope(I, J) Next I Next JEnd Sub计算坡向Private Sub Aspect(ZGrid() As Doubl

45、e, ZGridDir() As Double, nRows As Integer, nCols As Integer, dx As Double, dy As Double, Vmin As Double, Vmax As Double) Dim I As Integer, J As Integer Dim DYZx As Double, DXZy As Double Dim Zx As Double, Zy As Double ReDim ZGridDir(1 To nCols, 1 To nRows) '计算非边界点 For J = 2 To nRows - 1 For I =

46、2 To nCols - 1 Zx = ZGrid(I + 1, J - 1) - ZGrid(I - 1, J - 1) + ZGrid(I + 1, J) - ZGrid(I - 1, J) + ZGrid(I + 1, J + 1) - ZGrid(I - 1, J + 1) Zy = ZGrid(I - 1, J + 1) - ZGrid(I - 1, J - 1) + ZGrid(I, J + 1) - ZGrid(I, J - 1) + ZGrid(I + 1, J + 1) - ZGrid(I + 1, J - 1) DYZx = dy * Zx DXZy = dx * Zy I

47、f (Abs(DXZy) < 0.000000001 And Abs(DYZx) < 0.000000001) Then ZGridDir(I, J) = 0 ElseIf (Abs(DXZy) < 0.000000001) Then If (DYZx > 0) Then ZGridDir(I, J) = 270 Else ZGridDir(I, J) = 90 End If Else If (Abs(DYZx) < 0.000000001) Then If (DXZy > 0) Then ZGridDir(I, J) = 180 Else ZGridDir

48、(I, J) = 0 End If ElseIf (DYZx > 0 And DXZy > 0) Then ZGridDir(I, J) = 180 + Atn(DYZx / DXZy) * 180 / 3.14159265 ElseIf (DYZx < 0 And DXZy < 0) Then ZGridDir(I, J) = Atn(DYZx / DXZy) * 180 / 3.14159265 ElseIf (DYZx < 0 And DXZy > 0) Then ZGridDir(I, J) = 180 + Atn(DYZx / DXZy) * 180 / 3.14159265 ElseIf (DYZx > 0 And DXZy < 0) Then ZGridDir(I, J) = 360 + Atn(DYZ