Connecting the dots. Spanning trees have a wide variety of applications in GIS... or did. I have been finishing writings on distance and geometry and I almost forgot this one. The principle is simple. Connect the dots, so that no line (edge) overlaps, but only connects at a point. In practice, you want the tree to produce a minimal length/distance output. There are a variety of names forms, each with their subtle nuance and application, but my favorite is Prim's... only because I sortof understand it. I am not sure if this his implementation exactly, but it is close enough. So I post it here for those that want to try it. Besides, Prim's can be used to produce mazes, a far more useful application of dot connection.

My favorite code header to cover imports and array printing and some bare-bones graphing

`import sys`

import numpy as np

import matplotlib.pyplot as plt

from textwrap import dedent, indent

ft = {'bool': lambda x: repr(x.astype('int32')),

'float': '{: 0.1f}'.format}

np.set_printoptions(edgeitems=10, linewidth=100, precision=2,

suppress=True, threshold=120,

formatter=ft)

np.ma.masked_print_option.set_display('-')

script = sys.argv[0]

Now the ending of the script where the actual calls are performed and an explanation of what is happening occurs.

`# ---------------------------------------------------------------------`

if __name__ == "__main__":

"""Main section... """

#print("Script... {}".format(script))

# ---- Take a few points to get you started ----

a = np.array([[0, 0], [0,8], [10, 8], [10,0], [3, 4], [7,4]])

#

idx= np.lexsort((a[:,1], a[:,0])) # sort X, then Y

a_srt = a[idx,:] # slice the sorted array

d = _e_dist(a_srt) # determine the square form distances

pairs = mst(d) # get the orig-dest pairs for the mst

plot_mst(a_srt, pairs) # a little plot

o_d = connect(a_srt, d, pairs) # produce an o-d structured array

Now the rest is just filler. The code defs are given below.

`def _e_dist(a):`

"""Return a 2D square-form euclidean distance matrix. For other

: dimensions, use e_dist in ein_geom.py"""

b = a.reshape(np.prod(a.shape[:-1]), 1, a.shape[-1])

diff = a - b

d = np.sqrt(np.einsum('ijk,ijk->ij', diff, diff)).squeeze()

#d = np.triu(d)

return d

def mst(W, copy_W=True):

"""Determine the minimum spanning tree for a set of points represented

: by their inter-point distances... ie their 'W'eights

:Requires:

:--------

: W - edge weights (distance, time) for a set of points. W needs to be

: a square array or a np.triu perhaps

:Returns:

:-------

: pairs - the pair of nodes that form the edges

"""

if copy_W:

W = W.copy()

if W.shape[0] != W.shape[1]:

raise ValueError("W needs to be square matrix of edge weights")

Np = W.shape[0]

pairs = []

pnts_seen = [0] # Add the first point

n_seen = 1

# exclude self connections by assigning inf to the diagonal

diag = np.arange(Np)

W[diag, diag] = np.inf

#

while n_seen != Np:

new_edge = np.argmin(W[pnts_seen], axis=None)

new_edge = divmod(new_edge, Np)

new_edge = [pnts_seen[new_edge[0]], new_edge[1]]

pairs.append(new_edge)

pnts_seen.append(new_edge[1])

W[pnts_seen, new_edge[1]] = np.inf

W[new_edge[1], pnts_seen] = np.inf

n_seen += 1

return np.vstack(pairs)

def plot_mst(a, pairs):

"""plot minimum spanning tree test """

plt.scatter(a[:, 0], a[:, 1])

ax = plt.axes()

ax.set_aspect('equal')

for pair in pairs:

i, j = pair

plt.plot([a[i, 0], a[j, 0]], [a[i, 1], a[j, 1]], c='r')

lbl = np.arange(len(a))

for label, xpt, ypt in zip(lbl, a[:,0], a[:,1]):

plt.annotate(label, xy=(xpt, ypt), xytext=(2,2), size=8,

textcoords='offset points',

ha='left', va='bottom')

plt.show()

plt.close()

def connect(a, dist_arr, edges):

"""Return the full spanning tree, with points, connections and distance

: a - point array

: dist - distance array, from _e_dist

: edge - edges, from mst

"""

p_f = edges[:, 0]

p_t = edges[:, 1]

d = dist_arr[p_f, p_t]

n = p_f.shape[0]

dt = [('Orig', '<i4'), ('Dest', 'i4'), ('Dist', '<f8')]

out = np.zeros((n,), dtype=dt)

out['Orig'] = p_f

out['Dest'] = p_t

out['Dist'] = d

return out

The output from the sample points is hardly exciting. But you can see the possibilities for the other set.

This one is for a 100 points, with a minimum spacing of 3 within a 100x100 unit square. Sprightly solution even on an iThingy using python 3.5

Now on to maze creation .....

Thanks for sharing, Dan. It probably doesn't have the performance of numpy, but here's my try at simplifying it in terms of arcpy: