Reprojecting JPG with gcps that are not embedded in the file raises an error

[joostmeulenbeld]

I'm trying to reproject an aerial RGB image stored as a JPG, using ground control points. I thought this should be possible with rasterio.warp.reproject. Unfortunately, I didn't get this to work.

Below I posted a MWE of what I'm trying to do, including test image. When running this, a warning is shown during opening of the input JPG image: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix will be returned. However, as I'm using warp.reproject later with gcps, I would assume this warning should not be important.

Then, when calling warp.reproject, the following error is shown:rasterio._err.CPLE_AppDefinedError: Unable to compute a transformation between pixel/line and georeferenced coordinates for input.jpg. There is no affine transformation and no GCPs. Specify transformation option SRC_METHOD=NO_GEOTRANSFORM to bypass this check.

This error is weird because I'm supplying the GCPs with the reproject call.

Adding the suggested SRC_METHOD=NO_GEOTRANSFORM to the warp.reproject call, no error is raised, but the output becomes a fully 0-filled geotiff with the correct bounds and transform. (Also, I'd like to use SRC_METHOD=GCP_TPS)

Does anyone know how to do this?

MWE

The following code shows the error if it is run in the same folder as the attached JPG. Note that my real case is more complex with many GCPs per image, that cannot be represented with an Affine transformation.

import rasterio
from rasterio import warp

# Height and width of image in pixels
jpg_height = 323
jpg_width = 181

# Approximate ground coordinates in EPSG:28992 (RD-new)
left = 109881
top = 483106
right = 111960
bottom = 479345

gcps = [
    rasterio.control.GroundControlPoint(row=0, col=0, x=left, y=top),
    rasterio.control.GroundControlPoint(row=0, col=jpg_width - 1, x=right, y=top),
    rasterio.control.GroundControlPoint(row=jpg_height - 1, col=0, x=left, y=bottom),
    rasterio.control.GroundControlPoint(row=jpg_height - 1, col=jpg_width - 1, x=right, y=bottom),
]

transform, tif_width, tif_height = warp.calculate_default_transform(
    src_crs=rasterio.CRS.from_epsg(28992),  # Assuming this is the CRS of the GCP geographical coordinates
    dst_crs=rasterio.CRS.from_epsg(28992),
    width=jpg_width,
    height=jpg_height,
    gcps=gcps,
)
# Results:
# tif_height: 324
# tif_width: 179
# transform:
# | 11.65, 0.00, 109881.00|
# | 0.00,-11.65, 483106.00|
# | 0.00, 0.00, 1.00|

# So far so good: width and height approx. equal to input image, transform offsets match input, resolution seems about right.

output_profile = rasterio.profiles.DefaultGTiffProfile(
    crs=rasterio.CRS.from_epsg(28992),
    width=tif_width,
    height=tif_height,
    transform=transform,
    count=3,
)

# Opening the src jpg, the following warning is shown:
# NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix will be returned.
with rasterio.open("input.jpg") as src:
    with rasterio.open("output.tif", "w", **output_profile) as dst:
        # The warp.reproject call raises this error:
        # rasterio._err.CPLE_AppDefinedError: Unable to compute a transformation between pixel/line and georeferenced
        # coordinates for input.jpg. There is no affine transformation and no GCPs. Specify transformation option
        # SRC_METHOD=NO_GEOTRANSFORM to bypass this check.
        warp.reproject(
            rasterio.band(src, [1, 2, 3]),
            destination=rasterio.band(dst, [1, 2, 3]),
            src_crs=rasterio.CRS.from_epsg(28992),  # Again, assuming this CRS must match the GCP CRS
            gcps=gcps,
        )

[vincentsarago]

@joostmeulenbeld this might be of interest cogeotiff/rio-tiler#565

I wanted to do the same in rio-tiler/titiler, but I ended up having to create a custom VRT file https://github.com/developmentseed/titiler-image/blob/main/titiler/image/reader.py#L97
https://github.com/developmentseed/titiler-image/blob/main/docs/notebooks/titiler_image_GCPSReader.ipynb

I didn't want to have the file reprojected directly but only to have a virtual file which is why I used VRT

[joostmeulenbeld]

Thanks for the links! I think I would prefer calling out to gdalwarp and gdal_translate with a temporary file-based VRT instead of changing the VRT in-memory as XML. But if at some point I can use in-memory reprojection your links will offer a good solution.

[sgillies]

@joostmeulenbeld Rasterio supports GCPs for reprojection but only when the source is an array https://github.com/rasterio/rasterio/blob/main/tests/test_warp.py#L1647. GDAL's warper requires GCPs to be attached to a dataset and rasterio doesn't automatically wrap your input JPEG up in a VRT to do this.

[joostmeulenbeld]

I think I can get away with loading the entire input JPEG in memory in my use case, so using an array would work. (I do have to check memory usage of the reprojection in rasterio itself)

Implementing this in the MWE shows that indeed reprojection works as expected:

with rasterio.open("input.jpg") as src:
    data = src.read()

with rasterio.open("output.tif", "w", **output_profile) as dst:
    warp.reproject(
        data,
        destination=rasterio.band(dst, [1, 2, 3]),
        src_crs=rasterio.CRS.from_epsg(28992),  # Again, assuming this CRS must match the GCP CRS
        gcps=gcps,
    )

Thanks for the reply!