0

I have two point layers, one with 233 row, other with 1531. For every point in layer 1, I would like to get the closest point from layer 2, (ok simple enough, but here's a trick) that has the value of attribute 1 lower than the compared point in layer 2. This attribute being time, so it occurred before (I converted time to integer for easier comparison).

I wanted to go with distance matrix, but there's no filter option, so I just tried to do the complete matrix, but that didn't work since there s too much data. The solution doesn't need to include distance matrix, I can also do some simplescripting I guess.

Sorry if I didn't express myself clearly, I will edit, update if needed.

  • 1
    When you say "didn't work since there is too much data", do you mean the tool failed to run, or the result produced more than you wanted? Because it's easy enough to filter the results down to what you want with a couple of joins. There are other ways to do it (for example the SQL answer at gis.stackexchange.com/questions/126639). Do you want the nearest only one value lower, or nearest of any that came before? A couple of sample rows from each table might help in understanding the question, or a more generic description (they're gps points, or you're making tracks?). – Chris W May 13 '15 at 5:35
  • Sorry, I pressed enter to fast. The distance matrix didn't work means it failed to run, I tried limiting it to nearest 10, it worked, took it more than an hour, but sadly not all rows from layer 1 got a corresponding row from layer 2 (events occurred later). So I will try with SQL, table join would probably be best. I want the nearest point that comes before. Layer header (for both layers) looks like this: Date | lon | lat | int_time | some | other | attributes – Matos4df May 13 '15 at 16:16
  • Oh and point are just gps point, occurrence time span is the same for both layers. If the final result would be the map of red points(layer 1), each connected to the closest blue point that fits the criteria (layer 2), with some color coded (based on length) lines, that would be awesome. – Matos4df May 13 '15 at 16:31
  • I guess there is no straight answer for this. Anyway I've accomplished the task, the key was... well faster computer that could put the matrix together (well the closest 30 was enough for all but 2 cases). When I had the matrix done, I went on with filtering in excel (not best fun on the planet). To create the visual output I made a csv join on 1st layer and that's about it. So should I mark this as answered/stupid question/anything else, or leave it up for someone to come up with a great answer? – Matos4df May 29 '15 at 17:39
  • If you got a solution that works, I'd go ahead and make an answer out of it. You can answer your own question and accept that answer. Or if you want to accept the current answer you can do that (I'm not sure how different your ultimate solution was to said answer). Without some sample data to play with I was having a hard time coming up with a solution, and not really sure why the distance matrix tool would be failing - that shouldn't have been too many points (or anywhere near). I'd have gone for each of the 233 to all the 1531, then use sorting/selection to filter that result down. – Chris W May 29 '15 at 23:01
1

NetworkX (http://cheeseshop.python.org/pypi/networkx/) is a Python package with many functions for graph and network analysis. You can use it in QGIS analyze graphs without prior construction of topology.

Your problem of connecting every point of a set A with points of a set B can be modeled as a bipartite graph. I use the data structure and some methods for handling the connections from Graph object.

You may render the linesegments using the length (distance) or timediff (difference of time stamps of start and end point).

To use the following code in Python console you have to load the csv files as layer layer_a and layer_b.

from PyQt4.QtCore import QVariant
import networkx as nx
from math import sqrt

# for testing
layer_a = [(0, 4, 8, 1000), (1, 9, 2, 1200)]
layer_b = [(0, 1, 7, 958), (1, 2, 6, 1010), (2, 3, 5, 1002), (3, 4, 9, 1057), (4, 6, 10, 1005), (5, 6, 5, 959), (6, 7, 4, 1030), (7, 7, 2, 1205), (8, 9, 1, 1158), (9, 11, 0, 1210), (10, 10, 5, 1159), (11, 11, 5, 1155)]

# helper function to calculate euclidian distance
def length(p1, p2): return sqrt((p1[0]-p2[0])**2 + (p1[1]-p2[1])**2)

# define QGIS layer with some attributes
layer = QgsVectorLayer('LineString?crs=EPSG:4326', 'Connected', 'memory')
prov = layer.dataProvider()
prov.addAttributes([QgsField("distance", QVariant.Double), QgsField("timediff", QVariant.Int)])
QgsMapLayerRegistry.instance().addMapLayer(layer)

# lets go
layer.startEditing()

# create an empty graph, then add edges from every point in set A to every point in 
# set B assuming an attribute with lon at index 1, lat at index 2, and time at index 3.
# Please change indizes if this is not appropriate.
bigraph = nx.DiGraph()
bigraph.add_edges_from([((a[1], a[2]), (b[1], b[2]), {'length': length((a[1], a[2]), (b[1], b[2])), 'time': a[3]-b[3]}) for b in layer_b for a in layer_a])

# now remove from this completed graph all edges to points in B having a time difference less then 0 (timestamp of B is newer then A)
remove_list = []
for n, nbrs in bigraph.adjacency_iter():
    for nbr, eattr in nbrs.items():
        data = eattr['time']
        if data < 0: remove_list.append((n, nbr))

bigraph.remove_edges_from(remove_list)

# now iterate over every edge of every point in set A to find the one with minimal length. If the edge with minimal length is found, create a polyline from the edge
feats = []
for n in [(a[1], a[2]) for a in layer_a]:
    elen = 99999.0
    for e in bigraph.edges_iter([n], data=True):
        if e[2]['length'] < elen:
            elen = e[2]['length']
            gstr = e
    feat = QgsFeature()
    feat.setGeometry(QgsGeometry.fromPolyline([QgsPoint(gstr[0][0], gstr[0][1]), QgsPoint(gstr[1][0], gstr[1][1])]))
    feat.setAttributes([gstr[2]['length'], gstr[2]['time']])
    feats.append(feat)

# finally add all created features and save edits
prov.addFeatures(feats)
layer.updateExtents()
layer.commitChanges()
  • Thanks for answer, looks like something that would get the job done. I did the changes I thought were necessary, but I get an error: "DigitizingTools: Prolong Line is disabled because layer CRS and project CRS do not match!" I did change CRS of project to 4326, same for csv files, I even removed other layers from project, no difference :/ – Matos4df May 27 '15 at 19:23
  • Sounds like some process from DigitizingTools plugin is listening to edit events. Unfortunately I am not familiar with this very useful plugin. I would suggest to (temporarily) unload DigitizingTools and try the code again. But even with this plugin activated I couldn't reproduce your error (QGIS 2.4.0, WinXP32). – Detlev May 28 '15 at 20: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.