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