0

I have an interconnected weighted graph. The graph is a representation of a Multilinestring connecting points. See the visualisation in QGIS:

enter image description here

So from each point I know how much it costs to go to every other point. This information I compute from a Shapefile called "CostShape" to a graph in networkx. I use then this command to produce a numpy array:

    matrix = nx.to_numpy_matrix(graph)

The result looks as follows:

[[0.00000000e+00 2.36988656e+01 1.67140574e+04 ... 5.42805416e+02
5.66551990e+02 2.80353597e+04] 
[2.36988656e+01 0.00000000e+00 1.66903585e+04 ... 5.66504281e+02
5.90250855e+02 2.80590586e+04]
[1.67140574e+04 1.66903585e+04 0.00000000e+00 ... 1.72568628e+04
1.72806093e+04 4.47494171e+04]
...
[5.42805416e+02 5.66504281e+02 1.72568628e+04 ... 0.00000000e+00
2.37465741e+01 2.74925543e+04]
[5.66551990e+02 5.90250855e+02 1.72806093e+04 ... 2.37465741e+01
0.00000000e+00 2.74688077e+04]
[2.80353597e+04 2.80590586e+04 4.47494171e+04 ... 2.74925543e+04
2.74688077e+04 0.00000000e+00]]

What already is annoying that I am loosing the names of the node in the graph, but I assume that the order stays the same. I proceed then with the clustering of this Numpy matrix and feed the labels back into the "GraphShape":

labeler = KMeans()
labeler.fit(matrix)
for (row, label) in enumerate(labeler.labels_):
    GraphShape.at[row, 'costcluster'] = label
GraphShape.to_file('graphShape.geojson', driver='GeoJSON')

The result is the image below this text block and it is everything else but clusters. I tried many different cluster methods and none of them produce the things that I expect.

What am I doing wrong?

enter image description here

4
  • I would start by verifying if your assumption about the order in your matrix is actually correct? – Hornbydd May 11 '20 at 10:44
  • thats a good point. any idea how to do that? – Leo May 11 '20 at 11:44
  • What do you precisely expect? Do you want to cluster the point "by" what? And what exactly are your costs? Just to be sure I understood... Maybe what you want to compute is your graph centrality. After that, you may be able to interpolate along each of its edges to have "continuous" values, and finally get back this values to the points. – s.k May 11 '20 at 13:29
  • I extended my answer with another possible solution. – Bence Mélykúti May 16 '20 at 11:18
2
+25

I tried to replicate your setup as best as I could. I downloaded a juncture of two winding secondary highways from OpenStreetMap, loaded them into GeoPandas in Python. For every pair of points I computed the physical distance between the two points and used this as edge weight. (So I have a complete graph, a graph with an edge between any two nodes. I had 208 points. This calculation took some 15-20 sec.) Then I conducted the clustering according to your recipe. I loaded the result into QGIS and when I plot it, the clusters look how I expected them.

road juncture

This led me to suspect that perhaps you were displaying your data wrong in QGIS, although I couldn't easily replicate a random colour distribution. In the Layer Styling, you have to categorise according to the variable costcluster. If you categorise according to e.g. point id, then it can look very random.

QGIS Layer Styling

I don't know the edge or weight structure of your graph. But I realised that it can also be the source of the problem. I tested it with two other graphs derived from the same data. In the case when I only ever enter edges with the physical distance into the graph for pairs of points which are adjacent (consecutive) in any way, the output is rather random.

The reason is clear: all rows of the matrix are dissimilar. The clustering tries more or less to put every node into a separate class. As I observe, if I make n classes, then I get n-1 classes with one or at most two nodes in each and all remaining nodes are in one class. If I increase the number of classes a lot, then QGIS starts returning to colours that are very similar to previously used colours for other classes and the pattern is indiscernible.

If your plot shows 8 classes with 8 colours or even fewer, then I still can't replicate that. (In the case when I entered distances into the graph for any pair of points which are parts of the same way, it behaves more or less as expected.)

Here is a sketch of the code that I used when computing the physical distance between every pair of points and using this as edge weight. I loaded the points from OpenStreetMap into a Pandas DataFrame called nodes whose columns are id, lat, lon.

