Does GeoUtils reprojection handle GCPs and RPCs ?

[remi-braun]

Hi there!

Thanks for your open source library that seems very promising.
I have one question about your reprojection implementation (which leverages dask if I'm right): does it handle GCPs and RPCs ?

[rhugonnet]

Hi @remi-braun,

The short answer is not yet, but soon!

Dask support: We indeed finalized our Dask reprojection earlier this year, which is fully tested. It is available here: https://github.com/GlacioHack/geoutils/blob/main/geoutils/raster/delayed.py#L674. But it is not yet public, it will be once we finalize our Xarray rst accessor mirroring the Raster class (in #446) that is well on its way (probably done this spring). The Raster class itself uses NumPy masked-arrays, which is practical to save RAM when using integer-type data, but does not have underlying Dask support.

GCPs/RPCs: The GCPs/RPCs are only treated as keyword argument right now, they are not discriminated depending on chunks. So they would be passed to all chunks. I'm not so familiar with the use of GCPs/RCPs during reprojection, do you know if passing them as keyword argument to all chunks would be sufficient to work as expected? Or would we need to split them in the chunk metadata somehow? In any case, it feels like it wouldn't be a big task to add support for this.
Probably just passing and adjusting the proper additional kwargs here.

Running your delayed reprojection: In the meantime, if you want to use/adapt the delayed_reproject functionality before our Xarray accessor is finalized, you'd have to open your raster chunked via Xarray, extract the da.Array, pass it to delayed_reproject, than re-build the Xarray dataset.

This should look something like this:
(I did not test this code, but based it off different working routines either in the tests of the delayed reproject at https://github.com/GlacioHack/geoutils/blob/main/tests/test_raster/test_delayed.py#L391, or in the Xarray accessor PR mentioned above, so should require only minor tweaks)

import xarray as xr
import rioxarray
from geoutils.raster.delayed import delayed_reproject

# Open raster as delayed Dask array through Xarray
filename = "myexample.tif"
chunksizes = (200, 200)
# Or maybe rioxarray.open_rasterio here, not sure if it changes anything
ds = xr.open_dataset(filename, chunks={"x": chunksizes[0], "y": chunksizes[1])

# Specify your output CRS/transform/shape
dst_transform = 
dst_shape = 
dst_crs = 
dst_nodata = 

# Call delayed_reproject on dask array stored in ds.data
da_reproj = delayed_reproject(
        darr=ds.data,
        src_transform=ds.rio.transform(),
        src_crs=ds.rio.crs,
        dst_transform=dst_transform,
        dst_shape=dst_shape,
        dst_crs=dst_crs,
        resampling=Resampling.bilinear,
        src_nodata=ds.rio.nodata,
        dst_nodata=dst_nodata,
        dst_chunksizes=None,  # reproject will use the destination chunksizes if set to None.
    )

# Re-build the Xarray object with delayed output before saving

# First, get netCDF coordinates from transform and shape
from rioxarray.rioxarray import affine_to_coords
coords = affine_to_coords(affine=dst_transform, width=dst_shape[0], height=dst_shape[1])

# Copy raster tags/attributes if useful (GCPs, RCPs, potentially updated?)
tags = ds.attrs

# Build a data array
ds_reproj = xr.DataArray(
    data=da_reproj,
    coords={"x": coords["x"], "y": coords["y"]}, 
    attrs=tags,
)

# Set other attributes
ds_reproj.rio.write_transform(dst_transform)
ds_reproj.rio.set_crs(dst_crs)
ds_reproj.rio.set_nodata(dst_nodata)

then you can either perform more operations on the delayed Xarray object, or save to file in a delayed fashion:

write_delayed = ds_reproj.to_netcdf("myexample_reproj.tif", compute=False)
write_delayed.compute()

Once the rst accessor is out, it will be as simple as calling ds_reproj = ds.rst.reproject() after opening ds with chunks, all of this will be done in the background 😉.

We plan to explain clearly out-of-memory support in the documentation once the accessor is out (comparing the Raster class and the rst accessor), and add a table listing all functionality supported by Dask for the Xarray accessor.

[remi-braun]

Thanks for your prompt reply!
Thanks a lot for all your amazing work and the excellent effort you do in documenting what you do 🙏

However, I fear that RPC (I don't know for GCP) need more than just passing the keyword to all chunks, see errors I encountered here: opendatacube/odc-geo#193

I've also created a feature request in rioxarray to expose RPCs instead of using rasterio directly. (see corteva/rioxarray#837)

I see this is tricky, and if not doable with Dask easily I'll stay with current rioxarray implementation.

[rhugonnet]

Perfect! I agree Rioxarray is the right place to expose the RPCs.
Supporting RPCs during Dask-reprojection would have to live directly in ODC or GeoUtils (Rioxarray likely won't support a Dask implementation directly, as the mapping of chunks creates a dependency to Shapely).

After looking into it a bit more:

Adding GCP support looks quite straightforward, potentially just estimating the source transform ahead of the Dask reprojection (a couple of lines): https://rasterio.readthedocs.io/en/stable/api/rasterio.transform.html#rasterio.transform.from_gcps.

Adding RPCs indeed seems more complex, and I would need to inform myself more to be certain. But, from what I know, at least the RPCs are consistent on the raster extent. And I see only two steps to add support:
1/ Applying to chunks: This seems like it could be fairly easy. We'd just need to "trick" reproject in applying the RPCs for a shifted row/col numbering, which could be as simple as shifting the row/col offsets of the RPCs by the chunk size location.
2/ Mapping input to output chunks: This seems slightly more tricky. My intuition would be to build vector tiles from the 4 corner x/y of each input chunk using their row/col (using xy() in Rasterio with the RPCs). Then map the input/output chunks the same way it's already done.

The good news is that a Dask implementation is easy to test: We can just compare the output raster processed in chunks to the one processed all at once in-memory! 🙂
If you ever feel inspired to add this, feel free to open a PR, happy to help! I have some colleagues quite familiar with RPCs (as I work a lot with stereo-imagery), that could potentially provide input.

And thanks for all your work with EOReader! 😉

[remi-braun]

I'd definitely choose RPC solution 1 👍
I cannot promise any PR, but I'll do it if feasible with my schedule!

And thank you for your feedback 🥰. You looked at EOReader?

[rhugonnet]

To clarify: we'd have to add both steps 1/ and 2/ to support it.
I'm tagging @adehecq to see if he has more insights/ideas on RPCs (might take a bit of time to react, just back from fieldwork).

If you do feel motivated for a PR, take your time. I'm leaving for more than a month, so I will react at the earliest end-of-January.

Yes, curious to try EOReader once I get back to doing more remote sensing again! My work of the past years has been more methodological, but I normally use a lot of these sensors 😉