rioxarray.clip returns some nan values even if data value exist for these coordinates

[TheRed86]

Hi everyone,

I am using rioxarray.clip function to select a subset of a NetCDF dataset.

I get nan values in some coordinates after clipping, but I have verified that data exist in the original dataset (before clipping).
This problem happens when the geometry is not a square or a rectangle.

Here is an example of geometry that returns a few nan values:

"POLYGON((-31.486297087721766 60.87669368042079,-36.71370291227825 60.87669368042079,-36.52308966145061 58.335063073172215,-31.676910338549394 58.335063073172215,-31.486297087721763 60.87669368042079))".

Here is my simple code:

import rioxarray
import xarray as xr

data = xr.open_dataset('2023-10-18_00_t.nc')
ds = data.rio.write_crs("EPSG:4326")
geo = 'POLYGON((-31.486297087721766 60.87669368042079,-36.71370291227825 60.87669368042079,-36.52308966145061 58.335063073172215,-31.676910338549394 58.335063073172215,-31.486297087721763 60.87669368042079))'

coords = geo.split("((")[-1].split("))")[0].split(",")
coords_array = []
for cc in coords:
    point = cc.split(" ")
    point_to_float = [float(i) for i in point]
    coords_array.append(point_to_float)

geometries = [
    {
        'type': 'Polygon',
        'coordinates': [coords_array]
    }
]

data_clipped = ds.rio.clip(geometries, ds.rio.crs)

Is it possible to clip data over an irregular geometry? If yes, which can be the reason why I get nan values for a few coordinates?

Thanks!

[snowman2]

The output raster is always going to be a grid. The nan values are for regions of the grid that don't overlap with the shape.

It sounds like you may want to use the envelope of the shape to clip to get the behavior you want.

[TheRed86]

Many thanks @snowman2!

Your suggestion works perfectly.

I now have the following code, which returns no nan values:

import rioxarray
import xarray as xr

data = xr.open_dataset('2023-10-18_00_t.nc')
ds = data.rio.write_crs("EPSG:4326")
geo = 'POLYGON((-31.486297087721766 60.87669368042079,-36.71370291227825 60.87669368042079,-36.52308966145061 58.335063073172215,-31.676910338549394 58.335063073172215,-31.486297087721763 60.87669368042079))'

coords = geo.split("((")[-1].split("))")[0].split(",")
coords_array = []
for cc in coords:
    point = cc.split(" ")
    point_to_float = [float(i) for i in point]
    coords_array.append(point_to_float)

geometries = [
    {
        'type': 'Polygon',
        'coordinates': [coords_array]
    }
]
# Extract the envelope coordinates
        envelope_coords = [
            [min(p[0] for p in coords_array), min(p[1] for p in coords_array)],
            [max(p[0] for p in coords_array), min(p[1] for p in coords_array)],
            [max(p[0] for p in coords_array), max(p[1] for p in coords_array)],
            [min(p[0] for p in coords_array), max(p[1] for p in coords_array)],
        ]

        envelope = {
            'type': 'Polygon',
            'coordinates': [envelope_coords]
        }

        # Clip using the envelope
        data_clipped = ds.rio.clip([envelope], ds.rio.crs, all_touched=True, from_disk=True)