0

I finally got an efficient Python implementation of a minimum-spanning tree algorithm, which I want to use to compute MST's among geographic locations.

My code computes the MST from a list of longitude, latitude data and returns a list of edges joining the locations.

Here is the code and some test data:

import warnings
import numpy as np
from hcluster import pdist, squareform
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import shapefile
import simplekml

warnings.filterwarnings("ignore")

def minimum_spanning_tree(X, copy_X=True):
    """X are edge weights of fully connected graph"""
    if copy_X:
        X = X.copy()

    if X.shape[0] != X.shape[1]:
        raise ValueError("X needs to be square matrix of edge weights")
    n_vertices = X.shape[0]
    spanning_edges = []

    # initialize with node 0:                                                                                        
    visited_vertices = [0]                                                                                           
    num_visited = 1
    # exclude self connections:
    diag_indices = np.arange(n_vertices)
    X[diag_indices, diag_indices] = np.inf

    while num_visited != n_vertices:
        new_edge = np.argmin(X[visited_vertices], axis=None)
        # 2d encoding of new_edge from flat, get correct indices                                                     
        new_edge = divmod(new_edge, n_vertices)
        new_edge = [visited_vertices[new_edge[0]], new_edge[1]]                                                      
        # add edge to tree
        spanning_edges.append(new_edge)
        visited_vertices.append(new_edge[1])
        # remove all edges inside current tree
        X[visited_vertices, new_edge[1]] = np.inf
        X[new_edge[1], visited_vertices] = np.inf                                                                    
        num_visited += 1
    return np.vstack(spanning_edges)

def test_mst():
    P = np.loadtxt("testdata.csv", delimiter=',', skiprows=1, usecols=(1,2))

    X = squareform(pdist(P))

    edge_list = minimum_spanning_tree(X)

    map = Basemap()
    map.drawcoastlines(color='#777799')
    map.drawcountries(color='#ccccee')
    map.drawmapboundary()
    map.bluemarble()
    lon = P[:,0]
    lat = P[:,1]
    x,y = map(lon,lat)
    map.plot(x,y,'ro')

    edges = []
    for edge in edge_list:
        i, j = edge
        plt.plot([P[i, 0], P[j, 0]], [P[i, 1], P[j, 1]], c='r')
        edges.append([(P[i, 0], P[j, 0]), (P[i, 1], P[j, 1])])
    plt.show()

    # Save shapefile
    w = shapefile.Writer()
    w.line(parts=[edge_list])
    w.field("COMMON_ID", 'C')
    w.record("Point")
    w.save("mst")

    # Save KML file
    kml = simplekml.Kml()
    for i in range(len(edges)):
        line = kml.newlinestring(name="Track" + str(i + 1), coords=edges[i])
        line.style.linestyle.width = 3
        line.style.linestyle.color = simplekml.Color.red
    kml.save("mst.kml")

if __name__ == "__main__":
    test_mst()

And here is my test data:

Locality,Longitude,Latitude
Bush Bush Forest,-61.05,10.4
Cerro Azul,-79.4333333333,9.15
Dividive,-70.6666666667,9.53333333333
Hato El Frio,-63.1166666667,7.91666666667
Finca Vuelta Larga,-63.1166666667,10.55
Isla Cebaco,-81.1833333333,7.51666666667
Kayserberg Airstrip,-56.4833333333,3.1
Limao,-60.5,3.93333333333
Montijo Bay,-81.0166666667,7.66666666667
Parcela 200,-67.4333333333,8.93333333333
Rio Chico,-65.9666666667,10.3166666667
San Miguel Island,-78.9333333333,8.38333333333
Tukuko,-72.8666666667,9.83333333333
Urama,-68.4,10.6166666667
Valledup,-72.9833333333,10.6166666667

When I plot the tree on the map, using Matplotlib/Basemap, it displays nicely (see the figure below):

enter image description here

But when I attempt to save the same coordinates used to create the map plot into a KML file (using the simpleKML library) I got the (horrible) figure below:

enter image description here

The shapefile generated using the PyShp library is no less horrible:

enter image description here

This problem is simliar to one I myself presented some time ago (Storing a geographic minimum spanning tree with PyShp) but the solution offered there (interleaving of the points) seems to be not working in this case.

EDIT 1: Here is the list of edges generated by the MST function:

[[ 0  4]
 [ 4  3]
 [ 4 10]
 [10  9]
 [ 9 13]
 [13  2]
 [ 2 12]
 [12 14]
 [ 3  7]
 [ 7  6]
 [12 11]
 [11  1]
 [ 1  8]
 [ 8  5]]

EDIT 2: Here is the KML file output by the program:

