图幅面积.bas

根据西南角坐标及比例尺

Attribute VB_Name = "Module1"
Option Explicit
'54坐标系
Public Const PI As Double = 3.14159265358979

Public Const A0 As Double = 6367558.4969
Public Const B0 As Double = 32005.7798
Public Const C0 As Double = 133.9238
Public Const D0 As Double = 0.6973
Public Const E0 As Double = 0.0039

Public Const EE As Double = 0.006693421623
Public Const e2 As Double = 0.0067385254147
Public Const C As Double = 6399698.90178271
Public Const B As Double = 6356863.01877305

Public Const DB10 As Double = 1 / 3             '十万图纬差十进制度
Public Const DL10 As Double = 0.5               '十万图经差十进制度
Public Const DB5 As Double = 1 / 6
Public Const DL5 As Double = 0.25
Public Const DB25 As Double = 1 / 12
Public Const DL25 As Double = 0.125
Public Const DB1 As Double = 1 / 24
Public Const DL1 As Double = 1 / 16

Public Const D10 As Integer = 8             '纵向或横向一万图幅数,即10万图加密至个8*4图廓点
Public Const D5 As Integer = 4              '纵向或横向一万图幅数,即5万图加密至个4*4图廓点
Public Const D25 As Integer = 2             '纵向或横向一万图幅数,即2.5万图加密至2*4个图廓点
Public Const D1 As Integer = 1              '纵向或横向一万图幅数,即一万图1*4个图廓点
Sub Main()
  Load Form1
  Form1.Show
End Sub
Function BLTOX(ByVal B As Double, ByVal L As Double) As Double           'B-纬度,L-与中央子午经差(单位为弧度),返回值X,单位为m
  Dim X As Double, M As Double, T As Double, n2 As Double, N As Double, BB As Double
  X = A0 * B - (B0 * Sin(B) + C0 * Sin(B) ^ 3 + D0 * Sin(B) ^ 5 + E0 * Sin(B) ^ 7) * Cos(B)
  M = Cos(B) * L
  T = Tan(B)
  n2 = e2 * Cos(B) ^ 2
  N = C / Sqr(1 + n2)
  BB = 0.5 * M ^ 2 + ((5 - T ^ 2 + 9 * n2 + 4 * n2 ^ 2) * M ^ 4) / 24 + ((61 - 58 * T ^ 2 + T ^ 4) * M ^ 6) / 720
  BLTOX = X + N * T * BB
End Function
Function BLTOY(ByVal B As Double, ByVal L As Double) As Double            'B-纬度,L-与中央子午经差(单位为弧度),返回值Y,单位为m
  Dim M As Double, T As Double, n2 As Double, N As Double, BB As Double
  M = Cos(B) * L
  T = Tan(B)
  n2 = e2 * Cos(B) ^ 2
  N = C / Sqr(1 + n2)
  BB = (M + ((1 - T ^ 2 + n2) * M ^ 3) / 6 + ((5 - 18 * T ^ 2 + T ^ 4 + 14 * n2 - 58 * n2 * T ^ 2) * M ^ 5) / 120)
  BLTOY = N * BB
End Function
Sub PXY(ByVal B As Double, ByVal L As Double, ByVal DB As Double, ByVal DL As Double, ByVal DN As Integer, P() As Double)
    'PXY()-计算图廓点坐标数组.参数表:(左下角纬度,左下角经度[与中央子午线经差],纬差,经差,纵向或横向的一万图幅数[用于加密图廓点],
    '返回值[图廓点坐标数组P(X,Y)]从左下角起顺时针沿图廓回至起点)   [角度均为十进制度]
    '图廓线按1万图的经纬差(DB1,DL1)加密
     Dim I As Integer, K As Integer
     Dim BL As Double, LL As Double, BR As Double, LR As Double, BW As Double, LC As Double
     
     BL = B * PI / 180                          '左下角纬度
     LL = L * PI / 180                          '左下角经度
     BR = (B + DB) * PI / 180                        '右上角纬度
     LR = (L + DL) * PI / 180                        '右上角经度

     For I = 1 To DN                       '从左下角起分别加密西和南图廓线
       BW = (B + (I - 1) * DB1) * PI / 180       '西图廓线顺时针第I点图廓点纬度(按一万经纬差加密)总图廓点数为4*DN+1
       P(I, 1) = BLTOX(BW, LL)                   '从左下角起顺时针计算第I点X
       P(I, 2) = BLTOY(BW, LL) + 500000          '从左下角起顺时针计算第I点Y(加500公里)
 
       LC = (L + (I - 1) * DL1) * PI / 180              '南图廓线逆时针倒数第I点图廓点经度(按一万经纬差加密)
       P(DN * 4 + 2 - I, 1) = BLTOX(BL, LC)             '从左下角起逆时针计算倒数第I点X
       P(DN * 4 + 2 - I, 2) = BLTOY(BL, LC) + 500000    '从左下角起逆时针计算倒数第I点Y
     Next I
       K = DN * 3 + 2                                '顺时针已计算点数
     For I = 1 To DN                      '从左上角和右下角起分别加密北和东图廓线
       LC = (L + (I - 1) * DL1) * PI / 180           '北图廓线顺时针第I点图廓点纬度(按一万经纬差加密)
       P(DN + I, 1) = BLTOX(BR, LC)                  '从左上角起顺时针计算第I点X
       P(DN + I, 2) = BLTOY(BR, LC) + 500000         '从左上角起顺时针计算第I点Y
       
       BW = (B + (I - 1) * DB1) * PI / 180        '东图廓线逆时针倒数第I点图廓点纬度
       P(K - I, 1) = BLTOX(BW, LR)                '从右下角起逆时针计算倒数第I点X
       P(K - I, 2) = BLTOY(BW, LR) + 500000       '从右下角起逆时针计算倒数第I点Y
     Next I
     P(DN * 2 + 1, 1) = BLTOX(BR, LR)               '计算右上角X
     P(DN * 2 + 1, 2) = BLTOY(BR, LR) + 500000      '计算右上角Y
