slope.frm

地理信息系统：坡度、坡向程序

    Close (1)
    Call SortXY(Xcontou, Ycontou, Zcontou, nContou)
    
    '判断是否是网格数据
    Call CheckGrid(Xcontou, Ycontou, nContou, NX, NY)
    If (DataType = 0) Then
        ReDim ZGrid(1 To NX, 1 To NY), bZgrid(1 To NX, 1 To NY)
        II = -1
        For IX = 1 To NX
            For IY = 1 To NY
                II = II + 1
                ZGrid(IX, IY) = Zcontou(II)
            Next IY
        Next IX
    End If

    Xmin = Xcontou(0)
    Xmax = Xcontou(0)
    Ymin = Ycontou(0)
    Ymax = Ycontou(0)
    Vmin = Zcontou(0)
    Vmax = Zcontou(0)
    For I = 0 To nContou
        If (Xcontou(I) < Xmin) Then Xmin = Xcontou(I)
        If (Xcontou(I) > Xmax) Then Xmax = Xcontou(I)
            
        If (Ycontou(I) < Ymin) Then Ymin = Ycontou(I)
        If (Ycontou(I) > Ymax) Then Ymax = Ycontou(I)
            
        If (Zcontou(I) < Vmin) Then Vmin = Zcontou(I)
        If (Zcontou(I) > Vmax) Then Vmax = Zcontou(I)
    Next I
    
    'Begin判断经纬度是否颠倒
    If (Ymin >= -180 And Ymax <= 180 And Xmin >= -90 And Xmax <= 90) Then
        ReDim Xcontou(0 To nContou), Ycontou(0 To nContou), Zcontou(0 To nContou)
        Open TheContouPath For Input As #1
        For I = 0 To nContou
            Input #1, Xcontou(I), Ycontou(I), Zcontou(I)
        Next I
        Close (1)
        Call SortXY(Xcontou, Ycontou, Zcontou, nContou)
    
        '判断是否是网格数据
        Call CheckGrid(Xcontou, Ycontou, nContou, NX, NY)
        If (DataType = 0) Then
            ReDim ZGrid(1 To NX, 1 To NY), bZgrid(1 To NX, 1 To NY)
            II = -1
            For IX = 1 To NX
                For IY = 1 To NY
                    II = II + 1
                    ZGrid(IX, IY) = Zcontou(II)
                Next IY
            Next IX
        End If
    End If
    'End判断经纬度是否颠倒
    If (DataType <> 0) Then '非网格化数据
        NX = 0
        NY = 0
    End If
    dx = (Xmax - Xmin) / (NX - 1)
    dy = (Ymax - Ymin) / (NY - 1)
End If
End Sub
'排序子程序
Private Sub ShellSort(X() As Double, Y() As Double, Z() As Double, N As Long, N1 As Long, N2 As Long)
Dim B As Long, M As Long, L As Long, I As Long, J As Long, K As Long
Dim Xtemp As Double, Ytemp As Double, Ztemp As Double, jk As Long

