Extracting projection from grib files

Specific Domains Geo
1

Any recommendation on how to get the projection from grib files and to robustly create a proj string?

2

That might depend on the grib file and library you’re using to read it. One can use GDAL to read a grib file, that I downloaded from here

julia> using GeoArrays
julia> using Proj
julia> import GeoFormatTypes as GFT

julia> ga = GeoArrays.read("/Users/evetion/Downloads/era5-2mt-2019-03-uk.grib")
Warning: Inside GRIB2Inventory, Message # 745
ERROR: Ran out of file reading SECT0
49x33x744 Array{Float64, 3} with AffineMap([0.25 0.0; 0.0 -0.25], [-10.125, 58.125]) and CRS GEOGCS["Coordinate System imported from GRIB file",DATUM["unnamed",SPHEROID["Sphere",6367470,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST]]

julia> crs(ga)
GeoFormatTypes.WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"Coordinate System imported from GRIB file\",DATUM[\"unnamed\",SPHEROID[\"Sphere\",6367470,0]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST]]")

julia> convert(GFT.ProjString, crs(ga))
GeoFormatTypes.ProjString("+proj=longlat +R=6367470 +no_defs")
2 Likes
3

Another way is with GMT.jl (that also uses GDAL for this)

using GMT

getproj("era5-2mt-2019-03-uk.grib")
"GEOGCS[\"Coordinate System imported from GRIB file\",DATUM[\"unnamed\",SPHEROID[\"Sphere\",6367470,0]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST]]"

prj = wkt2proj(getproj("era5-2mt-2019-03-uk.grib"))
"+proj=longlat +R=6367470 +no_defs"
2 Likes
4

Thanks for the suggestions. I was planning to use GRIBDatasets.jl where I can extract the needed info from a Dict and then create a proj string.

Here is a simple test using GMT.jl with a grib file (186 MB!) relevant to me.

using GMT

file = "MEPS_20230816_0600_h_9.grib2"
julia> wkt2proj(getproj(file))
Warning 3: Cannot find grib2_table_4_2_192_201.csv (GDAL_DATA is not defined)
"+proj=lcc +lat_0=63.3 +lon_0=15 +lat_1=63.3 +lat_2=63.3 +x_0=0 +y_0=0 +R=6371229 +units=m +no_defs"
5

Hmm, that’s more than I know on grib files but you got your wished result, I think.
Funny, that works for me on Windows. Also complains (twice) of not finding the grib2_table_4_2_192_201.csv (GDAL_DATA is not defined) but works, whilst on WSL it doesn’t work

ERROR 1: Cannot find grib2_table_4_5.csv
getdrstemplate: DRS Template 5.42 not defined.
ERROR 4: /mnt/c/v/MEPS_20230816_0600_h_9.grib2 is a grib file, but no raster dataset was successfully identified.

Those grib2_table_4*.csv come from the GDAL installations but the grib2_table_4_2_192_201.csv is not in there (and I have GDAL master build). I wonder why GMT.jl in my Linux installation does not have them.