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).
To get the appropriate cell coordinates from your latitude and longitude you need to know the coverage of the netCDF file (for example, the coordinate of the upper left cell and the size of the cells in the x and y directions). @Luke 's solution below works for a srtaight x/y raster with no cell rotation. If you don't have that, you could look at the affine library to transform latitude and longitude to cell coordinates. Given GDAL GeoTransform in the form:
- Upper left cell
x coordinate
x cell size in projected unitsx directional rotation (usually 0)- Upper left cell
y coordinate
y directional rotation (usually 0)- Negative
y cell size in projected units
So for example the -237481.5, 425.0, 0.0, 237536.4, 0.0, -425.0 would translate as an upper left cell of the raster with coordinates -237481.5, 237536.4, and a 425.0 unit square with no rotation.
Using the affine library you can transform this into a transformation object like so:
from affine import Affine
aff = Affine.from_gdal(-237481.5, 425.0, 0.0, 237536.4, 0.0, -425.0)
Which can transform cell coordinates to projected coordinates:
xs = np.array([0, 1, 2])
ys = np.array([3, 4, 5])
x_coords, y_coords = aff * (xs, ys)
Or (in your case) from coordinates back to cell coordinates:
xs, ys = ~aff * (np.array(lons), np.array(lats))
These values are floats, so you'll need to transform them into integers to get cell coordinates you can use.
xs = np.round(xs).astype(np.int)
ys = np.round(ys).astype(np.int)
You can then use these as indexes to your netCDF4 array (use the latest version of netCDF4 - 1.2.1 as this means you don't need to sort or remove duplicates).
variable = fnc.variables['variable'][xs, ys] # be careful of your index order
This will return a square array due to the slightly different netCDF indexing, but you can get the actual values that you're after as the diagnoal:
values = variable.diagonal()
Best Answer
I don't know much at all about writing NetCDF files, but the following allowed me to write a NetCDF file that GDAL recognised the CRS.
gdalinfooutput: