Reproject a 3-D DataArray with Nodata

[QRW555]

Following this thread, when using reproject for a 2-D DataArray, after explicitly setting NoData value using write_nodata, Resampling.average can correctly calculate the average of all non-Nodata pixels.

However, this seems not to work for a 3-D DataArray. reproject will propagate nan value, instead of calculating the average of all non-Nodata pixels, even after setting the NoData value. Here is an example:

import xarray as xr
import rioxarray as rxr
import numpy as np
from rasterio.enums import Resampling

data = [[[1, 2, 3, 4], [5, np.nan, 7, 8], [9, 10, 11, 12], [13, 14, 15, 16]],
        [[17, 18, 19, 20], [21, 22, 23, 24], [25, np.nan, 27, 28], [29, 30, 31, np.nan]]]
x = np.arange(100, 116, 4)
y = np.arange(36, 20, -4)
time = [1, 2]

foo = xr.DataArray(data, coords=[time, y, x], dims=['time', 'y', 'x'])
foo = foo.rio.write_crs('EPSG:4326')
foo = foo.rio.write_nodata(np.nan)

foo_resampled = foo.rio.reproject(
    foo.rio.crs,
    resolution=8,
    resampling=Resampling.average
)

print(f'foo: \n {foo}')
print(f'foo_resampled: \n {foo_resampled}')

foo:

 <xarray.DataArray (time: 2, y: 4, x: 4)> Size: 256B
array([[[ 1.,  2.,  3.,  4.],
        [ 5., nan,  7.,  8.],
        [ 9., 10., 11., 12.],
        [13., 14., 15., 16.]],

       [[17., 18., 19., 20.],
        [21., 22., 23., 24.],
        [25., nan, 27., 28.],
        [29., 30., 31., nan]]])
Coordinates:
  * time         (time) int64 16B 1 2
  * y            (y) int64 32B 36 32 28 24
  * x            (x) int64 32B 100 104 108 112
    spatial_ref  int64 8B 0
Attributes:
    _FillValue:  nan

foo_resampled:

 <xarray.DataArray (time: 2, y: 2, x: 2)> Size: 64B
array([[[ nan,  5.5],
        [11.5, 13.5]],

       [[19.5, 21.5],
        [ nan,  nan]]])
Coordinates:
  * x            (x) float64 16B 102.0 110.0
  * y            (y) float64 16B 34.0 26.0
  * time         (time) int64 16B 1 2
    spatial_ref  int64 8B 0
Attributes:
    _FillValue:  nan

What I expected for foo_resampled:

<xarray.DataArray (time: 2, y: 2, x: 2)> Size: 64B
array([[[ 2.66666667,  5.5       ],
        [11.5       , 13.5       ]],

       [[19.5       , 21.5       ],
        [28.        , 28.66666667]]])
Coordinates:
  * x            (x) float64 16B 102.0 110.0
  * y            (y) float64 16B 34.0 26.0
  * time         (time) int64 16B 1 2
    spatial_ref  int64 8B 0
Attributes:
    _FillValue:  nan

Does that imply reproject can only correctly handle Nan value in 2-D case? I appreciate any suggestions or helps.

[VolodymyrLinuxovich]

This comes from the warp option UNIFIED_SRC_NODATA. rasterio's reproject sets it to YES by default (same as gdalwarp), which makes GDAL treat a source pixel as nodata only when it is nodata in all bands. In your example each NaN is only in one time slice, so GDAL treats those pixels as valid and the NaN gets averaged in. A 2-D array has only one band, so it doesn't happen there.

rio.reproject forwards extra kwargs as warp options, so you can make nodata per-band:

foo_resampled = foo.rio.reproject(
    foo.rio.crs,
    resolution=8,
    resampling=Resampling.average,
    UNIFIED_SRC_NODATA="NO",   # "PARTIAL" also works here
)

which gives

[[[ 2.667,  5.5  ], [11.5, 13.5  ]],
 [[19.5  , 21.5  ], [28. , 28.667]]]

i.e. the same result as reprojecting each time slice separately (checked with rioxarray main, rasterio 1.5.1, GDAL 3.12).