'xr_reproject' takes user-given resampling into account only in case of upscaling, not in case of downscaling

[remi-braun]

Hello,

When trying to replace some rioxarray code by odc.geo code to make it lazy and dask-compliant, I was surprised to see that xr_reproject seems to take into account the resampling value given by the user only in case of upscaling. In case of downscaling, I'm suspecting odc.geo performs only a simple decimation.

---

Here is my code:

import xarray as xr
import rioxarray
from affine import Affine
from odc.geo import xr as odc_geo_xr  # noqa: F401
from odc.geo.geobox import GeoBox
def test_reproj(path: str, new_width, new_height, resampling):
   # Open with rioxarray
    with (
        xr.set_options(keep_attrs=True),
        rioxarray.set_options(export_grid_mapping=False),
        rioxarray.open_rasterio(
            path,
            chunks="auto",
            masked=True,
        ) as xda,
    ):
        # Reproject with rioxarray
        xda1 = xda.rio.reproject(
            xda.rio.crs,
            shape=(new_height, new_width),
            resampling=resampling,
        )
        xda1.rio.to_raster(fr"rioxarray_{round(xda1.rio.resolution()[0])}m_{resampling}.tif")

        # Reproject with dask + odc-geo
        src_affine = xda.rio.transform()
        xda = xda.odc.reproject(
            how=GeoBox(
                (new_height, new_width),
                src_affine * Affine.scale(xda.rio.width / new_height, xda.rio.height / new_height),
                xda.rio.crs,
            ),
            resampling=resampling,
            dst_nodata=xda.rio.encoded_nodata,
        )
        # Set nodata in rioxr's way and remove odc.geo nodata in attributes
        xda.attrs.pop("nodata", None)
        xda.rio.write_nodata(xda.rio.encoded_nodata, encoded=True, inplace=True)
        xda.rio.to_raster(rf"odc_geo_{round(xda.rio.resolution()[0])}m_{resampling}.tif")

Maybe there is something I haven't correctly understood regarding the way to use xr_reproject?

# Around 1500m of resolution (downscaling) + nearest
test_reproj(  
    "LC09_L2SP_152041_20220828_20220830_02_T1_SR_B4.TIF",
    resampling=0,
    new_width=153, new_height=156,
)


# Around 15m of resolution (upscaling) + nearest
test_reproj(      
    "LC09_L2SP_152041_20220828_20220830_02_T1_SR_B4.TIF",
    resampling=0,
    new_width=15342, new_height=15622,
)


# Around 1500m of resolution (downscaling) + bilinear
test_reproj(  
    "LC09_L2SP_152041_20220828_20220830_02_T1_SR_B4.TIF",
    resampling=1,
    new_width=153, new_height=156,
)


# Around 15m of resolution (upscaling) + bilinear
test_reproj(      
    "LC09_L2SP_152041_20220828_20220830_02_T1_SR_B4.TIF",
    resampling=1,
    new_width=15342, new_height=15622,
)

For the upsampling, the result is quite convincing, the images are almost the same.

	Rioxarray	odc.geo
Nearest + upsampling ✔️
Bilinear + upsampling ✔️

In the case of nearest + downscaling, odc.geo gives sth that looks like nearest (but might be simple decimation), when rioxarray gives a results with correct resampling.

	Rioxarray	odc.geo
Nearest + downsampling ❓

In the case of bilinear+ downscaling, odc.geo gives sth completely wrong, when rioxarray gives again a results with correct resampling.

	Rioxarray	odc.geo
Bilinear + downsampling ❌

Here is the link to a zip containing all the results: https://seafile.unistra.fr/f/772545d7b8bc4b15a7c4/?dl=1

---

PS: there are also other edge effects, I can open dedicated issues if you want:

• odc.geo doesn't compute properly the right corner (see screenshots)
• the size of odc.geo files is almost doubled on disk
• for a given size when downsampling, the resolution is not exactly the same as GDAL's

---

Tested with Python 3.9.19

rioxarray (0.15.0) deps:
  rasterio: 1.3.10
    xarray: 2024.3.0
      GDAL: 3.8.4
      GEOS: 3.11.3
      PROJ: 9.3.1

odc-geo version: 0.4.8

And Python 3.11.13

rioxarray (0.19.0) deps:
  rasterio: 1.4.3
    xarray: 2025.7.0
      GDAL: 3.6.2
      GEOS: 3.11.1
      PROJ: 9.1.1

odc-geo version: 0.4.10

All with the same result.

[robbibt]

A long time ago I posted this issue on datacube-core - I suspect it's the same problem: opendatacube/datacube-core#779

I agree that that this is really non-intuitive and in many cases undesirable functionality. Having the ability to turn off decimation would be really valuable - I've had to use non-ODC tooling for this kind of task in the past to get around this.

[Kirill888]

