Rasterio: Georeference metadata not read from some GRIB files

[nawendt]

Expected behavior and actual behavior.

I have been working with GRIB (version 2) data and have run into an issue where the georeference metadata are not extracted (i.e., the CRS is missing, affine parameters are incorrect, etc.). However, if I simply use the standard GDAL Python bindings, the metadata are extracted. The data I am working with is from the NCEP SREF weather forecast model.

Steps to reproduce the problem.

I use rasterio as follows:

import rasterio
dat = rasterio.open('my_grib_file.grib2')
dat.profile

which produces the following output

>>> {'affine': Affine(1.0, 0.0, -0.5, 0.0, 1.0, -0.5),
 'count': 523,
 'crs': CRS({}),
 'driver': 'GRIB',
 'dtype': 'float64',
 'height': 129,
 'nodata': None,
 'tiled': False,
 'transform': (-0.5, 1.0, 0.0, -0.5, 0.0, 1.0),
 'width': 185}

With osgeo.gdal (in a separate session)

from osgeo import gdal
dat = gdal.Open('my_grib_file.grib2')
dat.GetProjectionRef()
dat.GetGeoTransform()

which produces

>>> 'PROJCS["unnamed",GEOGCS["Coordinate System imported from GRIB file",DATUM["unknown",SPHEROID["Sphere",6371229,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["standard_parallel_1",25],PARAMETER["standard_parallel_2",25],PARAMETER["latitude_of_origin",25],PARAMETER["central_meridian",265],PARAMETER["false_easting",0],PARAMETER["false_northing",0]]'
>>> (-4246424.496915466, 40635.0, 0.0, 4388899.238982436, 0.0, -40635.0)

Operating systems

RHEL 6

Rasterio version and provenance

rasterio 0.36.0 via conda-forge
gdal 2.1.3 via conda-forge
python 3.6.4 via conda-forge

Other information

I have verified that both osgeo.gdal and rasterio are using the same GDAL C library. When I use the same methods as above on another GRIB (this time version 1) from a different weather forecast model, the georeference metadata are read correctly. This would indicate it could be an issue with individual files, but, as I show above, the same GDAL library reads the correct georeference data when not using rasterio.

I determined that the issue occurs when GDALGetProjectionRef is called in read_crs from _base.pyx. What I did then was to create a small C program that uses GDALGetProjectionRef to access the georeference information and print it to STDOUT. That worked just fine and produced the same information as when using osgeo.gdal. I also built rasterio from source in a separate conda environment and still had the same problem. This is what ultimately led me to believe that the issue has something to do with what is getting passed around in rasterio versus the GDAL C library itself or something with my conda environments.

[sgillies]

I think it may be that rasterio 0.36 is in error. There are more recent 1.0 pre-releases in conda-forge. Could you try the latest one?

Also, can you link a specific GRIB file here that will help me reproduce the problem? I'm not sure how to navigate the NCEP FTP site.

[jdmcbr]

@nawendt To install the pre-release version, use conda install -c conda-forge/label/dev rasterio.

[nawendt]

@sgillies Here is a link to a SREF grib file. They are only kept on the server 6 days or so. I had no better luck using rasterio 1.0a12.

[snowman2]

This is a very strange issue indeed. I confirmed the same behavior with the file sref_arw.t15z.pgrb132.ctl.f00.grib2 when opening with rasterio (1.0.2) & with the GDAL python bindings.

However, when I do rio info sref_arw.t15z.pgrb132.ctl.f00.grib2, I get:

...
"transform": [16232.0, 0.0, -5640783.863757296, 0.0, -16232.0, 4355613.353550392, 0.0, 0.0, 1.0]
...
"crs": "+a=6371229 +b=6371229 +lat_0=50 +lat_1=50 +lat_2=50 +lon_0=253 +no_defs +proj=lcc +units=m +x_0=0 +y_0=0"
...

Which matches gdalinfo:

PROJCS["unnamed",
    GEOGCS["Coordinate System imported from GRIB file",
        DATUM["unknown",
            SPHEROID["Sphere",6371229,0]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Lambert_Conformal_Conic_2SP"],
    PARAMETER["standard_parallel_1",50],
    PARAMETER["standard_parallel_2",50],
    PARAMETER["latitude_of_origin",50],
    PARAMETER["central_meridian",253],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0]]
Origin = (-5640783.863757296465337,4355613.353550392203033)
Pixel Size = (16232.000000000000000,-16232.000000000000000)

[snowman2]

Looks like it is needing the environment for some reason. It works like this:

with rasterio.Env(), rasterio.open('sref_arw.t15z.pgrb132.ctl.f00.grib2') as dat:
    print(dat.profile)

Output:

{'driver': 'GRIB', 'dtype': 'float64', 'nodata': None, 'width': 697, 'height': 553, 'count': 551, 'crs': CRS({'proj': 'lcc', 'lat_1': 50, 'lat_2': 50, 'lat_0': 50, 'lon_0': 253, 'x_0': 0, 'y_0': 0, 'a': 6371229, 'b': 6371229, 'units': 'm', 'no_defs': True}), 'transform': Affine(16232.0, 0.0, -5640783.863757296,
       0.0, -16232.0, 4355613.353550392), 'tiled': False}

[sgillies]

Confirmed.

>>> with rasterio.open('/Users/sean/Downloads/nam.t00z.awip3d00.tm02.grib2') as src:
...     print(src.crs)
...
None
>>> with rasterio.Env(), rasterio.open('/Users/sean/Downloads/nam.t00z.awip3d00.tm02.grib2') as src:
...     print(src.crs)
...
+a=6371229 +b=6371229 +lat_0=25 +lat_1=25 +lat_2=25 +lon_0=265 +no_defs +proj=lcc +units=m +x_0=0 +y_0=0

Interesting! 🤔

[drnextgis]

It is still persists.

[sgillies]

This is a matter of GDAL being able to find PROJ data. Setting PROJ_LIB in our environments is the only means we have of configuring GDAL 2.x.

With PROJ_LIB set, this is what I get when I open one of GDAL's test files.

>>> import rasterio
>>> with rasterio.open("/home/sean/projects/gdal/autotest/gdrivers/data/grib/one_one.grib2") as dataset:
...     print(dataset.crs)
...
GEOGCS["Coordinate System imported from GRIB file",DATUM["unnamed",SPHEROID["Sphere",6371229,0]],PRIMEM["Greenwich",0],UNIT["degree",
0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST]]

With an installed rasterio 1.0.28 wheel and no PROJ_LIB set I get

>>> import rasterio
>>> with rasterio.open("/home/sean/projects/gdal/autotest/gdrivers/data/grib/one_one.grib2") as dataset:
...     print(dataset.crs)
...
GEOGCS["Coordinate System imported from GRIB file",DATUM["unknown",SPHEROID["Sphere",6371229,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]]

It's not empty like you report, but is missing data.

Calling rasterio.open() does import rasterio.env and that is supposed to patch our environments appropriately in the wheel case, setting PROJ_LIB. It looks like that is not working. I'll investigate.

[snowman2]

Just ran into this again recently: https://gis.stackexchange.com/a/357686/144357

[sgillies]

@snowman2 it's not a rasterio regression. If PROJ_LIB is set, there won't be a problem.

[snowman2]

I just tested that as I was pretty sure PROJ_LIB was set:

[snowman2]

Just did some debugging and it seems to have something to do with defaults used in Env.from_defaults():

>>> import rasterio
>>> with rasterio.Env.from_defaults(), rasterio.open("CMC_caps_HGT_ISBL_0050_ps3km_2018090900_P005.grib2") as rds:
...     print(rds.crs)
... 
None

>>> import rasterio
>>> with rasterio.Env(), rasterio.open("CMC_caps_HGT_ISBL_0050_ps3km_2018090900_P005.grib2") as rds:
...     print(rds.crs)
... 
CRS.from_wkt('PROJCS["unnamed",GEOGCS["Coordinate System imported from GRIB file",DATUM["unnamed",SPHEROID["Sphere",6371229,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]]],PROJECTION["Polar_Stereographic"],PARAMETER["latitude_of_origin",60],PARAMETER["central_meridian",264],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Metre",1],AXIS["Easting",SOUTH],AXIS["Northing",SOUTH]]')

I modified this part:
https://github.com/mapbox/rasterio/blob/5af075c1c0f260100e831b281909e8e5dee7f571/rasterio/env.py#L226-L228
To return just the Env object without the default options and it works:

>>> import rasterio
>>> rds = rasterio.open("CMC_caps_HGT_ISBL_0050_ps3km_2018090900_P005.grib2")
>>> rds.crs
CRS.from_wkt('PROJCS["unnamed",GEOGCS["Coordinate System imported from GRIB file",DATUM["unnamed",SPHEROID["Sphere",6371229,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]]],PROJECTION["Polar_Stereographic"],PARAMETER["latitude_of_origin",60],PARAMETER["central_meridian",264],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Metre",1],AXIS["Easting",SOUTH],AXIS["Northing",SOUTH]]')

[snowman2]

I think it has to do with this option:
https://github.com/mapbox/rasterio/blob/5af075c1c0f260100e831b281909e8e5dee7f571/rasterio/env.py#L101

>>> import rasterio
>>> rds = rasterio.open("CMC_caps_HGT_ISBL_0050_ps3km_2018090900_P005.grib2")
>>> rds.crs
CRS.from_wkt('PROJCS["unnamed",GEOGCS["Coordinate System imported from GRIB file",DATUM["unnamed",SPHEROID["Sphere",6371229,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]]],PROJECTION["Polar_Stereographic"],PARAMETER["latitude_of_origin",60],PARAMETER["central_meridian",264],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Metre",1],AXIS["Easting",SOUTH],AXIS["Northing",SOUTH]]')

[snowman2]

More debugging:

>>> import rasterio
>>> with rasterio.Env(CHECK_WITH_INVERT_PROJ=False), rasterio.open("CMC_caps_HGT_ISBL_0050_ps3km_2018090900_P005.grib2") as rds:
...     print(rds.crs)
... 
PROJCS["unnamed",GEOGCS["Coordinate System imported from GRIB file",DATUM["unnamed",SPHEROID["Sphere",6371229,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]]],PROJECTION["Polar_Stereographic"],PARAMETER["latitude_of_origin",60],PARAMETER["central_meridian",264],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Metre",1],AXIS["Easting",SOUTH],AXIS["Northing",SOUTH]]
>>> with rasterio.Env(CHECK_WITH_INVERT_PROJ=True), rasterio.open("CMC_caps_HGT_ISBL_0050_ps3km_2018090900_P005.grib2") as rds:
...     print(rds.crs)
... 
None

[snowman2]

Using GDAL Python, we get a useful error message:

>>> from osgeo import gdal
>>> gdal.SetConfigOption("CHECK_WITH_INVERT_PROJ", "YES")
>>> dat = gdal.Open("CMC_caps_HGT_ISBL_0050_ps3km_2018090900_P005.grib2")
Warning 1: Unable to perform coordinate transformations, so the correct projected geotransform could not be deduced from the lat/long control points.  Defaulting to ungeoreferenced.
>>> dat.GetProjectionRef()
''

[sgillies]

Interesting :thinking: By the way, you see an error message with the Python bindings because it prints to stderr by default regardless of error level. Rasterio sends this to the logger, which has to be configured.

[sgillies]

@snowman2 without setting CHECK_WITH_INVERT_PROJ in that script, do you get the results you expect?

And what is up with this GRIB file that PROJ functions are called when it is opened? Generally, Open does very little work at all, but the GRIB driver must be an exception.

[snowman2]

without setting CHECK_WITH_INVERT_PROJ in that script, do you get the results you expect?

I just checked that the CRS and geotransform were provided.

And what is up with this GRIB file that PROJ functions are called when it is opened?

:man_shrugging:. I guess GRIB stores things in a different way.

[sgillies]

Resolved in #1942. Thanks for your persistence and patience @snowman2 !