From 9de23efb4b8fb4b8305c6e3dc09d2957433e5939 Mon Sep 17 00:00:00 2001 From: snowman2 Date: Thu, 1 Jul 2021 12:13:45 -0500 Subject: [PATCH 1/3] ENH: Allow passing in kwargs to reproject & overide nodata. Provide default nodata based on dtype. --- docs/history.rst | 4 +- rioxarray/raster_array.py | 49 ++++++++++++++++--- .../integration/test_integration_rioxarray.py | 7 +-- 3 files changed, 49 insertions(+), 11 deletions(-) diff --git a/docs/history.rst b/docs/history.rst index 75761f5e..c230441e 100644 --- a/docs/history.rst +++ b/docs/history.rst @@ -3,7 +3,9 @@ History Latest ------ -- BUG: Remove duplicate acquire in open_rasterio +- ENH: Allow passing in kwargs to reproject & overide nodata. Provide +default nodata based on dtype. (issue ##369) +- BUG: Remove duplicate acquire in open_rasterio (pull #364) 0.4.3 ------ diff --git a/rioxarray/raster_array.py b/rioxarray/raster_array.py index 75f7024f..881d506d 100644 --- a/rioxarray/raster_array.py +++ b/rioxarray/raster_array.py @@ -18,6 +18,7 @@ import rasterio.mask import rasterio.warp import xarray +from rasterio.dtypes import dtype_rev from rasterio.enums import Resampling from rasterio.features import geometry_mask from scipy.interpolate import griddata @@ -38,6 +39,24 @@ ) from rioxarray.rioxarray import XRasterBase, _get_data_var_message, _make_coords +# DTYPE TO NODATA MAP +# Based on: https://github.com/OSGeo/gdal/blob/ +# cde27dc7641964a872efdc6bbcf5e3d3f7ab9cfd/gdal/ +# swig/python/gdal-utils/osgeo_utils/gdal_calc.py#L62 +_NODATA_DTYPE_MAP = { + 1: 255, # GDT_Byte + 2: 65535, # GDT_UInt16 + 3: -32768, # GDT_Int16 + 4: 4294967293, # GDT_UInt32 + 5: -2147483647, # GDT_Int32 + 6: 3.402823466e38, # GDT_Float32 + 7: 1.7976931348623158e308, # GDT_Float64 + 8: -32768, # GDT_CInt16 + 9: -2147483647, # GDT_CInt32 + 10: 3.402823466e38, # GDT_CFloat32 + 11: 1.7976931348623158e308, # GDT_CFloat64 +} + def _generate_attrs(src_data_array, dst_nodata): # add original attributes @@ -318,7 +337,9 @@ def reproject( resolution=None, shape=None, transform=None, + nodata=None, resampling=Resampling.nearest, + **kwargs, ): """ Reproject :obj:`xarray.DataArray` objects @@ -332,6 +353,7 @@ def reproject( .. versionadded:: 0.0.27 shape .. versionadded:: 0.0.28 transform + .. versionadded:: 0.5.0 nodata, kwargs Parameters ---------- @@ -343,10 +365,21 @@ def reproject( shape: tuple(int, int), optional Shape of the destination in pixels (dst_height, dst_width). Cannot be used together with resolution. - transform: optional + transform: Affine, optional The destination transform. resampling: rasterio.enums.Resampling, optional See :func:`rasterio.warp.reproject` for more details. + nodata: float, optional + The nodata value used to initialize the destination; + it will remain in all areas not covered by the reprojected source. + Defaults to the nodata value of the source image if none provided + and exists or attempts to find an appropriate value by dtype. + **kwargs: dict + Additional keyword arguments to pass into :func:`rasterio.warp.reproject`. + To override: + - src_transform: `rio.write_transform` + - src_crs: `rio.write_crs` + - src_nodata: `rio.write_nodata` Returns @@ -382,22 +415,24 @@ def reproject( else: dst_data = np.zeros((dst_height, dst_width), dtype=self._obj.dtype.type) - dst_nodata = self._obj.dtype.type( - self.nodata if self.nodata is not None else -9999 - ) - src_nodata = self._obj.dtype.type( - self.nodata if self.nodata is not None else dst_nodata + default_nodata = ( + _NODATA_DTYPE_MAP[dtype_rev[self._obj.dtype.name]] + if self.nodata is None + else self.nodata ) + dst_nodata = default_nodata if nodata is None else nodata + rasterio.warp.reproject( source=self._obj.values, destination=dst_data, src_transform=src_affine, src_crs=self.crs, - src_nodata=src_nodata, + src_nodata=self.nodata, dst_transform=dst_affine, dst_crs=dst_crs, dst_nodata=dst_nodata, resampling=resampling, + **kwargs, ) # add necessary attributes new_attrs = _generate_attrs(self._obj, dst_nodata) diff --git a/test/integration/test_integration_rioxarray.py b/test/integration/test_integration_rioxarray.py index bced49b7..16d74031 100644 --- a/test/integration/test_integration_rioxarray.py +++ b/test/integration/test_integration_rioxarray.py @@ -721,10 +721,10 @@ def test_reproject__no_nodata(modis_reproject): # replace -9999 with original _FillValue for testing if hasattr(mds_repr, "variables"): for var in mds_repr.rio.vars: - mds_repr[var].values[mds_repr[var].values == -9999] = orig_fill + mds_repr[var].values[mds_repr[var].values == -32768] = orig_fill else: - mds_repr.values[mds_repr.values == -9999] = orig_fill - _mod_attr(mdc, "_FillValue", val=-9999) + mds_repr.values[mds_repr.values == -32768] = orig_fill + _mod_attr(mdc, "_FillValue", val=-32768) # test _assert_xarrays_equal(mds_repr, mdc) @@ -1069,6 +1069,7 @@ def test_geographic_reproject__missing_nodata(): mds_repr = mda.rio.reproject("epsg:32721") # mds_repr.to_netcdf(sentinel_2_utm) # test + _mod_attr(mdc, "_FillValue", val=65535) _assert_xarrays_equal(mds_repr, mdc, precision=4) From 5c30e91b5cfc6c1a7fa3f13ee9b8a19ea6c30d26 Mon Sep 17 00:00:00 2001 From: snowman2 Date: Tue, 6 Jul 2021 10:39:20 -0500 Subject: [PATCH 2/3] ENH: Add support for passing in gcps to rio.reproject --- docs/history.rst | 5 +- rioxarray/raster_array.py | 29 +++++----- rioxarray/raster_dataset.py | 17 +++++- .../integration/test_integration_rioxarray.py | 54 +++++++++++++++++-- 4 files changed, 82 insertions(+), 23 deletions(-) diff --git a/docs/history.rst b/docs/history.rst index c230441e..1b6eced8 100644 --- a/docs/history.rst +++ b/docs/history.rst @@ -3,8 +3,9 @@ History Latest ------ -- ENH: Allow passing in kwargs to reproject & overide nodata. Provide -default nodata based on dtype. (issue ##369) +- ENH: Allow passing in kwargs to `rio.reproject` (issue #369; pull #370) +- ENH: Allow nodata override and provide default nodata based on dtype in `rio.reproject` (pull #370) +- ENH: Add support for passing in gcps to rio.reproject (issue #339; pull #370) - BUG: Remove duplicate acquire in open_rasterio (pull #364) 0.4.3 diff --git a/rioxarray/raster_array.py b/rioxarray/raster_array.py index 881d506d..de0460e5 100644 --- a/rioxarray/raster_array.py +++ b/rioxarray/raster_array.py @@ -106,10 +106,10 @@ def _add_attrs_proj(new_data_array, src_data_array): def _make_dst_affine( - src_data_array, src_crs, dst_crs, dst_resolution=None, dst_shape=None + src_data_array, src_crs, dst_crs, dst_resolution=None, dst_shape=None, **kwargs ): """Determine the affine of the new projected `xarray.DataArray`""" - src_bounds = src_data_array.rio.bounds() + src_bounds = () if "gcps" in kwargs else src_data_array.rio.bounds() src_height, src_width = src_data_array.rio.shape dst_height, dst_width = dst_shape if dst_shape is not None else (None, None) # pylint: disable=isinstance-second-argument-not-valid-type @@ -117,22 +117,21 @@ def _make_dst_affine( dst_resolution = tuple(abs(res_val) for res_val in dst_resolution) elif dst_resolution is not None: dst_resolution = abs(dst_resolution) - resolution_or_width_height = { - k: v - for k, v in [ - ("resolution", dst_resolution), - ("dst_height", dst_height), - ("dst_width", dst_width), - ] - if v is not None - } + + for key, value in ( + ("resolution", dst_resolution), + ("dst_height", dst_height), + ("dst_width", dst_width), + ): + if value is not None: + kwargs[key] = value dst_affine, dst_width, dst_height = rasterio.warp.calculate_default_transform( src_crs, dst_crs, src_width, src_height, *src_bounds, - **resolution_or_width_height, + **kwargs, ) return dst_affine, dst_width, dst_height @@ -337,8 +336,8 @@ def reproject( resolution=None, shape=None, transform=None, - nodata=None, resampling=Resampling.nearest, + nodata=None, **kwargs, ): """ @@ -394,10 +393,10 @@ def reproject( "CRS not found. Please set the CRS with 'rio.write_crs()'." f"{_get_data_var_message(self._obj)}" ) - src_affine = self.transform(recalc=True) + src_affine = None if "gcps" in kwargs else self.transform(recalc=True) if transform is None: dst_affine, dst_width, dst_height = _make_dst_affine( - self._obj, self.crs, dst_crs, resolution, shape + self._obj, self.crs, dst_crs, resolution, shape, **kwargs ) else: dst_affine = transform diff --git a/rioxarray/raster_dataset.py b/rioxarray/raster_dataset.py index 6cf80a0b..e6f772bf 100644 --- a/rioxarray/raster_dataset.py +++ b/rioxarray/raster_dataset.py @@ -58,6 +58,8 @@ def reproject( shape=None, transform=None, resampling=Resampling.nearest, + nodata=None, + **kwargs, ): """ Reproject :class:`xarray.Dataset` objects @@ -69,6 +71,7 @@ def reproject( .. versionadded:: 0.0.27 shape .. versionadded:: 0.0.28 transform + .. versionadded:: 0.5.0 nodata, kwargs Parameters ---------- @@ -84,7 +87,17 @@ def reproject( The destination transform. resampling: rasterio.enums.Resampling, optional See :func:`rasterio.warp.reproject` for more details. - + nodata: float, optional + The nodata value used to initialize the destination; + it will remain in all areas not covered by the reprojected source. + Defaults to the nodata value of the source image if none provided + and exists or attempts to find an appropriate value by dtype. + **kwargs: dict + Additional keyword arguments to pass into :func:`rasterio.warp.reproject`. + To override: + - src_transform: `rio.write_transform` + - src_crs: `rio.write_crs` + - src_nodata: `rio.write_nodata` Returns -------- @@ -102,6 +115,8 @@ def reproject( shape=shape, transform=transform, resampling=resampling, + nodata=nodata, + **kwargs, ) ) return resampled_dataset diff --git a/test/integration/test_integration_rioxarray.py b/test/integration/test_integration_rioxarray.py index 16d74031..baf04e87 100644 --- a/test/integration/test_integration_rioxarray.py +++ b/test/integration/test_integration_rioxarray.py @@ -15,6 +15,7 @@ from dask.delayed import Delayed from numpy.testing import assert_almost_equal, assert_array_equal from pyproj import CRS as pCRS +from rasterio.control import GroundControlPoint from rasterio.crs import CRS from rasterio.windows import Window @@ -699,7 +700,8 @@ def test_reproject__no_transform(modis_reproject): _assert_xarrays_equal(mds_repr, mdc) -def test_reproject__no_nodata(modis_reproject): +@pytest.mark.parametrize("nodata", [None, -9999]) +def test_reproject__no_nodata(nodata, modis_reproject): mask_args = ( dict(masked=False, mask_and_scale=False) if "rasterio" in str(modis_reproject["open"]) @@ -712,19 +714,20 @@ def test_reproject__no_nodata(modis_reproject): _del_attr(mda, "_FillValue") _del_attr(mda, "nodata") # reproject - mds_repr = mda.rio.reproject(modis_reproject["to_proj"]) + mds_repr = mda.rio.reproject(modis_reproject["to_proj"], nodata=nodata) # overwrite test dataset # if isinstance(modis_reproject['open'], xarray.DataArray): # mds_repr.to_netcdf(modis_reproject['compare']) # replace -9999 with original _FillValue for testing + fill_nodata = -32768 if nodata is None else nodata if hasattr(mds_repr, "variables"): for var in mds_repr.rio.vars: - mds_repr[var].values[mds_repr[var].values == -32768] = orig_fill + mds_repr[var].values[mds_repr[var].values == fill_nodata] = orig_fill else: - mds_repr.values[mds_repr.values == -32768] = orig_fill - _mod_attr(mdc, "_FillValue", val=-32768) + mds_repr.values[mds_repr.values == fill_nodata] = orig_fill + _mod_attr(mdc, "_FillValue", val=fill_nodata) # test _assert_xarrays_equal(mds_repr, mdc) @@ -750,6 +753,47 @@ def test_reproject__no_nodata_masked(modis_reproject): _assert_xarrays_equal(mds_repr, mdc) +def test_reproject__gcps_kwargs(tmp_path): + tiffname = tmp_path / "test.tif" + src_gcps = [ + GroundControlPoint(row=0, col=0, x=156113, y=2818720, z=0), + GroundControlPoint(row=0, col=800, x=338353, y=2785790, z=0), + GroundControlPoint(row=800, col=800, x=297939, y=2618518, z=0), + GroundControlPoint(row=800, col=0, x=115698, y=2651448, z=0), + ] + crs = CRS.from_epsg(32618) + with rasterio.open( + tiffname, + mode="w", + height=800, + width=800, + count=3, + dtype=numpy.uint8, + driver="GTiff", + ) as source: + source.gcps = (src_gcps, crs) + + rds = rioxarray.open_rasterio(tiffname) + rds.rio.write_crs(crs, inplace=True) + rds = rds.rio.reproject( + crs, + gcps=src_gcps, + ) + assert rds.rio.height == 923 + assert rds.rio.width == 1027 + assert rds.rio.crs == crs + assert rds.rio.transform().almost_equals( + Affine( + 216.8587081056465, + 0.0, + 115698.25, + 0.0, + -216.8587081056465, + 2818720.0, + ) + ) + + def test_reproject_match(modis_reproject_match): mask_args = ( dict(masked=False, mask_and_scale=False) From 17472b6f13ffb525829c2bfd8bfa9e3e30a24444 Mon Sep 17 00:00:00 2001 From: snowman2 Date: Tue, 6 Jul 2021 11:14:04 -0500 Subject: [PATCH 3/3] CI: Pin libgdal<3.3 --- .github/workflows/tests.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/tests.yaml b/.github/workflows/tests.yaml index f5f2196c..73455bf8 100644 --- a/.github/workflows/tests.yaml +++ b/.github/workflows/tests.yaml @@ -60,7 +60,7 @@ jobs: - name: Install Env shell: bash run: | - conda create -n test python=${{ matrix.python-version }} rasterio=${{ matrix.rasterio-version }} xarray=${{ matrix.xarray-version }} scipy pyproj netcdf4 dask pandoc + conda create -n test python=${{ matrix.python-version }} rasterio=${{ matrix.rasterio-version }} xarray=${{ matrix.xarray-version }} 'libgdal<3.3' scipy pyproj netcdf4 dask pandoc source activate test python -m pip install -e .[all]