excel vba中大数据集中lat / long之间的最近距离

问题描述 投票:4回答:2

初学者在这里徘徊......我正在研究这个井间距项目,它看着纬度/长度并确定下一个最近的井。我想我可能正在创建一个无限循环,或者程序只是永远运行(它循环遍历15,000行)。我的主要努力一直在努力确保将每个位置与数据集中的每个位置进行比较。从那里我采取第二个最低距离(因为当它与自身相比时,最低距离将为零)。

Sub WellSpacing()
Dim r As Integer, c As Integer, L As Integer, lastrow As Integer
Dim lat1 As Double, lat2 As Double, long1 As Double, long2 As Double
Dim distance As Double, d1 As Double, d2 As Double, d3 As Double
Dim PI As Double

PI = Application.WorksheetFunction.PI()
L = 2
r = 3
c = 10
lastrow = Sheets("Test").Cells(Rows.Count, "J").End(xlUp).Row

For L = 2 To lastrow
    For r = 2 To lastrow
        lat1 = Sheets("Test").Cells(L, c)
        long1 = Sheets("Test").Cells(L, c + 1)
        lat2 = Sheets("Test").Cells(r, c)
        long2 = Sheets("Test").Cells(r, c + 1)
        d1 = Sin((Abs((lat2 - lat1)) * PI / 180 / 2)) ^ 2 + Cos(lat1 * PI / 180) * Cos(lat2 * PI / 180) * Sin(Abs(long2 - long1) * PI / 180 / 2) ^ 2
        d2 = 2 * Application.WorksheetFunction.Atan2(Sqr(1 - d1), Sqr(d1))
        d3 = 6371 * d2 * 3280.84
        Sheets("Working").Cells(r - 1, c - 9) = d3
    Next r

    Sheet2.Activate
    Range("A:A").Sort Key1:=Range("A1"), Order1:=xlAscending
    distance = Sheet2.Range("A2")
    Sheets("Test").Cells(L, c + 2) = distance
    Sheet2.Range("A:A").Clear
    Sheet1.Activate

Next L
End Sub
excel vba excel-vba
2个回答
1
投票

我建议使用数组,如@Qharr所说。我还希望通过包含一些逻辑步骤来加速这个过程,避免在每组点上进行复杂的数学运算。

我的意思是你可以先做一个粗略的估计,看是否打算做实际的计算。我一边看着当前位置的纬度或长度是否比最后一个点更接近,但你可以做任何你想做的事情。

我会将您的代码更改为:

Sub WellSpacing()
    Dim R As Integer, C As Integer, L As Integer, LastRow As Integer, Shortest() As Integer
    Dim Lats() As Double, Longs() As Double, Distances() As Double
    Dim Distance As Double, D1 As Double, D2 As Double, D3 As Double
    Dim PI As Double

    On Error Resume Next
    PI = Application.WorksheetFunction.PI()
    L = 2
    R = 3
    C = 10
    LastRow = Sheets("Test").Cells(Rows.Count, 10).End(xlUp).Row
    ReDim Lats(1 To (LastRow - 1)) As Double
    ReDim Longs(1 To (LastRow - 1)) As Double
    ReDim Distances(1 To (LastRow - 1)) As Double
    ReDim Shortest(1 To (LastRow - 1)) As Integer

    For L = 2 To LastRow
        Lats(L - 1) = Sheets("Test").Range("J" & L).Value
        Longs(L - 1) = Sheets("Test").Range("K" & L).Value
    Next L

    For L = 1 To (LastRow - 1)
        'This is a method of setting an initial value that can't be obtained  through the caclucations (so you will know if any calcs have been done or not).
        Distances(L) = -1
        For R = 1 To (LastRow - 1)
            'This minimises your calculations by 15,000 to begin with
            If R = L Then GoTo Skip_This_R
            'This skips checking the previous distances if it is the first calculation being checked.
            If Distances(L) = -1 Then GoTo Skip_Check
            'If there has already been a distance calculated, this does a rough check of whether the Lat or Long is closer. If neither
            'the Lat or Long are closer than the current closest, then it will skip it. This reduces the code by 7 lines for most pairs.
            If Abs(Lats(L) - Lats(R)) < Abs(Lats(L) - Lats(Shortest(L))) Or Abs(Longs(L) - Longs(R)) < Abs(Longs(L) - Longs(Shortest(L))) Then
Skip_Check:
                    D1 = Sin((Abs((Lats(R) - Lats(L))) * PI / 180 / 2)) ^ 2 + Cos(Lats(L) * PI / 180) * Cos(Lats(R) * PI / 180) * Sin(Abs(Longs(R) - Longs(L)) * PI / 180 / 2) ^ 2
                    D2 = 2 * Application.WorksheetFunction.Atan2(Sqr(1 - D1), Sqr(D1))
                    D3 = 6371 * D2 * 3280.84
                    If D3 < Distances(L) Or Distances(L) = -1 Then
                            Distances(L) = D3
                            'This stores the index value in the array of the closest Lat/Long point so far.
                            Shortest(L) = R
                    End If
            End If
