Skip navigation
All People > Dan_Patterson > Py... blog > 2017 > January
2017
Dan_Patterson

Spanning trees

Posted by Dan_Patterson Champion Jan 31, 2017

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

Dan_Patterson

Vincenty

Posted by Dan_Patterson Champion Jan 24, 2017

Distance calculations using longitude and latitude on the ellipsoid.   Haversine or Vincenty?

The choice depends on the resolution needed.  This thread on Computing Distances reminded me that I had the latter kicking around and Bruce Harold's Reverse Geocoding uses it.  

 

So I thought I would throw the reference in the Py..Links so I/we don't forget.  

The code can be found on my GitHub link but I have attached a copy to this post.

 

Ispared no effort in code verbosity and I tried to remain faithful to Chris Veness's documentation and presentation of the original code implementation.  His page contains an interactive calculator implemented in Java if you just need a few calculations.  I have included a few sample calculations here should you be interested

 

:Examples:
: long0  lat0  long1  lat1   dist       initial    final  head to
: -75.0, 45.0, -75.0, 46.0   111141.548   0.000,   0.000   N
: -75.0, 46.0, -75.0, 45.0   111141.548 180.000, 180.000   S
: -76.0, 45.0, -75.0, 45.0    78846.334  89.646,  90.353   E
: -75.0, 45.0, -76.0, 45.0    78846.334 270.353, 269.646   W
: -76.0, 46.0, -75.0, 45.0   135869.091 144.526, 145.239   SE
: -75.0, 46.0, -76.0, 45.0   135869.091 215.473, 214.760   SW
: -76.0, 45.0, -75.0, 46.0   135869.091  34.760,  35.473   NE
: -75.0, 45.0, -76.0, 46.0   135869.091 325.239, 324.526   NW
: -90.0,  0.0    0.0   0.0 10018754.171  90.000   90.000   1/4 equator
: -75.0   0.0  -75.0  90.0 10001965.729   0.000    0.000   to N pole
:

 

So use Bruce's if you need to geocode... use Chris's if you just need a few points.... or play with this incarnation should you need to incorporate a function in your own code. I will get around to converting it to NumPy eventually so one can process large sets of origin-destinations that need a tad more than a spherical estimate.  The current version of the program uses iteration which makes it a poor candidate for vectorization on arrays, but there are other implementations one can use.

 

The most recent version of vincenty.py can be found on my github repository numpy_samples/vincenty.py at master · Dan-Patterson/numpy_samples · GitHub 

Generate closest features by distance


From near.py  See references


Emulating Generate Near Table from ArcMap

Let us begin with finding the closest 3 points to every point in a point set.
Well, that is easy... we just use the 'Generate Near Table tool'.

You have been spoiled. This tool is only available at the Advanced license level. You are now working in a job that uses a Standard license... what to do!?

Of course!.... roll out your own.
We will begin with a simple call to 'n_near' in near.py.

We can step through the process...

Begin with array 'a'. Since we are going to use einsum to perform the distance calculations, we need to clone and reshape the array to facilitate the process.

The following array, 'a', represents the 4 corners of a 2x2 unit square, with a centre point. The points are arranged in clockwise order.

>>> a # a.shape => (5, 2)
array([[0, 0],
       [0, 2],
       [2, 2],
       [2, 0],
       [1, 1]], dtype=int32)
The array reshaping is needed in order subtract the arrays.
>>> b = a.reshape(np.prod(a.shape[:-1]), 1, a.shape[-1])
>>> b # b.shape => (5, 1, 2)
array([[[0, 0]],
       [[0, 2]],
       [[2, 2]],
       [[2, 0]],
       [[1, 1]]], dtype=int32)

I have documented the details of the array construction and einsum notation elsewhere. Suffice to say, we can now subtract the two arrays, perform the einsum product summation and finish with the euclidean distance calculation.

The difference array produces 5 blocks of 5x2 values. The summation of the products of these arrays essentially yields the squared distance, from which, euclidean distance is derived. There are other ways of doing this, such as dot product calculations. I prefer einsum methods since it can be scaled up from 1D to n-D unlike most other approaches.

 

>>> diff = b - a # diff.shape => (5, 5, 2)

 

