reproject_match() output is identical regardless of requested resampling method.

trchudley commented 1 year ago

Code Sample, a copy-pastable example if possible

import rioxarray as rxr
from rasterio.enums import Resampling

original_tif_fpath = './original.tif'
target_tif_fpath = './target.tif'

original_ds = rxr.open_rasterio(original_tif_fpath)
target_ds = rxr.open_rasterio(target_tif_fpath)

print(f'Original: {original_ds.rio.crs}, resolution {original_ds.rio.resolution()}, shape {original_ds.rio.shape}')
print(f'Target: {target_ds.rio.crs}, resolution {target_ds.rio.resolution()}, shape {target_ds.rio.shape}')

# Original: EPSG:32622, resolution (60.0, -60.0), shape (811, 1100)
# Target: EPSG:3413, resolution (15.0, -15.0), shape (2840, 4140)

nearest = original_ds.rio.reproject_match(reference_ds, resampling=Resampling.nearest)
cubic = original_ds.rio.reproject_match(reference_ds, resampling=Resampling.cubic)

print(f'Resampled (nearest): {nearest.rio.crs}, resolution {nearest.rio.resolution()}, shape {nearest.rio.shape}')
print(f'Resampled (cubic):   {cubic.rio.crs}, resolution {cubic.rio.resolution()}, shape {cubic.rio.shape}')

# Resampled (nearest): EPSG:3413, resolution (15.0, -15.0), shape (2840, 4140)
# Resampled (cubic):   EPSG:3413, resolution (15.0, -15.0), shape (2840, 4140)

# Outputs now match the original extent/resolution/crs but both are identical:

diff = cubic - nearest

print(f'Min: {diff.max().item(0)}, Max: {diff.min().item(0)}')

# Min: 0, Max: 0

Problem description

I am attempting to resample older Landsat 60 m data to match newer Landsat panchromatic (15 m) data using the reproject_match() function. Whilst the resampling is apparently successful, the output is identical regardless of the rasterio.enums.Resampling method I choose (I would like cubic, but I achieve identical results when specifying nearest, average, max, cubic_spline, etc.).

Based upon visual inspection of the output data, I suspect that the resampling method is reverting to 'nearest'.

Expected Output

reproject_match() output to differ according to requested resampling method.

Environment Information

rioxarray 0.12.0 rasterio 1.2.10 gdal 3.4.2 python 3.10.6 xarray 2022.6.0

Installation method

Mamba/Conda

Conda environment information (if you installed with conda):

Environment (conda list):

``` $ conda list | grep -E "rasterio|xarray|gdal" gdal 3.4.2 pypi_0 pypi intake-xarray 0.6.1 pyhd8ed1ab_0 conda-forge libgdal 3.4.2 he830e9d_4 conda-forge rasterio 1.2.10 py310hdadb6e9_5 conda-forge rioxarray 0.12.0 pyhd8ed1ab_0 conda-forge xarray 2022.6.0 pyhd8ed1ab_1 conda-forge xarray-spatial 0.3.5 pyhd8ed1ab_0 conda-forge ```

Details about conda and system ( conda info ):

``` $ mamba info active environment : geomamba active env location : /home/chudley.1/software/miniconda3/envs/geomamba shell level : 2 user config file : /home/chudley.1/.condarc populated config files : /home/chudley.1/.condarc conda version : 4.12.0 conda-build version : not installed python version : 3.9.9.final.0 virtual packages : __linux=3.10.0=0 __glibc=2.17=0 __unix=0=0 __archspec=1=x86_64 base environment : /home/chudley.1/software/miniconda3 (writable) conda av data dir : /home/chudley.1/software/miniconda3/etc/conda conda av metadata url : None channel URLs : https://conda.anaconda.org/conda-forge/linux-64 https://conda.anaconda.org/conda-forge/noarch https://repo.anaconda.com/pkgs/main/linux-64 https://repo.anaconda.com/pkgs/main/noarch https://repo.anaconda.com/pkgs/r/linux-64 https://repo.anaconda.com/pkgs/r/noarch package cache : /home/chudley.1/software/miniconda3/pkgs /home/chudley.1/.conda/pkgs envs directories : /home/chudley.1/software/miniconda3/envs /home/chudley.1/.conda/envs platform : linux-64 user-agent : conda/4.12.0 requests/2.25.1 CPython/3.9.9 Linux/3.10.0-1160.76.1.el7.x86_64 rhel/7.9 glibc/2.17 UID:GID : 80204617:11003741 netrc file : /home/chudley.1/.netrc offline mode : False ```

snowman2 commented 1 year ago

I just inspected the code and don't see anything that appears to be incorrect.

scottyhq commented 1 year ago

For what it's worth I've noticed this behavior on microsoft planetary computer (https://planetarycomputer.microsoft.com/compute, JUPYTER_IMAGE=pcccr.azurecr.io/public/planetary-computer/python:2022.05.11.0)

rasterio                  1.2.10           py38h667dea4_5    conda-forge
gdal                      3.4.2            py38hdd69c9e_6    conda-forge
rioxarray                 0.11.1             pyhd8ed1ab_0    conda-forge

Figured it was something to do with caching but never looked into it. Upgrading to rasterio>1.3 / gdal>3.5 resampling behaves as expected.

trchudley commented 1 year ago

Thanks Scott, good to see it's reproducible. I was ultimately able to move to a new environment with an updated rasterio/gdal that worked as desired. This is a something that's easily missable if you're not paying attention, so I'm not sure whether it might be desirable to have some sort of check or warning on the versioning - will leave it for others to decide.

Cheers.

corteva / rioxarray