Translate and Rotate features (without skew)

Question

I want to transform features from one CRS to another, but without the skew that is normally associated with it, as I need to retain the feature shapes ( = rectangles) within specific boundary regions. A minor loss in accuracy is to be expected and acceptable.

My approach:

1. Calculate translate and rotate parameters for the centroid of each boundary polygon using regular QGis transformations
2. Apply the specific translate+rotate parameters to each feature that lies within the boundary.

I wrote the following custom QGIS3 processing script. I include only the interesting part for brevity. If more info about the other stuff is needed, please ask. This is based on the QGIS3 processing script example. It is not optimized for performance.

def processAlgorithm(self, parameters, context, feedback):
    source = self.parameterAsSource(parameters, self.INPUT, context)
    target_crs = self.parameterAsCrs(parameters, 'crs_target', context)
    src_extents = self.parameterAsSource(parameters, self.EXTENT_POLY, context)
    (sink, dest_id) = self.parameterAsSink(parameters, self.OUTPUT,
        context, source.fields(), source.wkbType(), target_crs)

    # Compute the number of steps to display within the progress bar and
    # get features from source
    total = 100.0 / source.featureCount() if source.featureCount() else 0


    features = source.getFeatures()
    feat_extents = src_extents.getFeatures()

    translat = []
    rotate =  []
    for feat_ext in feat_extents:
        # Get the centroid of each extent
        centroid = feat_ext.geometry().centroid().asPoint()
        tr = QgsCoordinateTransform(source.sourceCrs(), target_crs, QgsCoordinateTransformContext ())
        # Transform the centroid to the target CRS
        centroid_trg = QgsGeometry.fromPointXY(centroid)
        centroid_trg.transform(tr)
        centroid_trg = centroid_trg.asPoint()
        # Generate a dummy point to calculate the angle, distance doesn't matter
        pt_rotate_src = QgsPointXY(centroid.x()+ 10, centroid.y())
        # Do the same for the target
        pt_rotate_trg = QgsGeometry.fromPointXY(pt_rotate_src)
        pt_rotate_trg.transform(tr)
        pt_rotate_trg = pt_rotate_trg.asPoint()

        # Actually calculate the angles (in radians!)
        angle_rotate_src = math.atan2(centroid.y() - pt_rotate_src.y(), centroid.x() - pt_rotate_src.x())
        angle_rotate_trg = math.atan2(centroid_trg.y() - pt_rotate_trg.y(), centroid_trg.x() - pt_rotate_trg.x())

        # Get the difference, and convert it to degrees, as expected by QGIS rotate(). Also, gather the rotation center to be used.
        angle_rotate = angle_rotate_src - angle_rotate_trg
        rotate.append((math.degrees(angle_rotate),centroid))

        # The translation parameters for this extent.
        translat.append((
            centroid_trg.x() - centroid.x(),
            centroid_trg.y() - centroid.y()
        ))

    # Now we treat the actual features
    for current, feature in enumerate(features):
         #Stop the algorithm if cancel button has been clicked
        if feedback.isCanceled():
            break
        # Must reload this every time, because the iterator doesn't reset otherwise
        feat_extents = src_extents.getFeatures()
        for feat_ext, translat_ext, rotate_ext in zip(feat_extents, translat, rotate):
            # Figure out in which boundary this feature lies
            if feature.geometry().within(feat_ext.geometry()):
                # Get the feature, rotate it first ( based on centroid)
                # Then translate it and write it to output
                geom = feature.geometry()
                geom.rotate(*rotate_ext)
                geom.translate(*translat_ext)
                feat_translated = QgsFeature(feature)
                feat_translated.setGeometry(geom)
                sink.addFeature(feat_translated, QgsFeatureSink.FastInsert)
                break

        # Update the progress bar
        feedback.setProgress(int(current * total))

    # Return the results of the algorithm. In this case our only result is
    # the feature sink which contains the processed features, but some
    # algorithms may return multiple feature sinks, calculated numeric
    # statistics, etc. These should all be included in the returned
    # dictionary, with keys matching the feature corresponding parameter
    # or output names.
    return {self.OUTPUT: dest_id}

This script works well, except all the features are shifted slightly off-target. I "transform" from EPSG:31466 to EPSG:25832 (DHDN to UTM, very common task in Germany/EU right now).

The rotation is near perfect, as expected, but the translation is slightly off. Which I interpret as either the translation or the rotation center being off, but I can't figure out how and why that would be the case. Note that the magnitude of the translation is very large, as these two CRS are entirely different. Translation sets (x,y) look like this:

-2211678.9598990344, 2175.2763667702675

The error is different for each boundary (but uniform within), typically between 0.5 and 2 meters horizontally and 0.1-0.5 vertically.

The following image shows an example. Purple is the boundary, yellow are the original features, and cyan are the translated&rotated features. The cyan layer is in UTM32N, while the main map&project are in DHDN, and QGIS uses On-the-Fly-Reprojection to show them matching up.

EDIT: Further investigation shows the offset ranging from about ~1.75m to 1.85m over all the boundaries, always eastwards and very slightly to the north.

Answer 1

After further research, I found the culprit. In summary, the issue was caused by my unfamiliarity with all the QGIS3 novelties.

The error lies in this line:

tr = QgsCoordinateTransform(source.sourceCrs(), target_crs, QgsCoordinateTransformContext ())

The API docs for the QgsCoordinateTransform constructors clearly say, even in a big warning label:

Warning

Do NOT use an empty/default constructed QgsCoordinateTransformContext() object when creating QgsCoordinateTransform objects. This prevents correct datum transform handling and may result in inaccurate transformations. Always ensure that the QgsCoordinateTransformContext object is correctly retrieved based on the current code context, or use the constructor variant which accepts a QgsProject argument instead.

Hence, the solution is to alter said line to:

tr = QgsCoordinateTransform(source.sourceCrs(), target_crs, QgsProject.instance())

Now the offsets are in the sub-centimeter range, as I would have expected.

I still don't fully understand why the project instance suddenly makes everything work (where does it take the correct transform rules from? why aren't those set for the default context?), but at least it works. If anyone can enlighten on this, I'll gladly switch the accepted answer around :) .