Skip_This_R:
        Next R
        'This puts the resulting closest distance into the corresponding cell.
        Sheets("Test").Range("L" & (L + 1)).Value = Distances(L)
        'This clears any previous comments on the cell.
        Sheets("Test").Range("L" & (L + 1)).Comments.Delete
        'This adds a nice comment to let you know which Lat/Long position it is closest to.
        Sheets("Test").Range("L" & (L + 1)).AddComment "Matched to Row " & (Shortest(L) + 1)
    Next L
End Sub

2
投票

我最近一直在使用地理位置数学(aka,coordinate geometry)并且编写了一个sub来做你正在寻找的相同的东西。

您的代码可能不是创建无限循环,但计算数千个坐标之间的距离可能是处理器密集型的,即使对代码进行微小更改也会对处理时间产生巨大影响。


计算最近坐标对:强力方法

有许多用于确定最近点的算法,但最简单的编码(因此可能最适合一次性使用)被称为暴力方法。

For p1 = 1 to numPoints
    For p2 = p1 + 1 to numPoints
        ...calculate {distance}
        ...if {distance} < minDistance then minDist = {distance}
    Next p2
Next p1

使用此方法,将在x * ( n - 1 ) / 2点对之间计算距离。

例如,5个点的列表需要10个比较:

  1. 第1点↔第2点
  2. 第1点↔第3点
  3. 第1点↔第4点
  4. 第1点↔第5点
  5. 第2点↔第3点
  6. 第2点↔第4点
  7. 第2点↔第5点
  8. 第3点↔第4点
  9. 第3点↔第5点
  10. 第4点↔第5点

由于额外的点会以指数方式增加执行时间,因此该方法可能会产生一些冗长的处理时间,尤其是在较慢的机器上或具有过多的点数时。

我用来计算点之间的距离和比较点列表之间的距离的方法远远不是[代码更重]最有效的替代方案,但它们适用于我的“一次性”需求。

根据我的目的,我将在Excel和Access之间切换(几乎完全相同的代码),但Access更快,因此您可能希望将列表移动到表中并按此方式执行。

我比较的一个点列表有252个项目,这需要使用这种“简单代码”方法进行31,628次单独比较。在Excel中,该过程在1.12秒内完成,即Access只需0.16秒。

在我们开始使用更长的点列表之前,这似乎不是一个很大的区别:我的另一个列表(更接近你的大小)有大约12,000点,这需要使用Brute Force方法进行71,994,000次计算。在Access中,该过程在8.6分钟内完成,因此我估计在Excel中需要大约一个小时。

当然,所有这些时间都基于我的操作系统,处理能力,Office版本等.VBA并不适合这种级别的计算,你可以采取的一切措施来减少代码长度将产生很大的不同,包括评论状态栏更新,即时窗口输出,关闭屏幕更新等。

这段代码有点乱,没有评论,因为我为了自己的目的把它打了一遍,但它对我有用。如果您对其工作原理有任何疑问,请与我们联系。所有计算均以公制为单位,但可以轻松转换。

Sub findShortestDist_Excel()

    Const colLatitude = "C" ' Col.C = Lat, Col.D = Lon
    Dim pointList As Range, pointCount As Long, c As Range, _
        arrCoords(), x As Long, y As Long
    Dim thisDist As Double, minDist As Double, minDist_txt As String
    Dim cntCurr As Long, cntTotal As Long, timerStart As Single

    timerStart = Timer
    Set pointList = Sheets("Stops").UsedRange.Columns(colLatitude)
    pointCount = WorksheetFunction.Count(pointList.Columns(1))

    'build array of numbers found in Column C/D
    ReDim arrCoords(1 To 3, 1 To pointCount)
    For Each c In pointList.Columns(1).Cells
        If IsNumeric(c.Value) And Not IsEmpty(c.Value) Then
            x = x + 1
            arrCoords(1, x) = c.Value
            arrCoords(2, x) = c.Offset(0, 1).Value
        End If
    Next c

    minDist = -1
    cntTotal = pointCount * (pointCount + 1) / 2

    'loop through array
    For x = 1 To pointCount
        For y = x + 1 To pointCount
            If (arrCoords(1, x) & arrCoords(2, x)) <> (arrCoords(1, y) & arrCoords(2, y)) Then
                cntCurr = cntCurr + 1
                thisDist = Distance(arrCoords(1, x), arrCoords(2, x), _
                    arrCoords(1, y), arrCoords(2, y))
                'check if this distance is the smallest yet
                If ((thisDist < minDist) Or (minDist = -1)) And thisDist > 0 Then
                    minDist = thisDist
                    'minDist_txt = arrCoords(1, x) & "," & arrCoords(2, x) & " -> " & arrCoords(1, y) & "," & arrCoords(2, y)
                End If
                'Application.StatusBar = "Calculating Distances: " & Format(cntCurr / cntTotal, "0.0%")
            End If
        Next y
        'DoEvents
    Next x

    Debug.Print "Minimum distance: " & minDist_txt & " = " & minDist & " meters"
    Debug.Print "(" & Round(Timer - timerStart, 2) & "sec)"
    Application.StatusBar = "Finished.  Minimum distance: " & minDist_txt & " = " & minDist & "m"

End Sub

请注意,上述过程取决于以下过程(Access与Excel的版本略有不同):

Excel:计算点之间的距离

Public Function Distance(ByVal lat1 As Double, ByVal lon1 As Double, _
    ByVal lat2 As Double, ByVal lon2 As Double) As Double
'returns Meters distance in Excel (straight-line)
    Dim theta As Double: theta = lon1 - lon2
    Dim Dist As Double: Dist = Math.Sin(deg2rad(lat1)) * Math.Sin(deg2rad(lat2)) + Math.Cos(deg2rad(lat1)) * Math.Cos(deg2rad(lat2)) * Math.Cos(deg2rad(theta))
    Dist = rad2deg(WorksheetFunction.Acos(Dist))
    Distance = Dist * 60 * 1.1515 * 1.609344 * 1000
End Function

Function deg2rad(ByVal deg As Double) As Double
    deg2rad = (deg * WorksheetFunction.PI / 180#)
End Function

Function rad2deg(ByVal rad As Double) As Double
    rad2deg = rad / WorksheetFunction.PI * 180#
End Function

...和Microsoft Access的替代代码:

访问:最短距离

Sub findShortestDist_Access()

    Const tableName = "Stops"
    Dim pointCount As Long, arrCoords(), x As Long, y As Long
    Dim thisDist As Double, minDist As Double
    Dim cntCurr As Long, cntTotal As Long, timerStart As Single
    Dim rs As Recordset

    timerStart = Timer

    Set rs = CurrentDb.OpenRecordset("SELECT * FROM " & tableName)
    With rs
        .MoveLast
        .MoveFirst
        pointCount = .RecordCount

        'build array of numbers found in Column C/D
        ReDim arrCoords(1 To 2, 1 To pointCount)
        Do While Not .EOF
            x = x + 1
            arrCoords(1, x) = !stop_lat
            arrCoords(2, x) = !stop_lon
            .MoveNext
        Loop
        .Close
    End With

    minDist = -1
    cntTotal = pointCount * (pointCount + 1) / 2
    SysCmd acSysCmdInitMeter, "Calculating Distances:", cntTotal

    'loop through array
    For x = 1 To pointCount
        For y = x + 1 To pointCount
                cntCurr = cntCurr + 1
                thisDist = Distance(arrCoords(1, x), arrCoords(2, x), _
                    arrCoords(1, y), arrCoords(2, y))
                'check if this distance is the smallest yet
                If ((thisDist < minDist) Or (minDist = -1)) And thisDist > 0 Then
                    minDist = thisDist
                End If
                SysCmd acSysCmdUpdateMeter, cntCurr
        Next y
        DoEvents
    Next x
    SysCmd acSysCmdRemoveMeter
    Debug.Print "Minimum distance: " & minDist_txt & " = " & minDist & " meters"
    Debug.Print "(" & Round(Timer - timerStart, 2) & "sec)"

End Sub

请注意,上面的过程取决于以下...(Access可以更快地处理质量计算,但我们必须自己构建一些内置于Excel的函数)

访问:计算点之间的距离

Const pi As Double = 3.14159265358979

Public Function Distance(ByVal lat1 As Double, ByVal lon1 As Double, _
    ByVal lat2 As Double, ByVal lon2 As Double) As Double
'returns Meters distance in Access (straight-line)
    Dim theta As Double: theta = lon1 - lon2
    Dim dist As Double
    dist = Math.Sin(deg2rad(lat1)) * Math.Sin(deg2rad(lat2)) + Math.Cos(deg2rad(lat1)) _
        * Math.Cos(deg2rad(lat2)) * Math.Cos(deg2rad(theta))
    dist = rad2deg(aCos(dist))
    Distance = dist * 60 * 1.1515 * 1.609344 * 1000
End Function

Function deg2rad(ByVal deg As Double) As Double
    deg2rad = (deg * pi / 180#)
End Function

Function rad2deg(ByVal rad As Double) As Double
    rad2deg = rad / pi * 180#
End Function

Function aTan2(x As Double, y As Double) As Double
    aTan2 = Atn(y / x)
End Function

Function aCos(x As Double) As Double
    On Error GoTo aErr
    If x = 0 Or Abs(x) = 1 Then
        aCos = 0
    Else
        aCos = Atn(-x / Sqr(-x * x + 1)) + 2 * Atn(1)
    End If
    Exit Function
aErr:
    aCos = 0
End Function

平面案例

另一种计算近点的方法称为平面情况。我还没有看到任何现成的代码示例,我不需要它就足以打扰它编码,但它的要点是:

Planar Case

阅读关于维基百科上的Closest pair of points problem以及更多内容。

© www.soinside.com 2019 - 2024. All rights reserved.