Seems to be a shapely problem: even (a := LineString(((0,0), (0,1), (1,1)))).intersection(a) results in a MultiLineString.

Thank you, this was indeed almost the solution. In order to also obtain the linework of the outer boundary, I added an additional row with geometry df.union_all().boundary to the boundary dataframe. However, many of the resulting geometries are actually MultiLineStrings. This can be fixed with a simple df.geometry = df.geometry.line_merge(), but I don't know why it is happening in the first place.

Be aware that rasterio.transform.from_gcps tries to fit an affine linear transformation based on the GCP information. The precise transformation however, will usually not be linear. Instead of using the transform attribute of the dataset, you can store your list of GCPs in the gcps attribute to save the actual GCP information to the GeoTiff. Opening such a file with QGIS for example will then perform a more accurate reprojection based on thin plate splines.