1

I was trying to use rio.reproject_match() to match the resolution of two xarray datasets. However, when I used resample methods like Resampling.bilinear or Resampling.average, a large area would become Nan value after reprojection and a lot of information was lost.

I guess the reason is that rio.reproject() will not skip the Nan value. A simple 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]]
x = np.arange(100, 116, 4)
y = np.arange(36, 20, -4)

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

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

foo:

<xarray.DataArray (y: 4, x: 4)>
array([[ 1.,  2.,  3.,  4.],
       [ 5., nan,  7.,  8.],
       [ 9., 10., 11., 12.],
       [13., 14., 15., 16.]])
Coordinates:
  * y            (y) int32 36 32 28 24
  * x            (x) int32 100 104 108 112
    spatial_ref  int32 0

The result foo_resampled is:

<xarray.DataArray (y: 2, x: 2)>
array([[ nan,  5.5],
       [11.5, 13.5]])
Coordinates:
  * x            (x) float64 102.0 110.0
  * y            (y) float64 34.0 26.0
    spatial_ref  int32 0
Attributes:
    _FillValue:  1.7976931348623157e+308

However, what I want is:

<xarray.DataArray (y: 2, x: 2)>
array([[ 2.66666667,  5.5       ],
       [11.5       , 13.5       ]])
Coordinates:
  * y        (y) float64 34.0 26.0
  * x        (x) float64 102.0 110.0

The last one will retain more information. From rasterio document, it says:

Average resampling, computes the weighted average of all non-NODATA contributing pixels.

I think it should have skipped the Nan value when doing calculations? How can I get the result I expected?

1 Answer 1

4

You haven't actually set a NoData value. Explicitly set nan as the NoData value.

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]]
x = np.arange(100, 116, 4)
y = np.arange(36, 20, -4)

foo = xr.DataArray(data, coords=[y, x], dims=['y', 'x'])
foo.rio.write_crs('EPSG:4326', inplace=True)

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

# Set NaN as NoData
foo.rio.write_nodata(np.nan, inplace=True)
foo_resampled_nodata = foo.rio.reproject(
    foo.rio.crs,
    resolution=8,
    resampling=Resampling.average,
    nodata=np.nan # Set NaN as NoData
)

print(foo_resampled)
print(foo_resampled_nodata)

Output:

<xarray.DataArray (y: 2, x: 2)>
array([[ nan,  5.5],
       [11.5, 13.5]])
Coordinates:
  * x            (x) float64 102.0 110.0
  * y            (y) float64 34.0 26.0
    spatial_ref  int64 0
Attributes:
    _FillValue:  1.7976931348623157e+308

<xarray.DataArray (y: 2, x: 2)>
array([[ 2.66666667,  5.5       ],
       [11.5       , 13.5       ]])
Coordinates:
  * x            (x) float64 102.0 110.0
  * y            (y) float64 34.0 26.0
    spatial_ref  int64 0
Attributes:
    _FillValue:  nan
1
  • This is actually what I want. I had no idea that I needed to set the NoData value. Hope others could check the Nodata Management document first. Thanks a lot!
    – QRW
    Jul 17, 2022 at 3:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.