Closed kevinrosa closed 1 year ago
Dear @kevinrosa Can you share with me your two GRIB files (both grib1 and grib2)? Just one message will be enough as long as it shows this issue.
This will hopefully allow me to reproduce what you are seeing. Many thanks
Thanks @shahramn for following up. 2 GRIBs attached in zip, 418 KB each. sample0.zip
I have reproduced the problem using the ecCodes tools (grib_get_data and grib_filter). The problem is the direction increment in GRIB1 cannot be defined with enough precision: its value should be 0.0015 but it is 0.002 According to the formula: ( lat2 - lat1 ) / (Nj -1 )
ecCodes is not detecting this and returning incorrect latitudes. It should either print a warning or fail. To fix the grib1 file, you can do:
grib_set -s ijDirectionIncrementGiven=0 sample0_grib1.grb fixed.grib1
Once this key is set, ecCodes will compute the increment rather than use what is encoded. To see this we can use grib_dump:
# GRIB edition 1: invalid increment
grib_dump -wcount=1 -Da sample0_grib1.grb | grep jDirectionIncrementInDegrees
88-88 latlon_increment jDirectionIncrementInDegrees = 0.002 [DjInDegrees, DyInDegrees]
# GRIB edition 2: Increment is correct
grib_dump -wcount=1 -Da sample0_grib2.grb | grep jDirectionIncrementInDegrees
109-109 g2latlon jDirectionIncrementInDegrees = 0.001501 [DjInDegrees, DyInDegrees]
# GRIB 1 fixed
grib_dump -wcount=1 -Da fixed.grib | grep jDirectionIncrementInDegrees
88-88 latlon_increment jDirectionIncrementInDegrees = 0.0015 [DjInDegrees, DyInDegrees]
Ah thank you, the increment precision was the issue!
For completeness in closing out this issue, here is the updated version of my initial python code example that now works for GRIB1 and GRIB2:
gid = eccodes.codes_grib_new_from_samples('GRIB1')
eccodes.codes_set(gid, 'dataDate', int(date))
eccodes.codes_set(gid, 'dataTime', int(time))
eccodes.codes_set(gid, 'Ni', len(lon_values))
eccodes.codes_set(gid, 'Nj', len(lat_values))
eccodes.codes_set(gid, 'paramId', 3049)
eccodes.codes_set(gid, 'gridType', 'regular_ll')
eccodes.codes_set(gid, 'latitudeOfFirstGridPointInDegrees', lat_values.max())
eccodes.codes_set(gid, 'latitudeOfLastGridPointInDegrees', lat_values.min())
eccodes.codes_set(gid, 'longitudeOfFirstGridPointInDegrees', lon_values.min())
eccodes.codes_set(gid, 'longitudeOfLastGridPointInDegrees', lon_values.max())
eccodes.codes_set(gid, 'ijDirectionIncrementGiven', 0)
eccodes.codes_set(gid, 'shapeOfTheEarth', 5) # WGS-84 https://codes.ecmwf.int/grib/format/grib2/ctables/3/2/
eccodes.codes_set_array(gid, 'values', np.flipud(u_values[t, :, :]).flatten())
eccodes.codes_write(gid, f)
eccodes.codes_release(gid)
Thanks for letting me know. Have a nice day
I am creating a GRIB file to store data defined on a grid that has uniform spacing of 0.0015 degrees in the i and j directions. The latitude ranges from 30.17 to 30.50 and the longitude ranges from 272.53 to 273.0.
When I create this as a GRIB2 (via
eccodes.codes_grib_new_from_samples('GRIB2')
) and then read the file (viads = xr.open_dataset(fname, engine='cfgrib')
or pygrib), the latitude and longitude values match my initial dataset.But when I create this as a GRIB1, the latitudes get remapped based on what appears to be a uniform i and j direction increment in distance, not degrees. The upper latitude of the first point is still 30.5, but the latitude of the last point becomes 30.06 (instead of 30.17 like it was in the initial dataset. The result is that when I plot the data, the points are out of alignment (in my case, this is ocean model data so it's clear that the coastline is wrong).
Is this the expected behavior for GRIB1? Is there a way to preserve my initial mapping? I can't just switch to GRIB2 due to a compatibility issue.