import shapely
from shapely.geometry import Point
import geopandas

from pyproj import Geod
geod = Geod(ellps="WGS84")

import numpy as np
import networkx as nx
from sklearn.cluster import KMeans

'''
[...]
Here I load the points from my OSM data download into the Pandas DataFrame nodes.
[...]
'''

G = nx.Graph()
G.add_nodes_from(nodes['id'])

edges = list()
#nr_nodes = 50
nr_nodes = len(G.nodes)
for i in range(nr_nodes):
    for j in range(i, nr_nodes):
        lons = [nodes.iloc[k]['lon'] for k in [i, j]]
        lats = [nodes.iloc[k]['lat'] for k in [i, j]]
        dist = geod.line_length(lons, lats)
        edges.append((nodes.iloc[i]['id'], nodes.iloc[j]['id'], {'weight': dist}))

G.add_edges_from(edges)
matrix = nx.to_numpy_matrix(G)
labeler = KMeans(n_clusters=8)
labeler.fit(matrix)

p_list = nodes.apply(lambda x: Point(x[2], x[1]), axis=1) # The list of points with lon, lat.
GraphShape = geopandas.GeoDataFrame(nodes['id'], geometry=p_list)

for (row, label) in enumerate(labeler.labels_):
    GraphShape.at[row, 'costcluster'] = label
GraphShape.to_file('graphShape8.geojson', driver='GeoJSON')

My graphShape8.geojson looks like this:

{
"type": "FeatureCollection",
"features": [
{ "type": "Feature", "properties": { "id": "29587036", "costcluster": 2.0 }, "geometry": { "type": "Point", "coordinates": [ 7.9154778, 47.8702779 ] } },
{ "type": "Feature", "properties": { "id": "264508857", "costcluster": 2.0 }, "geometry": { "type": "Point", "coordinates": [ 7.9151344, 47.8706669 ] } },
{ "type": "Feature", "properties": { "id": "301435315", "costcluster": 2.0 }, "geometry": { "type": "Point", "coordinates": [ 7.9162008, 47.8698175 ] } },
{ "type": "Feature", "properties": { "id": "669212106", "costcluster": 2.0 }, "geometry": { "type": "Point", "coordinates": [ 7.9164004, 47.8696363 ] } },
...
1
  • Hey Bence thank you very much for that effort. It is good to see that my recipe works with your simpler data. And I think you are quite right that the problem comes from the weights on my edges. I meanwhile found a solution for my problem and will post it as the answer. – Leo May 18 '20 at 6:26
0

So I found a solution to cluster the points with a fixed threshold:

clusterThreshold = 1000
clusterNumber = 0
cluster = {}

for pt in range(len(GraphShape)):
    firstpt = list(GraphShape.at[pt, 'geometry'].coords)[0]
    notInClusterYet = True
    # check if point is already somewhere in the dict. if so this cluster value is taken
    for key in cluster.keys():
        if firstpt in cluster[key]:
            notInClusterYet = False
            clusterOfFirstPt = key
            break
    # if not already in the clusterdict, a new clusternumber is registered
    if notInClusterYet:
        clusterNumber += 1
        cluster[clusterNumber] = [firstpt]
        clusterOfFirstPt = clusterNumber
    for spt in range(pt + 1, len(GraphShape)):
        secondpt = list(GraphShape.at[spt, 'geometry'].coords)[0]
        weight, shortestpath = nx.single_source_dijkstra(graph, firstpt, secondpt)
        print(weight)
        if weight < clusterThreshold:
            coordinatelist = cluster[clusterOfFirstPt]
            coordinatelist.append(secondpt)
            cluster[clusterOfFirstPt] = coordinatelist
finished = False
while not finished:
    finished = True
    for key in cluster.keys():
        for coordinate in cluster[key]:
            for subkey in cluster.keys():
                if subkey == key:
                    continue
                for subCoordinate in cluster[subkey]:
                    if subCoordinate == coordinate:
                        newcoordinates = cluster[subkey] + cluster[key]
                        cluster[key] = newcoordinates
                        del cluster[subkey]
                        finished = False
                        break
                if not finished:
                    break
            if not finished:
                break
        if not finished:
            break

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.