Find nearest neighbours among boundary sets

General Usage
,
21

I’m sorry but I still can’t translate your MWE into my code:

function mindist(g1, g2, thold) # Are two constituencies closer than `thold` at their closest point?
    let dist
        g1_v = vertices(g1)
        g2_v = vertices(g2)

        ball = MetricBall(10u"km")

        searcher = BallSearch(PointSet(g1_v), ball)

        for p2 in g2_v
            _, dists = searchdists(p2, searcher) # <-- This method still produces an error
            dist = min(dists)                        # <--

            if dist < thold
                break
            end
        end
        return dist
    end
end

I’ve tried defining

        searcher = BallSearch(PointSet(g1_v), ball)
        searcher = BallSearch(g1_v, ball)
        searcher = BallSearch(g1), ball)

but they all fail one way or another.

the example above gives:

ERROR: LoadError: MethodError: no method matching searchdists(::Point{𝔼{2}, TransverseMercator{0.9996012717, Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Uni
tful.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}}}}, ::BallSearch{PointSet{𝔼{2}, TransverseMercator{0.9996012717, Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimen 
sions{()}(), 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}}}}, MetricBall{1, Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(3, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}, Nothing, Euclidean}, NearestNeighbors.KDTree{StaticArraysCore.SVector{2, Float64}, Euclidean, Float64, StaticArraysCore.SVector{2, Float64}}})
The function `searchdists` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  searchdists(::Point, ::BoundedNeighborSearchMethod; mask)
   @ Meshes C:\Users\TGebbels\.julia\packages\Meshes\R0qPK\src\neighborsearch.jl:76
22

The problem seems to be with searchdists().
Using search() works fine in the above code.

23

Thanks for reporting. We will take a closer look at searchdists to see if something is missing regarding unitful CRS.

1 Like
24

@TimG , the searchdists function is only available for BoundedNeighborSearchMethod subtypes such as KBallSearch and KNearestSearch. The BallSearch is unbounded in the sense that the number of neighbors can grow indefinitely inside the ball.

If you need the distances during neighbor search, use KBallSearch instead. Otherwise, stick to BallSearch to find all neighbors inside the ball, and then compute distances as a post-processing step.

1 Like
25

Got it! Thanks!

function doestouch(g1, g2, thold) # Are two constituencies closer than `thold` at their closest point?
    let touching

        searcher = BallSearch(PointSet(vertices(g1)), MetricBall(thold))

        for p2 in vertices(g2) # compare every point in g2 boundary...
            inds = search(p2, searcher)
            touching = (length(inds) > 1)
            if touching # Only one touch is enough!
                break
            end
        end
        
        return touching
    end
end

Renamed mindist to doestouch, returning a bool. This is now superfast (<8s)!

For comparison with above, here is my pathological constituency:

563 S14000027 0             ...   0.255018 seconds (804.57 k allocations: 209.722 MiB, 10.31% gc time)

What’s more, I can now see a way to limit the search to the “region of possibility” using a similar approach with a BallSearch() based on the centre and radius of the two disks. [Edit: This actually costs more than it saves in this use case, so not worth it here.]

1 Like
26

@TimG thank you for giving us the chance to improve the documentation and related software stack :juliaheartpulsing: We will submit a patch release with a couple of fixes soon.

1 Like
27

@TimG we’ve introduced a new function eachvertex that you can use in place of vertices to perform iteration on vertices of polytopes and meshes without memory allocation. This can be useful in your use case with nested for loops:

for v1 in eachvertex(g1)
  for v2 in eachvertex(g2)
    ...
  end
end

and can give you an additional speedup.

The function should be available soon:

2 Likes
28

That’s great - thank you!

1 Like
29

Whichever way adjacency is defined, it will be helpful to make sure the calculation would be symmetric. So there would be no case of (A near B) being true while (B near A) is false.
To achieve this consistency, thinking about rounding and floating point and tolerance is required.
It might seem a waste now, but it would make the resulting function more robust and composable.

30

In my code, I use a threshold (thold) to define proximity. If two neighbours are separated by less than this threshold I define them as adjacent. This works well to allow for boundaries on opposite banks of a river, for example. I set my threshold to 500m based on trial and error in London, where I want neighbours separated by the Thames to be adjacent.

Here is the resulting map for a part of central London (around Chelsea and Fulham constituency).

image
image618×512 62.6 KB

When creating my adjacency matrix, I only compute the half below the diagonal. I then copy to make it symmetrical to ensure adjacency is associative and also define each area to be adjacent to itself.

I still need to manually edit the adjacency matrix for island constituencies that are entirely more than 500m offshore, but there are only 4 of these in the UK.

2 Likes
31

I think you really need to consider the fine details of why you need an adjacency matrix. Are Thurrock and Gravesham really adjacent, in view of (non)existing bridges + ferries? To what do you edit the adjacency matrix of the island constituencies? The Humber bridge is longer than 500m, do you really want to consider the two banks non-adjacent?

This doesn’t sound like a job for software, it sounds like a job for human data entry + a paragraph discussing / justifying your choices.

32

If the underlying calculation is not symmetric, this would make the adjacency depend on the order of nodes. It would be a “pleasant” surprise when sorting the input table (say alphabetically) and having all the downstream data change.

33

My use case doesn’t really justify this. This is not an academic study but an expedient way to quickly visualise data. I want to share with a non-technical audience a comparison of area based data in their area and a few near neighbours. This is not a precision mapping exercise nor is it politically particularly sensitive, especially since these kind of area-based comparisons have not previously been done.

If I generate any interest and there are questions, I’ll refine my methods as necessary. If not, and there is no interest in this kind of view, I won’t have lost too much time pursuing a dead end.

2 Likes
34

I haven’t tested this. However, it is clear that all the constituencies surrounding Chelsea and Fulham are shown on the map I posted above. The same is true for Ceredigion and Preseli constituency in the map I posted before. You can see this by inspection. The same is true for the other maps generated.

To @foobar_lv2’s point, the choice of when to assert adjacency across a gap (eg a body of water) is arbitrary and down to my judgement. But for my use case this seems entirely satisfactory. I’m not complacent about this but, really, my maps are not being used for any high purpose!