Bug in WarpedVRT support of open_rasterio()

[jmichel-otb]

Code Sample, a copy-pastable example if possible

Using the following data:

https://gitlab.orfeo-toolbox.org/orfeotoolbox/otb/blob/develop/Data/Input/QB_Toulouse_Ortho_XS.tif

import rasterio as rio
from rasterio.crs import CRS
from rasterio.warp import calculate_default_transform,aligned_target
from rasterio.enums import Resampling
from rasterio.vrt import WarpedVRT

import xarray as xr

path = 'QB_Toulouse_Ortho_XS.tif'

# Read input metadata with rasterio                                                                                                                                                                                                           
with rio.open(path) as src:
    print('Input file CRS is {}'.format(src.crs))
    print('Input file shape is {}'.format(src.shape))
    print('Input file transform is {}'.format(src.transform))
    # Create a different CRS                                                                                                                                                                                                                  
    dst_crs = CRS.from_epsg(2154)
    # Compute a transform that will resample to dst_crs and change resolution to dst_crs                                                                                                                                                      
    transform, width, height = calculate_default_transform(
        src.crs, dst_crs, src.width, src.height,resolution=20, *src.bounds)
    # Fill vrt options as shown in https://rasterio.readthedocs.io/en/stable/topics/virtual-warping.html                                                                                                                                      
    vrt_options = {
    'resampling': Resampling.cubic,
    'crs': dst_crs,
    'transform': transform,
    'height': height,
    'width': width
    }
    # Create a WarpedVRT using the vrt_options                                                                                                                                                                                                
    with WarpedVRT(src,**vrt_options) as vrt:
        print('VRT shape is {}'.format(vrt.shape))
        # Open VRT with xarray                                                                                                                                                                                                                
        ds = xr.open_rasterio(vrt)
        # Shape does not match vrt shape:                                                                                                                                                                                                     
        print(ds)

Output:

$ python test_rio_vrt.py
Input file CRS is EPSG:32631
Input file shape is (500, 500)
Input file transform is | 0.60, 0.00, 374149.98|
| 0.00,-0.60, 4829183.99|
| 0.00, 0.00, 1.00|
VRT shape is (16, 16)
<xarray.DataArray (band: 4, y: 500, x: 500)>
[1000000 values with dtype=int16]
Coordinates:
  * band     (band) int64 1 2 3 4
  * y        (y) float64 6.28e+06 6.28e+06 6.28e+06 ... 6.279e+06 6.279e+06
  * x        (x) float64 5.741e+05 5.741e+05 5.741e+05 ... 5.744e+05 5.744e+05
Attributes:
    transform:   (0.6003151072155879, 0.0, 574068.2261249251, 0.0, -0.6003151...
    crs:         EPSG:2154
    res:         (0.6003151072155879, 0.6003151072155879)
    is_tiled:    0
    nodatavals:  (nan, nan, nan, nan)

Problem description

In the above example, xarray.open_rasterio() is asked to read a WarpedVRT created with rasterio. This WarpedVRT has custom transform, width and height, which is a very common use case of WarpedVRT where you upsample or downsample your data on the fly during read. Las, open_rasterio() ignores those custom attributes, which result in a wrong DataArray shape:

VRT shape is (16, 16)
<xarray.DataArray (band: 4, y: 500, x: 500)>

Expected Output

Correct output should be:

VRT shape is (16, 16)
<xarray.DataArray (band: 4, y: 16, x: 16)>

Which is fairly easy to obtain by modifying the following lines in open_rasterio():

xarray/xarray/backends/rasterio_.py

Line 222 in 0c73a38

vrt_params = dict(crs=vrt.crs.to_string(),

With the following:

      vrt_params = dict(crs=vrt.crs.to_string(),
                          resampling=vrt.resampling,
                          src_nodata=vrt.src_nodata,
                          dst_nodata=vrt.dst_nodata,
                          tolerance=vrt.tolerance,
                          # Edit                                                                                                                                                                                                              
                          transform=vrt.transform,
                          width=vrt.width,
                          height=vrt.height,
                          # End edit                                                                                                                                                                                                          
                          warp_extras=vrt.warp_extras)

I can provide this patch in a pull request if needed.

Output of xr.show_versions()

>>> xr.show_versions()

INSTALLED VERSIONS

commit: None
python: 3.6.7 | packaged by conda-forge | (default, Nov 21 2018, 03:09:43)
[GCC 7.3.0]
python-bits: 64
OS: Linux
OS-release: 3.10.0-862.2.3.el7.x86_64
machine: x86_64
processor: x86_64
byteorder: little
LC_ALL: None
LANG: fr_FR.utf8
LOCALE: fr_FR.UTF-8
libhdf5: 1.10.4
libnetcdf: 4.6.2

xarray: 0.11.3
pandas: 0.24.1
numpy: 1.16.1
scipy: 1.2.0
netCDF4: 1.4.2
pydap: None
h5netcdf: None
h5py: None
Nio: None
zarr: None
cftime: 1.0.3.4
PseudonetCDF: None
rasterio: 1.0.17
cfgrib: None
iris: None
bottleneck: None
cyordereddict: None
dask: 1.1.1
distributed: 1.25.3
matplotlib: 3.0.2
cartopy: 0.17.0
seaborn: 0.9.0
setuptools: 40.7.3
pip: 19.0.1
conda: None
pytest: 4.2.0
IPython: 7.1.1
sphinx: None

[fmaussion]

Thanks for the great report!

I can provide this patch in a pull request if needed.

Please, do ;-)

[jhamman]

cc @scottyhq

[jmichel-otb]

@fmaussion done.