<?xml version="1.0" encoding="UTF-8"?>
<kml xmlns="http://www.opengis.net/kml/2.2" xmlns:gx="http://www.google.com/kml/ext/2.2">
    <Document id="feat_1">
        <Style id="stylesel_0">
            <LineStyle id="substyle_0">
                <color>ff0000ff</color>
                <colorMode>normal</colorMode>
                <width>3</width>
            </LineStyle>
        </Style>
        <Style id="stylesel_1">
            <LineStyle id="substyle_1">
                <color>ff0000ff</color>
                <colorMode>normal</colorMode>
                <width>3</width>
            </LineStyle>
        </Style>
        <Style id="stylesel_2">
            <LineStyle id="substyle_2">
                <color>ff0000ff</color>
                <colorMode>normal</colorMode>
                <width>3</width>
            </LineStyle>
        </Style>
        <Style id="stylesel_3">
            <LineStyle id="substyle_3">
                <color>ff0000ff</color>
                <colorMode>normal</colorMode>
                <width>3</width>
            </LineStyle>
        </Style>
        <Style id="stylesel_4">
            <LineStyle id="substyle_4">
                <color>ff0000ff</color>
                <colorMode>normal</colorMode>
                <width>3</width>
            </LineStyle>
        </Style>
        <Style id="stylesel_5">
            <LineStyle id="substyle_5">
                <color>ff0000ff</color>
                <colorMode>normal</colorMode>
                <width>3</width>
            </LineStyle>
        </Style>
        <Style id="stylesel_6">
            <LineStyle id="substyle_6">
                <color>ff0000ff</color>
                <colorMode>normal</colorMode>
                <width>3</width>
            </LineStyle>
        </Style>
        <Style id="stylesel_7">
            <LineStyle id="substyle_7">
                <color>ff0000ff</color>
                <colorMode>normal</colorMode>
                <width>3</width>
            </LineStyle>
        </Style>
        <Style id="stylesel_8">
            <LineStyle id="substyle_8">
                <color>ff0000ff</color>
                <colorMode>normal</colorMode>
                <width>3</width>
            </LineStyle>
        </Style>
        <Style id="stylesel_9">
            <LineStyle id="substyle_9">
                <color>ff0000ff</color>
                <colorMode>normal</colorMode>
                <width>3</width>
            </LineStyle>
        </Style>
        <Style id="stylesel_10">
            <LineStyle id="substyle_10">
                <color>ff0000ff</color>
                <colorMode>normal</colorMode>
                <width>3</width>
            </LineStyle>
        </Style>
        <Style id="stylesel_11">
            <LineStyle id="substyle_11">
                <color>ff0000ff</color>
                <colorMode>normal</colorMode>
                <width>3</width>
            </LineStyle>
        </Style>
        <Style id="stylesel_12">
            <LineStyle id="substyle_12">
                <color>ff0000ff</color>
                <colorMode>normal</colorMode>
                <width>3</width>
            </LineStyle>
        </Style>
        <Style id="stylesel_13">
            <LineStyle id="substyle_13">
                <color>ff0000ff</color>
                <colorMode>normal</colorMode>
                <width>3</width>
            </LineStyle>
        </Style>
        <Placemark id="feat_2">
            <name>Track1</name>
            <styleUrl>#stylesel_0</styleUrl>
            <LineString id="geom_0">
                <coordinates>-61.05,-63.1166666667,0.0 10.4,10.55,0.0</coordinates>
            </LineString>
        </Placemark>
        <Placemark id="feat_3">
            <name>Track2</name>
            <styleUrl>#stylesel_1</styleUrl>
            <LineString id="geom_1">
                <coordinates>-63.1166666667,-63.1166666667,0.0 10.55,7.91666666667,0.0</coordinates>
            </LineString>
        </Placemark>
        <Placemark id="feat_4">
            <name>Track3</name>
            <styleUrl>#stylesel_2</styleUrl>
            <LineString id="geom_2">
                <coordinates>-63.1166666667,-65.9666666667,0.0 10.55,10.3166666667,0.0</coordinates>
            </LineString>
        </Placemark>
        <Placemark id="feat_5">
            <name>Track4</name>
            <styleUrl>#stylesel_3</styleUrl>
            <LineString id="geom_3">
                <coordinates>-65.9666666667,-67.4333333333,0.0 10.3166666667,8.93333333333,0.0</coordinates>
            </LineString>
        </Placemark>
        <Placemark id="feat_6">
            <name>Track5</name>
            <styleUrl>#stylesel_4</styleUrl>
            <LineString id="geom_4">
                <coordinates>-67.4333333333,-68.4,0.0 8.93333333333,10.6166666667,0.0</coordinates>
            </LineString>
        </Placemark>
        <Placemark id="feat_7">
            <name>Track6</name>
            <styleUrl>#stylesel_5</styleUrl>
            <LineString id="geom_5">
                <coordinates>-68.4,-70.6666666667,0.0 10.6166666667,9.53333333333,0.0</coordinates>
            </LineString>
        </Placemark>
        <Placemark id="feat_8">
            <name>Track7</name>
            <styleUrl>#stylesel_6</styleUrl>
            <LineString id="geom_6">
                <coordinates>-70.6666666667,-72.8666666667,0.0 9.53333333333,9.83333333333,0.0</coordinates>
            </LineString>
        </Placemark>
        <Placemark id="feat_9">
            <name>Track8</name>
            <styleUrl>#stylesel_7</styleUrl>
            <LineString id="geom_7">
                <coordinates>-72.8666666667,-72.9833333333,0.0 9.83333333333,10.6166666667,0.0</coordinates>
            </LineString>
        </Placemark>
        <Placemark id="feat_10">
            <name>Track9</name>
            <styleUrl>#stylesel_8</styleUrl>
            <LineString id="geom_8">
                <coordinates>-63.1166666667,-60.5,0.0 7.91666666667,3.93333333333,0.0</coordinates>
            </LineString>
        </Placemark>
        <Placemark id="feat_11">
            <name>Track10</name>
            <styleUrl>#stylesel_9</styleUrl>
            <LineString id="geom_9">
                <coordinates>-60.5,-56.4833333333,0.0 3.93333333333,3.1,0.0</coordinates>
            </LineString>
        </Placemark>
        <Placemark id="feat_12">
            <name>Track11</name>
            <styleUrl>#stylesel_10</styleUrl>
            <LineString id="geom_10">
                <coordinates>-72.8666666667,-78.9333333333,0.0 9.83333333333,8.38333333333,0.0</coordinates>
            </LineString>
        </Placemark>
        <Placemark id="feat_13">
            <name>Track12</name>
            <styleUrl>#stylesel_11</styleUrl>
            <LineString id="geom_11">
                <coordinates>-78.9333333333,-79.4333333333,0.0 8.38333333333,9.15,0.0</coordinates>
            </LineString>
        </Placemark>
        <Placemark id="feat_14">
            <name>Track13</name>
            <styleUrl>#stylesel_12</styleUrl>
            <LineString id="geom_12">
                <coordinates>-79.4333333333,-81.0166666667,0.0 9.15,7.66666666667,0.0</coordinates>
            </LineString>
        </Placemark>
        <Placemark id="feat_15">
            <name>Track14</name>
            <styleUrl>#stylesel_13</styleUrl>
            <LineString id="geom_13">
                <coordinates>-81.0166666667,-81.1833333333,0.0 7.66666666667,7.51666666667,0.0</coordinates>
            </LineString>
        </Placemark>
    </Document>