End Sub
Function PMJ(ByVal DN As Double, P() As Double) As Double   '计算平面面积.参数表:(纵向或横向一万图幅数即加密点数,坐标数组)
   Dim S As Double, I As Integer
   S = 0
   For I = 1 To DN * 4
       S = S + P(I, 1) * P(I + 1, 2) - P(I, 2) * P(I + 1, 1)
   Next I
   PMJ = 0.5 * Abs(S)
End Function
Function QMJ(ByVal B1 As Double, ByVal B2 As Double, ByVal DL As Double)       '计算曲面面积.参数表:(左下角纬度,左上角纬度,经差)  [角度为十进制度]
   Dim S As Double
   B1 = B1 * PI / 180
   B2 = B2 * PI / 180
   DL = DL * PI / 180
   S = Sin(B2) - Sin(B1) + EE * (Sin(B2) ^ 3 - Sin(B1) ^ 3) * 2 / 3 + EE ^ 2 * (Sin(B2) ^ 5 - Sin(B1) ^ 5) * 3 / 5 + EE ^ 3 * (Sin(B2) ^ 7 - Sin(B1) ^ 7) * 4 / 7
   QMJ = B ^ 2 * DL * S
End Function
Function DEG(ByVal A As String) As Double              '度.分秒制转换为十进制度,参数表:(度.分秒制角度字符串)
  Dim N As Integer, N0 As Integer, N1 As Integer, K As Integer
  Dim DS As String, MS As String, SS As String, A0 As String, TS As String, SSS As String
  Dim D As Double, M As Double, S As Double
  K = 0
  If A = "0" Then                       '输入参数为0
     DEG = 0
  ElseIf InStr(1, A, ".") = 0 Then      '没有小数点,输入参数为整度数
     DEG = Val(A)
  Else                                  '输入参数非整度数
     If Val(A) < 0 Then                 '负数,去 - 号
        K = 1
        A = Str(Abs(Val(A)))
     End If
     N0 = InStr(1, A, ".")                '小数点的位置
     TS = Mid(A, N0 + 1)                  '小数点右侧字符串
     N1 = Len(TS)                         '小数点右侧字符串的长度,若长度>=4,则秒的个位有数,不需补0
     If N1 < 4 Then   '小数点右侧字符串的长度<4,秒的个位左侧还有空位,应补0至秒位,以避免分的十位或秒的十位以后无数字的情况而产生的如30处理为3错误,
        A0 = A & String(4 - N1, "0")
     Else                                 'TS的长度>=4,输入参数的秒的个位(或秒的小数位)有数,不需补0
        A0 = A
     End If
     DS = Left(A0, N0 - 1)              '小数点左侧整数位
     MS = Mid(A0, N0 + 1, 2)            '小数点右侧第1,2位(分位)
     SS = Mid(A0, N0 + 3, 2)            '小数点右侧第3,4位(秒位)
     If N1 > 4 Then               '秒位有小数
        SSS = Mid(A0, N0 + 5)     '秒的小数位
        SS = SS & "." & SSS       '连接秒的整数位和小数位
     End If
     D = Val(DS)
     M = Val(MS)
     S = Val(SS)
     DEG = D + M / 60 + S / 3600
     DEG = IIf(K = 1, -DEG, DEG)        '输入参数为负值
  End If
End Function