Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Convert gribapi to eccodes #325

Merged
merged 6 commits into from
Mar 29, 2023
Merged

Convert gribapi to eccodes #325

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
d0cb87f
Convert gribapi references/calls to eccodes.
s-boardman Mar 16, 2023
9dfbe75
pep8 changes
s-boardman Mar 16, 2023
e137886
Correct gribapi references in comments.
s-boardman Mar 16, 2023
4d73df5
Swap references to gribapi to ecCodes in docstrings.
s-boardman Mar 21, 2023
6aca8c8
Correct typo in comment.
s-boardman Mar 28, 2023
012d43f
Swap from GribInternalError to CodesInternalError for consistency
s-boardman Mar 28, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
59 changes: 33 additions & 26 deletions iris_grib/__init__.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@

import cartopy.crs as ccrs
import cf_units
import gribapi
import eccodes
import numpy as np
import numpy.ma as ma

Expand Down Expand Up @@ -98,9 +98,11 @@ def ndim(self):
def __getitem__(self, keys):
with open(self.path, 'rb') as grib_fh:
grib_fh.seek(self.offset)
grib_message = gribapi.grib_new_from_file(grib_fh)
grib_message = eccodes.codes_new_from_file(
grib_fh, eccodes.CODES_PRODUCT_GRIB
)
data = _message_values(grib_message, self.shape)
gribapi.grib_release(grib_message)
eccodes.codes_release(grib_message)

result = data.__getitem__(keys)

Expand Down Expand Up @@ -128,6 +130,7 @@ class GribWrapper:
provides alternative means of working with GRIB message instances.

