Clipping rotated pole grib2 file with polygon

[haileye]

Hi All, My end goal is to read a grib2 file using rioxarray and clip it with a polygon. I'm running into issues with the grib2 projection definition. I'm fairly new to grib2 file formats, any help is appreciated.

import rioxarray
import xarray as xr
import cfgrib

grid_file_path = 'CMC_hrdps_west_PRATE_SFC_0_rotated_latlon0.009x0.009_20221018T00Z_P005-00.grib2'
wrf = xr.open_dataset(grid_file_path, engine = 'cfgrib')

wrf
<xarray.Dataset>
Dimensions:      (y: 1180, x: 1330)
Coordinates:
    time         datetime64[ns] ...
    step         timedelta64[ns] ...
    surface      float64 ...
    latitude     (y, x) float64 ...
    longitude    (y, x) float64 ...
    valid_time   datetime64[ns] ...
    spatial_ref  int64 0
Dimensions without coordinates: y, x
Data variables:
    prate        (y, x) float32 ...
Attributes:
    GRIB_edition:            2
    GRIB_centre:             cwao
    GRIB_centreDescription:  Canadian Meteorological Service - Montreal
    GRIB_subCentre:          0
    Conventions:             CF-1.7
    institution:             Canadian Meteorological Service - Montreal
    history:                 2022-10-24T15:55 GRIB to CDM+CF via cfgrib-0.9.1...


## clip by geometry
## coordinates in the grib2 file's rotated pole projection
## I've overlain the grib2 file and the geometry below in QGIS to confirm that they overlap
## the geometry for this post is a simple rectangle. My real use case is a more complicated polygon, not a simple rectangle.

clip_geom = [
			{"type":"Polygon",
			"coordinates":[[
				[-20.053465032,4.855709248],
				[-19.468843208,4.521462474],
				[-20.068351733,3.435063352],
				[-20.705115362,3.786607157],
				[-20.053465032,4.855709248]
				]]
			}
]

clipped = wrf.rio.clip(clip_geom)

