Enable setting gcps in Reader option

[vincentsarago]

Let's say you have an Non-geo image (as COG already) and some external GCP (e.g from a Georeferencing UI), it would be nice to pass those GCPS to rio_tiler.io.Reader directly to allow virtual georeferencing

rio-tiler/rio_tiler/io/rasterio.py

Lines 98 to 110 in 93cf68e

	def __attrs_post_init__(self):
	"""Define _kwargs, open dataset and get info."""
	if not self.dataset:
	dataset = self._ctx_stack.enter_context(rasterio.open(self.input))
	if dataset.gcps[0]:
	self.dataset = self._ctx_stack.enter_context(
	WarpedVRT(
	dataset,
	src_crs=dataset.gcps[1],
	src_transform=transform.from_gcps(dataset.gcps[0]),
	)
	)
	else:

[vincentsarago]

After couple hours digging, it seems there is no easy solution.
Notes:

• WarpedVRT will no work if the dataset do not have internal GCPS
• I was hopping we could pass a list of gcps to WarpedVRT but it's not possible (the underneath GDAL function won't get it either)

The only solution I found was to create a VRT

import xml.etree.ElementTree as ET
from rasterio.vrt import _boundless_vrt_doc
from rasterio.warp import WarpedVRT

with rasterio.open("image.jpeg") as dataset:
    # Use rasterio boundless_vrt_doc to create a VRT xml document
    vrt_xml = _boundless_vrt_doc(dataset)

    # Update the VRT to add the GCPS
    vrtdataset = ET.fromstring(vrt_xml)
    gcp_list = ET.SubElement(vrtdataset, 'GCPList')
    gcp_list.attrib['Projection'] = str(gcps_crs)
    for gcp in self.gcps:
        g = ET.SubElement(gcp_list, 'GCP')
        g.attrib["Id"] = gcp.id
        g.attrib['Pixel'] = str(gcp.col)
        g.attrib['Line'] = str(gcp.row)
        g.attrib['X'] = str(gcp.x)
        g.attrib['Y'] = str(gcp.y)
    vrt_xml = ET.tostring(vrtdataset)

# Open the VRT Dataset and pass it to the WarpedVRT
with rasterio.open(vrt_xml.decode()) as dataset:
    with WarpedVRT(
        dataset,
        src_crs=dataset.gcps[1],
        src_transform=transform.from_gcps(dataset.gcps[0]),
    ) as vrt:
        ...

Notes:

• I'm not sure what _boundless_vrt_doc does so we might not to check if there are any specifics
• in GDAL if we wanted to allow passing GCPS I thing we would need to update https://github.com/OSGeo/gdal/blob/e641b1356a3ac2ef9f0676fe75f5dc95162afc36/frmts/vrt/vrtwarped.cpp#L139-L298 to allow GDALAutoCreateWarpedVRTEx to set the projection from external GCPS instead of https://github.com/OSGeo/gdal/blob/e641b1356a3ac2ef9f0676fe75f5dc95162afc36/frmts/vrt/vrtwarped.cpp#L292

Ps: I'm going to convert the issue to discussion because there Is no real actionable right now

[mradamcox]

I'd like to add one idea here (a complication, not a solution!), which is that I think the best implementation would also allow passing in a region extent, such that the GCPs could be applied only to a clipped subset of the original image. This would allow various pieces of the same image (think, a scanned atlas sheet with multiple different maps on it, e.g. https://www.loc.gov/resource/g4014nm.g03376188501/?sp=3) to be virtually georeferenced without having to do one of the following:

1. cut the original file into two separate files, virtually warp each separately, or
2. virtually warp the entire image multiple times with different sets of GCPs, and later apply a polygon mask to obscure the irrelevant part of the sheet.

I'm not sure if this makes much sense in terms of how COGs are structured and read, but I wanted to put it down here for consideration.

[vincentsarago]

FYI, I think I have something which work quite ok. I ended up creating a temporary VRT file (using a custom version of rasterio's _boundless_vrt_doc: https://github.com/developmentseed/titiler-image/blob/a898fa0dc7f237e70ce2758772730df6a927760b/titiler/image/reader.py#L86-L216)

I think there are still some edge cases about nodata/alpha/mask that needs more tests but you can see how it works in https://github.com/developmentseed/titiler-image/blob/main/docs/notebooks/titiler_image_GCPSReader.ipynb

[vincentsarago]

I also added a cutline option to be applied before warping (in the local image space)

[mradamcox]

Woah, incredible, I'm really looking forward to trying this out! Somewhat related to this, the IIIF-Maps committee has been working on a IIIF spec extension for holding georeferencing/GCP information: https://preview.iiif.io/api/georef/extension/georef/... I see now you used resourceCoords as a property name so I have a feeling you have seen this document already :)

[mradamcox]

Haha yes, I hope to improve that PR a little bit, but the topic of GCPs + projection information brings a lot of interesting ideas. I realize you can get much better warped results if you first chose a projection that is close to the original historical map's projection, and then reproject that to 3857 for web tiling. Seems like it could be done with just a chain of VRTs (I already chain two together, so it would just be third I guess??), or perhaps all in the same function.

[vincentsarago]

👀

[mradamcox]

I never followed up on this, and now it is May... but here is a demo gif I made a little while ago.

To explain: First, three GCPs are created using Web Mercator (the basemap projection), and a 1st order polynomial transformation is used to georeference the image. Distortion is obvious.

Second, a slightly more complex process is used behind the scenes:

• 1. the same 3 GCPs are projected to Lambert North America, a projection that closely fits the original input image,
• 2. 1st order polynomial transformation is applied, resulting in a VRT with CRS of Lambert, and then
• 3. this VRT is reprojected to Web Mercator, in order to properly display on a web map.

You can see, this is a much better way to handle maps that aren't originally in a projection that is close to Web Mercator.

It's very slick that this can all be done with a chain of VRTs, no extra GeoTIFFs are made (in fact, no GeoTiffs at all, in this specific case).

[vincentsarago]

🤯 that's really cool!

are you using https://github.com/developmentseed/titiler-image/blob/a898fa0dc7f237e70ce2758772730df6a927760b/titiler/image/reader.py or you ended up writing your own?

ps I edited the format of your comment to make it a bit more clear 😄

[mradamcox]

No, in the illustration here, the preview is served with MapServer. I plan to switch that over to Titiler eventually, to match what I now use in the rest of the app, but MapServer has worked like a dream for this live preview (alter/add a GCP, push them into the VRTs, refresh the layer and that's it).

Here is the commit where I added the projection switcher you see above: ohmg-dev/OldInsuranceMaps@996b4cc