Skip to content
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

rasterio.vrt.WarpedVRT warp_extras ignored by rioxarray #598

Closed
scottyhq opened this issue Oct 31, 2022 · 5 comments · Fixed by #599
Closed

rasterio.vrt.WarpedVRT warp_extras ignored by rioxarray #598

scottyhq opened this issue Oct 31, 2022 · 5 comments · Fixed by #599
Labels
bug Something isn't working
Milestone

Comments

@scottyhq
Copy link
Contributor

#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
@scottyhq scottyhq added the bug Something isn't working label Oct 31, 2022
@snowman2
Copy link
Member

Looks like it is being passed in:

warp_extras=vrt.warp_extras,

@snowman2
Copy link
Member

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

@scottyhq
Copy link
Contributor Author

@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'}}

@snowman2
Copy link
Member

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?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug Something isn't working
Projects
None yet
Development

Successfully merging a pull request may close this issue.

2 participants