I run into this error:

  File "/opt/venvs/flux/lib/python3.8/site-packages/rioxarray/raster_dataset.py", line 380, in clip
    self._obj[var]
  File "/opt/venvs/flux/lib/python3.8/site-packages/rioxarray/raster_array.py", line 906, in clip
    raise MissingCRS(
rioxarray.exceptions.MissingCRS: CRS not found. Please set the CRS with 'rio.write_crs()'. Data variable: prate

Next, I define the grib2 projection using the WKT projection string found when I open the file in QGIS:

hrpds_proj_wkt = """GEOGCRS["Coordinate System imported from GRIB file",
	BASEGEOGCRS["Coordinate System imported from GRIB file",
		DATUM["unnamed",
			ELLIPSOID["Sphere",6371229,0,
				LENGTHUNIT["metre",1,
					ID["EPSG",9001]]]],
		PRIMEM["Greenwich",0,
			ANGLEUNIT["degree",0.0174532925199433,
				ID["EPSG",9122]]]],
	DERIVINGCONVERSION["Pole rotation (GRIB convention)",
		METHOD["Pole rotation (GRIB convention)"],
		PARAMETER["Latitude of the southern pole (GRIB convention)",-33.443381,
			ANGLEUNIT["degree",0.0174532925199433,
				ID["EPSG",9122]]],
		PARAMETER["Longitude of the southern pole (GRIB convention)",-93.536426,
			ANGLEUNIT["degree",0.0174532925199433,
				ID["EPSG",9122]]],
		PARAMETER["Axis rotation (GRIB convention)",0,
			ANGLEUNIT["degree",0.0174532925199433,
				ID["EPSG",9122]]]],
	CS[ellipsoidal,2],
		AXIS["latitude",north,
			ORDER[1],
			ANGLEUNIT["degree",0.0174532925199433,
				ID["EPSG",9122]]],
		AXIS["longitude",east,
			ORDER[2],
			ANGLEUNIT["degree",0.0174532925199433,
				ID["EPSG",9122]]]]"""

wrf.rio.write_crs(hrpds_proj_wkt,  inplace=True)

Then I run into a KeyError using the projection:

  File "/opt/venvs/flux/lib/python3.8/site-packages/rioxarray/rioxarray.py", line 493, in write_crs
    grid_map_attrs = pyproj.CRS.from_user_input(data_obj.rio.crs).to_cf()
  File "/opt/venvs/flux/lib/python3.8/site-packages/pyproj/crs/crs.py", line 700, in to_cf
    _INVERSE_GEOGRAPHIC_GRID_MAPPING_NAME_MAP[
KeyError: 'pole rotation (grib convention)'

As the error is a KeyError, I try to define the projection using a proj4string instead:

hrdps_proj4string = "+proj=ob_tran +o_proj=longlat +o_lon_p=-0 +o_lat_p=33.443381 +lon_0=-93.536426 +R=6371229 +no_defs +type=crs"

wrf.rio.write_crs(hrdps_proj4string, inplace=True)
clipped = wrf.rio.clip(clip_geom)

And then I run into an error stating that no data found in bounds. I have confirmed that data is within the bounds using QGIS.

  File "<stdin>", line 1, in <module>
  File "/opt/venvs/flux/lib/python3.8/site-packages/rioxarray/raster_dataset.py", line 380, in clip
    self._obj[var]
  File "/opt/venvs/flux/lib/python3.8/site-packages/rioxarray/raster_array.py", line 941, in clip
    raise NoDataInBounds(
rioxarray.exceptions.NoDataInBounds: No data found in bounds. Data variable: prate

I'm fairly stumped. Any help is greatly appreciated.

Please find the grib2 file attached.

grib2_file.zip

[snowman2]

Your latitude and longitude coordinates are 2D. rioxarray only works with 1D coordinates.

Here are issues for reference for how to deal with this:

• DOC: 2D Coordinate Clipping Examples #202
• DOC: Example of converting 2D coordinates to regular grid #209

[haileye]

Thanks for your response. Through your links and a stack exchange post, I have a workflow that is working. Here is what worked.

import rasterio
import rioxarray
from rioxarray.rioxarray import affine_to_coords
from pyproj.crs import CRS
import xarray as xr
import cfgrib

grid_file_path = 'CMC_hrdps_west_PRATE_SFC_0_rotated_latlon0.009x0.009_20221026T00Z_P003-00.grib2'
hrdps_proj4string = "+proj=ob_tran +o_proj=longlat +o_lon_p=-0 +o_lat_p=33.443381 +lon_0=-93.536426 +R=6371229 +no_defs +type=crs"
clip_geom = [
			{"type":"Polygon",
			"coordinates":[[
				[-20.053465032,4.855709248],
				[-19.468843208,4.521462474],
				[-20.068351733,3.435063352],
				[-20.705115362,3.786607157],
				[-20.053465032,4.855709248]
				]]
			}
]


with rasterio.Env(), rasterio.open(grid_file_path) as rrr:
	riocrs = rrr.crs
	coords = affine_to_coords(rrr.transform, rrr.width, rrr.height)

## flip the y coordinates so it aligns with how xarray has assigned y coords (asc)
coords['y'].sort()

wrf = xr.open_dataset(grid_file_path, engine = 'cfgrib')
wrf_w_coords = wrf.assign_coords(coords)
wrf_w_coords.rio.write_crs(hrdps_proj4string, inplace=True)

clipped = wrf_w_coords.rio.clip(clip_geom)

[snowman2]

Found a simpler solution. The issue you are encountering is occurring when the CRS is exporting to CF.
This can be disabled using set_options like so:

import rioxarray

clip_geom = [
    {
        "type":"Polygon",
        "coordinates":[[
            [-20.053465032,4.855709248],
            [-19.468843208,4.521462474],
            [-20.068351733,3.435063352],
            [-20.705115362,3.786607157],
            [-20.053465032,4.855709248]
        ]]
    }
]

with rioxarray.set_options(export_grid_mapping=False):
    rds = rioxarray.open_rasterio(
        "grib2_file/CMC_hrdps_west_PRATE_SFC_0_rotated_latlon0.009x0.009_20221018T00Z_P005-00.grib2"
    )
    clipped = rds.rio.clip(clip_geom)

Side note: This will no longer be an issue with the next release of pyproj (pyproj4/pyproj#1167).

[snowman2]

rioxarray fix: #595

[winnerhash]

I want to process GFS's grib2 files using rio.clip to extract data for a specific area using shp files, but I keep getting an error message saying 'No data found in bounds' and it also indicates that the grib2 file has no CRS. I've tried many methods but none of them worked, until I came across your answer, which was the only correct and useful solution. Thank you.