The 'diff' array looks like the following. I took the liberty of using a function in arr_tools (on github) to rearrange the array into a more readable form ( https://github.com/Dan-Patterson/numpy_samples/blob/master/formatting/arr_frmts.py )

>>> import arr_tools as art
>>> art.frmt_(diff)
Array...
-shape (5, 5, 2), ndim 3
. 0 0 0 -2 -2 -2 -2 0 -1 -1
. 0 2 0 0 -2 0 -2 2 -1 1
. 2 2 2 0 0 0 0 2 1 1
. 2 0 2 -2 0 -2 0 0 1 -1
. 1 1 1 -1 -1 -1 -1 1 0 0



The distance calculation is pretty simple, just a bit of einsum notation, get rid of some extraneous dimensions if present and there you have it...

>>> dist = np.einsum('ijk,ijk->ij', diff, diff) # the magic happens...
>>> d = np.sqrt(dist).squeeze() # get rid of extra 'stuff'
>>> d # the distance array...
array([[ 0.0, 2.0, 2.8, 2.0, 1.4],
       [ 2.0, 0.0, 2.0, 2.8, 1.4],
       [ 2.8, 2.0, 0.0, 2.0, 1.4],
       [ 2.0, 2.8, 2.0, 0.0, 1.4],
       [ 1.4, 1.4, 1.4, 1.4, 0.0]])

The result as you can see from the above is a row-column structure much like that derived from scipy's cdist function. Each row and column represents a point, resulting in the diagonal having a distance of zero.

The next step is to get a sorted list of the distances. This is where np.argsort comes into play, since it returns a list of indices that represent the sorted distance values. The indices are used to pull out the coordinates in the appropriate order.

>>> kv = np.argsort(d, axis=1) # sort 'd' on last axis to get keys
>>> kv
array([[0, 4, 1, 3, 2],
       [1, 4, 0, 2, 3],
       [2, 4, 1, 3, 0],
       [3, 4, 0, 2, 1],
       [4, 0, 1, 2, 3]])

>>> coords = a[kv] # for each point, pull out the points in closest order
>>> a[kv].shape # the shape is still not ready for use...
(5, 5, 2)

The coordinate array (coords) needs to be reshaped so that the X, Y pair values can be laid out in row format for final presentation. Each point calculates the distances to itself and the other points, so the array has 5 groupings of 5 pairs of coordinates. This can be reshaped, to produce 5 rows of x, y values using the following.

>>> s0, s1, s2 = coords.shape
>>> coords = coords.reshape((s0, s1*s2)) # the result will be a 2D array...
>>> coords
array([[0, 0, 1, 1, 0, 2, 2, 0, 2, 2],
       [0, 2, 1, 1, 0, 0, 2, 2, 2, 0],
       [2, 2, 1, 1, 0, 2, 2, 0, 0, 0],
       [2, 0, 1, 1, 0, 0, 2, 2, 0, 2],
       [1, 1, 0, 0, 0, 2, 2, 2, 2, 0]], dtype=int32)

Each row represents an input point in the order they were input. Compare input array 'a' with the first two columns of the 'coords' array to confirm. The remaining columns are pairs of the x, y values arranged by their distance sorted order (more about this later).

The distance values are then sorted in ascending order. Obviously, the first value in each list will be the distance of each point to itself (0.0) so it is sliced off leaving the remaining distances.

>>> dist = np.sort(d)[:,1:] # slice sorted distances, skip 1st
>>> dist
array([[ 1.4, 2.0, 2.0, 2.8],
       [ 1.4, 2.0, 2.0, 2.8],
       [ 1.4, 2.0, 2.0, 2.8],
       [ 1.4, 2.0, 2.0, 2.8],
       [ 1.4, 1.4, 1.4, 1.4]])

If you examine the points that were used as input, they formed a rectangle with a point in the middle. It should come as no surprise that the first column represents the distance of each point to the center point (the last row). The next two columns are the distance of each point to its adjacent neighbour while the last column is the distance of each point to its diagonal. The exception is of course the center point (last row) which is equidistant to the other 4 points.

The rest of the code is nothing more that a fancy assemblage of the resultant data into a structure that can be used to output a structured array of coordinates and distances, which can be brought in to ArcMap to form various points or polyline assemblages.


Here are the results from the script...

:-----------------------------------------------------------------
:Closest 2 points for points in an array. Results returned as
: a structured array with coordinates and distance values.
Demonstrate n_near function ....
:Input points... array 'a'
[[0 0]
[0 2]
[2 2]
[2 0]
[1 1]]
:output array
ID     Xo   Yo C0_X C0_Y C1_X C1_Y Dist0 Dist1
0.00 0.00 0.00 1.00 1.00 0.00 2.00 1.41 2.00
1.00 0.00 2.00 1.00 1.00 0.00 0.00 1.41 2.00
2.00 1.00 1.00 0.00 0.00 0.00 2.00 1.41 1.41
3.00 2.00 2.00 1.00 1.00 0.00 2.00 1.41 2.00
4.00 2.00 0.00 1.00 1.00 0.00 0.00 1.41 2.00
:------------------------------------------------------------------
This is the final form of the array.
array([(0, 0.0, 0.0, 1.0, 1.0, 0.0, 2.0, 1.4142135..., 2.0),
       (1, 0.0, 2.0, 1.0, 1.0, 0.0, 0.0, 1.4142135..., 2.0),
       (2, 1.0, 1.0, 0.0, 0.0, 0.0, 2.0, 1.4142135..., 1.4142135...),
       (3, 2.0, 2.0, 1.0, 1.0, 0.0, 2.0, 1.4142135..., 2.0),
       (4, 2.0, 0.0, 1.0, 1.0, 0.0, 0.0, 1.4142135..., 2.0)],
      dtype=[('ID', '<i4'),
             ('Xo', '<f8'), ('Yo', '<f8'),
             ('C0_X', '<f8'), ('C0_Y', '<f8'),
             ('C1_X', '<f8'), ('C1_Y', '<f8'),
             ('Dist0', '<f8'), ('Dist1', '<f8')])

I took the liberty of doing some fiddling with the format to make it easier to read. It should be readily apparent that this array could be used as input to NumPyArrayToFeatureClass so that you can produce a featureclass or shapefile of the data.

That is about all now. There are a variety of ways to perform the same thing... hope this adds to your arsenal of tools.


References:
----------
http://desktop.arcgis.com/en/arcmap/latest/tools/analysis-toolbox/generate-near-table.htm

https://github.com/Dan-Patterson/numpy_samples/blob/master/geometry/scripts/near.py


----------------------------------------------------------------------------

Code posted on my GitHub repository... called circle_make.py, perhaps a bit of a misnomer since it works to create all things associated with 'circular' features, which would include ellipses, triangles, squares, rectangles, sectors, arcs, pentagons, hexagons, octagons and n-gons... anything whose points can be placed on a circle.

 

Here are some pictures, you can examine the code at your leisure.  The functions (def) can be used in scripts to work with arcpy, numpy and with some stretching... the field calculator.  If you have any useful examples, pass them on.

 

With donut holes in the middle... the radiating line is a matplotlib artifact...I didn't want to waste time removing it.  

Note that the ring widths are not equal... they need not be, you just set the threshold distances you want.

These last two have a rotation set.  And the ellipse is the result of scaling the y-values and rotating the coordinates.

 

As a simple example of the internal structure of the inputs, the following is an example of two triangles (3 points on a circle) with holes.  The data input is simply an input array as shown by the coordinates and the plotting routine handles the output.

a = buffer_ring(outer=10, inner=8, theta=120, xc=0.0, yc=0.0)
b = buffer_ring(outer=10, inner=8, theta=120, xc=10.0, yc=0.0)
a0 = Polygon(a, closed=False)
b0 = Polygon(b, closed=False)
#plot_([a,b])
a0p = a0.get_xy()
b0p = b0.get_xy()
props = b0.properties()

 

Header 1
>>> a0p
array([[-10.000, 0.000],
       [ 5.000, 8.660],
       [ 5.000, -8.660],
       [-10.000, -0.000],
       [-8.000, -0.000],
       [ 4.000, -6.928],
       [ 4.000, 6.928],
       [-8.000, 0.000]])
>>> b0p
array([[ 0.000, 0.000],
       [ 15.000, 8.660],
       [ 15.000, -8.660],
       [ 0.000, -0.000],
       [ 2.000, -0.000],
       [ 14.000, -6.928],
       [ 14.000, 6.928],
       [ 2.000, 0.000]])

       The triangles are as follows:

     

 

 

The coordinates are above.

The outer rings go clockwise and the inner rings are counter clockwise.  You will notice that the first and last point of each ring are identical.