diff --git a/doc/whats-new.rst b/doc/whats-new.rst index 472cd47adf1..74bc3ee5ba9 100644 --- a/doc/whats-new.rst +++ b/doc/whats-new.rst @@ -20,6 +20,8 @@ v0.10.1 (unreleased) Enhancements ~~~~~~~~~~~~ +- Added nodatavals attribute to DataArray when using :py:func:`~xarray.open_rasterio`. (:issue:`1736`). + By `Alan Snow `_. Bug fixes ~~~~~~~~~ diff --git a/xarray/backends/rasterio_.py b/xarray/backends/rasterio_.py index 13f9c0a3de6..3f7dc3992b9 100644 --- a/xarray/backends/rasterio_.py +++ b/xarray/backends/rasterio_.py @@ -157,6 +157,10 @@ def open_rasterio(filename, chunks=None, cache=None, lock=None): # Affine transformation matrix (tuple of floats) # Describes coefficients mapping pixel coordinates to CRS attrs['transform'] = tuple(riods.transform) + if hasattr(riods, 'nodatavals'): + # The nodata values for the raster bands + attrs['nodatavals'] = tuple([np.nan if nodataval is None else nodataval + for nodataval in riods.nodatavals]) data = indexing.LazilyIndexedArray(RasterioArrayWrapper(riods)) diff --git a/xarray/tests/test_backends.py b/xarray/tests/test_backends.py index 017e290558a..a60b5b5e622 100644 --- a/xarray/tests/test_backends.py +++ b/xarray/tests/test_backends.py @@ -1876,6 +1876,68 @@ def test_serialization(self): with xr.open_dataarray(tmp_nc_file) as ncds: assert_identical(rioda, ncds) + @requires_scipy_or_netCDF4 + def test_nodata(self): + import rasterio + from rasterio.transform import from_origin + + # Create a geotiff file in utm proj + with create_tmp_file(suffix='.tif') as tmp_file: + # data + nx, ny, nz = 4, 3, 3 + data = np.arange(nx*ny*nz, + dtype=rasterio.float32).reshape(nz, ny, nx) + transform = from_origin(5000, 80000, 1000, 2000.) + with rasterio.open( + tmp_file, 'w', + driver='GTiff', height=ny, width=nx, count=nz, + crs={'units': 'm', 'no_defs': True, 'ellps': 'WGS84', + 'proj': 'utm', 'zone': 18}, + transform=transform, + nodata=-9765, + dtype=rasterio.float32) as s: + s.write(data) + expected_nodatavals = [-9765, -9765, -9765] + with xr.open_rasterio(tmp_file) as rioda: + np.testing.assert_array_equal(rioda.attrs['nodatavals'], + expected_nodatavals) + with create_tmp_file(suffix='.nc') as tmp_nc_file: + rioda.to_netcdf(tmp_nc_file) + with xr.open_dataarray(tmp_nc_file) as ncds: + np.testing.assert_array_equal(ncds.attrs['nodatavals'], + expected_nodatavals) + + @requires_scipy_or_netCDF4 + def test_nodata_missing(self): + import rasterio + from rasterio.transform import from_origin + + # Create a geotiff file in utm proj + with create_tmp_file(suffix='.tif') as tmp_file: + # data + nx, ny, nz = 4, 3, 3 + data = np.arange(nx*ny*nz, + dtype=rasterio.float32).reshape(nz, ny, nx) + transform = from_origin(5000, 80000, 1000, 2000.) + with rasterio.open( + tmp_file, 'w', + driver='GTiff', height=ny, width=nx, count=nz, + crs={'units': 'm', 'no_defs': True, 'ellps': 'WGS84', + 'proj': 'utm', 'zone': 18}, + transform=transform, + dtype=rasterio.float32) as s: + s.write(data) + + expected_nodatavals = [np.nan, np.nan, np.nan] + with xr.open_rasterio(tmp_file) as rioda: + np.testing.assert_array_equal(rioda.attrs['nodatavals'], + expected_nodatavals) + with create_tmp_file(suffix='.nc') as tmp_nc_file: + rioda.to_netcdf(tmp_nc_file) + with xr.open_dataarray(tmp_nc_file) as ncds: + np.testing.assert_array_equal(ncds.attrs['nodatavals'], + expected_nodatavals) + def test_utm(self): import rasterio from rasterio.transform import from_origin