I have written a custom processing tool, which I use within a large graphical model in QGIS. I am looking for a way to speed it up.

First step of this processing tool is to connect every point of one input layer to every point of another input layer and measure its line length. This takes quite a while, because usually both layers have more than 10000 points. But I am only interested in point pairs within a given distance to each other. So currently I am going to delete the lines afterwards, which are too long. I think this is quite ineffective, and therefore wanted to ask if there is an easy way to speed up this step, by using a spatial index or something like that. I could create one by using index = QgsSpatialIndex(layer.getFeatures()), but how would I proceed to only connect points I am interested in?

This is a minimal part of my current code:

Memorylayer_VL = QgsVectorLayer(str('LineString?crs=25832')), "tmp_lines", "memory")
Memorylayer_PR = Memorylayer_VL.dataProvider()
maxdistance = 150

# iterate over stop features and create lines
for stop_feat in stop_layer.getFeatures():
    point1 = QgsPoint(stop_feat.geometry().asPoint())
    for source_feat in possibility_layer.getFeatures():
        point2 = QgsPoint(source_feat.geometry().asPoint())
        temp_feat = QgsFeature(fields)
        temp_feat.setGeometry(QgsGeometry.fromPolyline([point1, point2]))
        temp_feat["line_length"] = temp_feat.geometry().length()    
        # Only add lines shorter than maxdistance and longer than 1m (Already simplified like this thanks to @BERAs suggestion in comments) 
        # However, couldn't using a spatial index even have more effect?
        if (temp_feat.geometry().length() < maxdistance) and (tmp_line_feat.geometry().length() > 1):       

# Doing more processing stuff on the just created lines layer....
  • 1
    Dont create a line if the distance is > maxdistance instead of removing afterwards. Are you comparing each point twice? p1-p2 and then p2-p1? – BERA Oct 9 '20 at 11:31
  • 1
    @BERA Good and simple point, wondering why I did not think about this. Still I think using a spatial index would have even more effect wouldnt it? No, only once: p1-p2. – MrXsquared Oct 9 '20 at 11:38

Below the code for select all solution lines between 0 (excluded) and maxdistance (excluded too). All POSSIBILITY features are inserted inside a spatial index (QgsSpatialIndex). A filter, thanks to the spatial index, with a maxdistance square around the feature of the layer STOP is applied for speed the process and select (intersects) only POSSIBILITY features in the square.

If the layer STOP and POSSIBILITY are the same, the code below makes duplicates and isn't not very efficient because of the duplicate.

qgis_project = QgsProject.instance()
layers = {lyr.name(): lyr for lyr in qgis_project.mapLayers().values()}
lyr_from = layers["STOP"]
lyr_to = layers["POSSIBILITY"]
maxdistance = 150

spatial_idx_to = QgsSpatialIndex()
pr_from = lyr_from.dataProvider()
pr_to = lyr_to.dataProvider()
feats_from = pr_from.getFeatures()
feats_to = pr_to.getFeatures()

# insert destination features to the spatial index

result_lines = []
# iterate over origin features
for ft_from in feats_from:
    # find all the nearest feature thanks to the spatial index
    geom_from = ft_from.geometry().asPoint()
    x_from = geom_from.x()
    y_from = geom_from.y()
    # create a rectangle at maxdistance of the from point
    rect = QgsRectangle(
        QgsPointXY(x_from - maxdistance, y_from - maxdistance),
        QgsPointXY(x_from + maxdistance, y_from + maxdistance)
    # get all points in that intersects the rectangle
    distance_ids = spatial_idx_to.intersects(rect)
    for fid in distance_ids:
        geom_to = lyr_to.getFeature(fid).geometry().asPoint()
        # if the distance between origin and destination < maxdistance
        if 0 < geom_from.distance(geom_to) < maxdistance:
            line_feat = QgsFeature()
            # create a PolylineXY from origin and destination PointXY
            line_feat.setGeometry(QgsGeometry.fromPolylineXY([geom_from, geom_to]))
            # append the line feature to the result list

# create the memory layer
memlyr = QgsVectorLayer("LineString?crs=epsg:25832", "tmp_lines", "memory")
memlyr_pr = memlyr.dataProvider()
# add all resulting lines features
# finally add the line memory layer to the project

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.