search

ecmwf / cfgrib

A Python interface to map GRIB files to the NetCDF Common Data Model following the CF Convention using ecCodes
Apache License 2.0
407 stars 78 forks source link

No latitudes/longitudes provided by ecCodes for gridType #28

Open alexamici opened 6 years ago

alexamici commented 6 years ago

cfgrib relies on the ECMWF ecCodes C-library for encoding and decoding the GRIB format, including everything related to the coordinate systems. GRIB files encoded in a gridType not supported by the installed ecCodes version will print the warning: No latitudes/longitudes provided by ecCodes for gridType.

The GRIB file will be opened but the geographic portion of the data will be represented by as single dimension without coordinate named values.

For example:

>>> cfgrib.open_dataset('../ST4.2018080204.01h')
No latitudes/longitudes provided by ecCodes for gridType = 'polar_stereographic'
<xarray.Dataset>
Dimensions:     (values: 987601)
Coordinates:
    time        datetime64[ns] ...
    step        timedelta64[ns] ...
    surface     int64 ...
    valid_time  datetime64[ns] ...
Dimensions without coordinates: values
Data variables:
    tp          (values) float32 ...
Attributes:
    GRIB_edition:            1
    GRIB_centre:             kwbc
    GRIB_centreDescription:  US National Weather Service - NCEP 
    GRIB_subCentre:          4
    history:                 GRIB to CDM+CF via cfgrib-0.9.../ecCodes-2.8...

The list of known supported gridTypes is:

see http://xarray.pydata.org/en/stable/data-structures.html#coordinates for details on xarray naming conventions.

Initially noted in #27.

alexamici commented 6 years ago

The reason the issue is labeled wontfix is that cfgrib is not designed to handle coordinate systems and we can just wait for ecCodes to add support for new gridTypes.

marcowurth commented 4 years ago

When opening unstructured_grid files from DWD's ICON model with parallel=False and engine='cfgrib' there is a warning message "ecCodes provides no latitudes/longitudes for gridType='unstructured_grid'" for every opened grib file. Since opening, reading and e.g. saving as netcdf4 works fine and one can add unstructured coordinates for the 'values' dimension afterwards to the Dataset, could you suppress the warning message?

alexamici commented 4 years ago

@marcowurth you can add backend_kwargs={'errors': 'ignore'} to your xr.open_dataset call to have all errore silenced, e.g.:

>>> xr.open_dataset('dwd.grib', engine='cfgrib', backend_kwargs={'errors': 'ignore'})
matteodefelice commented 4 years ago

I have this issue when opening EFAS files, so is any workarounds available? The NetCDFs are much bigger and I wanted to stick to GRIB format. I got ecCodes provides no latitudes/longitudes for gridType='lambert_azimuthal_equal_area'

shahramn commented 4 years ago

Dear Matteo, ecCodes does support the Lambert Azimuthal Equal Area projection but only on a perfect sphere i.e. not on the oblate spheroid

alexamici commented 4 years ago

@matteodefelice I added 'lambert_azimuthal_equal_area' and 'albers' to the supported gridTypes.

winash12 commented 3 years ago

When opening unstructured_grid files from DWD's ICON model with parallel=False and engine='cfgrib' there is a warning message "ecCodes provides no latitudes/longitudes for gridType='unstructured_grid'" for every opened grib file. Since opening, reading and e.g. saving as netcdf4 works fine and one can add unstructured coordinates for the 'values' dimension afterwards to the Dataset, could you suppress the warning message?

I am also trying to read an ICON grib file with eccodes and cfgrib I get this message - is there a way to get the latitudes and longitudes ?

