测量计算方位角万能通用公式及VB、VBA源代码

使用该公式不用判断象限，直接得出方位角值

设有两点A、B，坐标分别为A（Xa，Ya）、B（Xb、Yb），则
$\Delta X=X_b-X_a$
$\Delta Y=Y_b-Y_a$
$\Delta Y=\Delta Y+1^{-10}$ 为了使除数不为零而加一个很小的数(注：在参与计算前加一个极小值后参与sign函数计算，这样sign($\Delta Y$)就只会出现1或-1两种情况)

则方位角：$\alpha=\pi \left(1-\frac {sign(\Delta Y)}2 \right)-arctan \left( \frac {\Delta X}{\Delta Y} \right)$，计算值单位为弧度，用公式 $\alpha=\frac {180\cdot \alpha}\pi$ 将角度单位换算为十进制度，进而换算为度、分、秒格式。

其中，sign()为求符号函数，有些软件该函数名为sgn()，该函数在计算时若参数$\Delta X$<0时其值为-1，$\Delta X$>0时值为1，$\Delta X$=0时取值为0。使用此公式不用判断所在象限，直接将坐标增量代入即可求出方位角值，在用计算器编程时若没有sign()函数可自行判断并用一个变量代替!

VBA代码：

'方位角计算函数 Azimuth()
'Sx为起点X，Sy为起点Y
'Ex为终点X，Ey为终点Y
'Style指明返回值格式
'Style=-1为弧度格式
'Style=0为“DD MM SS”格式
'Style=1为“DD-MM-SS”格式
'Style=2为“DD°MMˊSS""”格式
'Style=其它值时返回十进制度值
Function Azimuth(Sx As Double, Sy As Double, Ex As Double, Ey As Double, Style As Integer)
Dim DltX As Double, DltY As Double, A_tmp As Double, Pi As Double
Pi = Atn(1) * 4 '定义PI值
DltX = Ex - Sx
DltY = Ey - Sy + 1E-20
A_tmp = Pi * (1 - Sgn(DltY) / 2) - Atn(DltX / DltY) '计算方位角
A_tmp = A_tmp * 180 / Pi '转换为360进制角度
Azimuth = Deg2DMS(A_tmp, Style)
End Function

'转换角度为度分秒
'Style=-1为弧度格式
'Style=0为“DD MM SS”格式
'Style=1为“DD-MM-SS”格式
'Style=2为“DD°MMˊSS""”格式
'Style=其它值时返回十进制度值
Function Deg2DMS(DegValue As Double, Style As Integer)
Dim tD As Integer, tM As Integer, Ts As Double, tmp As Double
tD = Int(DegValue)
tmp = (DegValue - tD) * 60
tM = Int(tmp)
tmp = (tmp - tM) * 60
Ts = Round(tmp, 1)
Select Case Style
Case -1 '返回弧度
Deg2DMS = DegValue * Atn(1) * 4 / 180
Case 0
Deg2DMS = tD & " " & Format(tM, "00") & " " & Format(Ts, "00.0")
Case 1
Deg2DMS = tD & "-" & Format(tM, "00") & "-" & Format(Ts, "00.0")
Case 2
Deg2DMS = tD & "°" & Format(tM, "00") & "ˊ" & Format(Ts, "00.0") & """"
Case Else
Deg2DMS = DegValue
End Select
End Function

Function aa(area1 As Double, area2 As Double) As Double
Dim rat As Double
rat = area1 / area2
If (rat < 0.6 Or rat > (1 / 0.6)) And area1 <> 0 And area2 <> 0 Then
aa = (area1 + area2 + sqrt(area1 * area2)) / 3
Else
aa = (area1 + area2) / 2
End If
End Function

Function Distance(Sx As Double, Sy As Double, Ex As Double, Ey As Double, Precision As Integer) As Double
Dim DltX As Double, DltY As Double
DltX = Ex - Sx
DltY = Ey - Sy
Distance = Round(Sqr(DltX * DltX + DltY * DltY), Precision)
End Function

Function inValue(stgA As Double, Av As Double, stgB As Double, Bv As Double, stgC As Double) As Double
If stgB <> stgA Then
inValue = Av + (Bv - Av) / (stgB - stgA) * (stgC - stgA)
Else
inValue = -9.9999999
End If
End Function

Function pol(AX As Double, AY As Double, Bx As Double, By As Double) As String
pol = Azimuth(AX, AY, Bx, By, 2) & " " & Distance(AX, AY, Bx, By, 3)
End Function

Function rec(alpha As String, dist As Double) As String
Dim Alpha_Rad As Double
Alpha_Rad = StringToRad(alpha)
rec = "dx:" & Round(Cos(Alpha_Rad) * dist, 3) & " dy:" & Round(Sin(Alpha_Rad) * dist, 3)
End Function

Function StringToRad(strAz) '将字符串格式方位角转换成弧度格式
Dim azSubStr
If strAz <> "" Then
azSubStr = Split(strAz, "-")
If UBound(azSubStr) = 2 Then
StringToRad = (azSubStr(0) + azSubStr(1) / 60 + azSubStr(2) / 3600) * Atn(1) * 4 / 180
Else
StringToRad = 0
End If
Else
StringToRad = 0
End If
End Function

'竹山龙背湾 2010-09-08
'判断是否存在坐标系定义表
Function CoSysTableExist() As Boolean
Dim i As Long
CoSysTableExist = False
For i = 1 To Sheets.Count
If Sheets(i).Name = "CoSys" Then
CoSysTableExist = True
Exit For
End If
Next
'If Not CoSysTableExist Then
'Dim NewTable As Sheets
'End If
End Function

'查找坐标系名称并返回参数
Function CoSysFndPara(CoSysName As String) As String
Dim FndIndex As Long
If CoSysTableExist Then
    For FndIndex = 1 To 100
        If Trim(Sheets("CoSys").Range("A" & FndIndex).Text) = Trim(CoSysName) Then
            CoSysFndPara = Trim(Sheets("CoSys").Range("B" & FndIndex).Text)                      'AX
            CoSysFndPara = CoSysFndPara & "," & Trim(Sheets("CoSys").Range("C" & FndIndex).Text) 'AY
            CoSysFndPara = CoSysFndPara & "," & Trim(Sheets("CoSys").Range("D" & FndIndex).Text) 'Ax
            CoSysFndPara = CoSysFndPara & "," & Trim(Sheets("CoSys").Range("E" & FndIndex).Text) 'Ay
            If InStr(Trim(Sheets("CoSys").Range("F" & FndIndex).Text), "-") <> 0 Then
            CoSysFndPara = CoSysFndPara & "," & Trim(Sheets("CoSys").Range("F" & FndIndex).Text) 'az
            Else
            CoSysFndPara = CoSysFndPara & "," & Azimuth(Trim(Sheets("CoSys").Range("B" & FndIndex).Text), Trim(Sheets("CoSys").Range("C" & FndIndex).Text), Trim(Sheets("CoSys").Range("F" & FndIndex).Text), Trim(Sheets("CoSys").Range("G" & FndIndex).Text), 1) 'BY or Type
            End If
            Exit For
        End If
    Next
Else
    CoSysFndPara = ""
End If
End Function

'测图坐标转施工坐标
Function NE2SO_STG(CoSysName As String, P_N As Double, P_E As Double) As Double
Dim coSysParaStr As String
Dim coSysPara
Dim O_X As Double, O_Y As Double, O_Stage As Double, O_Offset As Double, X_Line_Azimuth_Str As Double
'读取坐标系参数
coSysParaStr = CoSysFndPara(CoSysName)
coSysPara = Split(coSysParaStr, ",")
O_X = coSysPara(0)         '基点测图坐标
O_Y = coSysPara(1)
O_Stage = coSysPara(2)     '基点施工坐标
O_Offset = coSysPara(3)
X_Line_Azimuth_Str = StringToRad(coSysPara(4)) '施工坐标系X轴方位角,弧度
NE2SO_STG = Round((P_N - O_X) * Cos(X_Line_Azimuth_Str) + (P_E - O_Y) * Sin(X_Line_Azimuth_Str) + O_Stage, 3)
End Function

'测图坐标转施工坐标
Function NE2SO_OFF(CoSysName As String, P_N As Double, P_E As Double) As Double
Dim coSysParaStr As String
Dim coSysPara
Dim O_X As Double, O_Y As Double, O_Stage As Double, O_Offset As Double, X_Line_Azimuth_Str As Double
'读取坐标系参数
coSysParaStr = CoSysFndPara(CoSysName)
coSysPara = Split(coSysParaStr, ",")
O_X = coSysPara(0)         '基点测图坐标
O_Y = coSysPara(1)
O_Stage = coSysPara(2)     '基点施工坐标
O_Offset = coSysPara(3)
X_Line_Azimuth_Str = StringToRad(coSysPara(4)) '施工坐标系X轴方位角,弧度
NE2SO_OFF = Round(-(P_N - O_X) * Sin(X_Line_Azimuth_Str) + (P_E - O_Y) * Cos(X_Line_Azimuth_Str) + O_Offset, 3)
End Function

'测图坐标转施工坐标
Function SO2NE_N(CoSysName As String, P_x As Double, P_y As Double) As Double
Dim coSysParaStr As String
Dim coSysPara
Dim O_X As Double, O_Y As Double, O_Stage As Double, O_Offset As Double, X_Line_Azimuth_Str As Double
'读取坐标系参数
coSysParaStr = CoSysFndPara(CoSysName)
coSysPara = Split(coSysParaStr, ",")
O_X = coSysPara(0)         '基点测图坐标
O_Y = coSysPara(1)
O_Stage = coSysPara(2)     '基点施工坐标
O_Offset = coSysPara(3)
X_Line_Azimuth_Str = StringToRad(coSysPara(4)) '施工坐标系X轴方位角,弧度
SO2NE_N = Round(O_X + (P_x - O_Stage) * Cos(X_Line_Azimuth_Str) - (P_y - O_Offset) * Sin(X_Line_Azimuth_Str), 3)
End Function

'测图坐标转施工坐标
Function SO2NE_E(CoSysName As String, P_x As Double, P_y As Double) As Double
Dim coSysParaStr As String
Dim coSysPara
Dim O_X As Double, O_Y As Double, O_Stage As Double, O_Offset As Double, X_Line_Azimuth_Str As Double
'读取坐标系参数
coSysParaStr = CoSysFndPara(CoSysName)
coSysPara = Split(coSysParaStr, ",")
O_X = coSysPara(0)         '基点测图坐标
O_Y = coSysPara(1)
O_Stage = coSysPara(2)     '基点施工坐标
O_Offset = coSysPara(3)
X_Line_Azimuth_Str = StringToRad(coSysPara(4)) '施工坐标系X轴方位角,弧度
SO2NE_E = Round(O_Y + (P_x - O_Stage) * Sin(X_Line_Azimuth_Str) + (P_y - O_Offset) * Cos(X_Line_Azimuth_Str), 3)
End Function

转载自：https://blog.csdn.net/QinDongZ/article/details/80388113