opening file with forced CRS


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, I'm unable to derive the native crs from rasterio in Python:

>>> import rasterio
>>> ds ='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)
>>> is None

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)?



Alan Snow

This should be fixed in 1.1.5 IIRC: #1248