` Dimensions: (depthBelowLand: 5, depthBelowLandLayer: 4, generalVertical: 71, generalVerticalLayer: 70, values: 76376) Coordinates: time datetime64[ns] ... step timedelta64[ns] ... surface int64 ... valid_time datetime64[ns] ...

shahramn commented 3 years ago

For unstructured grids, the latitudes and longitudes are stored in separate GRIB messages as their field values. So for example let's imagine we have 4 messages: The field values (data section) of message 1 will store the latitudes The field values of message 2 will store the longitudes The field values of messages 3 and 4 could be temperature or pressure To get the lat/lon for the temperature field, you read the 1st and 2nd messages. Of course these messages must all have the same number of values i.e., the number of grid points

shahramn commented 3 years ago

Re: Unstructured grids Here is a worked example:

Let's say you download some DWD ICON model data. For example http://icon-downloads.zmaw.de/dwd_grids.xml#grid2

First look for the message whose shortName is "CLON" (center longitude), then a message whose shortName is "CLAT" (center latitude). The lat/lon params may be different e.g. discipline=0,parameterCategory=191,parameterNumber=1 (Geographical latitude)

Store their values in separate files. Then get one of the other fields e.g. HSURF (Geometric Height of the earths surface above sea level) and store its values. % echo 'if(shortName is "CLAT") { print "[values!1]"; }' | grib_filter - dwd.grib2 > lats % echo 'if(shortName is "CLON") { print "[values!1]"; }' | grib_filter - dwd.grib2 > lons % echo 'if(shortName is "HSURF"){ print "[values!1]"; }' | grib_filter - dwd.grib2 > vals

Now put these columns together: % paste lons lats vals > dwd.geopts

winash12 commented 3 years ago

@shahramn Thanks for the detailed response. I am sorry I did not understand a word of your second message. I get a subset of ICON data for running COSMO directly from DWD. I uploaded a sample here -https://gofile.io/d/r0I6EN and can you tell me whether I can print the lats and lons ?

shahramn commented 3 years ago

@winash12 Apologies. I am using the ecCodes command line tools to decode these ICON files.

The sample file you mentioned (igfff00000000) doesn't seem to have any lat/lon info. None of its parameters are using the parameter category 191 which gives you "Geographical latitude" and "Geographical longitude". In the ICON files I had seen before, there were two GRIB messages providing those values.

So for this one I'm afraid you would have to contact the data provider

kdl0013 commented 1 year ago

I'm working with EMC GEFSv12 data and I'm having the same issues, but I found a weird workaround (very slow and I'm not a huge fan of it). I can download the grib and then open up each file in python, manipulate it (select variables), then save it as a netcdf. After I can then remap data and re-open the files in python.

But if I take the original grib file, perform cdo operators, and save it as a netcdf file, then when I try to open it in python I get ECCODE errors. This may be useful or it just a weird workaround and it also may just be specific to my datasets. Here was how I opened and extracted my variables in python

grib_o = xr.open_dataset(file_o,engine='cfgrib',filter_by_keys={'typeOfLevel': 'isobaricInhPa', 'dataType': 'cf'})

meteoDaniel commented 1 year ago

Dear @alexamici @iainrussell , I running the latest eccodes and cfgrib version atm and having trouble to read ICON data transformed into regular latitude longitude grid following this guide: CDO remap

In my test case I use this grib file . 

import subprocess
from pathlib import Path
import xarray 

in_file = Path('icon_global_icosahedral_single-level_2023091400_000_T_2M.grib2')
    intermediate_file = Path(
        str(in_file).replace("icon_global_", "icon_global_to_transform_")
    )
in_file.replace(intermediate_file)
subprocess.run(
        [
            "cdo",
            "-f",
            "grb2",
            f"remap,{GRID_FIXTURE_BASE_PATH}/target_grid_world_0125.txt,"
            f"{GRID_FIXTURE_BASE_PATH}/weights_icogl2world_0125.nc",
            str(intermediate_file),
            str(in_file),
        ]
)
  subprocess.run(
        [
            "cdo",
            "-f",
            "grb2",
            "remap,target_grid_world_0125.txt,"
            "weights_icogl2world_0125.nc",
            str(intermediate_file),
            str(in_file),
        ]
    )
    intermediate_file.unlink()
import xarray
xarray.open_dataset(in_file, engine='cfgrib')

This yields to: 

ECCODES ERROR   :  Lat/Lon Geoiterator: First and last latitudes are inconsistent with scanning order: lat1=-90, lat2=90 jScansPositively=0
ECCODES ERROR   :  Geoiterator factory: Error instantiating iterator latlon (Grid description is wrong or inconsistent)
ECCODES ERROR   :  latitudes: Unable to create iterator
---------------------------------------------------------------------------
WrongGridError                            Traceback (most recent call last)
Cell In[13], line 1
----> 1 xarray.open_dataset(in_file, engine='cfgrib')

File /usr/local/lib/python3.10/site-packages/xarray/backends/api.py:570, in open_dataset(filename_or_obj, engine, chunks, cache, decode_cf, mask_and_scale, decode_times, decode_timedelta, use_cftime, concat_characters, decode_coords, drop_variables, inline_array, chunked_array_type, from_array_kwargs, backend_kwargs, **kwargs)
    558 decoders = _resolve_decoders_kwargs(
    559     decode_cf,
    560     open_backend_dataset_parameters=backend.open_dataset_parameters,
   (...)
    566     decode_coords=decode_coords,
    567 )
    569 overwrite_encoded_chunks = kwargs.pop("overwrite_encoded_chunks", None)
--> 570 backend_ds = backend.open_dataset(
    571     filename_or_obj,
    572     drop_variables=drop_variables,
    573     **decoders,
    574     **kwargs,
    575 )
    576 ds = _dataset_from_backend_dataset(
    577     backend_ds,
    578     filename_or_obj,
   (...)
    588     **kwargs,
    589 )
    590 return ds

File /usr/local/lib/python3.10/site-packages/cfgrib/xarray_plugin.py:108, in CfGribBackend.open_dataset(self, filename_or_obj, mask_and_scale, decode_times, concat_characters, decode_coords, drop_variables, use_cftime, decode_timedelta, lock, indexpath, filter_by_keys, read_keys, encode_cf, squeeze, time_dims, errors, extra_coords)
     87 def open_dataset(
     88     self,
     89     filename_or_obj: T.Union[str, abc.MappingFieldset[T.Any, abc.Field]],
   (...)
    106     extra_coords: T.Dict[str, str] = {},
    107 ) -> xr.Dataset:
--> 108     store = CfGribDataStore(
    109         filename_or_obj,
    110         indexpath=indexpath,
    111         filter_by_keys=filter_by_keys,
    112         read_keys=read_keys,
    113         encode_cf=encode_cf,
    114         squeeze=squeeze,
    115         time_dims=time_dims,
    116         lock=lock,
    117         errors=errors,
    118         extra_coords=extra_coords,
    119     )
    120     with xr.core.utils.close_on_error(store):
    121         vars, attrs = store.load()  # type: ignore

File /usr/local/lib/python3.10/site-packages/cfgrib/xarray_plugin.py:40, in CfGribDataStore.__init__(self, filename, lock, **backend_kwargs)
     38 else:
     39     opener = dataset.open_fieldset
---> 40 self.ds = opener(filename, **backend_kwargs)

File /usr/local/lib/python3.10/site-packages/cfgrib/dataset.py:784, in open_file(path, grib_errors, indexpath, filter_by_keys, read_keys, time_dims, extra_coords, **kwargs)
    781 index_keys = compute_index_keys(time_dims, extra_coords)
    782 index = open_fileindex(stream, indexpath, index_keys, filter_by_keys=filter_by_keys)
--> 784 return open_from_index(index, read_keys, time_dims, extra_coords, **kwargs)

File /usr/local/lib/python3.10/site-packages/cfgrib/dataset.py:726, in open_from_index(index, read_keys, time_dims, extra_coords, **kwargs)
    719 def open_from_index(
    720     index: abc.Index[T.Any, abc.Field],
    721     read_keys: T.Sequence[str] = (),
   (...)
    724     **kwargs: T.Any,
    725 ) -> Dataset:
--> 726     dimensions, variables, attributes, encoding = build_dataset_components(
    727         index, read_keys=read_keys, time_dims=time_dims, extra_coords=extra_coords, **kwargs
    728     )
    729     return Dataset(dimensions, variables, attributes, encoding)

File /usr/local/lib/python3.10/site-packages/cfgrib/dataset.py:653, in build_dataset_components(index, errors, encode_cf, squeeze, log, read_keys, time_dims, extra_coords)
    651 var_index = index.subindex(paramId=param_id)
    652 try:
--> 653     dims, data_var, coord_vars = build_variable_components(
    654         var_index,
    655         encode_cf,
    656         filter_by_keys,
    657         errors=errors,
    658         squeeze=squeeze,
    659         read_keys=read_keys,
    660         time_dims=time_dims,
    661         extra_coords=extra_coords,
    662     )
    663 except DatasetBuildError as ex:
    664     # NOTE: When a variable has more than one value for an attribute we need to raise all
    665     #   the values in the file, not just the ones associated with that variable. See #54.
    666     key = ex.args[1]

File /usr/local/lib/python3.10/site-packages/cfgrib/dataset.py:529, in build_variable_components(index, encode_cf, filter_by_keys, log, errors, squeeze, read_keys, time_dims, extra_coords)
    526 header_dimensions = tuple(d for d, c in coord_vars.items() if not squeeze or c.data.size > 1)
    527 header_shape = tuple(coord_vars[d].data.size for d in header_dimensions)
--> 529 geo_dims, geo_shape, geo_coord_vars = build_geography_coordinates(first, encode_cf, errors)
    530 dimensions = header_dimensions + geo_dims
    531 shape = header_shape + geo_shape

File /usr/local/lib/python3.10/site-packages/cfgrib/dataset.py:392, in build_geography_coordinates(first, encode_cf, errors, log)
    390 geo_dims = ("latitude", "longitude")  # type: T.Tuple[str, ...]
    391 geo_shape = (first["Ny"], first["Nx"])  # type: T.Tuple[int, ...]
--> 392 latitudes = np.array(first["distinctLatitudes"], ndmin=1)
    393 geo_coord_vars["latitude"] = Variable(
    394     dimensions=("latitude",), data=latitudes, attributes=COORD_ATTRS["latitude"].copy()
    395 )
    397 if latitudes[0] > latitudes[-1]:

File /usr/local/lib/python3.10/site-packages/cfgrib/messages.py:246, in ComputedKeysAdapter.__getitem__(self, item)
    244     return getter(self)
    245 else:
--> 246     return self.context[item]

File /usr/local/lib/python3.10/site-packages/cfgrib/messages.py:168, in Message.__getitem__(self, item)
    166     raise ValueError("key type not supported %r" % key_type_text)
    167 key_type = KEY_TYPES[key_type_text]
--> 168 return self.message_get(key, key_type=key_type)

File /usr/local/lib/python3.10/site-packages/cfgrib/messages.py:130, in Message.message_get(self, item, key_type, default)
    128 """Get value of a given key as its native or specified type."""
    129 try:
--> 130     if eccodes.codes_get_size(self.codes_id, item) > 1:
    131         values = eccodes.codes_get_array(self.codes_id, item, key_type)
    132     else:

File /usr/local/lib/python3.10/site-packages/gribapi/gribapi.py:595, in grib_get_size(msgid, key)
    593 size_p = ffi.new("size_t*")
    594 err = lib.grib_get_size(h, key.encode(ENC), size_p)
--> 595 GRIB_CHECK(err)
    596 return size_p[0]

File /usr/local/lib/python3.10/site-packages/gribapi/gribapi.py:226, in GRIB_CHECK(errid)
    218 """
    219 Utility function checking the ecCodes error code and raising
    220 an error if that was set.
   (...)
    223 @exception CodesInternalError
    224 """
    225 if errid:
--> 226     errors.raise_grib_error(errid)

File /usr/local/lib/python3.10/site-packages/gribapi/errors.py:381, in raise_grib_error(errid)
    377 def raise_grib_error(errid):
    378     """
    379     Raise the GribInternalError corresponding to ``errid``.
    380     """
--> 381     raise ERROR_MAP[errid](errid)

WrongGridError: Grid description is wrong or inconsistent

I tried to set another gridType with grib_set but this did not worked out. I hope you have some hints to get this working again.

Using ECCODES 2.27. and cfgrib 0.9.5.5 worked for me. (my old setup)