rasterio.vrt.WarpedVRT warp_extras ignored by rioxarray

[scottyhq]

339 illustrated passing transform options to rio.reproject(), and in theory I think you should be able to pass those to rasterio.vrt.WarpedVRT() for example to open a file with GCPs and have the reprojection and resampling handled by GDAL upon reading.

But it seems options like SRC_METHOD: GCP_TPS are currently ignored:

Code Sample, a copy-pastable example if possible

import rasterio
import xarray as xr
import rioxarray
import pystac
import planetary_computer

# Get a test file

STAC = "https://planetarycomputer.microsoft.com/api/stac/v1/collections/sentinel-1-grd/items"
item_url = f"{STAC}/S1A_IW_GRDH_1SDV_20220320T230514_20220320T230548_042411_050E99"
stac_item = pystac.read_file(item_url)
item = planetary_computer.sign(stac_item)
url = item.assets["vh"].href
url # https://sentinel1euwest.blob.core.windows.net/s1-grd/GRD/2022/3/20/IW/DV/S1A_IW_GRDH_1SDV_20220320T230514_20220320T230548_042411_050E99_EDB0/measurement/iw-vh.tiff?st=2022-10-30T19%3A45%3A46Z&se=2022-11-07T19%3A45%3A47Z&sp=rl&sv=2021-06-08&sr=c&skoid=c85c15d6-d1ae-42d4-af60-e2ca0f81359b&sktid=72f988bf-86f1-41af-91ab-2d7cd011db47&skt=2022-10-31T19%3A45%3A45Z&ske=2022-11-07T19%3A45%3A45Z&sks=b&skv=2021-06-08&sig=xjX3K/OzXixf1MCx9TLhPcFwBsAKaxnN6XLTv80LkhQ%3D

# !gdalwarp -overwrite /vsicurl/"{url}" test-gcp.vrt
# (Creating output file that is 29343P x 28276L.)

# # -tps ==> -to SRC_METHOD=GCP_TPS
# !gdalwarp -overwrite -tps /vsicurl/"{url}" test-tps.vrt
# Creating output file that is 29342P x 28283L.

with rasterio.open(url) as src:
    print('original shape', src.shape)
    with rasterio.vrt.WarpedVRT(src,
                                **{'SRC_METHOD': 'GCP_TPS'}
                               ) as vrt:

        print('TPS shape', vrt.shape) # matches GDAL (28283, 29342) 
        assert vrt.warp_extras['SRC_METHOD'] == 'GCP_TPS'

        da = xr.open_dataarray(vrt, engine='rasterio')
        print(da.rio.shape) # matches default polynomial order 2 (ignores GCP_TPS)
        assert vrt.shape == da.rio.shape # AssertionError  

Expected Output

vrt.shape == da.rio.shape == (28283, 29342)

Environment Information

rioxarray (0.12.1) deps:
  rasterio: 1.3.2
    xarray: 2022.6.0
      GDAL: 3.5.2
      GEOS: 3.11.0
      PROJ: 9.0.1
 PROJ DATA: /srv/conda/envs/notebook/share/proj
 GDAL DATA: /srv/conda/envs/notebook/share/gdal

Other python deps:
     scipy: 1.9.1
    pyproj: 3.4.0

System:
    python: 3.10.6 | packaged by conda-forge | (main, Aug 22 2022, 20:35:26) [GCC 10.4.0]
executable: /srv/conda/envs/notebook/bin/python
   machine: Linux-5.4.0-1091-azure-x86_64-with-glibc2.35

[snowman2]

Looks like it is being passed in:

https://github.com/corteva/rioxarray/blob/041e6bdd81bf9378ff03d70380095df9cbb3c9df/rioxarray/_io.py#L897

[snowman2]

@scottyhq are you able to upload the file to this issue?

[scottyhq]

@snowman2 yes! should that just be **vrt.warp_extras ? otherwise I think the dictionary is passed and not expanded/recognized when the IO eventually happens:

'transform': None, 'dtype': None, 'warp_extras': {'SRC_METHOD': 'GCP_TPS', 'init_dest': 'NO_DATA'}}

[scottyhq]

i think the signed url above is publicly accessible for a week?

https://sentinel1euwest.blob.core.windows.net/s1-grd/GRD/2022/3/20/IW/DV/S1A_IW_GRDH_1SDV_20220320T230514_20220320T230548_042411_050E99_EDB0/measurement/iw-vh.tiff?st=2022-10-30T19%3A45%3A46Z&se=2022-11-07T19%3A45%3A47Z&sp=rl&sv=2021-06-08&sr=c&skoid=c85c15d6-d1ae-42d4-af60-e2ca0f81359b&sktid=72f988bf-86f1-41af-91ab-2d7cd011db47&skt=2022-10-31T19%3A45%3A45Z&ske=2022-11-07T19%3A45%3A45Z&sks=b&skv=2021-06-08&sig=xjX3K/OzXixf1MCx9TLhPcFwBsAKaxnN6XLTv80LkhQ%3D

[snowman2]

should that just be **vrt.warp_extras

Good catch! Yes, I think that is it. Did you want to submit a PR with the fix?