home / github / issues

Menu
  • Search all tables
  • GraphQL API

issues: 428300345

This data as json

id node_id number title user state locked assignee milestone comments created_at updated_at closed_at author_association active_lock_reason draft pull_request body reactions performed_via_github_app state_reason repo type
428300345 MDU6SXNzdWU0MjgzMDAzNDU= 2864 Bug in WarpedVRT support of open_rasterio() 10595679 closed 0     3 2019-04-02T15:37:08Z 2019-04-11T16:24:13Z 2019-04-11T16:24:13Z CONTRIBUTOR      

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

```python 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(): https://github.com/pydata/xarray/blob/0c73a380745c4792ab440eb020f78f203897abe5/xarray/backends/rasterio_.py#L222

With the following: python 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
{
    "url": "https://api.github.com/repos/pydata/xarray/issues/2864/reactions",
    "total_count": 0,
    "+1": 0,
    "-1": 0,
    "laugh": 0,
    "hooray": 0,
    "confused": 0,
    "heart": 0,
    "rocket": 0,
    "eyes": 0
}
  completed 13221727 issue

Links from other tables

  • 0 rows from issues_id in issues_labels
  • 3 rows from issue in issue_comments
Powered by Datasette · Queries took 0.723ms · About: xarray-datasette