</kml>

Can anyone give me a hand?

  • Have you considered using networkx to generate a nicer MST? – barrycarter Sep 15 '17 at 15:58
  • Can't understand nothing of NetworkX data structures, frankly. But why is the algorithm above not nice? It works fine, in approximately linear time. My problem is not with the algorithm, but with using its output to create shapefiles and KML files. – maurobio Sep 15 '17 at 16:02
  • My mistake, sorry. However, if you have the points and the edges of the MST, isn't there only one possible output plot? Can you look at the generated KML files, compare them to the output of your program, and see what went wrong? – barrycarter Sep 15 '17 at 16:04
  • That is just my problem. The plot is absolutely correct, why is it that the shapefile and KML file using the same data are wrong? – maurobio Sep 15 '17 at 16:07
  • Could you provide us with the output of your program and the KML and SHP outputs (use some sort of pastebin or dropbox or something, since you can't attach non-image files here). – barrycarter Sep 15 '17 at 16:09
1

Your KML still isn't interleaved correctly. Your original data is approximately -70 longitude and 10 latitude, you have a cluster of points and around (-70, -70), and another cluster around (10, 10).

<coordinates>-67.4333333333,-68.4,0.0 8.93333333333,10.6166666667,0.0</coordinates>

Looking at the table of of the node values, you can see that 67.43 should be with 8.93, and -68.4 should be with 10.61. So you list that should be (in JSON-like format):

[[x1,y1], [x2,y2]

is still stuck in

[[x1,x2], [y1,y2]]

so that suggests

edges.append([(P[i, 0], P[j, 0]), (P[i, 1], P[j, 1])])

should be

edges.append([(P[i, 0], P[i, 1]), (P[j, 0], P[j, 1])])

Now for the crazy shapefile. I wasn't even sure where it was at first until I spotted the problem. Then I realized that it was the west coast of Africa, and since that's close to (0,0) this is a big hint.

w.line(parts=[edge_list])

edge list is your node numbers, which is why your shapefile is jogging around between (1,1) and (14,14). You just got the wrong variable:

w.line(parts=[edges])

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.