I’m trying to spatially subset a GFS 0.5 degree forecast read form a GRIB file downloaded from NOAA website (https://www.ncei.noaa.gov/data/global-forecast-system/access/grid-004-0.5-degree/analysis/)
I can load the file using rasterio with cfgrib backend, but when it comes to clipping the dataset to a selected region, rasterio crashes with “No valid geometry objects found for rasterize” error message. Identical clipping operation works on MODIS files, so so I don’t expect problems here.

My end goal is re-rasterization and co-registering of the forecast results with satellite observations, and clipping the dataset is the first step in the process.

Here is an MWE reproducing the problem:

import cfgrib
import xarray as xr

data_source="gfs_4_20240201_0000_000.grb2"
variable_name = "t"

south_america_bbox = {
    "coordinates": [[
        [-90.246917, -57.266272], 
        [-90.246917, 17.482479], 
        [-33.603376, 17.482479], 
        [-33.603376, -57.266272], 
        [-90.246917, -57.266272]]], 
    "type": "Polygon"}

ds = xr.open_dataset(data_source, engine="cfgrib", filter_by_keys={'typeOfLevel': 'surface'})
dataset = ds[variable_name]
dataset.rio.write_crs("epsg:4326", inplace=True)

ds_clipped = dataset.rio.clip(geometries=south_america_bbox, crs=4326)

And here is a slightly abbreviated error message:

venv/lib64/python3.12/site-packages/rasterio/features.py:328: ShapeSkipWarning: Invalid or empty shape coordinates at index 0 will not be rasterized.
  warnings.warn('Invalid or empty shape {} at index {} will not be rasterized.'.format(geom, index), ShapeSkipWarning)
venv/lib64/python3.12/site-packages/rasterio/features.py:328: ShapeSkipWarning: Invalid or empty shape type at index 1 will not be rasterized.
  warnings.warn('Invalid or empty shape {} at index {} will not be rasterized.'.format(geom, index), ShapeSkipWarning)

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[1], line 23
     20 dataset = ds[variable_name]
     21 dataset.rio.write_crs("epsg:4326", inplace=True)
---> 23 ds_clipped = dataset.rio.clip(geometries=south_america_bbox, crs=4326)

File venv/lib64/python3.12/site-packages/rioxarray/raster_array.py:911, in RasterArray.clip(self, geometries, crs, all_touched, drop, invert, from_disk)
    903     cropped_ds = _clip_from_disk(
    904         self._obj,
    905         geometries=geometries,
   (...)
    908         invert=invert,
    909     )
    910 if cropped_ds is None:
--> 911     cropped_ds = _clip_xarray(
    912         self._obj,


[...]


File venv/lib64/python3.12/site-packages/rasterio/features.py:331, in rasterize(shapes, out_shape, fill, out, transform, all_touched, merge_alg, default_value, dtype)
    328         warnings.warn('Invalid or empty shape {} at index {} will not be rasterized.'.format(geom, index), ShapeSkipWarning)
    330 if not valid_shapes:
--> 331     raise ValueError('No valid geometry objects found for rasterize')
    333 shape_values = np.array(shape_values)
    335 if not validate_dtype(shape_values, valid_dtypes):

ValueError: No valid geometry objects found for rasterize

At this point I’m stuck trying to make rasterio understand geospatial information embedded in the GRIB file, and I will appreciate your help in getting around the problem.