Skip to content

This issue was moved to a discussion.

You can continue the conversation there. Go to discussion →

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Enable setting gcps in Reader option #564

Closed
vincentsarago opened this issue Jan 4, 2023 · 1 comment
Closed

Enable setting gcps in Reader option #564

vincentsarago opened this issue Jan 4, 2023 · 1 comment

Comments

@vincentsarago
Copy link
Member

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

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
Copy link
Member Author

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:

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

@cogeotiff cogeotiff locked and limited conversation to collaborators Jan 9, 2023
@vincentsarago vincentsarago converted this issue into discussion #565 Jan 9, 2023

This issue was moved to a discussion.

You can continue the conversation there. Go to discussion →

Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant