Rasterio reproject using from_gcps moves raster boundaries but not pixels

Question

Raster was misaligned with ground truth (as observed in google earth) so I used rasterio.warp.reproject with some manually obtained GCPs to try to align it. This is the code to do so (somewhat abridged):

import rasterio
import rasterio.crs
from rasterio.control import GroundControlPoint
from rasterio.transform import from_gcps
from rasterio.warp import reproject, Resampling
from rasterio.crs import CRS

if __name__ == '__main__':
    # Open the input TIFF file
    with rasterio.open(input_img_path) as src:
        gcps = [
            GroundControlPoint(col=(...),row=(...),x=(...),y=(...)),
            (...)
        ] # x and y are lon and lat extracted from Google Earth
        transform = from_gcps(gcps)
        gcp_crs = CRS.from_epsg(4326)  # This is the crs of the GCPs
        # Warp the image
        kwargs = src.profile.copy()
        kwargs.update(transform=transform, crs=gcp_crs)
        with rasterio.open(output_img_path, "w", **kwargs) as dst:
            for i in range(1, src.count + 1):
                reproject(
                    source=rasterio.band(src, i),
                    destination=rasterio.band(dst, i),
                    src_transform=src.transform,
                    src_crs=src.crs,
                    dst_transform=dst.transform,
                    dst_crs=dst.crs,
                    resampling=Resampling.nearest
                )

However, what happened was that the image bounds were realigned correctly while the image pixels were kept in the same positions. These images in QGIS illustrate what I mean:

Before any changes Reprojected, black square is in the correct placement Original and reprojected, shows that pixels were kept in the same positions, only bounding box changed

I have tried including the gcps directly in the reproject function. I have tried keeping the same crs as source image. Tried different combinations of gcps.

Appreciate any suggestions, thanks!

Answer 1

In the end this is what worked for anyone who has the same issue in the future:

for i in range(1, src.count + 1):
source = src.read(i)
dest = np.zeros_like(source)
reproject(
    source=source,
    destination=dest,
    gcps=gcps,
    src_crs=gcp_crs,
    dst_transform=dst.transform,
    dst_crs=dst.crs,
    resampling=Resampling.nearest
)
dst.write(dest, indexes=i)

The idea was to not provide it the source transform since that is what it seems it was using to place the pixels at the destination.