# Saving geographic minimum spanning tree in KML and shapefile using Python?

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):

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:

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

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?
– user1462
Commented Sep 15, 2017 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. Commented Sep 15, 2017 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?
– user1462
Commented Sep 15, 2017 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? Commented Sep 15, 2017 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).
– user1462
Commented Sep 15, 2017 at 16:09

## 1 Answer

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])
``````