getFeatures() is inaccurate with memory layers

Question

I was writing a plugin for QGIS when I wanted to copy all my data to memory layers and edit them instead of actual shapefiles to speed up the process. In a QgsMapTool-derived class, I am using a method call like this

def canvasPressEvent(self, e):
    mapCoords = self.toMapCoordinates(e.pos())
    features = self.__layer.getFeatures(QgsFeatureRequest(QgsRectangle(mapCoords, mapCoords)))
    #...

This worked just fine for the old, file-based layers, however, for the new in-memory ones, I'm getting horribly inaccurate results.

My code for copying the layer is pretty straightfoward:

with fiona.open(path) as f:
    geometry = f.schema["geometry"] # yeah, this is a bit hacky, sry

originalLayer = QgsVectorLayer(path, name + "_temp_layer", "ogr")

# build the creation string
props = geometry + "?"
props += "crs=WKT:" + originalLayer.crs().toWkt()
props += "&"
props += "index=yes"

# create layer and start filling it
layer = QgsVectorLayer(props, name, "memory")
#layer.setCrs(originalLayer.crs()) #this doesn't seem to change anything
layer.startEditing()

# copy fields
for field in originalLayer.fields():
    layer.addAttribute(QgsField(field.name(), field.type()))

# copy features
newFeatures = []
for feature in originalLayer.getFeatures():
    newFeature = QgsFeature()
    newFeature.setGeometry(feature.geometry())
    newFeature.setAttributes(feature.attributes())
    newFeatures.append(newFeature)
result = layer.dataProvider().addFeatures(newFeatures)

# commit
layer.commitChanges()
layer.updateExtents()

Now, this code seems to work just fine, the result looks just as it looked before, but when using the tool, I can click on a spot a few centimetres away from a feature and the tool detects it as if I clicked the feature itself. This also creates massive overlap between the features. These things worked just fine with file-based layers.

It feels like something is wrong with the CRS... Do I somehow have to apply the CRS to be taken into account during spatial querys?

I already tried manually creating spatial indices, but that yielded the same results.

Answer 1

As no one could answer me, I opened this bug report: https://issues.qgis.org/issues/16068

It looks like you need to set the ExactIntersect flag in the QgsFeatureRequest. I still hope that this is getting fixed (or at least documented).

The proposed fix does not work, but you can refine your results by checking the geometry intersection directly:

features = [feature for feature in features 
    if feature.geometry().intersects(QgsRectangle(mapCoords, mapCoords))]

Note that feature.geometry().contains(mapCoords) method does not work. I'll probably open another bug report for this.