XL 2016 Affichage villes proches à partir coordonnées GPS

    If Abs(Lat - t(j, 3)) * RT < dmax Then 'test sur les latitudes
      da = sinLat * Sin(t(j, 3)) + cosLat * Cos(t(j, 3)) * Cos(Lon - t(j, 4)) 'cosinus de la distance angulaire
      da = Atn(Sqr(Abs(1 - da ^ 2)) / da) 'distance angulaire en radian
      If da * RT < dmax Then

Sub Classer()
Dim dur#, Ni&, i&, a, d!, vmax!, ville$, vmin&
dur = Timer
With Feuil2.[D2:D36467] 'plage à adapter éventuellement
  Ni = .Count
  For i = 1 To Ni
    If .Cells(i) <> "" Then
      a = Split(.Cells(i), "#")
      tri a, 0, UBound(a)
      .Cells(i, 0).Resize(, UBound(a) + 1) = a
      d = Val(Replace(.Cells(i), ",", "."))
      If d > vmax Then vmax = d: ville = Feuil1.Cells(i + 1, 1)
    End If
    If i Mod 100 = 1 Then Application.StatusBar = "Temps écoulé " & Format((Timer - dur) / 86400, "hh:mm:ss") & _
      "   Réalisé " & Format(i / Ni, "0.0 %") 'affichage de la progression dans la barre d'état
  Next i
  .Columns(0) = "=COUNTA(" & .Cells(1).Resize(, Columns.Count - 3).Address(0, 0) & ")" 'Nb villes proches
  vmin = Application.Min(.Columns(0))
End With
Application.StatusBar = ""
dur = (Timer - dur) / 86400
MsgBox "Durée du classement " & Minute(dur) & " min " & Second(dur) & " s" & vbLf & _
  "Nombre minimum de villes proches " & vmin & vbLf & _
  "Distance maximum de la ville la plus proche " & Format(vmax, "0.0") & " km pour " & ville
End Sub

Sub tri(a, gauc, droi)  ' Quick sort
Dim ref, g, d, temp
ref = a((gauc + droi) \ 2)
g = gauc: d = droi
Do
    Do While a(g) < ref: g = g + 1: Loop
    Do While ref < a(d): d = d - 1: Loop
    If g <= d Then
      temp = a(g): a(g) = a(d): a(d) = temp
      g = g + 1: d = d - 1
    End If
Loop While g <= d
If g < droi Then Call tri(a, g, droi)
If gauc < d Then Call tri(a, gauc, d)
End Sub

Dim flag As Boolean 'mémorise la variable (bloque la MsgBox de la macro RAZ)

Sub Calcul()
Dim dur#, RT#, dmax&, t, ub&, Ncombi&, resu$(), f$, i&, Lat#, Lon#, sinLat#, cosLat#, Latmax#, j&, combi&, da#, x$, n&
dur = Timer
RT = 6378.137 'rayon terrestre en km
dmax = Int(Val(Feuil1.[P1])) 'distance maximum retenue en km, cellule à adapter éventuellement
With Feuil1.[A1].CurrentRegion.Offset(1)
  If .Rows.Count = 1 Then Exit Sub
  .Sort .Columns(3), xlAscending, Header:=xlNo 'tri croissant sur les latitudes, pour accélérer
  t = .Resize(.Rows.Count - 1)
  ub = UBound(t)
  Ncombi = Application.Combin(ub, 2) 'nombre de distances à calculer
  '---tableau des résultats---
  ReDim resu(1 To ub, 1 To 1)
  f = String(Len(CStr(dmax)), "0") & ".0 k\m " 'format des distances
  For i = 1 To ub - 1
    Lat = t(i, 3): Lon = t(i, 4): sinLat = Sin(Lat): cosLat = Cos(Lat)
    Latmax = Lat + dmax / RT
    For j = i + 1 To ub
      If t(j, 3) > Latmax Then combi = combi + ub - j + 1: Exit For 'test sur les latitudes
      combi = combi + 1
      da = sinLat * Sin(t(j, 3)) + cosLat * Cos(t(j, 3)) * Cos(Lon - t(j, 4)) 'cosinus de la distance angulaire
      da = Atn(Sqr(Abs(1 - da ^ 2)) / da) 'distance angulaire en radian
      If da * RT < dmax Then
        x = Format(da * RT, f)
        resu(i, 1) = resu(i, 1) & x & t(j, 1) & "#"
        resu(j, 1) = resu(j, 1) & x & t(i, 1) & "#"
        If n Mod 10000 = 0 Then Application.StatusBar = "Temps écoulé " & Format((Timer - dur) / 86400, "hh:mm:ss") & _
          "   Réalisé " & Format(combi / Ncombi, "0.0 %") 'affichage de la progression dans la barre d'état
        n = n + 1
      End If
  Next j, i
  '---restitution---
  flag = True: Call RAZ: flag = False
  .Cells(1).Resize(ub, 2).Copy Feuil2.[A2]
  Feuil2.[D2].Resize(ub) = resu
  Feuil2.[A2].Resize(ub, 4).Sort Feuil2.[A2], xlAscending, Header:=xlNo, Orientation:=xlTopToBottom 'tri alphabétique
  Feuil2.Activate
  .Sort .Cells(1), xlAscending, Header:=xlNo 'tri alphabétique
