# How to quickly identify nearby stations?

**URL:** <https://discourse.julialang.org/t/how-to-quickly-identify-nearby-stations/88124>\
**Category:** General Usage\
**Tags:** question\
**Created:** [October 2, 2022, 2:36pm UTC](https://discourse.julialang.org/t/how-to-quickly-identify-nearby-stations/88124 "2022-10-02T14:36:54Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![Leon6](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/leon6/32/205483_2.png) [@Leon6](https://discourse.julialang.org/u/Leon6)\
**Post date:** [October 2, 2022, 2:36pm UTC](https://discourse.julialang.org/t/how-to-quickly-identify-nearby-stations/88124/1 "2022-10-02T14:36:55Z")

</div>

I have a total of 1 million stations numbered from 1 to 1 million. For each of them, how do I quickly identify the rest of the stations (Station numbers) that are within 2 km of distance?

Any tips would be greatly appreciated. Many thanks.

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [October 2, 2022, 2:40pm UTC](https://discourse.julialang.org/t/how-to-quickly-identify-nearby-stations/88124/2 "2022-10-02T14:40:59Z")

</div>

How are the station positions defined? Which type of distance you want (euclidean, or something else)?

(The solution is probably in the [NearestNeighbors](https://github.com/KristofferC/NearestNeighbors.jl) package, or [similar alternatives](https://m3g.github.io/CellListMap.jl/stable/neighborlists/)).

---

<div class="post-metadata">

**Author:** ![Leon6](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/leon6/32/205483_2.png) [@Leon6](https://discourse.julialang.org/u/Leon6)\
**Post date:** [October 2, 2022, 3:47pm UTC](https://discourse.julialang.org/t/how-to-quickly-identify-nearby-stations/88124/3 "2022-10-02T15:47:38Z")

</div>

Many thanks for the reply.

These stations are from the global surface ocean. Each of them have its own longitude and latitude information. The distance should be the closest point-to-point distance on a globe.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [October 2, 2022, 4:04pm UTC](https://discourse.julialang.org/t/how-to-quickly-identify-nearby-stations/88124/4 "2022-10-02T16:04:04Z")

</div>

step 1 is to sort by laditude. then divide the points into 2km laditude bins and sort each bin by longitude. then for each point, you can get the list of possible close points by looking at 3 bins and doing a binary search to see which have close enough longitudes (this is latitude dependent). then for each candidate identified, you can perform a distance test to find the actual distance. Note that you will need a minor modification of this for the northern and southern most bins, but it is still relatively straightforward.

As described, this algorithm is `O(n*log(n)+close_pairs)` but with a little care, could be dropped to `O(n+close_pairs)` although as described it should be pretty fast already.

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [October 3, 2022, 12:12pm UTC](https://discourse.julialang.org/t/how-to-quickly-identify-nearby-stations/88124/5 "2022-10-03T12:12:06Z")

</div>

using the `NearestNeighbors.jl` package you can do this: first convert the spherical coordinates into cartesian coordinates, adjust the cutoff to a 3D cutoff, and do:

```julia
using NearestNeighbors
using StaticArrays

function points_in_sphere(r,N)
    p = SVector{3,Float64}[]
    for ip in 1:N
        θ = π*rand()
        ϕ = 2π*rand()
        x = r * sin(ϕ) * cos(θ)
        y = r * sin(ϕ) * sin(θ) 
        z = r * cos(ϕ)
        push!(p, SVector(x,y,z))
    end
    return p
end

function cartesian_cutoff(r, arch_cutoff)
    angle = arch_cutoff / r 
    cutoff = r * sqrt(2 * (1-cos(angle)))
    return cutoff
end

function main(;N=10^6)
    r = 6371 # km - earth radius
    arch_cutoff = 2 # km
    p = points_in_sphere(r, N)
    cutoff = cartesian_cutoff(r, arch_cutoff)
    tree = BallTree(p)
    inrange(tree, p, cutoff)
end

```

This takes 2s in my computer, thus probably solves the problem, although there may be smarter ways that do not involve converting the points to 3D space.

---

<div class="post-metadata">

**Author:** ![Leon6](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/leon6/32/205483_2.png) [@Leon6](https://discourse.julialang.org/u/Leon6)\
**Post date:** [October 3, 2022, 2:52pm UTC](https://discourse.julialang.org/t/how-to-quickly-identify-nearby-stations/88124/6 "2022-10-03T14:52:16Z")

</div>

This is awesome! Thank you so much.