# Find nearest neighbours among boundary sets

**URL:** <https://discourse.julialang.org/t/find-nearest-neighbours-among-boundary-sets/122342>\
**Category:** General Usage\
**Tags:** question, geo\
**Created:** [November 6, 2024, 5:22pm UTC](https://discourse.julialang.org/t/find-nearest-neighbours-among-boundary-sets/122342 "2024-11-06T17:22:10Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![TimG](https://avatars.discourse-cdn.com/v4/letter/t/82dd89/32.png) [@TimG](https://discourse.julialang.org/u/TimG)\
**Post date:** [November 6, 2024, 5:22pm UTC](https://discourse.julialang.org/t/find-nearest-neighbours-among-boundary-sets/122342/1 "2024-11-06T17:22:10Z")

</div>

I have a set of boundaries defining UK electoral constituencies. I can pick one and look for its nearest neighbours like this:

```julia
con_bounds = GeoIO.load("Westminster_Parliamentary_Constituencies_July_2024_Boundaries_UK_BFC") |> Rename("PCON24CD" => "con_code")
inds = search(centroid(con_bounds[con_bounds.con_code.=="W07000093", :geometry]), KNearestSearch(con_bounds.geometry, 8))
gt=con_bounds[inds, :]

```

Or I could use

```julia
distance=evaluate(Euclidean(), centroid(g1), centroid(g2))
dist = distance.(Ref(con_bounds[con_bounds.con_code.=="W07000093", :geometry]), con_bounds.geometry)

```

Either case will produce something like this:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/4/c/4c337ead9ab40f2b4c88d1207b7cd84b01588009.png)

The constituency I’m interested in is the very dark one on the north west. It is on the coast. However, at the north east of this constituency, it has another immediately adjacent neighbour. This isn’t found in the nearest neighbour search (presumably) because the centroids are too far apart. Both constituencies are relatively large compared to others nearby.

However, I want to find nearest neighbours preferring immediately adjacent constituencies over those whose centroids are closer but that are not immediately adjacent.

How can I do this?

I have created a naive distance function based on the vertices in the boundary:

```julia
function distance(g1, g2)
    first=true
    let dist
        for p in vertices(g2)
            if first
                dist = evaluate(Euclidean(), centroid(g1), p)
                first=false
            else
                dist = min(dist, evaluate(Euclidean(), centroid(g1), p))
            end
        end
        return dist
    end
end

```

I’m hoping to find the nearest point on the neighbouring boundary to the centroid of the constituency of interest. But there are many points in each boundary and 650 constituencies, so this doesn’t run fast.  
There must be a better way!  
Thanks!

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [November 6, 2024, 6:20pm UTC](https://discourse.julialang.org/t/find-nearest-neighbours-among-boundary-sets/122342/2 "2024-11-06T18:20:07Z")

</div>

That is a very interesting problem. A couple of years ago I implemented a `ModifiedHausdorff` distance in ImageDistances.jl that may have the definition you need. The Wikipedia article is a bit heavy on the math notation, but the idea is very simple. Check their illustration and search for more didactic material online. I can help with writing efficient code if that distance solves your problem:

> **[Hausdorff distance](https://en.wikipedia.org/wiki/Hausdorff_distance)**
>
> In mathematics, the Hausdorff distance, or Hausdorff metric, also called Pompeiu–Hausdorff distance, measures how far two subsets of a metric space are from each other. It turns the set of non-empty compact subsets of a metric space into a metric space in its own right. It is named after Felix Hausdorff and Dimitrie Pompeiu.
> Informally, two sets are close in the Hausdorff distance if every point of either set is close to some point of the other set. The Hausdorff distance is the longest distance...

This distance would take two geometries as input and would produce a single distance value, which you could then use to find the neighbor of interest.

---

<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:** [November 6, 2024, 6:40pm UTC](https://discourse.julialang.org/t/find-nearest-neighbours-among-boundary-sets/122342/3 "2024-11-06T18:40:59Z")

</div>

If I understand the problem correctly (significant chances I didn’t), you might get an acceptable performance if you (with, for example, the NearestNeighbors.jl package):

1. Construct a BallTree for each constituencies.
2. Loop over the centroids of the other constituencies searching for the nearest point in the KDTree above to the given centroid, and keeping (as I understood) the closest one among the other constituencies.

Step 1 should be ok because it is O(N) for each constituency and must be performed only once. Step 2 should be really quick for each centroid.

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [November 6, 2024, 6:53pm UTC](https://discourse.julialang.org/t/find-nearest-neighbours-among-boundary-sets/122342/4 "2024-11-06T18:53:27Z")

</div>

Length of the shared boundary?

---

<div class="post-metadata">

**Author:** ![TimG](https://avatars.discourse-cdn.com/v4/letter/t/82dd89/32.png) [@TimG](https://discourse.julialang.org/u/TimG)\
**Post date:** [November 7, 2024, 9:11am UTC](https://discourse.julialang.org/t/find-nearest-neighbours-among-boundary-sets/122342/5 "2024-11-07T09:11:59Z")

</div>

Thank you for the suggestions.

I want to be able to arbitrarily pick any one of the 650 UK constituencies and produce a map like the one above showing it and its immediate neighbours. It doesn’t matter if I have a couple of extra constituencies shown but I must show all immediate neighbours.

In my example, the neighbours to the north are the 10th and 16th nearest by centroid difference.

A simple way to speed up my naive method is therefore to pick the nearest (say) 40 constituencies based on centroid distance and only look for the closest point on the boundaries of these 40. This is obviously vastly quicker, and produces this (for 12 nearest constituencies), which is what I need:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/f/4/f4a6c2e7c472df7011a6e81b22957c14e66aeb34.png)  
However, even this took 1253 seconds, which still seems too long to me.

@lmiq I think I will try your method first to look for a speed-up. It seems easiest since I’m already using NearestNeighbors.jl.  
@juliohm I’m a bit intimidated by your suggestion. I’ll rashly suggest that I could tackle it successfully, but I think I’ll keep it in reserve while I try the other suggestion first.  
@PetrKryslUCSD Some of the constituencies are islands and so do not have any common boundaries with other constituencies. Also, I’m not sure if the boundaries are sufficiently accurately drawn to guarantee that they share vertices along common boundaries (although they probably do).

Thanks again for your suggestions.

---

<div class="post-metadata">

**Author:** ![cshen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cshen/32/217287_2.png) [@cshen](https://discourse.julialang.org/u/cshen)\
**Post date:** [November 7, 2024, 9:57am UTC](https://discourse.julialang.org/t/find-nearest-neighbours-among-boundary-sets/122342/6 "2024-11-07T09:57:02Z")

</div>

Have you tried computing the nearest neighbour directly from the points on the border of your desired constituency instead?  
You have a geometry, you can most probably just define points along the perimeter, and then for each of those points compute the nearest neighbour between the centroids of all other geometries. There will be a lot of duplication but then you can just start with an empty `set` and it will take care of itself.  
It should also be quite stable in picking up various islands since there is no need for the borders to intersect.

Same step as Imiq suggested, you build the BallTree for the centroids (exclude the centroid of your constituency otherwise that would add some complications in trying to filter that one out of the results), and then you iterate on the border points instead.

I am not very sure if you can do this in Julia (I mean, you can, but I don’t know how easily), but from my experience it is something that GIS libraries have at hand quite easily from my experience with QGIS.  
For example some kind of `points along geometry` function or a multistep `polygon to lines -> points along lines ` kind of thing.  
This depends a lot on what kind of data you have at hand though.

---

<div class="post-metadata">

**Author:** ![TimG](https://avatars.discourse-cdn.com/v4/letter/t/82dd89/32.png) [@TimG](https://discourse.julialang.org/u/TimG)\
**Post date:** [November 7, 2024, 10:58am UTC](https://discourse.julialang.org/t/find-nearest-neighbours-among-boundary-sets/122342/7 "2024-11-07T10:58:53Z")

</div>

I can’t persuade `BallTree` to accept a vector of `unitful Points`.

BallTree won’t accept a vector of centroids:

```julia
ERROR: LoadError: MethodError: no method matching vertices(::Vector{Point{𝔼{2}, TransverseMercator{0.9996012717, Quantity{Float64, Unitful.Dimensions{()}(), 
Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}(49.0), OSGB36, CoordRefSystems.Shift{Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}, Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}}(Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}(-2.0), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}(400000.0), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}(-100000.0)), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}}}})
The function `vertices` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  vertices(::SimpleMesh)
   @ Meshes C:\Users\TGebbels\.julia\packages\Meshes\R0qPK\src\domains\meshes\simplemesh.jl:48
  vertices(::Mesh)
   @ Meshes C:\Users\TGebbels\.julia\packages\Meshes\R0qPK\src\domains\meshes.jl:26
  vertices(::Topology)
   @ Meshes C:\Users\TGebbels\.julia\packages\Meshes\R0qPK\src\topologies.jl:30
  ...

Stacktrace:
  [1] balltree(g1::Vector{Point{𝔼{2}, TransverseMercator{0.9996012717, Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1 
),), Unitful.Dimensions{()}(), nothing}}(49.0), OSGB36, CoordRefSystems.Shift{Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}, Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}}(Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}(-2.0), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}(400000.0), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}(-100000.0)), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}}}})
    @ Main c:\Users\TGebbels\...\Documents\Public Affairs\LAinRegions.jl:135
  [2] macro expansion
    @ .\timing.jl:581 [inlined]
  [3] top-level scope
    @ c:\Users\TGebbels\...\Documents\Public Affairs\LAinRegions.jl:315
  [4] include(fname::String)
    @ Main .\sysimg.jl:38
  [5] run(debug_session::VSCodeDebugger.DebugAdapter.DebugSession, error_handler::VSCodeDebugger.var"#3#4"{String})
    @ VSCodeDebugger.DebugAdapter c:\Users\TGebbels\.vscode\extensions\julialang.language-julia-1.127.2\scripts\packages\DebugAdapter\src\packagedef.jl:122
  [6] startdebugger()
    @ VSCodeDebugger c:\Users\TGebbels\.vscode\extensions\julialang.language-julia-1.127.2\scripts\packages\VSCodeDebugger\src\VSCodeDebugger.jl:45
  [7] top-level scope
    @ c:\Users\TGebbels\.vscode\extensions\julialang.language-julia-1.127.2\scripts\debugger\run_debugger.jl:12
  [8] include(mod::Module, _path::String)
    @ Base .\Base.jl:557
  [9] exec_options(opts::Base.JLOptions)
    @ Base .\client.jl:323
 [10] _start()
    @ Base .\client.jl:531
in expression starting at c:\Users\TGebbels\...\Documents\Public Affairs\LAinRegions.jl:179

```

And it won’t accept a circular vector of vertices from a boundary set:

```julia
ERROR: LoadError: MethodError: no method matching BallTree(::CircularArrays.CircularVector{Point{𝔼{2}, TransverseMercator{0.9996012717, Quantity{Float64, Uni
tful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}(49.0), OSGB36, CoordRefSystems.Shift{Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}, Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}}(Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}(-2.0), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}(400000.0), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}(-100000.0)), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}}}, Vector{Point{𝔼{2}, TransverseMercator{0.9996012717, Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensi 
ons{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}(49.0), OSGB36, CoordRefSystems.Shift{Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}, Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}}(Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}(-2.0), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}(400000.0), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}(-100000.0)), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}}}}}, ::Euclidean)
The type `BallTree` exists, but no method is defined for this combination of argument types when trying to construct it.

Closest candidates are:
  BallTree(::AbstractVector{V}, ::Metric; leafsize, reorder, storedata, reorderbuffer) where V<:AbstractArray
   @ NearestNeighbors C:\Users\TGebbels\.julia\packages\NearestNeighbors\cgW1F\src\ball_tree.jl:21
  BallTree(::AbstractVecOrMat{T}, ::Metric; leafsize, storedata, reorder, reorderbuffer) where T<:AbstractFloat
   @ NearestNeighbors C:\Users\TGebbels\.julia\packages\NearestNeighbors\cgW1F\src\ball_tree.jl:71
  BallTree(::Vector{V}, ::Array{NearestNeighbors.HyperSphere{N, T}, 1}, ::Vector{Int64}, ::M, ::NearestNeighbors.TreeData, ::Bool) where {V<:(AbstractVector), N, T, M<:Metric}
   @ NearestNeighbors C:\Users\TGebbels\.julia\packages\NearestNeighbors\cgW1F\src\ball_tree.jl:7
  ...

Stacktrace:
  [1] balltree(g1::PolyArea{𝔼{2}, TransverseMercator{0.9996012717, Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1),), 
 Unitful.Dimensions{()}(), nothing}}(49.0), OSGB36, CoordRefSystems.Shift{Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}, Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}}(Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}(-2.0), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}(400000.0), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}(-100000.0)), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}}, Ring{𝔼{2}, TransverseMercator{0.9996012717, Quantity{Float64, Un 
itful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}(49.0), OSGB36, CoordRefSystems.Shift{Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}, Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}}(Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}(-2.0), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}(400000.0), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}(-100000.0)), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}}, CircularArrays.CircularVector{Point{𝔼{2}, TransverseMercator{0.9996012717, Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{: 
Degree, Unitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}(49.0), OSGB36, CoordRefSystems.Shift{Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}, Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}}(Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}(-2.0), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}(400000.0), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}(-100000.0)), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}}}, Vector{Point{𝔼{2}, 
 TransverseMercator{0.9996012717, Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}(49.0), OSGB36, CoordRefSystems.Shift{Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}, Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}}(Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}(-2.0), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}(400000.0), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}(-100000.0)), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}}}}}}, Vector{Ring{𝔼{2}, TransverseMercator{0.9996012717, Quantity{Float64, Unitful.Dimensions{()}(), 
 Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}(49.0), OSGB36, CoordRefSystems.Shift{Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}, Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}}(Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}(-2.0), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}(400000.0), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}(-100000.0)), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}}, CircularArrays.CircularVector{Point{𝔼{2}, TransverseMercator{0.9996012717, Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensi 
ons{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}(49.0), OSGB36, CoordRefSystems.Shift{Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}, Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}}(Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}(-2.0), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}(400000.0), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}(-100000.0)), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}}}, Vector{Point{𝔼{2}, TransverseMercator{0.9 
996012717, Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}(49.0), OSGB36, CoordRefSystems.Shift{Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}, Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}}(Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}(-2.0), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}(400000.0), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}(-100000.0)), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}}}}}}}})
    @ Main c:\Users\TGebbels\...\Documents\Public Affairs\LAinRegions.jl:135
  [2] _broadcast_getindex_evalf
    @ .\broadcast.jl:673 [inlined]
  [3] _broadcast_getindex
    @ .\broadcast.jl:646 [inlined]
  [4] getindex
    @ .\broadcast.jl:605 [inlined]
  [5] copy
    @ .\broadcast.jl:906 [inlined]
  [6] materialize(bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Nothing, typeof(balltree), Tuple{Vector{Geometry{𝔼{2}, TransverseMercator{0.9996012717, Quantity{Float6
4, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}(49.0), OSGB36, CoordRefSystems.Shift{Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}, Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}}(Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}(-2.0), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}(400000.0), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}(-100000.0)), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}}}}}})
    @ Base.Broadcast .\broadcast.jl:867
  [7] macro expansion
    @ .\timing.jl:581 [inlined]
  [8] top-level scope
    @ c:\Users\TGebbels\...\Documents\Public Affairs\LAinRegions.jl:315
  [9] include(fname::String)
    @ Main .\sysimg.jl:38
 [10] run(debug_session::VSCodeDebugger.DebugAdapter.DebugSession, error_handler::VSCodeDebugger.var"#3#4"{String})
    @ VSCodeDebugger.DebugAdapter c:\Users\TGebbels\.vscode\extensions\julialang.language-julia-1.127.2\scripts\packages\DebugAdapter\src\packagedef.jl:122
 [11] startdebugger()
    @ VSCodeDebugger c:\Users\TGebbels\.vscode\extensions\julialang.language-julia-1.127.2\scripts\packages\VSCodeDebugger\src\VSCodeDebugger.jl:45
 [12] top-level scope
    @ c:\Users\TGebbels\.vscode\extensions\julialang.language-julia-1.127.2\scripts\debugger\run_debugger.jl:12
 [13] include(mod::Module, _path::String)
    @ Base .\Base.jl:557
 [14] exec_options(opts::Base.JLOptions)
    @ Base .\client.jl:323
 [15] _start()
    @ Base .\client.jl:531
in expression starting at c:\Users\TGebbels\...\Documents\Public Affairs\LAinRegions.jl:179

```

I guess I need to make a conversion but I’m not sure how…

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [November 7, 2024, 11:05am UTC](https://discourse.julialang.org/t/find-nearest-neighbours-among-boundary-sets/122342/8 "2024-11-07T11:05:25Z")

</div>

> [@TimG](#):
>
> I can’t persuade `BallTree` to accept a vector of `unitful Points`.

You don’t need to. The `search` methods already do the necessary internal steps for you with unitful coordinates. The `BallTree` from NearestNeighbors.jl is constructed when you use the `BallSearch` method and the distance is not a Minkowski metric:

```julia
search(point, BallSearch(...))

```

---

<div class="post-metadata">

**Author:** ![pitsianis](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pitsianis/32/26588_2.png) [@pitsianis](https://discourse.julialang.org/u/pitsianis)\
**Post date:** [November 7, 2024, 11:07am UTC](https://discourse.julialang.org/t/find-nearest-neighbours-among-boundary-sets/122342/9 "2024-11-07T11:07:07Z")

</div>

> [@TimG](#):
>
> constituencies

You have a sequence of points that define the boundary of each constituency.

What happens when two constituencies share a common boundary? Are these points identical or do they randomly interpolate the common border? Most probably the second.

What I would do is the following:

Compute the centroid for each constituency, and the maximum distance from it. This will define a disk representing the constituency.

If two disks do not overlap, it should be impossible to share a common border. This will avoid the all-to-all points comparisons.

Bin the centroids in tiles with side length the max of the radii. Therefore, you only need to check for disk overlap if two tiles are same or adjacent.

If two disks overlap, find the minimum distance of the polygons they represent (the boundary of each constituency). Whenever the distance is zero, place a one into an sparse adjacency matrix.

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [November 7, 2024, 11:07am UTC](https://discourse.julialang.org/t/find-nearest-neighbours-among-boundary-sets/122342/10 "2024-11-07T11:07:20Z")

</div>

> [@cshen](#):
>
> For example some kind of `points along geometry` function or a multistep `polygon to lines -> points along lines ` kind of thing.

You can sample points on the geometry or on the boundary with

```julia
sample(geometry, method)
sample(boundary(geometry), method)

```

The list of methods is documented in the Meshes.jl documentation. `HomogeneousSampling` is probably a good choice.

---

<div class="post-metadata">

**Author:** ![TimG](https://avatars.discourse-cdn.com/v4/letter/t/82dd89/32.png) [@TimG](https://discourse.julialang.org/u/TimG)\
**Post date:** [November 7, 2024, 11:17am UTC](https://discourse.julialang.org/t/find-nearest-neighbours-among-boundary-sets/122342/11 "2024-11-07T11:17:50Z")

</div>

> [@juliohm](#):
>
> The `BallTree` from [NearestNeighbors.jl](https://juliahub.com/ui/Packages/General/NearestNeighbors) is constructed when you use the `BallSearch` method

I knew this to be true but “scoured” the docs inadequately and couldn’t find it!  
Thanks!

---

<div class="post-metadata">

**Author:** ![foobar\_lv2](https://avatars.discourse-cdn.com/v4/letter/f/ee59a6/32.png) [@foobar\_lv2](https://discourse.julialang.org/u/foobar_lv2)\
**Post date:** [November 7, 2024, 12:07pm UTC](https://discourse.julialang.org/t/find-nearest-neighbours-among-boundary-sets/122342/12 "2024-11-07T12:07:22Z")

</div>

> [@TimG](#):
>
> However, I want to find nearest neighbours preferring immediately adjacent constituencies over those whose centroids are closer but that are not immediately adjacent.

If you care about adjacent districts then you should probably start with a graph of which constituencies are adjacent.

Ideally you don’t compute that (in e.g. the sense of “share a boundary”), but instead take it as input – politics and maps are weird, and adjacency may be more of a legal / historical category than a geometric one.

“These two constituencies share a geometric boundary on the map, but the boundary is over a river and there is neither bridge nor ferry crossing along it, so they are not adjacent for your purpose.”

“These two districts may not share a geometric boundary, but there is a direct river ferry, so they count as adjacent for the sake of postal delivery.”

“These two communities used to share a boundary pre 16xy, and their legal status as neighbors has been grandfathered in”

---

<div class="post-metadata">

**Author:** ![TimG](https://avatars.discourse-cdn.com/v4/letter/t/82dd89/32.png) [@TimG](https://discourse.julialang.org/u/TimG)\
**Post date:** [November 11, 2024, 11:55am UTC](https://discourse.julialang.org/t/find-nearest-neighbours-among-boundary-sets/122342/13 "2024-11-11T11:55:43Z")

</div>

Thank you all again for your various suggestions. The value of this forum is immense! ❤

Here is an update on how I’ve got on with your various suggestions.

I found @pitsianis’ suggestions both clear and seemingly straightforward. So I started by creating the disks and tiles he recommended. However, electoral constituencies aim for similar populations, and population density varies considerably in the UK (as it does everywhere, I guess). The constituencies vary considerably in size as a result, spanning almost two orders of magnitude (from 2.3km at the smallest to 167km at the largest).

I used a tile size of 150km, resulting in a 5x8 grid covering the whole UK. The number of constituencies in each bin looks like this:

```julia
  0 0 0 15 0 0 1 0
 17 10 4 4 28 2 1 0
 23 84 136 35 19 5 0 1
 32 168 38 7 0 0 0 0
  3 11 6 0 0 0 0 0

```

(North to the right, west to the top!)

So binning, though relatively fast, didn’t reduce the scope of the point-by-point boundary searches very much.

However, with the disks calculated, it becomes straightforward and relatively fast to check constituencies pairwise for possible adjacency based on the sum of the radii in comparison with the distance between centroids. Such a simple comparison eliminates the vast majority of constituency pairs from more detailed checking and is at least as fast as the binning method.

Since the point-by-point checking of boundaries is so costly (see below), I also added a second check for possible adjacency by checking the closest point of one boundary was inside the radius of the other. This is slightly more stringent and more costly than the first check, but was marginally worth it, saving 10% in overall time.

The final test of adjacency I used was to make a point-by-point search to find the nearest pair of points in any two boundary sets that passed the previous tests.

The boundary sets I’ve been using are from [here](https://geoportal.statistics.gov.uk/search?sort=Title%7Ctitle%7Casc&tags=July_2024_PCON_Boundaries). ONS publish several versions at different resolutions. Using the “ultra-low resolution” BUC files, the entire search runs fast (1.7s) but the resolution of the maps is too poor for small, urban constituencies. On the other hand, using the “full resolution” BFC files, the resolution is great, even for the smallest constituency, but the timing is terrible (18,373s). Allocations are also huge with this dataset. Both timing and allocations are completely dominated by a few, difficult constituencies - namely those with long, complex boundaries - usually coastal. Here are two examples of the time taken to complete this point-by-point search:

```julia
 ID Constituency Adj_count
 51 E14001113 3 ... 0.226660 seconds (3.34 k allocations: 53.766 KiB)
563 S14000027 0 ... 5305.311048 seconds (31.42 M allocations: 11.984 TiB, 13.00% gc time)

```

The first constituency, Bootle, is a small urban constituency in Liverpool. The second is the collection of islands off the west coast of Scotland known as Na h-Eileanan an Iar (the Outer Hebrides).

The [GeoStats book](https://juliaearth.github.io/geospatial-data-science-with-julia/01-geodata.html#remarks-1) says:

> - A `GeoTable` provides all the data access patterns of a `DataFrame`. In particular, it supports the syntax `geotable[rows,cols]` without expensive copies of geometries.

Given the scale of the allocations I’m seeing, I wonder if I am doing something wrong? I tried using `@view` for example, but couldn’t get it to work on a `geotable`. I tried using `DataFrame`-style slicing, but switched to using `domain()` for type stability (as recommended in the book) but it didn’t make any difference.

I have two ideas for speeding things up:

- Identify the point in each boundary that is closest to the centroid of the other and start the point-by-point search with this pair, moving alternately clockwise and anti-clockwise out from these starting points. Presumably, the point of contact between the boundaries will be close to this pair, so ought to be found quicker than just starting at what ever the first defined point in the boundary set is. This won’t help with any constituencies that might touch but actually don’t.
- Only check points that lie in the region where the two constituency disks overlap. This is the “region of possibility” outside of which it is not possible for the constituencies to touch. I don’t know yet how to identify these points, nor whether such a test would be quick enough to be worthwhile.

Finally, I create an adjacency matrix, as @pitsianis suggested, which I save to a file. This has the overwhelming advantage of removing the need to run this process more than once. I can simply use any column of the adjacency matrix as a boolean index to select boundaries to put on a map. Here is my example from the OP, showing only those constituencies directly adjacent to the constituency of interest:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/c/0/c0c523cfb9105a3a42a28a690af745ee31e367b5.png)

So, in conclusion, I now have something that works and that I could apply more generally. I will need to do something similar for local authority areas, and for the constituency areas in the devolved assemblies of Northern Ireland, Scotland and Wales.

> **Full code**
>
> ```julia
> using CSV
> using DataFrames
> using GeoIO
> using GeoStats
> using LinearAlgebra
> 
> function finddisk(g) # Find disk to represent a geometry
> centre = centroid(g) # disk centre is at the centroid
> radius = 0.0u"m"
> for p in vertices(g) # radius of the disk is maximum distance from centre to any point on boundary
> radius = max(radius, evaluate(Euclidean(), centre, p))
> end
> return centre, radius
> end
> 
> function closest_distance(c, g) # find the closest point in a geometry (g) to a given point (c)
> first = true
> let dist
> for p in vertices(g)
> if first
> dist = evaluate(Euclidean(), c, p)
> first = false
> else
> dist = min(dist, evaluate(Euclidean(), c, p))
> end
> end
> return dist
> end
> end
> 
> function mighttouch(c1, r1, c2, r2, g1, g2, thold) # Check if it is possible the two boundaries might touch (within threshold)?
> if r1 + r2 + thold >= evaluate(Euclidean(), c1, c2) # Is the distance between disk centres closer than the sum of the disk radii?
> return r1 + thold >= closest_distance(c1, g2) || r2 + thold >= closest_distance(c2, g1) # Is closest point on other boundary withing radius?
> else
> return false
> end
> end
> 
> function mindist(g1, g2, thold) # Are two constituencies closer than `thold` at their closest point?
> let dist
> first=true
> for p1 in vertices(g1) # compare every point in one boundary...
> # inds = search(p1, KNearestSearch(vertices(g2), 1)) # <-- This method produced a stack overflow error
> # p2 = vertices(g2)[inds][1] # <--
> for p2 in vertices(g2) # ... with every point in the other, until ...
> if first
> dist = evaluate(Euclidean(), p1, p2)
> first=false
> else
> dist = min(dist, evaluate(Euclidean(), p1, p2))
> end
> if dist < thold # ... two points closer than `thold` are found.
> break
> end
> end
> if dist < thold # Need to break out of the outer loop too(?).
> break
> end
> end
> return dist
> end
> end
> 
> function createUKtiles(con_bounds) # assign each constituency to a grid of tiles/bins
> 
> disks=DataFrame(idx=Int[], concode=String[], centre=Point[], radius=Quantity[], binx=Int[], biny=Int[])
> 
> # using 150km tiles, assign each disk to a tile/bin
> for id in 1:length(con_bounds.con_code)
> centre, radius = finddisk(con_bounds[id, :geometry]) # obtain disk for each constituency
> binx = div(coords(centre).x, 150000u"m")+1 
> biny = div(coords(centre).y, 150000u"m")+1
> push!(disks, [con_bounds[id, :idx], con_bounds[id, :con_code], centre, radius, binx, biny])
> end
> 
> # How many constituencies are in each bin?
> bins = zeros(Int, maximum(disks.binx), maximum(disks.biny))
> for b in eachrow(disks)
> bins[b.binx, b.biny] += 1
> end
> display(bins)
> 
> # Five biggest and five smallest radii
> println(sort(disks, :radius, rev=true)[1:5, [:radius]]) # actual max radius is ~168km
> println(sort(disks, :radius)[1:5, [:radius]]) # actual min radius is ~2.3km
> 
> return disks
> end
> 
> function makeadjmatrix(con_bounds, disks, thold) # Create an adjacency matrix
> adjmat = falses(nrow(disks), nrow(disks)) 
> for (i, subj) in enumerate(eachrow(@view disks[begin:end-1, :])) # compare each constituency ... 
> print(rpad(i, 4), subj.concode, " ")
> # g1 = con_bounds[i, :geometry]
> g1 = domain(con_bounds)[i]
> adjmat[i, i] = true # each constituency considered adjacent to itself
> @time begin
> nbrs=0
> for (j, obj) in enumerate(eachrow(@view disks[i+1:end, :])) # ... with each subsequest one in the set
> # g2 = con_bounds[i+j, :geometry]
> g2 = domain(con_bounds)[i+j]
> adjmat[j+i, i] = (mighttouch(subj.centre, subj.radius, obj.centre, obj.radius, g1, g2, thold) && # Are they close enough to touch?
> mindist(g1, g2, thold) < thold # If so, check point by point to see if they are adjacent according to specified distance threshold.
> )
> if adjmat[j+i, i]
> nbrs += 1 # an indicative count of neighbours. Will not include count of any previous constituencies that already countied this one as their neighbour.
> end
> adjmat[i, i+j] = copy(adjmat[j+i, i]) # make the matrix symmetrical - neighbour is associative.
> end
> print(rpad(nbrs, 3), " ... ")
> end
> end
> adjmat[end, end] = true
> println("adjmat matrix symmetric: ",issymmetric(adjmat)) # Check to make sure!
> return adjmat
> end
> 
> function readconbounds() # read boundary sets for UK constituencies
> fname = "Westminster_Parliamentary_Constituencies_July_2024_Boundaries_UK_BFC"
> # Shape file from ONS
> print("Reading shapefile ... ")
> @time con_bounds = GeoIO.load(fname) |> Rename("PCON24CD" => "con_code")
> # can't assign using con_bounds.idx = 1:length(con_bounds.con_code)
> con = con_bounds.con_code
> ndf = DataFrame(:con_code => con, :idx => 1:length(con)) # idx is a precaution in case I need to re-order
> return tablejoin(con_bounds, ndf, on=:con_code)
> end
> 
> function processConNeighbours()
> con_bounds = readconbounds() # read boundary sets for UK constituencies
> print("binning constituencies ... ")
> @time disks = createUKtiles(con_bounds) # assign each constituency to a grid of tiles/bins
> println("creating adjacency matrix ... ")
> @time adjmat = makeadjmatrix(con_bounds, disks, 100u"m") # Create an adjacency matrix. I use the disks but not the bins.
> 
> mat = DataFrame(adjmat, disks.concode)
> insertcols!(mat, 1, :idx => disks.idx, :concode => disks.concode)
> CSV.write("dummy_mat.csv", mat)
> return nothing
> end
> processConNeighbours()
> 
> ```

If you have any thoughts on how I can make further improvements, please let me know!

---

<div class="post-metadata">

**Author:** ![cshen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cshen/32/217287_2.png) [@cshen](https://discourse.julialang.org/u/cshen)\
**Post date:** [November 11, 2024, 12:50pm UTC](https://discourse.julialang.org/t/find-nearest-neighbours-among-boundary-sets/122342/14 "2024-11-11T12:50:54Z")

</div>

That

```julia
for p1 in vertices(g1)
    for p2 in vertices(g2) 
        ...
    end
end

```

in the `mindist` function really smells wrt to allocations. I am not sure exactly how the `vertices` function is implemented or if it allocates (I suspect it does and returns an array of ids), but you are computing the same ( ~~very expensive~~ edit: see juliohm reply) thing `length(vertices(g1))` times.

```julia
vertices_g1 = vertices(g1)
vertices_g2 = vertices(g2)
for p1 in vertices_g1
    for p2 in vertices_g2
        ...
    end
end

```

may be better

Also, for NN-search, using the NearestNeighbours.jl package you would only need to instantiate a `BallTree` for the vertices of `g1` only once, and perform the search over all the vertices of `g2`.

```julia
tree = BallTree(vertices(g1), Euclidean())
# find the closest vertex of g1 to every 
# vertex of g2 in a single pass.
inds, dists = knn(tree, vertices(g2), 1, true) 
if minimum(dists) < thold
    ... do something ...
end

```

(assuming that `vertices()` returns a matrix/vector of points)

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [November 11, 2024, 1:06pm UTC](https://discourse.julialang.org/t/find-nearest-neighbours-among-boundary-sets/122342/15 "2024-11-11T13:06:22Z")

</div>

> [@TimG](#):
>
> Given the scale of the allocations I’m seeing, I wonder if I am doing something wrong?

If you can pin point the exact line of code that is allocating, we can investigate it further. Please feel free to share a specific function from the GeoStats.jl stack that you wish allocated less memory than it does currently, and we will take a look.

> [@cshen](#):
>
> I am not sure exactly how the `vertices` function is implemented or if it allocates

It doesn’t allocate in general. It returns the already allocated vertices of the geometry. However, the implementation for `Multi` is allocating inadequately. We will fix it soon:

> <https://github.com/JuliaGeometry/Meshes.jl/issues/1127>
>
> The current implementation currently allocates:
> 
> https://github.com/JuliaGeome…try/Meshes.jl/blob/42f8543991b0da3442cbe7accebc89a2ac558b0e/src/geometries/multigeom.jl#L53
> 
> We need a quick fix for that.

---

<div class="post-metadata">

**Author:** ![TimG](https://avatars.discourse-cdn.com/v4/letter/t/82dd89/32.png) [@TimG](https://discourse.julialang.org/u/TimG)\
**Post date:** [November 11, 2024, 2:45pm UTC](https://discourse.julialang.org/t/find-nearest-neighbours-among-boundary-sets/122342/16 "2024-11-11T14:45:09Z")

</div>

> [@cshen](#):
>
> ```julia
> vertices_g1 = vertices(g1)
> vertices_g2 = vertices(g2)
> for p1 in vertices_g1
> for p2 in vertices_g2
> ...
> end
> end
> 
> ```
> 
> may be better

I put this change in, and have restarted the task. It hasn’t finished yet (of course), but here are some comparisons from before and after:  
Before:

```julia
91 E14001153 5 ... 1.591338 seconds (3.19 k allocations: 53.078 KiB)
92 E14001154 3 ... 29.146777 seconds (1.25 M allocations: 57.522 GiB, 8.46% gc time)
93 E14001155 8 ... 203.594242 seconds (3.79 M allocations: 344.139 GiB, 9.65% gc time)
94 E14001156 4 ... 20.421569 seconds (603.53 k allocations: 32.684 GiB, 6.73% gc time)
95 E14001157 6 ... 3.199758 seconds (166.26 k allocations: 3.979 GiB, 5.58% gc time)

```

After:

```julia
91 E14001153 5 ... 1.778265 seconds (3.67 k allocations: 57.406 KiB)
92 E14001154 3 ... 15.964418 seconds (3.93 k allocations: 11.558 MiB, 0.01% gc time)
93 E14001155 8 ... 129.740393 seconds (3.91 k allocations: 28.531 MiB, 0.00% gc time)
94 E14001156 4 ... 6.814589 seconds (3.81 k allocations: 10.223 MiB, 0.04% gc time)
95 E14001157 6 ... 2.603089 seconds (3.87 k allocations: 5.190 MiB)

```

Before:

```julia
323 E14001385 3 ... 562.346507 seconds (4.90 M allocations: 910.051 GiB, 9.03% gc time)
324 E14001386 5 ... 55.998153 seconds (2.02 k allocations: 34.875 KiB)
325 E14001387 2 ... 23.120379 seconds (433.84 k allocations: 37.172 GiB, 4.59% gc time)
326 E14001388 4 ... 10.667701 seconds (341.59 k allocations: 14.829 GiB, 3.84% gc time)
327 E14001389 1 ... 2.421155 seconds (1.93 k allocations: 273.922 KiB)

```

After:

```julia
323 E14001385 3 ... 345.167239 seconds (2.64 k allocations: 50.923 MiB, 0.00% gc time)
324 E14001386 5 ... 46.530364 seconds (2.27 k allocations: 35.516 KiB)
325 E14001387 2 ... 16.632158 seconds (2.28 k allocations: 3.207 MiB)
326 E14001388 4 ... 5.148056 seconds (2.31 k allocations: 3.682 MiB, 0.12% gc time)
327 E14001389 1 ... 2.856778 seconds (2.26 k allocations: 521.594 KiB)

```

Allocations and GC time are drastically reduced!

Thanks,  
Tim

---

<div class="post-metadata">

**Author:** ![pitsianis](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pitsianis/32/26588_2.png) [@pitsianis](https://discourse.julialang.org/u/pitsianis)\
**Post date:** [November 11, 2024, 3:00pm UTC](https://discourse.julialang.org/t/find-nearest-neighbours-among-boundary-sets/122342/17 "2024-11-11T15:00:26Z")

</div>

> [@TimG](#):
>
> If you have any thoughts on how I can make further improvements, please let me know!

There should be some efficient implementation of the distance of two polygons. See if something like this [Readme · PolygonInbounds.jl](https://docs.juliahub.com/General/PolygonInbounds/stable/) helps.

When you are dealing with a large number of points, binning is your only choice. We have a hierarchical adaptive regular binning package ([GitHub - pitsianis/AdaptiveHierarchicalRegularBinning.jl](https://github.com/pitsianis/AdaptiveHierarchicalRegularBinning.jl)) but it will be an overkill for something that you need to run only once (in a while). You are not redistricting with every election, I hope.

Simply whenever two disks overlap, check the points of the one district that fall into the disk of the other district against the corresponding set of points of the second district that fall into the first district disk.

---

<div class="post-metadata">

**Author:** ![TimG](https://avatars.discourse-cdn.com/v4/letter/t/82dd89/32.png) [@TimG](https://discourse.julialang.org/u/TimG)\
**Post date:** [November 11, 2024, 3:14pm UTC](https://discourse.julialang.org/t/find-nearest-neighbours-among-boundary-sets/122342/18 "2024-11-11T15:14:31Z")

</div>

> [@cshen](#):
>
> ```julia
> tree = BallTree(vertices(g1), Euclidean())
> # find the closest vertex of g1 to every 
> # vertex of g2 in a single pass.
> inds, dists = knn(tree, vertices(g2), 1, true) 
> if minimum(dists) < thold
> ... do something ...
> end
> 
> ```
> 
> (assuming that `vertices()` returns a matrix/vector of points)

I did try this before. I couldn’t use NearestNeighbors.jl directly because the geometries are all unitful. As @juliohm says above, the `search()` method handles this internally, but I can’t get it to work. I tried again, but still fail. (As an aside, I find the unitful stack traces pretty intractable - too long even to post here in full!).

 ![image](https://global.discourse-cdn.com/julialang/original/3X/b/a/baa3ecfe1e0611fa92d384c70596812af4ae334c.png)

```julia
function mindist(g1, g2, thold) # Are two constituencies closer than `thold` at their closest point?
    let dist
# first=true
        g1_v = vertices(g1)
        g2_v = vertices(g2)
        searcher = KNearestSearch(g1_v, 1)
        for p2 in g2_v # compare every point in one boundary...
            _, dists = searchdists(p2, searcher) # <-- This method produced a stack overflow error
            dist = min(dists)

            if dist < thold
                break
            end
        end
        return dist
    end
end

```

gives me:

```julia
ERROR: LoadError: StackOverflowError:
Stacktrace:
  [1] _stable_typeof
    @ .\operators.jl:929 [inlined]
  [2] Base.Fix2(f::typeof(isequal), x::Type)
    @ Base .\operators.jl:1140
  [3] isequal
    @ .\operators.jl:1155 [inlined]
  [4] allequal
    @ .\set.jl:650 [inlined]
  [5] GeometrySet(geoms::CircularArrays.CircularVector{Point{𝔼{2}, TransverseMercator{0.9996012717, Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unit 
ful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}(49.0), OSGB36, CoordRefSystems.Shift{Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}, Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}}(Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}(-2.0), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}(400000.0), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}(-100000.0)), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}}}, Vector{Point{𝔼{2}, TransverseM 
ercator{0.9996012717, Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}(49.0), OSGB36, CoordRefSystems.Shift{Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}, Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}}(Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}(-2.0), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}(400000.0), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}(-100000.0)), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}}}}})
    @ Meshes C:\Users\TGebbels\.julia\packages\Meshes\R0qPK\src\domains\sets.jl:29
  [6] GeometrySet(geoms::CircularArrays.CircularVector{Point{𝔼{2}, TransverseMercator{0.9996012717, Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unit
ful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}(49.0), OSGB36, CoordRefSystems.Shift{Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}, Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}}(Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}(-2.0), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}(400000.0), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}(-100000.0)), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}}}, Vector{Point{𝔼{2}, TransverseM 
ercator{0.9996012717, Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}(49.0), OSGB36, CoordRefSystems.Shift{Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}, Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}}(Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}(-2.0), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}(400000.0), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}(-100000.0)), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}}}}}) (repeats 52172 times)
    @ Meshes C:\Users\TGebbels\.julia\packages\Meshes\R0qPK\src\domains\sets.jl:35
  [7] KNearestSearch(geoms::CircularArrays.CircularVector{Point{𝔼{2}, TransverseMercator{0.9996012717, Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, U 
nitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}(49.0), OSGB36, CoordRefSystems.Shift{Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}, Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}}(Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}(-2.0), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}(400000.0), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}(-100000.0)), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}}}, Vector{Point{𝔼{2}, Transver 
seMercator{0.9996012717, Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}(49.0), OSGB36, CoordRefSystems.Shift{Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}, Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}}(Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}(-2.0), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}(400000.0), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}(-100000.0)), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}}}}}, k::Int64; metric::Euclidean)
    @ Meshes C:\Users\TGebbels\.julia\packages\Meshes\R0qPK\src\neighborsearch\knearest.jl:26
  [8] mindist(g1::PolyArea{𝔼{2}, **[... stuff deleted here]**
    @ Main c:\Users\TGebbels\...\Documents\Public Affairs\LAinRegions.jl:193
  [9] macro expansion
    @ c:\Users\TGebbels\...\Documents\Public Affairs\LAinRegions.jl:254 [inlined]
 [10] macro expansion
    @ .\timing.jl:581 [inlined]
 [11] makeadjmatrix(con_bounds::GeoTable{GeometrySet{𝔼{2}, **[... stuff deleted here]**
    @ Main c:\Users\TGebbels\...\Documents\Public Affairs\LAinRegions.jl:249
 [12] macro expansion
    @ .\timing.jl:581 [inlined]
 [13] processConNeighbours()
    @ Main c:\Users\TGebbels\...\Documents\Public Affairs\LAinRegions.jl:286
 [14] macro expansion
    @ .\timing.jl:581 [inlined]
 [15] top-level scope
    @ c:\Users\TGebbels\...\Documents\Public Affairs\LAinRegions.jl:315
 [16] include(fname::String)
    @ Main .\sysimg.jl:38
 [17] run(debug_session::VSCodeDebugger.DebugAdapter.DebugSession, error_handler::VSCodeDebugger.var"#3#4"{String})
    @ VSCodeDebugger.DebugAdapter c:\Users\TGebbels\.vscode\extensions\julialang.language-julia-1.127.2\scripts\packages\DebugAdapter\src\packagedef.jl:122
 [18] startdebugger()
    @ VSCodeDebugger c:\Users\TGebbels\.vscode\extensions\julialang.language-julia-1.127.2\scripts\packages\VSCodeDebugger\src\VSCodeDebugger.jl:45
 [19] include(mod::Module, _path::String)
    @ Base .\Base.jl:557
 [20] exec_options(opts::Base.JLOptions)
    @ Base .\client.jl:323
 [21] _start()
    @ Base .\client.jl:531
in expression starting at c:\Users\TGebbels\...\Documents\Public Affairs\LAinRegions.jl:351
 * Terminal will be reused by tasks, press any key to close it. 

```

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [November 11, 2024, 3:55pm UTC](https://discourse.julialang.org/t/find-nearest-neighbours-among-boundary-sets/122342/19 "2024-11-11T15:55:17Z")

</div>

> [@TimG](#):
>
> As @juliohm says above, the `search()` method handles this internally, but I can’t get it to work.

Below is a MWE:

```julia
ball = MetricBall(10km)

dom = PointSet(rand(Point, 100))

searcher = BallSearch(dom, ball)

search(rand(Point), searcher)

```

We will improve the documentation to make sure this `MetricBall` is not confused with the geometric `Ball`.

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [November 11, 2024, 4:03pm UTC](https://discourse.julialang.org/t/find-nearest-neighbours-among-boundary-sets/122342/20 "2024-11-11T16:03:46Z")

</div>

Done here:

> <https://github.com/JuliaGeometry/Meshes.jl/commit/6630ef71554d8cf85a9321c02de3f58ada770c73>

Also, check the examples of search here:

> **[Neighbor search · Meshes.jl](https://juliageometry.github.io/MeshesDocs/stable/algorithms/neighborsearch/)**
>
> Documentation for Meshes.jl.

[Next page](https://discourse.julialang.org/t/find-nearest-neighbours-among-boundary-sets/122342.md?page=2)