3

I have a list of geographic coordinates (longitudes, latitudes in decimal degrees) as follows:

-72.9833,10.6167
-63.1167,10.55
-63.1167,7.91667
-65.9667,10.3167
-67.4333,8.93333
-68.4,10.6167
-70.6667,9.53333
-72.8667,9.83333
-72.9833,10.6167
-60.5,3.93333
-56.4833,3.1
-78.9333,8.38333
-79.4333,9.15
-81.0167,7.66667
-81.1833,7.51667

and a list of edges connecting these points with a minimum spanning tree (computed using Prim's algorithm) as shown below:

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

I am trying to save the minimum spanning tree as a shapefile using the PyShp library (so that the MST can be displayed on a GIS).

Here is my code so far:

    import numpy as np
    from mpl_toolkits.basemap import Basemap
    import matplotlib.pyplot as plt
    import shapefile

    # Load coordinates
    coords = np.loadtxt('coords.dat')
    lon = coords[:,0]
    lat = coords[:,1]

    # Load graph edges
    edges = np.loadtxt('edges.dat')
    edges = np.array(edges, dtype=np.int32)

    # Create map
    map = Basemap()

    map.drawcoastlines()
    map.drawcountries()
    map.drawmapboundary()
    x,y = map(lon,lat)
    map.plot(x,y,'ro')
    plt.title('Minimum Spanning Tree')

    # Ths MST plotting algorithm works!
    edge_list = []
    n = len(edges)
    for i in range(n):
        t = edges[i,0]-1
        u = edges[i,1]-1
        x = [coords[t,0], coords[u,0]]
        y = [coords[t,1], coords[u,1]]
        #edge_list.append(x)
        #edge_list.append(y)
        #edge_list.append((x,y))
        edge_list.append([(x[0],y[0]), (x[1], y[1])])
        plt.plot(x,y,'-r')
    plt.show()

# Sort the list of edges
edge_list = list(sorted(edge_list))

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

The code above correctly plots the points and draws the minimum spanning tree connecting them, as shown in the figure below:

enter image description here

But the tree is not being correctly saved in the shapefile, as shown in the figure below (displayed by QGIS):

enter image description here

An alternative would be to have the tree saved as a KML file (compatible Google Earth), which could then be converted into a shapefile using the ogr2ogr tool of GDAL. I have not tried this as yet.

Any hints?

UPDATE 1: Following the suggestion of Marc Pfister (see comments below), I stored the coordinates as lists of tuples in edge_list: edge_list.append((x, y)). The results are shown below:

[([-61.049999999999997, -63.116666666699999], [10.4, 10.550000000000001]), ([-63.116666666699999, -63.116666666699999], [7.9166666666700003, 10.550000000000001]), ([-63.116666666699999, -65.966666666699993], [10.550000000000001, 10.3166666667]), ([-67.433333333299998, -65.966666666699993], [8.9333333333299993, 10.3166666667]), ([-67.433333333299998, -68.400000000000006], [8.9333333333299993, 10.6166666667]), ([-70.666666666699996, -68.400000000000006], [9.5333333333300008, 10.6166666667]), ([-70.666666666699996, -72.866666666699999], [9.5333333333300008, 9.8333333333299997]), ([-72.866666666699999, -72.983333333299996], [9.8333333333299997, 10.6166666667]), ([-63.116666666699999, -60.5], [7.9166666666700003, 3.9333333333299998]), ([-56.483333333300003, -60.5], [3.1000000000000001, 3.9333333333299998]), ([-78.933333333299998, -72.866666666699999], [8.3833333333300004, 9.8333333333299997]), ([-79.433333333299998, -78.933333333299998], [9.1500000000000004, 8.3833333333300004]), ([-79.433333333299998, -81.016666666700004], [9.1500000000000004, 7.6666666666700003]), ([-81.183333333299998, -81.016666666700004], [7.5166666666699999, 7.6666666666700003])]

The shapefile created by storing these data is still uglier than the previous one.

enter image description here

UPDATE 2: Fixed the storing of the (x,y) tuples in edge_list. It now works fine! Have updated the code above to include this fix.

  • 1
    Your updated edge list is ((x,x), (y,y), (x,x), (y,y)...) – Marc Pfister Aug 31 '17 at 21:22
  • Sure, in a rush I haven't noticed. Now I fixed it. It works! Thanks! – maurobio Aug 31 '17 at 21:42
2

Looking at the second graphic it looks like edge_list is all the coordinate values mixed together and sorted in order, not coordinate pairs. For the shapefile geometry you should be appending (x,y) tuples to a list.

Pyplt wants two lists: [x1,x2,x3] and [y1,y2,y3]

Pyshp wants one list of a list of coordinate pairs: [[[x1,y1], [x2,y2], [y3,y3]]]

Essentially you want to interleave your plotting data.

  • Doing this edge_list.append((x,y)) does not results in a different shapefile (cannot include images in comments, unfortunately). – maurobio Aug 31 '17 at 21:06
  • Just update the original post with the contents of edge_list – Marc Pfister Aug 31 '17 at 21:09
  • Interleave, that's it. – maurobio Aug 31 '17 at 21:46

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.