iris-grib
iris-grib copied to clipboard
SciTools
Grid Type 3.12 Transverse Mercator not loading
Am trying to load a grib2 file and am getting the error 'iris.exceptions.TranslationError: Grid definition section 3 contains unsupported source of grid definition [12]' Grib header for record 1 listed below..
Grib header is:
GRIB {
# Meteorological products (grib2/tables/4/0.0.table)
discipline = 0;
editionNumber = 2;
# U.K. Met Office - Exeter (common/c-11.table)
centre = 74;
subCentre = 0;
# Start of forecast (grib2/tables/4/1.2.table)
significanceOfReferenceTime = 1;
dataDate = 20190225;
dataTime = 1200;
# Operational products (grib2/tables/4/1.3.table)
productionStatusOfProcessedData = 0;
# Analysis and forecast products (grib2/tables/4/1.4.table)
typeOfProcessedData = 2;
numberOfDataPoints = 385792;
# There is no appended list (grib2/tables/4/3.11.table)
interpretationOfNumberOfPoints = 0;
# Unknown code table entry (grib2/tables/4/3.1.table)
gridDefinitionTemplateNumber = 12;
# Earth assumed oblate spheroid with major and minor axes specified by data producer (grib2/tables/4/3.2.table)
shapeOfTheEarth = 3;
Ni = 548;
Nj = 704;
latitudeOfReferencePointInDegrees = 4.9e-05;
longitudeOfReferencePointInDegrees = -2e-06;
XRInMetres = 400000;
YRInMetres = -100000;
iScansNegatively = 0;
jScansPositively = 1;
jPointsAreConsecutive = 0;
alternativeRowScanning = 0;
DiInMetres = 2000;
DjInMetres = 2000;
X1InGridLengths = -238000;
Y1InGridLengths = 1.222e+06;
X2InGridLengths = 856000;
Y2InGridLengths = -184000;
gridType = transverse_mercator;
NV = 0;
# Analysis or forecast at a horizontal level or in a horizontal layer at a point in time (grib2/tables/4/4.0.table)
productDefinitionTemplateNumber = 0;
# Temperature (grib2/tables/4/4.1.0.table)
parameterCategory = 0;
# Temperature (K) (grib2/tables/4/4.2.0.0.table)
parameterNumber = 0;
#-READ ONLY- parameterUnits = K;
#-READ ONLY- parameterName = Temperature ;
# Analysis (grib2/tables/4/4.3.table)
typeOfGeneratingProcess = 0;
generatingProcessIdentifier = 1;
# Minute (grib2/tables/4/4.4.table)
indicatorOfUnitOfTimeRange = 0;
# Hour (stepUnits.table)
stepUnits = 1;
forecastTime = 0;
stepRange = 0;
# Specified height level above ground (m) (grib2/tables/4/4.5.table)
typeOfFirstFixedSurface = 103;
#-READ ONLY- unitsOfFirstFixedSurface = m;
scaleFactorOfFirstFixedSurface = 0;
scaledValueOfFirstFixedSurface = 1;
# Missing (grib2/tables/4/4.5.table)
typeOfSecondFixedSurface = 255;
#-READ ONLY- unitsOfSecondFixedSurface = unknown;
#-READ ONLY- nameOfSecondFixedSurface = Missing;
scaleFactorOfSecondFixedSurface = MISSING;
scaledValueOfSecondFixedSurface = MISSING;
level = 1;
shortNameECMF = t;
shortName = t;
nameECMF = Temperature;
name = Temperature;
cfNameECMF = air_temperature;
cfName = air_temperature;
cfVarNameECMF = t;
cfVarName = t;
#-READ ONLY- modelName = unknown;
numberOfValues = 385792;
packingType = grid_simple;
This seems to be a missing feature. If you want to give this a shot yourself, have a look at #95 to get an idea what might be needed.
thanks for your help - not sure if I'm going to be able to do this but will have a look
In any case, thanks for taking the time to report this! Maybe @pelson or @DPeterK have some more words of wisdom?
So GDT 12 (Transverse Mercator) is definitely supported by iris-grib – however that's not what the problem is here, despite initial appearances. Instead it's the sourceOfGridDefinition key that's causing issue. Only a value of 0 is supported by iris-grib for this flag.
Digging in the GRIB spec code tables suggests that this key should not be set to 12, as 12 is a reserved value. It's only 0 and 1 that should typically be used for this key, where 0 means use the existing grid definition templates and 1 means use a separate set of grid definition templates. In fact, given the coincidence of this value being the same as the GDT template number for Transverse Mercator, it would be my guess that this key has been set by accident as well as the GDT Template Number.
Given that, I reckon this GRIB message is slightly malformed. I'd recommend editing the GRIB message using one of the commandline GRIB operators (probably grib-set is the one you want) to edit the value of this key and set it back to 0. At that point you should be able to load this message with iris-grib.
Hi, thanks for your help - I think I have a pretty bad grib file (from a major supplier) I've changed that flag and now get: >>> import iris
levels = iris.load("test2.grib") Traceback (most recent call last): File "
", line 1, in File "/usr/local/lib/python2.7/site-packages/iris/init.py", line 343, in load return _load_collection(uris, constraints, callback).merged().cubes() File "/usr/local/lib/python2.7/site-packages/iris/init.py", line 313, in _load_collection result = iris.cube._CubeFilterCollection.from_cubes(cubes, constraints) File "/usr/local/lib/python2.7/site-packages/iris/cube.py", line 145, in from_cubes for cube in cubes: File "/usr/local/lib/python2.7/site-packages/iris/init.py", line 300, in _generate_cubes for cube in iris.io.load_files(part_names, callback, constraints): File "/usr/local/lib/python2.7/site-packages/iris/io/init.py", line 210, in load_files for cube in handling_format_spec.handler(fnames, callback): File "/usr/local/lib/python2.7/site-packages/iris/fileformats/rules.py", line 429, in load_cubes user_callback_wrapper=loadcubes_user_callback_wrapper): File "/usr/local/lib/python2.7/site-packages/iris/fileformats/rules.py", line 350, in _load_pairs_from_fields_and_filenames cube, factories, references = _make_cube(field, converter) File "/usr/local/lib/python2.7/site-packages/iris/fileformats/rules.py", line 293, in _make_cube metadata = converter(field) File "/usr/local/lib/python2.7/site-packages/iris_grib/_load_convert.py", line 2515, in convert grib2_convert(field, metadata) File "/usr/local/lib/python2.7/site-packages/iris_grib/_load_convert.py", line 2464, in grib2_convert grid_definition_section(field.sections[3], metadata) File "/usr/local/lib/python2.7/site-packages/iris_grib/_load_convert.py", line 1269, in grid_definition_section grid_definition_template_12(section, metadata) File "/usr/local/lib/python2.7/site-packages/iris_grib/_load_convert.py", line 820, in grid_definition_template_12 check_range(y1, y2, section['Nj'], section['Dj']) File "/usr/local/lib/python2.7/site-packages/iris_grib/_load_convert.py", line 818, in check_range raise TranslationError('Inconsistent grid definition') iris.exceptions.TranslationError: Inconsistent grid definition
Looks like it doesn't like the y values in the file in that order
Hi! Thanks loads to the iris-grib developers for building this great tool :)
I'm getting identical errors to @karldevinenetweather when using iris to load 201701010000_u1096_ng_umqv_Wholesale1.grib from JASMIN
Like @karldevinenetweather, I tried using grib_set -s sourceOfGridDefinition=0 201701010000_u1096_ng_umqv_Wholesale1.grib 201701010000_u1096_ng_umqv_Wholesale1_mod.grib
Loading the modified file with iris results in TranslationError: Inconsistent grid definition.
I poked around using the debugger and found these values:
%debug
> /Users/JackKelly/miniconda3/envs/nwp/lib/python3.7/site-packages/iris_grib/_load_convert.py(837)check_range()
835 max_last = v1 + (n - 1) * (d + 1)
836 if not (min_last < v2 < max_last):
--> 837 raise TranslationError('Inconsistent grid definition')
838 check_range(x1, x2, section['Ni'], section['Di'])
839 check_range(y1, y2, section['Nj'], section['Dj'])
ipdb> min_last
262799297
ipdb> max_last
262800703
ipdb> v2
-18400000
I tried commenting out line 837 raise TranslationError('Inconsistent grid definition') but then I got a third error: TranslationError: Product definition template [5] is not supported.
I also tried using cfgrib. cfgrib loads the grib file, but the loaded data has no geospatial coordinates: cfgrib complains that ecCodes provides no latitudes/longitudes for gridType='transverse_mercator'
What's the best way forwards?
update I get the same problem with more recent UKV GRIB files. I tried 201908290000_u1096_ng_umqv_Wholesale1.grib and got the same issues as above.
Hi @karldevinenetweather and @JackKelly,
I don't suppose you have had any further luck with this? I'm also working with Met Office files and running into the same issue.
I'm currently using the following workaround (snippet below) to reshape the data into the right structure, but the values used for the longitude and latitude indexes are made up. Do you know how the indexes should be reconstructed?
Best, Ayrton
longitude_size = 704
latitude_size = 548
ds_reshaped = (
ds.assign_coords(
longitude=-np.arange(longitude_size),
latitude=np.arange(latitude_size)
)
.stack(aux_dim=('longitude', 'latitude'))
.rename(aux_dim='values')
.unstack('values')
)
Hey @AyrtonB, great to bump into you here :slightly_smiling_face:
I must admit, I couldn't remember off the top of my head. So I had a look at the script that we currently use to convert .grib files from the UK Met Office to Zarr. Our script uses ECMWF's cfgrib Python library instead of iris-grib (sorry iris-grib!). Here's the line of code where we load each .grib file.
Likewise, hope you're well :)
That script is perfect, thank you!
For anyone in future this code snippet should achieve what you need
import xarray as xr
import numpy as np
def add_spatial_coords_to_ds(
ds: xr.core.dataset.Dataset,
y_label: str = 'northing',
x_label: str = 'easting',
dy_step: float = 2_000,
dx_step: float = 2_000,
y_max: float = 1223_000,
y_min: float = -185_000,
x_min: float = -239_000,
x_max: float = 857_000
):
return (
ds
.assign_coords(
northing=np.arange(start=y_max, stop=y_min, step=-dy_step, dtype=np.int32),
easting=np.arange(start=x_min, stop=x_max, step=dx_step, dtype=np.int32)
)
.stack(aux_dim=(y_label, x_label))
.rename(aux_dim='values')
.unstack('values')
)
@JackKelly I have been able to load a wholesale file with a bit of faffing.
Firstly, splitting the file up into a grib file per parameter with the following command:
$ grib_copy 201908290000_u1096_ng_umqv_Wholesale1.grib '201908290000_[shortName].grib2'
This creates a grib file per parameter from the wholesale file. Now modify incorrect grib field in section 3 as you had done previously via
$ grib_set -s sourceOfGridDefinition=0 201908290000_r.grib2 201908290000_r_mod.grib2
Iris's Inconsistent Grid Definition error was bounded in the fact that the wholesale file has a reversed set of Y coordinates, and as such Y1 is the larger value out of Y1 and Y2 when passed to the check_range function in _load_convert.py#L858. Modifying that function to no longer assume that v1 is smaller than v2, such as the below
def check_range(v1, v2, n, d):
min_last = min(v1, v2) + (n - 1) * (d - 1)
max_last = min(v1, v2) + (n - 1) * (d + 1)
if not (min_last < max(v1, v2) < max_last):
raise TranslationError('Inconsistent grid definition')
results in a successful opening of the parameter file as a cube, with x and y projection values.
import iris_grib
cubes = iris_grib.load_cubes('201908290000_r_mod.grib2')
print(list(cubes)[0])
# Returns:
#
# relative_humidity / (%) (projection_y_coordinate: 704; projection_x_coordinate: 548)
# Dimension coordinates:
# projection_y_coordinate x -
# projection_x_coordinate - x
# Scalar coordinates:
# forecast_period 0.0 hours
# forecast_reference_time 2019-08-29 00:00:00
# height 1.0 m
# time 2019-08-29 00:00:00
# Attributes:
# GRIB_PARAM GRIB2:d000c001n001
It actually seems like the necessary change to the check_range function in _load_convert.py has already been made last year: https://github.com/SciTools/iris-grib/pull/296/files - but there hasn't been a new release of the package since that was merged in. Installing the 0.19.dev0 build via
$ pip install iris_grib==0.19.dev0
enables a successful load of the parameter grib file as a cube, but also returns a warning: UserWarning: File grid Y definition inconsistent: scanningMode = positive, actual grid point direction is negative. warnings.warn(message).
Hi JackKelly
I installed and tried your script to read the UKV grib files, e.g. 202201010000_u1096_ng_umqv_Wholesale1.grib but I'm getting lots of python errors, e.g.:
Traceback (most recent call last): File "miniconda3/lib/python3.10/site-packages/xarray/core/dataarray.py", line 854, in _getitem_coord var = self._coords[key] KeyError: 'variables'
Any ideas how to fix
Cheers