@remi-braun instead of computing geobox you think rioxarray used for the down-sampled image, the test should use the exact output geobox, i.e. how = xda1.odc.geobox, for the comparison to make sense. Issue with the corner is probably due to that also.

@robbibt issue you have linked is about reproject performed as part of IO read operation and has to do with GDAL using overview images stored in the file, this is extremely helpful in most cases where you need a "thumbnail" of a large image. Not sure about datacube but in odc-stac there is an option to disable use of overviews with use_overviews=False, but I have just checked the docs and that option is not currently documented, and probably not properly tested, but it's there.

[robbibt]

in odc-stac there is an option to disable use of overviews with use_overviews=False, but I have just checked the docs and that option is not currently documented, and probably not properly tested, but it's there.

That's really cool! Yeah, that is exactly what I was looking for - will test it out and see if I can add some documention.

[Kirill888]

@remi-braun FYI: xr_reproject also understands non-geobox way of specifying destination pixel grid, i.e.

xda1 = xda.rio.reproject(
            xda.rio.crs,
            shape=(new_height, new_width),
            resampling=resampling,
        )
# same in ODC
xda2 = xda.odc.reproject(
            xda.rio.crs,
            shape=(new_height, new_width),
            resampling=resampling,
            anchor="floating", # << closest option to what rioxarray does
        )

although the grids produced differ somewhat (i.e. xda1.odc.geobox != xda2.odc.geobox)

[Kirill888]

Alright, it's due to GDAL's XSCALE/YSCALE work-around #100, #224 and affects bilinear,cubic,cubic_spline,lanczos resampling methods.

supplying extra parameters to .odc.reproject(..., XSCALE=None, YSCALE=None) makes bilinear output look very similar to what rioxarray produces. But, this will cause blurry output when using Dask some of the times, see rasterio/rasterio#2839 for why and how, discovered in datacube-core opendatacube/datacube-core#1448. Not sure if this was fixed in GDAL already, and if yes in what version of GDAL (maybe not? - OSGeo/gdal#7750, OSGeo/gdal#12595 )

In my opinion XSCALE=1, YSCALE=1, which is the current default, is a more reasonable default, as it always samples directly from the source image array. Default GDAL behaviour is to effectively first perform shrink via average sampling followed by actual resampling requested by the user afterwards, hence the image looks more "smooth" in the rio case, and "patchy" in the case of odc.

@emmaai, @robbibt we might need to rethink {X|Y}SCALE=1 default once more.

[remi-braun]

@remi-braun instead of computing geobox you think rioxarray used for the down-sampled image, the test should use the exact output geobox, i.e. how = xda1.odc.geobox, for the comparison to make sense. Issue with the corner is probably due to that also.

The thing is, in the new version of the code I won't be able to access the xda1 array as I want to go full odc.geo. This is why I tried to do a code sample without any help of the rioxarray reprojection

---

# same in ODC
xda2 = xda.odc.reproject(
            xda.rio.crs,
            shape=(new_height, new_width),
            resampling=resampling,
            anchor="floating", # << closest option to what rioxarray does
        )

Thanks for giving the simplified version, appreciate it ! 🙏
✔️ The corner is computed correctly here, but the size of the file is the same than previous version.

although the grids produced differ somewhat (i.e. xda1.odc.geobox != xda2.odc.geobox)

This different grid bothers me a bit, could you explain me why it happens (if you know)?

rioxarray	odc.geo + geobox	odc.geo without geobox

---

Moreover, the issue appears also in this new version, but I think you see that also since your last messages 😄

[Kirill888]

@remi-braun see my previous comment about real reason for the difference, see below for closest I can get to what rioxarray produces:

# same in ODC
xda2 = xda.odc.reproject(
            xda.rio.crs,
            shape=(new_height, new_width),
            resampling=resampling,
            anchor="floating", # << closest option to what rioxarray does
            XSCALE=None, YSCALE=None, # << turn off GDAL warp work-arounds (like rioxarray)
        )

difference in grid computation from rioxarray will need to be investigated, but I'm fairly confident that odc-geo version is computed "correctly". By correctly I mean output image has exactly the requested shape or resolution, and covers roughly the same area as the original. Footprint can't be exactly the same area as input due to discrete nature of pixels. In odc-geo output image footprint might overhang by some amount, smaller than the sampling distance of the output image.

The thing is, in the new version of the code I won't be able to access the xda1 array as I want to go full odc.geo. This is why I tried to do a code sample without any help of the rioxarray reprojection

I understand that, but this is "investigate difference in pixel sampling between two libs"...

[Kirill888]

As far as size of the file goes, it should become very similar once you add XSCALE= parameters. Images that are smooth compress better than images with speckle noise, nothing surprising there. Also I looked at the images in your zip file, and odc ones have no NODATA set, so I suspect the logic you have there for fiddling with nodata is not correct.