Hi all: as part of using rasterio for OGC API - Coverages work in pygeoapi, a client can pass a bbox parameter as a shortcut to spatially subset a coverage using lat/long coordinates. The idea here, then, is that lat/long coordinates would be reprojected into the data source's native coordinates and then the coverage would be subsetted accordingly.

Having said this, using example files in https://dd.weather.gc.ca/model_hrdps/continental/grib2/06/002, I'm unable to derive the native crs from rasterio in Python:

>>> import rasterio
>>> ds = rasterio.open('CMC_hrdps_continental_HGT_ISBL_0985_ps2.5km_2020102706_P002-00.grib2')
>>> ds.bounds
BoundingBox(left=-0.5, bottom=1455.5, right=2575.5, top=-0.5)
>>> ds.crs
>>> ds.crs is None
True
>>> 

However, when running through rio info, I am able to see bounds and CRS:

$ rio info CMC_hrdps_continental_HGT_ISBL_0985_ps2.5km_2020102706_P002-00.grib2
WARNING:rasterio._env:CPLE_AppDefined in Unable to perform coordinate transformations, so the correct projected geotransform could not be deduced from the lat/long control points.  Defaulting to ungeoreferenced.
{"bounds": [-2099127.4944969374, -5739388.521499627, 4340872.505503062, -2099388.5214996273], "colorinterp": ["undefined"], "count": 1, "crs": "PROJCS[\"unnamed\",GEOGCS[\"Coordinate System imported from GRIB file\",DATUM[\"unknown\",SPHEROID[\"Sphere\",6371229,0]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433]],PROJECTION[\"Polar_Stereographic\"],PARAMETER[\"latitude_of_origin\",60],PARAMETER[\"central_meridian\",252],PARAMETER[\"scale_factor\",1],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"Metre\",1]]", "descriptions": ["98500[Pa] ISBL=\"Isobaric surface\""], "driver": "GRIB", "dtype": "float64", "height": 1456, "indexes": [1], "lnglat": [-92.0404517295466, 52.147885433427845], "mask_flags": [["nodata"]], "nodata": 9999.0, "res": [2500.0, 2500.0], "shape": [1456, 2576], "tiled": false, "transform": [2500.0, 0.0, -2099127.4944969374, 0.0, -2500.0, -2099388.5214996273, 0.0, 0.0, 1.0], "units": [null], "width": 2576}

I need the crs value to be able to reproject the incoming lat/long bbox into native coordinates to subset accordingly. Is there a way to set a given dataset's CRS through code (in readonly mode)?

Thanks

..Tom