Using gcps with rio.reproject

[oarcher]

Hi, this is my first question here:

I'd like to be able to use gcps for reproject, like rasterio is able to do.

My goal is to plot a sentinel-1 image on a map, with geoviews.

My first guess was to use:

rasterize(gv.load_tiff(filename))

But this doesn't works, because the tiff has no transform. However, it has some gcps, and it is correctly mapped in qgis.

I was not able to find any solution to do the reprojection with rioxarray, so I've ended in a pure rasterio solution, converting to rioxarray at the end:

# 100% rasterio
src = rasterio.open(filename)
gcps, src_crs = src.get_gcps()
destination = np.zeros((5000,5000)) 
dest, transform = rasterio.warp.reproject(src.read(1), destination, gcps=gcps, src_crs=src_crs, dst_crs=src_crs)

# to rioxarray
dn = xr.DataArray(dest, dims=['y','x'])
dn = dn.rio.write_crs(src_crs)
dn = dn.rio.write_transform(transform)
dn = dn.assign_coords(x=(transform * (dn.x, dn.y[0]))[0] , y=(transform * (dn.x[0], dn.y))[1]) 

# plot with geoviews
rasterize(gv.Image(dn))

My input image is pretty big (25000 * 15000), so ideally, the reprojection should be done with dask (#119)

[snowman2]

gcps support could definitely be improved.

Using the transform.from_gcps method in open_rasterio might be a path forward on that.

[snowman2]

@mraspaud also had interest in using gcps with rioxarray.

[mraspaud]

Yes I do, also for sentinel 1 data

[mraspaud]

To be clear, I need the gcps both for reading and writing. The software I'm developing for processing the Sentinel1 SAR data (satpy) needs to be able to write the gcps to the resulting geotiff file. At the moment I we do that by passing the gcps to the rasterio.open call done for writing.

[snowman2]

@oarcher the image link gives a 404.

[snowman2]

needs to be able to write the gcps to the resulting geotiff file

Is there a reason you need the gcps versus the transform? I ask because it would be relatively simple to get the transform/crs from the gcps and then use that within rioxarray.

Portential method for rioxarray to get transform from gcps in the file in open_rasterio:

import rasterio

file_path = "s1a-iw-grd-vh-20160121t001321-20160121t001346-009585-00df26-002.tiff"
with rasterio.open(file_path) as rds:
    gcps, crs = rds.get_gcps()
    transform = rasterio.transform.from_gcps(gcps)

Current method for writing transform from gcps:

import rioxarray

rds = ...
gcps = [...]
transform = rasterio.transform.from_gcps(gcps)
rds.rio.write_transform(transform, inplace=True)

[oarcher]

@snowman2 , I've found that using rasterio.transform.from_gcps is not very accurate on big image with many gcps.

For example, we can have a 600 meters errors when applying the transform to the same gcps (for a 10m pixel, and a image 15000*25000). Probably because the transorm is affine, and very small error on coefficients gave large error on large image.

Using bi-linear interpolation from the gcps give more accurate results.

[oarcher]

@oarcher the image link gives a 404.

oops, fixed.

[snowman2]

That is good to know. I wonder how to accurately generate coordinates for a data array from gcps?

[scottyhq]

Generally gdalwarp with the thin plate spline option (-tps) will give better geocoding results for Sentinel1 images (but still approximate). I don't think -tps is exposed in rasterio... Another workaround is gdalwarp -of VRT tps -srcnodata 0 -dstnodata 0 INPUT.grd INPUT_wgs84.vrt and then open the VRT with rioxarray.

It would be nice though to have some sort of from_gcps keyword or convience function to simplify the process for GCP fit -> rioxarray transform.

[mraspaud]

IIRC, the Sentinel1 documentation recommends bilinear interpolation to expand the GCPs to full lon/lat arrays. We do this in satpy using a daskified version of scipy's interp2d function right now.

However, Sentinel 1 is not the only use case for GCPs and the GDAL documentation says:

The GDAL data model does not imply a transformation mechanism that must be generated from the GCPs … this is left to the application. However 1st to 5th order polynomials are common.

( https://gdal.org/user/raster_data_model.html#gcps )

So in my opinion, rioxarray should expose the GCPs but leave the interpolation scheme to the application. How to do the interpolation (what coordinate system, which order, which interpolation) is highly dependent both on the data at hand and the application the user has, so we would be opening a can of worms here if this was to be implemented in rioxarray. Just my 2 cents.

[snowman2]

Have y'all tried WarpedVRT for resampling datasets with gcps?

• https://corteva.github.io/rioxarray/stable/examples/reproject.html#Reproject-Large-Rasters-with-Virtual-Warping
• https://rasterio.readthedocs.io/en/latest/topics/virtual-warping.html

from rasterio.enums import Resampling
from rasterio.vrt import WarpedVRT
import rasterio

with rasterio.open("m_3907617_ne_18_1_20130924.tif") as src:
    with WarpedVRT(src, crs=src.crs, resampling=Resampling.bilinear) as vrt:
        rds = rioxarray.open_rasterio(vrt)

@oarcher this should enable loading in the image lazily with dask. Here is an example: ref.

[snowman2]

Nevermind. Seems like WarpedVRT doesn't support gcps. Might be worthwhile to open up an issue with rasterio about WarpedVRT and gcps/tps suppport. Then, we can update rioxarray to the changes made by rasterio.