End With
dur = (Timer - dur) / 86400
MsgBox "Nombre de distances retenues " & Format(n, "#,##0") & vbLf & "Durée du calcul " & Minute(dur) & " min " & Second(dur) & " s"
Application.StatusBar = ""
End Sub

      For j = 1 To IIf(ub1 < Nville, ub1, Nville)
        f = a(j)
        Do While Left(f, 1) = "0" 'suppression des zéros non significatifs
          If Not IsNumeric(Mid(f, 2, 1)) Then Exit Do
          f = Mid(f, 2)
        Loop
        resu(i, j) = f
      Next j

resu

Private Sub Worksheet_Change(ByVal Target As Range)
If Target.Address <> "$G$1" Then Exit Sub
Dim RT$, tablo, i As Variant, resu(), Lat, Lon, sinLat, cosLat, j&, da, n&
RT = 6378.137 'rayon terrestre en km
tablo = [A1].CurrentRegion.Resize(, 4)
i = Application.Match(Target, Columns(1), 0)
If IsNumeric(i) And UBound(tablo) > 2 Then
    ReDim resu(1 To UBound(tablo) - 2, 1 To 2)
    Lat = tablo(i, 3): Lon = tablo(i, 4): sinLat = Sin(Lat): cosLat = Cos(Lat)
    For j = 2 To UBound(tablo)
        If j <> i Then
            da = sinLat * Sin(tablo(j, 3)) + cosLat * Cos(tablo(j, 3)) * Cos(Lon - tablo(j, 4)) 'cosinus de la distance angulaire
            da = Atn(Sqr(Abs(1 - da ^ 2)) / da) 'distance angulaire en radian
            n = n + 1
            resu(n, 1) = tablo(j, 1)
            resu(n, 2) = RT * da
        End If
    Next j
End If
'---restitution---
Application.ScreenUpdating = False
If FilterMode Then ShowAllData 'si la feuille est filtrée
With [G2] '1ère cellule de destination, à adapter
    If n Then
        .Resize(n, 2) = resu
        .Resize(n, 2).Sort .Cells(1, 2), xlAscending, Header:=xlNo 'tri sur les distances
    End If
    .Offset(n).Resize(Rows.Count - n - .Row + 1, 2).ClearContents 'RAZ en dessous
End With
End Sub

Private Sub Worksheet_Change(ByVal Target As Range)
If Intersect(Target, [G1,K1]) Is Nothing Then Exit Sub
Dim RT$, daMax, tablo, i As Variant, resu(), Lat, Lon, sinLat, cosLat, j&, da, n&
RT = 6378.137 'rayon terrestre en km
daMax = [K1] / RT 'distance angulaire max
tablo = [A1].CurrentRegion.Resize(, 4)
i = Application.Match([G1], Columns(1), 0)
If IsNumeric(i) And UBound(tablo) > 2 Then
    ReDim resu(1 To UBound(tablo) - 2, 1 To 2)
    Lat = tablo(i, 3): Lon = tablo(i, 4): sinLat = Sin(Lat): cosLat = Cos(Lat)
    For j = 2 To UBound(tablo)
        If j <> i Then
            da = sinLat * Sin(tablo(j, 3)) + cosLat * Cos(tablo(j, 3)) * Cos(Lon - tablo(j, 4)) 'cosinus de la distance angulaire
            da = Atn(Sqr(Abs(1 - da ^ 2)) / da) 'distance angulaire en radian
            If da <= daMax Then
                n = n + 1
                resu(n, 1) = tablo(j, 1)
                resu(n, 2) = RT * da
            End If
        End If
    Next j
End If
'---restitution---
Application.ScreenUpdating = False
If FilterMode Then ShowAllData 'si la feuille est filtrée
With [G2] '1ère cellule de destination, à adapter
    If n Then
        .Resize(n, 2) = resu
        .Resize(n, 2).Sort .Cells(1, 2), xlAscending, Header:=xlNo 'tri sur les distances
    End If
    .Offset(n).Resize(Rows.Count - n - .Row + 1, 2).ClearContents 'RAZ en dessous
End With
End Sub

.Resize(n, 2).Sort .Cells(1, 2), xlAscending, Header:=xlNo 'tri sur les distances