Thursday 27 December 2018

python - Getting time dimension from GRIB file?


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



Answer



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



No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...