Description
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
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()
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