2

I'm working on this script in PyQGIS to take points from a layer called LEK_ALL and snap them to the closest feature of a line layer called LEID_filtered which fulfills a certain request.

In order to do so, I created the function snapPointToLine. This is the full code:

from qgis.PyQt.QtGui import *
from qgis.PyQt.QtWidgets import *
from qgis.PyQt.QtCore import *
from qgis.core import *
import os

def snapPointToLine(feature, lineLayer, request, AttrID):
        '''
        QgsFeature (point), QgsLayer (lines), QgsFeatureRequest, list -> QgsFeature
        Takes the feature feature and snaps it to the colest point of the feature in layer lineLayer that fullfills request.
        Returns the new feature with the same attributes as the old ones, plus the attributes of the snapped-to feature with indices in AttrID.
        '''
        pt = feature.geometry().asPoint()
        lineIter = lineLayer.getFeatures(request) #iterator of the features in lineLayer fulfilling the request
        spIndex = QgsSpatialIndex(lineIter) #create spatial index from selected feature
        nearestIds = spIndex.nearestNeighbor(pt,1) #list of 1 closest neighbour
        featureId = nearestIds[0] #item 0 (and only) of the list is the id
        nearestIterator = lineLayer.getFeatures(QgsFeatureRequest().setFilterFid(featureId)) #generate iterator containing only the closest feature
        nearFeature = QgsFeature() #generate an empty feature
        nearestIterator.nextFeature(nearFeature) #copy information from nearestIterator to nearFeature
        closeSegResult = nearFeature.geometry().closestSegmentWithContext(pt) #results in a description of the closest segment between the nearFeature and point
        closePoint = closeSegResult[1] #closest point on geometry to LEK
        #generate feature in new location
        newFeature = QgsFeature()
        newFeature.setGeometry(QgsGeometry.fromPoint(QgsPoint(closePoint[0],closePoint[1])))
        newFeature.setAttributes(feature.attributes() + [nearFeature.attributes()[i] for i in AttrID])
        return newFeature

registry = QgsMapLayerRegistry.instance()
LEK = registry.mapLayersByName('LEK_ALL')[0]
LEID = registry.mapLayersByName('LEID_filtered')[0]

LEKIter =  LEK.getFeatures()
Date_id = LEK.fieldNameIndex('Date') #gets id of field "Date"
LEID_id = LEID.fieldNameIndex('AssetNum') #gets id for field "AssetNum"

newFeat_list = [] #generate new empt list to save the new features

for feature in LEKIter:
        LEKDate = feature.attributes()[Date_id] #get failure date
        y,m,d = LEKDate.year(),LEKDate.month(),LEKDate.day() #split it into components
        request_str = u'"InstYear" <= to_date(\'%d-%02d-%02d\') AND "RemDate*" >= to_date(\'%d-%02d-%02d\')' % (y,m,d,y,m,d) #generate query pt.1
        request = QgsFeatureRequest().setFilterExpression(request_str) #generate query pt.2
        newFeature = snapPointToLine(feature, LEID, request, [LEID_id])
        newFeat_list.append(newFeature)

When I run it I get the following error:

NameError: global name 'QgsSpatialIndex' is not defined

Associated to the line

spIndex = QgsSpatialIndex(lineIter) #create spatial index from selected feature

Any suggestions?

1

It seems that you are using the QgsSpatialIndex() class in a wrong way.

Once you have defined it:

spIndex = QgsSpatialIndex()

you need to insert all the features by iterating over a layer (or a selection, as in your specific case):

for feat in lineIter:
    spIndex.insertFeature(feat)

That's all, so your function becomes:

def snapPointToLine(feature, lineLayer, request, AttrID):
    '''
    QgsFeature (point), QgsLayer (lines), QgsFeatureRequest, list -> QgsFeature
    Takes the feature feature and snaps it to the colest point of the feature in layer lineLayer that fullfills request.
    Returns the new feature with the same attributes as the old ones, plus the attributes of the snapped-to feature with indices in AttrID.
    '''
    pt = feature.geometry().asPoint()
    lineIter = lineLayer.getFeatures(request) #iterator of the features in lineLayer fulfilling the request
    spIndex = QgsSpatialIndex() #create spatial index from selected feature
    for feat in lineIter: #iterate over the selection
        spIndex.insertFeature(feat)
    nearestIds = spIndex.nearestNeighbor(pt,1) #list of 1 closest neighbour
    featureId = nearestIds[0] #item 0 (and only) of the list is the id
    nearestIterator = lineLayer.getFeatures(QgsFeatureRequest().setFilterFid(featureId)) #generate iterator containing only the closest feature
    nearFeature = QgsFeature() #generate an empty feature
    nearestIterator.nextFeature(nearFeature) #copy information from nearestIterator to nearFeature
    closeSegResult = nearFeature.geometry().closestSegmentWithContext(pt) #results in a description of the closest segment between the nearFeature and point
    closePoint = closeSegResult[1] #closest point on geometry to LEK
    #generate feature in new location
    newFeature = QgsFeature()
    newFeature.setGeometry(QgsGeometry.fromPoint(QgsPoint(closePoint[0],closePoint[1])))
    newFeature.setAttributes(feature.attributes() + [nearFeature.attributes()[i] for i in AttrID])
    return newFeature

BONUS I performed the reversed operation here, maybe it could help.

  • Thanks for your answer, unfortunately this does not seem to work. Also: spIndex = QgsSpatialIndex(lineIter) works fine in the console. – ltr Apr 13 '17 at 9:21
  • @ltr Maybe a simple import error? Do you have from qgis.core import * at the beginning of the script? – mgri Apr 13 '17 at 9:23
  • It's not in the code I copied, but I actually do have it there. The funny thing is that the code of the function works fine when I run it individually instead of defining a function with it. – ltr Apr 13 '17 at 9:27
  • Does it work if you use lineIter = lineLayer.getFeatures(request) without the request, so taking all the features? I mean, using this line lineIter = lineLayer.getFeatures() – mgri Apr 13 '17 at 9:32
  • For whatever reason, closing everything and starting over seems to have done the trick. I am always baffled by computers! Now the issue is that QGIS crashes running the script, I think it is the call to getFeatures(request) done for every point-feature (1300 or so times), but that's an issue for another question. – ltr Apr 13 '17 at 9:45

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.