1

I have a GRIB file from NOAA that contains several variables for a specific subregion in Norht-East Atlantic.

I'm using ECMWF grip-api Python library to parse the file. While I'm able to get the values for each {latitude, longitude} pair for each GRIB message with this code:

def iterate_latitude_longitude():
f = open(INPUT)

while 1:
    gid = grib_new_from_file(f)
    if gid is None: break

    iterid = grib_iterator_new(gid,0)

    #Get the missingValue for this GRIB message (e.g. NaN)
    missingValue = grib_get_double(gid,"missingValue")

    i=0
    while 1:
        #Get the value from the geo iterator obtained previously, value is a tuple with [latitude, longitude, value]
        result = grib_iterator_next(iterid)
        if not result: break
        [lat,lon,value] = result

        sys.stdout.write("-Point # %d - lat=%.6f lon=%.6f value=" % (i,lat,lon))
        #Check if missing value is present or not
        if value == missingValue:
            print "missing"
        else:
            print "%.6f" % value
        #increase iterator
        i += 1

    grib_iterator_delete(iterid)
    grib_release(gid)

f.close()

While this is great, this file contains forecast data, it was extracted from this place, more precisely, this file:

multi_1.nww3.t00z.grib2

Which contains world data. There are other files for different times: 00-60-12-18-00, thus, every 6 hours.

What I'm missing is the relation between the value obtained for each latitude,longitude and the time dimension. Does anybody know how to extract this info from a grib file, whether using Python or Java?

I know I can use things like wgrib2 or pygrib to convert to CSV and do the other way around...but it just looks so inefficient

2 Answers 2

1

You probably didn't want to just iterate over the whole file, but rather iterate over each of the (interesting) values, for each of the time steps you want.

Here is some code adapted from the sample code

import traceback
import sys

from gribapi import *

INPUT='multi_1.nww3.t00z.grib2'
VERBOSE=1 # verbose error reporting

def example():
    f = open(INPUT)

    keys = [
        'stepRange',
        'shortName',
        ]

    while 1:
        gid = grib_new_from_file(f)
        if gid is None: break

        for key in keys:
            if not grib_is_defined(gid, key): raise Exception("Key was not defined")
            print '%s=%s' % (key, grib_get(gid, key))

        missingValue = grib_get_double(gid,"missingValue")

        iterid = grib_iterator_new(gid,0)
        i=0
        while 1:
            result = grib_iterator_next(iterid)
            if not result: break

            [lat,lon,value] = result

            sys.stdout.write("- %d - lat=%.6f lon=%.6f value=" % (i,lat,lon))

            if value == missingValue:
                print "missing"
            else:
                print "%.6f" % value

            i += 1

        grib_iterator_delete(iterid) 
        grib_release(gid)

    f.close()

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

        return 1

if __name__ == "__main__":
    sys.exit(main())

And some sample output:

stepRange=0
shortName=unknown
....
- 288 - lat=77.000000 lon=0.000000 value=4.890000
- 289 - lat=77.000000 lon=1.250000 value=5.720000
- 290 - lat=77.000000 lon=2.500000 value=6.690000
- 291 - lat=77.000000 lon=3.750000 value=7.620000
- 292 - lat=77.000000 lon=5.000000 value=8.350000
- 293 - lat=77.000000 lon=6.250000 value=8.810000
- 294 - lat=77.000000 lon=7.500000 value=9.090000
- 295 - lat=77.000000 lon=8.750000 value=9.250000
....
stepRange=3
shortName=ws
....
- 288 - lat=77.000000 lon=0.000000 value=34.030000
- 289 - lat=77.000000 lon=1.250000 value=20.160000
- 290 - lat=77.000000 lon=2.500000 value=4.290000
- 291 - lat=77.000000 lon=3.750000 value=351.300000
- 292 - lat=77.000000 lon=5.000000 value=342.490000
- 293 - lat=77.000000 lon=6.250000 value=337.310000
- 294 - lat=77.000000 lon=7.500000 value=334.990000
- 295 - lat=77.000000 lon=8.750000 value=334.460000
- 296 - lat=77.000000 lon=10.000000 value=334.800000
- 297 - lat=77.000000 lon=11.250001 value=335.990000
- 298 - lat=77.000000 lon=12.500001 value=337.970000
- 299 - lat=77.000000 lon=13.750001 value=338.940000

Note that per your comments on the other answer, the results are in 3 hour blocks (not 6 hour as noted in the question).

1
  • That's a great answer, thanks for the detailed explanation Brad! Apr 11, 2014 at 6:49
1

I haven't worked very much with GRIB files, but I have worked with meteorological forecast data in other file types such as NetCDF. In my experience, there is often a time step associated with a forecast grid and thus, it will be recorded in the file as a forecast hour. The step must then be used along with the timestamp from the header to calculate the exact time of the forecast. An example of a GRIB1 file read using GrADS looks like this:

 1:0:d=04040200:UGRD:kpds5=33:kpds6=100:kpds7=1000:TR=10:P1=0:P2=0:TimeU=1:1000 mb:anl:NAve=0
 2:81534:d=04040200:VGRD:kpds5=34:kpds6=100:kpds7=1000:TR=10:P1=0:P2=0:TimeU=1:1000 mb:anl:NAve=0
 3:154922:d=04040200:UGRD:kpds5=33:kpds6=100:kpds7=850:TR=10:P1=0:P2=0:TimeU=1:850 mb:anl:NAve=0
 4:236456:d=04040200:VGRD:kpds5=34:kpds6=100:kpds7=850:TR=10:P1=0:P2=0:TimeU=1:850 mb:anl:NAve=0
 5:317990:d=04040200:UGRD:kpds5=33:kpds6=100:kpds7=500:TR=10:P1=0:P2=0:TimeU=1:500 mb:anl:NAve=0
 6:399524:d=04040200:VGRD:kpds5=34:kpds6=100:kpds7=500:TR=10:P1=0:P2=0:TimeU=1:500 mb:anl:NAve=0
 7:489202:d=04040200:UGRD:kpds5=33:kpds6=100:kpds7=200:TR=10:P1=0:P2=0:TimeU=1:200 mb:anl:NAve=0
 8:578880:d=04040200:VGRD:kpds5=34:kpds6=100:kpds7=200:TR=10:P1=0:P2=0:TimeU=1:200 mb:anl:NAve=0
 9:660414:d=04040200:UGRD:kpds5=33:kpds6=100:kpds7=1000:TR=10:P1=0:P2=6:TimeU=1:1000 mb:6hr fcst:NAve=0
 10:741948:d=04040200:VGRD:kpds5=34:kpds6=100:kpds7=1000:TR=10:P1=0:P2=6:TimeU=1:1000 mb:6hr fcst:NAve=0

and as you can see, the time dimension is listed as a forecast hour. I have no idea if the GRIB2 file you are working with would look anything like this. I would suggest perhaps contacting someone such as the person who wrote the iterator code on the ECMWF GRIB-API website, Daniel Varela Santoalla.

Here is a link to a good PDF on the basics of the GRIB-API http://www.ecmwf.int/services/computing/training/material/grib_api/grib_api_keys.pdf

2
  • Thanks for pointing to Daniel, jencarta, he has helped me quite a lot! Apr 9, 2014 at 14:37
  • Now for the sake of clarity, what I've done is to get the key named "stepRange", that determines the time step associated with the current GRIB message. Thus, if it is 3, it means that +3 hours have passed from starting time. This way is quite easy to map each value with the time value it belongs to. Apr 9, 2014 at 14:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.