I want to perform a QgsCoordinateTransform on a QgsGeometry object and do some operations on this afterwards. I am doing this in a processing script and am getting the original QgsGeometry object from a QgsSpatialIndex with flags=QgsSpatialIndex.FlagStoreFeatureGeometries set.
So basically I am doing
overlay_geom = QgsGeometry(overlay_layer_idx.geometry(id))
overlay_geom.transform(QgsCoordinateTransform(source_layer_crs, overlay_layer_crs, QgsProject.instance()))
nearest_point_geom = overlay_geom.nearestPoint(source_geom)
and it causes QGIS to crash with no informations provided. Tested on 3.16 and 3.28.
How can I transform my QgsGeometry object and do an operation on this afterwards?
I'd have a workaround by simply using the native reproject algorithm, but want to avoid it, because it makes the parameterAsSource input unusable and forces to use parameterAsLayer instead.
Here is the script I am actually using narrowed down to the essentials. Its working fine when both inputs have the same CRS, but will crash if not.
Edit: Just figured out it is not actually the transformation causing the crash but the next operation on the transformed geometry. And it is only crashing when reprojecting to 4326. But it does work with 4326 when both inputs have this crs and therefore no transformation is done.
# -*- coding: utf-8 -*-
import processing, math
from PyQt5.QtCore import QCoreApplication, QVariant
from qgis.core import (QgsProject, QgsFeature, QgsProcessing, QgsSpatialIndex, QgsGeometry, QgsCoordinateReferenceSystem, QgsCoordinateTransform,
QgsFeatureSink, QgsFeatureRequest, QgsProcessingAlgorithm, QgsProcessingParameterFeatureSink, QgsProcessingParameterFeatureSource)
class TransformTest(QgsProcessingAlgorithm):
SOURCE_LYR = 'SOURCE_LYR'
OVERLAY_LYR = 'OVERLAY_LYR'
OUTPUT = 'OUTPUT'
def initAlgorithm(self, config=None):
self.addParameter(
QgsProcessingParameterFeatureSource(
self.SOURCE_LYR, self.tr('Points'), [QgsProcessing.TypeVectorPoint]))
self.addParameter(
QgsProcessingParameterFeatureSource(
self.OVERLAY_LYR, self.tr('Lines'), [QgsProcessing.TypeVectorLine]))
self.addParameter(
QgsProcessingParameterFeatureSink(
self.OUTPUT, self.tr('Test')))
def processAlgorithm(self, parameters, context, feedback):
source_layer = self.parameterAsSource(parameters, self.SOURCE_LYR, context)
source_layer_vl = self.parameterAsLayer(parameters, self.SOURCE_LYR, context)
overlay_layer = self.parameterAsSource(parameters, self.OVERLAY_LYR, context)
overlay_layer_vl = self.parameterAsLayer(parameters, self.OVERLAY_LYR, context)
output_layer_fields = source_layer.fields()
(sink, dest_id) = self.parameterAsSink(parameters, self.OUTPUT, context,
output_layer_fields, 2, # LineString = 2
source_layer.sourceCrs())
total = 100.0 / source_layer.featureCount() if source_layer.featureCount() else 0
source_layer_crs = QgsCoordinateReferenceSystem(source_layer_vl.crs().authid())
overlay_layer_crs = QgsCoordinateReferenceSystem(overlay_layer_vl.crs().authid())
overlay_layer_idx = QgsSpatialIndex(overlay_layer.getFeatures(), flags=QgsSpatialIndex.FlagStoreFeatureGeometries)
for current, source_feat in enumerate(source_layer.getFeatures()):
if feedback.isCanceled():
break
nearest_lines = overlay_layer_idx.nearestNeighbor(source_feat.geometry().centroid().asPoint(), neighbors = 1, maxDistance = -1)
for nearest_line_id in nearest_lines:
new_feat = QgsFeature(output_layer_fields)
attridx = 0
for attr in source_feat.attributes():
new_feat[attridx] = attr
attridx += 1
nearest_line_geom = overlay_layer_idx.geometry(nearest_line_id)
if source_layer_vl.sourceCrs() != overlay_layer_vl.sourceCrs():
feedback.setProgressText('Doing transform now')
nearest_line_geom.transform(QgsCoordinateTransform(source_layer_crs, overlay_layer_crs, context.transformContext()))
feedback.setProgressText('You see me, so the transform did not cause the crash')
feedback.setProgressText('Will crash in 3.. 2.. 1.. go:')
point_on_nearest_line = nearest_line_geom.nearestPoint(source_feat.geometry().centroid())
feedback.setProgressText('You will not see me, because it already crashed')
sqrDist, minDistPoint, afterVertex, leftOf = nearest_line_geom.closestSegmentWithContext(point_on_nearest_line.asPoint(),1)
vertexOnSegment2 = nearest_line_geom.vertexAt(afterVertex)
vertexOnSegment1 = nearest_line_geom.vertexAt(afterVertex - 1)
segmentAngle = math.atan2(vertexOnSegment2.x() - vertexOnSegment1.x(), vertexOnSegment2.y() - vertexOnSegment1.y())
segmentAngleDegree = math.degrees(segmentAngle) if segmentAngle > 0 else math.degrees(segmentAngle) + 180
segmentgeom = QgsGeometry.fromPolyline([vertexOnSegment1,vertexOnSegment2])
perpendicularLinePoint1 = point_on_nearest_line.asPoint().project(10,segmentAngleDegree+90)
perpendicularLinePoint2 = point_on_nearest_line.asPoint().project(10,segmentAngleDegree-90)
perpendicularLineGeom = QgsGeometry.fromPolylineXY([perpendicularLinePoint1,perpendicularLinePoint2])
new_feat.setGeometry(perpendicularLineGeom)
sink.addFeature(new_feat, QgsFeatureSink.FastInsert)
feedback.setProgress(int(current * total))
return {self.OUTPUT: dest_id}
def tr(self, string):
return QCoreApplication.translate('Processing', string)
def createInstance(self):
return TransformTest()
def name(self):
return 'TransformTest'
def displayName(self):
return self.tr('Transform Test')
def group(self):
return self.tr(self.groupId())
def groupId(self):
return 'Test'
def shortHelpString(self):
return self.tr('Testing Transform')