B = N2 - N1 + 1
B = Log(B) / Log(2#)
M = B
L = 2 ^ M
For I = 1 To M
    K = L - 1
    L = L / 2
    For J = K + N1 To N2
        Xtemp = X(J)
        Ytemp = Y(J)
        Ztemp = Z(J)
        jk = J - K
        Do While (jk > N1 - 1 And X(jk) > Xtemp)
            X(jk + K) = X(jk)
            Y(jk + K) = Y(jk)
            Z(jk + K) = Z(jk)
            jk = jk - K
            If (jk < 0) Then Exit Do
        Loop
        X(jk + K) = Xtemp
        Y(jk + K) = Ytemp
        Z(jk + K) = Ztemp
    Next J
Next I
End Sub
'判断是否为网格数据
Private Sub CheckGrid(X() As Double, Y() As Double, N As Long, NX As Integer, NY As Integer)
Dim IX As Long, IY As Long
Dim RR() As Double, VV() As Double, V As Double, C As Double
Dim Xtemp As Double, Ytemp As Double, Error1 As Integer
Dim I As Long, dy As Double, dx As Double

DataType = 0
'判断是否是网格数据
Xtemp = X(0)
Ytemp = Y(0)
NX = 0
NY = 0
dy = Y(2) - Y(1)
For I = 0 To N
    If (X(I) <> Xtemp) Then Exit For
    NY = NY + 1
Next I
If (NY < 3) Then 'X方向数据小于3个,肯定非网格数据
    DataType = 1
    GoTo Error1
End If
NX = 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 Error1
        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 IY
Else '肯定不是网格数据
    DataType = 1
End If
Error1:
End Sub



'数据按X、Y排序
Private Sub SortXY(X() As Double, Y() As Double, Z() As Double, N As Long)
Dim N1 As Long, N2 As Long, X1 As Double
Dim I As Long, J As Long

'按X坐标排序
Call ShellSort(X, Y, Z, N, 0, N)
'按Y坐标排序
X1 = X(0)
N1 = 0
N2 = 0
For I = 1 To N
    If (X(I) = X1) Then
        N2 = N2 + 1
    Else
        If (N2 > N1) Then '相同X按Y坐标排序
            Call ShellSort(Y, X, Z, N, N1, N2)
        End If
        X1 = X(I)
        N1 = I
        N2 = I
    End If
Next I
If (N2 > N1) Then '相同X按Y坐标排序
    Call ShellSort(Y, X, Z, N, N1, N2)
End If

'平均重合点
J = -1
X1 = Z(0)
N1 = 1
For I = 1 To N
    If (Abs(X(I) - X(I - 1)) + Abs(Y(I) - Y(I - 1)) < 0.00001) Then
        X1 = X1 + Z(I)
        N1 = N1 + 1
    Else
        J = J + 1
        Z(J) = X1 / N1
        X(J) = X(I - 1)
        Y(J) = Y(I - 1)
        X1 = Z(I)
        N1 = 1
    End If
Next I
J = J + 1
Z(J) = X1 / N1
X(J) = X(N)
Y(J) = Y(N)
N = J
End Sub

Private Sub CommandExit_Click()
    Unload Me
End Sub
Private Sub CommandOK_Click()
    Dim ZGrid() As Double, ZGridDir() As Double, zGridSlope() As Double
    Dim nRows As Integer, nCols As Integer
    Dim nRow As Integer, nCol As Integer
    Dim dx As Double, dy As Double
    Dim Xmin As Double, Xmax As Double, Ymin As Double, Ymax As Double, Vmin As Double, Vmax As Double
    Dim I As Integer
    Dim TheInFile As String, TheOutFile As String

    Screen.MousePointer = 11

    ProgressBar1.Visible = True
    ProgressBar1.Min = 0
    ProgressBar1.Max = ListIn.ListCount - 1

    For I = 0 To ListIn.ListCount - 1

        If (ListIn.Selected(I) = True) Then
            ProgressBar1.Value = I
            TheInFile = TheInPath + FileInT.List(I)
            
            Call ReadContouFile(TheInFile, dx, dy, nCols, nRows, ZGrid, Xmin, Xmax, Ymin, Ymax, Vmin, Vmax)
            If (nRows > 0 And nCols > 0) Then
                Call Slope(ZGrid, zGridSlope, nRows, nCols, dx, dy, Vmin, Vmax)
                '保存坡度
                TheOutFile = Left(TheInFile, Len(TheInFile) - 4) + "_坡度.GRD"
                Open TheOutFile For Output As #1
                Print #1, "DSAA"
                Print #1, Format(nCols, "###0 "); Format(nRows, "###0")
                Print #1, Format(Xmin, "###0.000 "); Format(Xmax, "###0.000")
                Print #1, Format(Ymin, "###0.000 "); Format(Ymax, "###0.000")
                Print #1, Format(Vmin, "###0.0 "); Format(Vmax, "###0.0")
                For nRow = 1 To nRows
                    For nCol = 1 To nCols - 1
                        Print #1, Format(zGridSlope(nCol, nRow), "####0.0 ");
                    Next nCol
                    Print #1, Format(zGridSlope(nCol, nRow), "####0.0")
                Next nRow
                Close (1)
                
                Call Aspect(ZGrid, ZGridDir, nRows, nCols, dx, dy, Vmin, Vmax)
                '保存坡向
                TheOutFile = Left(TheInFile, Len(TheInFile) - 4) + "_坡向.GRD"
                Open TheOutFile For Output As #1
                Print #1, "DSAA"
                Print #1, Format(nCols, "###0 "); Format(nRows, "###0")
                Print #1, Format(Xmin, "###0.000 "); Format(Xmax, "###0.000")
                Print #1, Format(Ymin, "###0.000 "); Format(Ymax, "###0.000")
                Print #1, Format(Vmin, "###0.0 "); Format(Vmax, "###0.0")
                For nRow = 1 To nRows
                    For nCol = 1 To nCols - 1
                        Print #1, Format(ZGridDir(nCol, nRow), "####0.0 ");
                    Next nCol
                    Print #1, Format(ZGridDir(nCol, nRow), "####0.0")
                Next nRow
                Close (1)
            End If
        End If
    Next I
    Screen.MousePointer = 0

    Unload Me
End Sub
Private Sub dirin_Change()
    If Len(DirIn.Path) <= 3 Then
        TheInPath = DirIn.Path
    Else
        TheInPath = DirIn.Path + "\"
    End If
    
    FileInT.Path = DirIn.Path
End Sub

Private Sub dirin_LostFocus()
    DirIn.Path = DirIn.List(DirIn.ListIndex)

End Sub


Private Sub driveIn_Change()
    On Error GoTo DriveHandler
    DirIn.Path = DriveIn.Drive
    Exit Sub
DriveHandler:
    DriveIn.Drive = DirIn.Path
    Exit Sub

End Sub

Private Sub FileInT_PathChange()
Call SearchFile
End Sub


Private Sub FileInT_PatternChange()
Call SearchFile
End Sub


Private Sub Form_Load()
Dim TheInPathT As String

TheInPath = App.Path

If (Right(TheInPath, 1) <> "\") Then
    TheInPath = TheInPath + "\"
End If

bClick = False
InIndex.Text = "*.GRD"
FileInT.Pattern = InIndex.Text
TheInPathT = TheInPath
DriveIn.Drive = Left(TheInPathT, 2)
DirIn.Path = TheInPathT
bClick = True
TheInPath = TheInPathT
Call SearchFile
End Sub
Private Sub InIndex_KeyDown(KeyCode As Integer, Shift As Integer)
On Error Resume Next

If (KeyCode = 13) Then
    FileInT.Pattern = InIndex.Text
End If
End Sub


Private Sub ListIn_Click()
Dim I As Integer

If (bClick = False) Then Exit Sub

If (ListIn.SelCount > 0) Then
    CommandOK.Enabled = True
Else
    CommandOK.Enabled = False
End If

End Sub

Private Sub Slope(ZGrid() As Double, zGridSlope() 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
    Dim X As Double

    ReDim zGridSlope(1 To nCols, 1 To nRows)
    
    '计算非边界点
    For J = 2 To nRows - 1
        For I = 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
            
            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) = 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 J
End Sub
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)
    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 = 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
            
            If (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(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(DYZx / DXZy) * 180 / 3.14159265
                End If
                If (ZGridDir(I, J) >= 360# Or ZGridDir(I, J) < 0) Then
                    ZGridDir(I, J) = 0
                End If
            End If
        Next I
    Next J
    
    '边界点赋值
    For J = 1 To nRows
        ZGridDir(1, J) = ZGridDir(2, J)
        ZGridDir(nCols, J) = ZGridDir(nCols - 1, J)
    Next J
    For I = 1 To nCols
        ZGridDir(I, 1) = ZGridDir(I, 2)
        ZGridDir(I, nRows) = ZGridDir(I, nRows - 1)
    Next I
    
    Vmin = ZGridDir(1, 1)
    Vmax = ZGridDir(1, 1)
    For J = 1 To nRows
        For I = 1 To nCols
            If (ZGridDir(I, J) > Vmax) Then Vmax = ZGridDir(I, J)
            If (ZGridDir(I, J) < Vmin) Then Vmin = ZGridDir(I, J)
        Next I
    Next J
End Sub
'0～180
Private Function ArcCos(X As Double) As Double
    If (X >= 1) Then
        ArcCos = 0
    ElseIf (X <= -1) Then
        ArcCos = 4# * Atn(1#)
    Else
        ArcCos = 2 * Atn(1#) + Atn(-X / Sqr(1 - X * X))
    End If
End Function