Find nearest neighbours among boundary sets

TimG · November 6, 2024, 5:22pm

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

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

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:

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:

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!

juliohm · November 6, 2024, 6:20pm

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:

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.

lmiq · November 6, 2024, 6:40pm

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):

Construct a BallTree for each constituencies.
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.

PetrKryslUCSD · November 6, 2024, 6:53pm

Length of the shared boundary?

TimG · November 7, 2024, 9:11am

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:

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.

cshen · November 7, 2024, 9:57am

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.

TimG · November 7, 2024, 10:58am

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

BallTree won’t accept a vector of centroids:

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:

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…

juliohm · November 7, 2024, 11:05am

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:

search(point, BallSearch(...))

pitsianis · November 7, 2024, 11:07am

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.

juliohm · November 7, 2024, 11:07am

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

sample(geometry, method)
sample(boundary(geometry), method)

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

TimG · November 7, 2024, 11:17am

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

foobar_lv2 · November 7, 2024, 12:07pm

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”

TimG · November 11, 2024, 11:55am

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:

  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. 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:

 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 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:

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

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!

cshen · November 11, 2024, 12:50pm

That

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.

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.

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)

juliohm · November 11, 2024, 1:06pm

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.

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:

TimG · November 11, 2024, 2:45pm

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:

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:

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:

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:

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

pitsianis · November 11, 2024, 3:00pm

There should be some efficient implementation of the distance of two polygons. See if something like this Readme · PolygonInbounds.jl 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) 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.

TimG · November 11, 2024, 3:14pm

cshen:

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!).

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:

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.

juliohm · November 11, 2024, 3:55pm

Below is a MWE:

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.

juliohm · November 11, 2024, 4:03pm

Done here:

Also, check the examples of search here: