KD trees
Table of contents: (use browser back to return here)
Distance stuff, lots of tree types. Your distribution of Python comes with scipy which has a couple of implementations.
This is just a quick demo of its potential use. And an homage to Canadian university students, for which KD will always have a special meaning as a primary food group. Kraft Dinner for the non-student
Start with some points and you then want to calculate the closest 2 points to form point origin-destination pairs... because it can be done.
Steps.
- Just deal with the coordinates first, leave the attribute baggage to the side for now.
- Decide on the number of points you want to find the 'closest' of. Don't get ridiculous and ask for an origin-destination matrix with a couple of thousand points. Go back to the project design stage or look at scipy.distance.cdist and a new computer.
- Sorting the points by X, then Y coordinate is useful in some situations. An option to do so is provided.
- Building the KDTree is fairly straightforward using scipy.
- decide on the number of points to find
- the returned list of indices will include the origin point itself, so if you want the closest 2 points, then set your query to N = 3. This can be exploited to suck up the x,y values to form origin-destination pairs if you want to form lines, and/or polygons.
- Decide if you want to just pull out the indices of the closest pairs with their distance.
- Optionally, you can produce a structured array, which you can then bring into ArcGIS Pro as a table for use with a couple of ArcToolbox tools to create geometry
- You are done. Do the join thing if you really need the attributes.
The picture:
The code:
So this function just requires a point array of x,y pairs, the number of closest points (N), whether you want to do an x,y sort first and finally, whether you want an output table suitable for use in ArcGIS Pro.
From there, you simply use arcpy.NumPyArrayToTable to produce a gdb featureclass table.
You can them use... XY to Line … to produce line segments, connecting the various origins and destinations as you see fit, or just bask in the creation of an... XY event layer.
Note: lines 32 and 41 can use... cKDTree ...in place of... KDTree ..., if you just need speed for Euclidean calculations.
<SPAN class="keyword token">def</SPAN> <SPAN class="token function">nn_kdtree</SPAN><SPAN class="punctuation token">(</SPAN>a<SPAN class="punctuation token">,</SPAN> N<SPAN class="operator token">=</SPAN><SPAN class="number token">3</SPAN><SPAN class="punctuation token">,</SPAN> sorted<SPAN class="operator token">=</SPAN><SPAN class="token boolean">True</SPAN><SPAN class="punctuation token">,</SPAN> to_tbl<SPAN class="operator token">=</SPAN><SPAN class="token boolean">True</SPAN><SPAN class="punctuation token">)</SPAN><SPAN class="punctuation token">:</SPAN>
<SPAN class="string token">""</SPAN>"Produce the N closest neighbours array <SPAN class="keyword token">with</SPAN> their distances using
scipy<SPAN class="punctuation token">.</SPAN>spatial<SPAN class="punctuation token">.</SPAN>KDTree <SPAN class="keyword token">as</SPAN> an alternative to einsum<SPAN class="punctuation token">.</SPAN>
Parameters<SPAN class="punctuation token">:</SPAN>
<SPAN class="operator token">-</SPAN><SPAN class="operator token">-</SPAN><SPAN class="operator token">-</SPAN><SPAN class="operator token">-</SPAN><SPAN class="operator token">-</SPAN><SPAN class="operator token">-</SPAN><SPAN class="operator token">-</SPAN><SPAN class="operator token">-</SPAN><SPAN class="operator token">-</SPAN><SPAN class="operator token">-</SPAN><SPAN class="operator token">-</SPAN>
a <SPAN class="punctuation token">:</SPAN> array
Assumed to be an array of point objects <SPAN class="keyword token">for</SPAN> which `nearest` <SPAN class="keyword token">is</SPAN> needed<SPAN class="punctuation token">.</SPAN>
N <SPAN class="punctuation token">:</SPAN> integer
Number of neighbors to <SPAN class="keyword token">return</SPAN><SPAN class="punctuation token">.</SPAN> Note<SPAN class="punctuation token">:</SPAN> the point counts <SPAN class="keyword token">as</SPAN> <SPAN class="number token">1</SPAN><SPAN class="punctuation token">,</SPAN> so N<SPAN class="operator token">=</SPAN><SPAN class="number token">3</SPAN>
returns the closest <SPAN class="number token">2</SPAN> points<SPAN class="punctuation token">,</SPAN> plus itself<SPAN class="punctuation token">.</SPAN>
sorted <SPAN class="punctuation token">:</SPAN> boolean
A nice option to facilitate things<SPAN class="punctuation token">.</SPAN> See `xy_sort`<SPAN class="punctuation token">.</SPAN> Its mini<SPAN class="operator token">-</SPAN>version
<SPAN class="keyword token">is</SPAN> included <SPAN class="keyword token">in</SPAN> this function<SPAN class="punctuation token">.</SPAN>
References<SPAN class="punctuation token">:</SPAN>
<SPAN class="operator token">-</SPAN><SPAN class="operator token">-</SPAN><SPAN class="operator token">-</SPAN><SPAN class="operator token">-</SPAN><SPAN class="operator token">-</SPAN><SPAN class="operator token">-</SPAN><SPAN class="operator token">-</SPAN><SPAN class="operator token">-</SPAN><SPAN class="operator token">-</SPAN><SPAN class="operator token">-</SPAN><SPAN class="operator token">-</SPAN>
`<SPAN class="operator token"><</SPAN>https<SPAN class="punctuation token">:</SPAN><SPAN class="operator token">//</SPAN>stackoverflow<SPAN class="punctuation token">.</SPAN>com<SPAN class="operator token">/</SPAN>questions<SPAN class="operator token">/</SPAN><SPAN class="number token">52366421</SPAN><SPAN class="operator token">/</SPAN>how<SPAN class="operator token">-</SPAN>to<SPAN class="operator token">-</SPAN>do<SPAN class="operator token">-</SPAN>n<SPAN class="operator token">-</SPAN>d<SPAN class="operator token">-</SPAN>distance<SPAN class="operator token">-</SPAN>
<SPAN class="operator token">and</SPAN><SPAN class="operator token">-</SPAN>nearest<SPAN class="operator token">-</SPAN>neighbor<SPAN class="operator token">-</SPAN>calculations<SPAN class="operator token">-</SPAN>on<SPAN class="operator token">-</SPAN>numpy<SPAN class="operator token">-</SPAN>arrays<SPAN class="operator token">/</SPAN><SPAN class="number token">52366706</SPAN><SPAN class="comment token">#52366706>`_.</SPAN>
`<SPAN class="operator token"><</SPAN>https<SPAN class="punctuation token">:</SPAN><SPAN class="operator token">//</SPAN>stackoverflow<SPAN class="punctuation token">.</SPAN>com<SPAN class="operator token">/</SPAN>questions<SPAN class="operator token">/</SPAN><SPAN class="number token">6931209</SPAN><SPAN class="operator token">/</SPAN>difference<SPAN class="operator token">-</SPAN>between<SPAN class="operator token">-</SPAN>scipy<SPAN class="operator token">-</SPAN>
spatial<SPAN class="operator token">-</SPAN>kdtree<SPAN class="operator token">-</SPAN><SPAN class="operator token">and</SPAN><SPAN class="operator token">-</SPAN>scipy<SPAN class="operator token">-</SPAN>spatial<SPAN class="operator token">-</SPAN>ckdtree<SPAN class="operator token">/</SPAN><SPAN class="number token">6931317</SPAN><SPAN class="comment token">#6931317>`_.</SPAN>
<SPAN class="string token">"""
def _xy_sort_(a):
"""</SPAN>mini xy_sort<SPAN class="string token">""</SPAN>"
a_view <SPAN class="operator token">=</SPAN> a<SPAN class="punctuation token">.</SPAN>view<SPAN class="punctuation token">(</SPAN>a<SPAN class="punctuation token">.</SPAN>dtype<SPAN class="punctuation token">.</SPAN>descr <SPAN class="operator token">*</SPAN> a<SPAN class="punctuation token">.</SPAN>shape<SPAN class="punctuation token">[</SPAN><SPAN class="number token">1</SPAN><SPAN class="punctuation token">]</SPAN><SPAN class="punctuation token">)</SPAN>
idx <SPAN class="operator token">=</SPAN>np<SPAN class="punctuation token">.</SPAN>argsort<SPAN class="punctuation token">(</SPAN>a_view<SPAN class="punctuation token">,</SPAN> axis<SPAN class="operator token">=</SPAN><SPAN class="number token">0</SPAN><SPAN class="punctuation token">,</SPAN> order<SPAN class="operator token">=</SPAN><SPAN class="punctuation token">(</SPAN>a_view<SPAN class="punctuation token">.</SPAN>dtype<SPAN class="punctuation token">.</SPAN>names<SPAN class="punctuation token">)</SPAN><SPAN class="punctuation token">)</SPAN><SPAN class="punctuation token">.</SPAN>ravel<SPAN class="punctuation token">(</SPAN><SPAN class="punctuation token">)</SPAN>
a <SPAN class="operator token">=</SPAN> np<SPAN class="punctuation token">.</SPAN>ascontiguousarray<SPAN class="punctuation token">(</SPAN>a<SPAN class="punctuation token">[</SPAN>idx<SPAN class="punctuation token">]</SPAN><SPAN class="punctuation token">)</SPAN>
<SPAN class="keyword token">return</SPAN> a<SPAN class="punctuation token">,</SPAN> idx
<SPAN class="comment token">#</SPAN>
<SPAN class="keyword token">def</SPAN> <SPAN class="token function">xy_dist_headers</SPAN><SPAN class="punctuation token">(</SPAN>N<SPAN class="punctuation token">)</SPAN><SPAN class="punctuation token">:</SPAN>
<SPAN class="string token">"""Construct headers for the optional table output"""</SPAN>
vals <SPAN class="operator token">=</SPAN> np<SPAN class="punctuation token">.</SPAN>repeat<SPAN class="punctuation token">(</SPAN>np<SPAN class="punctuation token">.</SPAN>arange<SPAN class="punctuation token">(</SPAN>N<SPAN class="punctuation token">)</SPAN><SPAN class="punctuation token">,</SPAN> <SPAN class="number token">2</SPAN><SPAN class="punctuation token">)</SPAN>
names <SPAN class="operator token">=</SPAN> <SPAN class="punctuation token">[</SPAN><SPAN class="string token">'X_{}'</SPAN><SPAN class="punctuation token">,</SPAN> <SPAN class="string token">'Y_{}'</SPAN><SPAN class="punctuation token">]</SPAN><SPAN class="operator token">*</SPAN>N <SPAN class="operator token">+</SPAN> <SPAN class="punctuation token">[</SPAN><SPAN class="string token">'d_{}'</SPAN><SPAN class="punctuation token">]</SPAN><SPAN class="operator token">*</SPAN><SPAN class="punctuation token">(</SPAN>N<SPAN class="number token">-1</SPAN><SPAN class="punctuation token">)</SPAN>
vals <SPAN class="operator token">=</SPAN> <SPAN class="punctuation token">(</SPAN>np<SPAN class="punctuation token">.</SPAN>repeat<SPAN class="punctuation token">(</SPAN>np<SPAN class="punctuation token">.</SPAN>arange<SPAN class="punctuation token">(</SPAN>N<SPAN class="punctuation token">)</SPAN><SPAN class="punctuation token">,</SPAN> <SPAN class="number token">2</SPAN><SPAN class="punctuation token">)</SPAN><SPAN class="punctuation token">)</SPAN><SPAN class="punctuation token">.</SPAN>tolist<SPAN class="punctuation token">(</SPAN><SPAN class="punctuation token">)</SPAN> <SPAN class="operator token">+</SPAN> <SPAN class="punctuation token">[</SPAN>i <SPAN class="keyword token">for</SPAN> i <SPAN class="keyword token">in</SPAN> range<SPAN class="punctuation token">(</SPAN><SPAN class="number token">1</SPAN><SPAN class="punctuation token">,</SPAN> N<SPAN class="punctuation token">)</SPAN><SPAN class="punctuation token">]</SPAN>
n <SPAN class="operator token">=</SPAN> <SPAN class="punctuation token">[</SPAN>names<SPAN class="punctuation token">[</SPAN>i<SPAN class="punctuation token">]</SPAN><SPAN class="punctuation token">.</SPAN>format<SPAN class="punctuation token">(</SPAN>vals<SPAN class="punctuation token">[</SPAN>i<SPAN class="punctuation token">]</SPAN><SPAN class="punctuation token">)</SPAN> <SPAN class="keyword token">for</SPAN> i <SPAN class="keyword token">in</SPAN> range<SPAN class="punctuation token">(</SPAN>len<SPAN class="punctuation token">(</SPAN>vals<SPAN class="punctuation token">)</SPAN><SPAN class="punctuation token">)</SPAN><SPAN class="punctuation token">]</SPAN>
f <SPAN class="operator token">=</SPAN> <SPAN class="punctuation token">[</SPAN><SPAN class="string token">'<f8'</SPAN><SPAN class="punctuation token">]</SPAN><SPAN class="operator token">*</SPAN>N<SPAN class="operator token">*</SPAN><SPAN class="number token">2</SPAN> <SPAN class="operator token">+</SPAN> <SPAN class="punctuation token">[</SPAN><SPAN class="string token">'<f8'</SPAN><SPAN class="punctuation token">]</SPAN><SPAN class="operator token">*</SPAN><SPAN class="punctuation token">(</SPAN>N<SPAN class="number token">-1</SPAN><SPAN class="punctuation token">)</SPAN>
<SPAN class="keyword token">return</SPAN> list<SPAN class="punctuation token">(</SPAN>zip<SPAN class="punctuation token">(</SPAN>n<SPAN class="punctuation token">,</SPAN>f<SPAN class="punctuation token">)</SPAN><SPAN class="punctuation token">)</SPAN>
<SPAN class="comment token"># </SPAN>
<SPAN class="keyword token">from</SPAN> scipy<SPAN class="punctuation token">.</SPAN>spatial <SPAN class="keyword token">import</SPAN> cKDTree
<SPAN class="comment token">#</SPAN>
idx_orig <SPAN class="operator token">=</SPAN> <SPAN class="punctuation token">[</SPAN><SPAN class="punctuation token">]</SPAN>
<SPAN class="keyword token">if</SPAN> sorted<SPAN class="punctuation token">:</SPAN>
a<SPAN class="punctuation token">,</SPAN> idx_orig <SPAN class="operator token">=</SPAN> _xy_sort_<SPAN class="punctuation token">(</SPAN>a<SPAN class="punctuation token">)</SPAN>
<SPAN class="comment token"># ---- query the tree for the N nearest neighbors and their distance</SPAN>
t <SPAN class="operator token">=</SPAN> cKDTree<SPAN class="punctuation token">(</SPAN>a<SPAN class="punctuation token">)</SPAN>
dists<SPAN class="punctuation token">,</SPAN> indices <SPAN class="operator token">=</SPAN> t<SPAN class="punctuation token">.</SPAN>query<SPAN class="punctuation token">(</SPAN>a<SPAN class="punctuation token">,</SPAN> N<SPAN class="punctuation token">)</SPAN>
<SPAN class="keyword token">if</SPAN> to_tbl<SPAN class="punctuation token">:</SPAN>
dt <SPAN class="operator token">=</SPAN> xy_dist_headers<SPAN class="punctuation token">(</SPAN>N<SPAN class="punctuation token">)</SPAN> <SPAN class="comment token"># --- Format a structured array header</SPAN>
xys <SPAN class="operator token">=</SPAN> a<SPAN class="punctuation token">[</SPAN>indices<SPAN class="punctuation token">]</SPAN>
new_shp <SPAN class="operator token">=</SPAN> <SPAN class="punctuation token">(</SPAN>xys<SPAN class="punctuation token">.</SPAN>shape<SPAN class="punctuation token">[</SPAN><SPAN class="number token">0</SPAN><SPAN class="punctuation token">]</SPAN><SPAN class="punctuation token">,</SPAN> np<SPAN class="punctuation token">.</SPAN>prod<SPAN class="punctuation token">(</SPAN>xys<SPAN class="punctuation token">.</SPAN>shape<SPAN class="punctuation token">[</SPAN><SPAN class="number token">1</SPAN><SPAN class="punctuation token">:</SPAN><SPAN class="punctuation token">]</SPAN><SPAN class="punctuation token">)</SPAN><SPAN class="punctuation token">)</SPAN>
xys <SPAN class="operator token">=</SPAN> xys<SPAN class="punctuation token">.</SPAN>reshape<SPAN class="punctuation token">(</SPAN>new_shp<SPAN class="punctuation token">)</SPAN>
ds <SPAN class="operator token">=</SPAN> dists<SPAN class="punctuation token">[</SPAN><SPAN class="punctuation token">:</SPAN><SPAN class="punctuation token">,</SPAN> <SPAN class="number token">1</SPAN><SPAN class="punctuation token">:</SPAN><SPAN class="punctuation token">]</SPAN> <SPAN class="comment token">#[d[1:] for d in dists]</SPAN>
arr <SPAN class="operator token">=</SPAN> np<SPAN class="punctuation token">.</SPAN>concatenate<SPAN class="punctuation token">(</SPAN><SPAN class="punctuation token">(</SPAN>xys<SPAN class="punctuation token">,</SPAN> ds<SPAN class="punctuation token">)</SPAN><SPAN class="punctuation token">,</SPAN> axis<SPAN class="operator token">=</SPAN><SPAN class="number token">1</SPAN><SPAN class="punctuation token">)</SPAN>
arr <SPAN class="operator token">=</SPAN> arr<SPAN class="punctuation token">.</SPAN>view<SPAN class="punctuation token">(</SPAN>dtype<SPAN class="operator token">=</SPAN>dt<SPAN class="punctuation token">)</SPAN><SPAN class="punctuation token">.</SPAN>squeeze<SPAN class="punctuation token">(</SPAN><SPAN class="punctuation token">)</SPAN>
<SPAN class="keyword token">return</SPAN> arr
<SPAN class="keyword token">else</SPAN><SPAN class="punctuation token">:</SPAN>
<SPAN class="keyword token">return</SPAN> np<SPAN class="punctuation token">.</SPAN>array<SPAN class="punctuation token">(</SPAN>indices<SPAN class="punctuation token">)</SPAN><SPAN class="punctuation token">,</SPAN> idx_orig<SPAN class="line-numbers-rows"><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN></SPAN>The output
Just a slightly better formatting that you can get with one of my numpy functions... obviating the need for Pandas for table niceness.
id X_0 Y_0 X_1 Y_1 X_2 Y_2 d_1 d_2
--------------------------------------------------------------
000 3.00 98.00 10.00 94.00 23.00 94.00 8.06 20.40
001 10.00 94.00 3.00 98.00 23.00 94.00 8.06 13.00
002 13.00 18.00 19.00 22.00 34.00 16.00 7.21 21.10
003 19.00 22.00 13.00 18.00 34.00 16.00 7.21 16.16
004 23.00 94.00 10.00 94.00 3.00 98.00 13.00 20.40
005 34.00 16.00 19.00 22.00 43.00 1.00 16.16 17.49
006 37.00 64.00 43.00 89.00 56.00 84.00 25.71 27.59
007 43.00 1.00 34.00 16.00 66.00 6.00 17.49 23.54
008 43.00 89.00 56.00 84.00 61.00 87.00 13.93 18.11
009 56.00 84.00 61.00 87.00 43.00 89.00 5.83 13.93
010 61.00 87.00 56.00 84.00 43.00 89.00 5.83 18.11
011 66.00 6.00 76.00 20.00 43.00 1.00 17.20 23.54
012 67.00 41.00 78.00 50.00 76.00 20.00 14.21 22.85
013 76.00 20.00 66.00 6.00 67.00 41.00 17.20 22.85
014 78.00 50.00 67.00 41.00 80.00 67.00 14.21 17.12
015 80.00 67.00 91.00 66.00 78.00 50.00 11.05 17.12
016 82.00 91.00 94.00 95.00 61.00 87.00 12.65 21.38
017 91.00 66.00 80.00 67.00 78.00 50.00 11.05 20.62
018 94.00 95.00 82.00 91.00 91.00 66.00 12.65 29.15
019 96.00 40.00 78.00 50.00 91.00 66.00 20.59 26.48<SPAN class="line-numbers-rows"><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN></SPAN>
Summary
Behind the scenes, there should be some 'spatial cleanup'. Specifically, if you look at the image you have point pairs connected by a single segment, that is because they are the closest to one another. Rather than duplicating the segment with opposing directions, you can 'prune' the indices and remove those prior to producing the geometry.
There are lots of tools that you can produce/show geometric relationships. Use them to provide answers to your questions. This implementation will appear soon on the code sharing site. I will provide a link soon.