"""

def __init__(self, grib_message, grib_fh=None):
"""Store the grib message and compute our extra keys."""
self.grib_message = grib_message
Expand All @@ -145,7 +148,8 @@ def __init__(self, grib_message, grib_fh=None):
# Note that, the grib-api has already read this message and
# advanced the file pointer to the end of the message.
offset = grib_fh.tell()
message_length = gribapi.grib_get_long(grib_message, 'totalLength')
message_length = eccodes.codes_get_long(
grib_message, 'totalLength')

# Initialise the key-extension dictionary.
# NOTE: this attribute *must* exist, or the the __getattr__ overload
Expand All @@ -155,12 +159,12 @@ def __init__(self, grib_message, grib_fh=None):
self._compute_extra_keys()

# Calculate the data payload shape.
shape = (gribapi.grib_get_long(grib_message, 'numberOfValues'),)
shape = (eccodes.codes_get_long(grib_message, 'numberOfValues'),)

if not self.gridType.startswith('reduced'):
ni, nj = self.Ni, self.Nj
j_fast = gribapi.grib_get_long(grib_message,
'jPointsAreConsecutive')
j_fast = eccodes.codes_get_long(grib_message,
'jPointsAreConsecutive')
shape = (nj, ni) if j_fast == 0 else (ni, nj)

if deferred:
Expand Down Expand Up @@ -191,27 +195,30 @@ def __getattr__(self, key):
# we just get <type 'float'> as the type of the "values"
# array...special case here...
if key in ["values", "pv", "latitudes", "longitudes"]:
res = gribapi.grib_get_double_array(self.grib_message, key)
res = eccodes.codes_get_double_array(self.grib_message, key)
elif key in ('typeOfFirstFixedSurface',
'typeOfSecondFixedSurface'):
res = np.int32(gribapi.grib_get_long(self.grib_message, key))
res = np.int32(eccodes.codes_get_long(self.grib_message, key))
else:
key_type = gribapi.grib_get_native_type(self.grib_message, key)
key_type = eccodes.codes_get_native_type(
self.grib_message, key)
if key_type == int:
res = np.int32(gribapi.grib_get_long(self.grib_message,
key))
res = np.int32(eccodes.codes_get_long(self.grib_message,
key))
elif key_type == float:
# Because some computer keys are floats, like
# longitudeOfFirstGridPointInDegrees, a float32
# is not always enough...
res = np.float64(gribapi.grib_get_double(self.grib_message,
key))
res = np.float64(eccodes.codes_get_double(
self.grib_message, key
)
)
elif key_type == str:
res = gribapi.grib_get_string(self.grib_message, key)
res = eccodes.codes_get_string(self.grib_message, key)
else:
emsg = "Unknown type for {} : {}"
raise ValueError(emsg.format(key, str(key_type)))
except gribapi.errors.GribInternalError:
except eccodes.CodesInternalError:
res = None

# ...or is it in our list of extras?
Expand Down Expand Up @@ -258,7 +265,7 @@ def _compute_extra_keys(self):
longitudeOfSouthernPoleInDegrees = 0.0
latitudeOfSouthernPoleInDegrees = 90.0

centre = gribapi.grib_get_string(self.grib_message, "centre")
centre = eccodes.codes_get_string(self.grib_message, "centre")

# default values
self.extra_keys = {'_referenceDateTime': -1.0,
Expand Down Expand Up @@ -287,7 +294,7 @@ def _compute_extra_keys(self):

# cf phenomenon translation
# Get centre code (N.B. self.centre has default type = string)
centre_number = gribapi.grib_get_long(self.grib_message, "centre")
centre_number = eccodes.codes_get_long(self.grib_message, "centre")
# Look for a known grib1-to-cf translation (or None).
cf_data = gptx.grib1_phenom_to_cf_info(
table2_version=self.table2Version,
Expand Down Expand Up @@ -345,7 +352,7 @@ def _compute_extra_keys(self):
# Earth assumed spherical with radius 6367.47 km
geoid = coord_systems.GeogCS(semi_major_axis=6367470)

gridType = gribapi.grib_get_string(self.grib_message, "gridType")
gridType = eccodes.codes_get_string(self.grib_message, "gridType")

if gridType in ["regular_ll", "regular_gg", "reduced_ll",
"reduced_gg"]:
Expand Down Expand Up @@ -416,8 +423,8 @@ def _compute_extra_keys(self):
# get the distinct latitudes, and make sure they are sorted
# (south-to-north) and then put them in the right direction
# depending on the scan direction
latitude_points = gribapi.grib_get_double_array(
self.grib_message, 'distinctLatitudes').astype(np.float64)
latitude_points = eccodes.codes_get_double_array(
self.grib_message, 'distinctLatitudes').astype(np.float64)
latitude_points.sort()
if not self.jScansPositively:
# we require latitudes north-to-south
Expand Down Expand Up @@ -639,8 +646,8 @@ def _longitude_is_cyclic(points):


def _message_values(grib_message, shape):
gribapi.grib_set_double(grib_message, 'missingValue', np.nan)
data = gribapi.grib_get_double_array(grib_message, 'values')
eccodes.codes_set_double(grib_message, 'missingValue', np.nan)
data = eccodes.codes_get_double_array(grib_message, 'values')
data = data.reshape(shape)

# Handle missing values in a sensible way.
Expand Down Expand Up @@ -787,7 +794,7 @@ def save_pairs_from_cube(cube):

# Save each latlon slice2D in the cube
for slice2D in cube.slices([y_coords[0], x_coords[0]]):
grib_message = gribapi.grib_new_from_samples("GRIB2")
grib_message = eccodes.codes_grib_new_from_samples("GRIB2")
_save_rules.run(slice2D, grib_message, cube)
yield (slice2D, grib_message)

Expand Down Expand Up @@ -824,8 +831,8 @@ def save_messages(messages, target, append=False):

try:
for message in messages:
gribapi.grib_write(message, grib_file)
gribapi.grib_release(message)
eccodes.codes_write(message, grib_file)
eccodes.codes_release(message)
finally:
# (this bit is common to the pp and grib savers...)
if isinstance(target, str):
Expand Down
4 changes: 2 additions & 2 deletions iris_grib/_load_convert.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -204,7 +204,7 @@ def _masker(item):
return result


# Use ECCodes gribapi to recognise missing value
# Use ECCodes to recognise missing value
_MDI = None


Expand Down Expand Up @@ -1875,7 +1875,7 @@ def coord_timedelta(coord, value):
def time_coords(section, metadata, rt_coord):
if 'forecastTime' in section.keys():
forecast_time = section['forecastTime']
# The gribapi encodes the forecast time as 'startStep' for pdt 4.4x;
# ecCodes encodes the forecast time as 'startStep' for pdt 4.4x;
# product_definition_template_40 makes use of this function. The
# following will be removed once the suspected bug is fixed.
elif 'startStep' in section.keys():
Expand Down
Loading