How to extract specific information from GRIB files?

Question

How to extract fetched GRIB data with Python and gribapi package (link goes to official ECMWF Wiki)? I tried to follow few examples from their "documentation" and do it myself, but I just cannot figure out how to retrieve only specific parameter (e.g. surface temperature or wind) for given latitude and longitude. I don't need to plot the data or visualizations, I just need to pass few parameters when reading the file and get numbers back.

So far, I tried to download original GRIB2 file from www.ftp.ncep.noaa.gov and run few Python API example scripts from mentioned ECMWF Wiki page. Since Python API documentation doesn't exist, I have some hard time to understand which methods should I use in order to:

• read the file,
• select only specific "message" (e.g. surface temperature only) if there's more "messages" in a file,
• print out the value for a specific location (using grib_find_nearest method).

Example code which I have so far, mostly copied from ECMWF examples and this answer from Getting time dimension from GRIB file?:

import traceback
import sys

# hopefully the only external package I need
from gribapi import *

# file with all possible 'messages'
INPUT = '2015121406_gfs.t06z.pgrb2.0p25.f000'
VERBOSE=1 # verbose error reporting

def example():
    f = open(INPUT)

    # location we are interested in
    lat, lon = 64.1353, -21.8952

    while 1:
        # STEP 1:
        # open downloaded GRIB2 file
        gid = grib_new_from_file(f)
        if gid is None: break

        # define the iterator (which is throughout the program the same?)
        iterid = grib_iterator_new(gid,0)

        # get the result for the nearest location
        nearest = grib_find_nearest(gid, lat, lon)[0]

        while 1:
            # STEP 2:
            # ???
            #
            # the loop goes through the whole file
            # instead just selecting messages we need beforehand...
            # fictional function which selects
            # 'TMP' (temperature) on 'surface' messages only
            # without need to iterate all of them:
            # 
            # result = grib_select_specific_message(gid, 'TMP', 'surface')

            result = grib_iterator_next(iterid)
            if not result: break

        # STEP 3:
        # profit!
        # result variable returns numbers which I process in any way I need, yay!

        # more undocumented stuff,
        # put here just because examples do it the same way
        grib_iterator_delete(iterid)
        grib_release(gid)

    f.close()

# main program function
def main():
    try:
        example()
    except GribInternalError,err:
        if VERBOSE:
            traceback.print_exc(file=sys.stderr)
        else:
            print >>sys.stderr,err.msg

        return 1

# run the program
if __name__ == "__main__":
    sys.exit(main())

The problem with this code is that it iterates trough the whole file, all the messages and all possible coordinates (whole planet, 0.25 degrees precision). If file is 300+ MB it takes a while to read it. I feel like using "iterator" here (whatever that is) is a bit wrong since I need to "cherry-pick" only specific information (few weather parameters for a particular location only) and grib_find_nearest() function doesn't need an iterator at all to select only part of the data.
This could be half-solved with filtering the data I need before downloading the file (results in much smaller file) but I'm still not sure how to do it (part of Downloading GRIB GFS files with specific filters?), but still I'd like to figure this out since it might happen that I will occasionally have to download "full" files.

Answer 1

Another idea would be to use pygrib module:

    import pygrib

    grbs= pygrib.open("my_file.grb")

    # use grbs.select to select the grids you are interested in (shortName, typeOfLevel, level=500, marsParam, dataDate, dataTime, ...)

    DATA=np.array(grbs.select(marsParam=my_param,dataDate=my_date,dataTime=my_time))

    # DATA will contain 3 arrays
    # DATA[0] for values
    # DATA[1] for longitudes
    # DATA[2] for latitudes

    # from the "values" array, extract in lon and lat
    DATA1=DATA[0].data(lat1=my_y1,lat2=my_y2,lon1=my_x1,lon2=my_x2))

Answer 2

I eventually ended up with ditching gribapi and switching to Met Office's Iris Python package which solves my problem in a very elegant way. Although it was a bit of a pain to install it (dependencies could be really tricky) at least the documentation is good and it is really easy to use it.

Answer 3

Using pygrib was eventually easy once I figured out the keys() of the .grib file I had. One thing I noticed is that the keys which have specific data are the same keys I used in the API query to download the weather data.

I downloaded my data from the Climate Data Store and the ECMWF parameter database is a good source to see parameter names (the names of the parameters are slightly different than that used to pass through the API query).

wash_2014_grib = "Washington2014Data.grib" # Local downloaded grib file
grbs = pygrib.open(wash_2014_grib)

grb = grbs.select(name="2 metre temperature",year=2014,month=5,day=1,time=500)[0]
temporal_data = grb.values

The above will extract 2 m temperature data for May 1, 2014 at 5:00. (Note that the .select method of pygrib does not permit leading zeros in decimal integer literals. Without calling the zeroeth index ([0]), it returns a list of pygribmessages, which can be accessed iteratively or simply by calling the index.

In the above output, grb is a pygribmessage object and temporal_data is an array of values of temperatures at every coordinate.

To extract specific coordinate data,

y1,y2,x1,x2 = 46,46.5,-120.5,-120
local_data,lats,lons = grb.data(lat1=y1,lat2=y2,lon1=x1,lon2=x2)

which gives

>>> local_data
array([[287.25024414, 289.18969727, 292.21118164],
       [290.43579102, 293.13500977, 292.69360352],
       [285.51391602, 288.09790039, 291.69165039]])
>>> lats
array([[46.5 , 46.5 , 46.5 ],
       [46.25, 46.25, 46.25],
       [46.  , 46.  , 46.  ]])
>>> lons
array([[-120.5 , -120.25, -120.  ],
       [-120.5 , -120.25, -120.  ],
       [-120.5 , -120.25, -120.  ]])

All set for further analysis!

Answer 4

Try using xarray with cfgrib (these can be installed via pip or conda).

See the GRIB Data Example to view and extract data from an example file.

import xarray as xr

ds = xr.open_dataset("path/to/era5-2mt-2019-03-uk.grib", engine="cfgrib")

# extract a timeseries from a coordinate
ds.t2m.sel(longitude=0, latitude=51.5)