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.

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

How to add date and step info when writing grib file with to_grib #156

Closed
blaylockbk opened this issue Aug 19, 2020 · 10 comments
Closed

How to add date and step info when writing grib file with to_grib #156

blaylockbk opened this issue Aug 19, 2020 · 10 comments

Comments

@blaylockbk
Copy link

blaylockbk commented Aug 19, 2020

I followed the example in the README for writing a GRIB file from a simple Dataset

import xarray as xr
import cfgrib

ds2 = xr.Dataset({'skin_temperature': (('latitude', 'longitude'), np.zeros((5, 6)) + 300.)})
ds2.coords['latitude'] = np.linspace(90., -90., 5)
ds2.coords['longitude'] = np.linspace(0., 360., 6, endpoint=False)
ds2.skin_temperature.attrs['GRIB_shortName'] = 'skt'

cfgrib.to_grib(ds2, 'out.grib2', grib_keys={'edition': 2})

That was successful.

Now, I would like to specify the datetime to the grib2 file, but when I add the coordinate time, valid_time, or step I get an error.

import xarray as xr
import cfgrib
import pandas as pd

ds2 = xr.Dataset({'skin_temperature': (('latitude', 'longitude'), np.zeros((5, 6)) + 300.)})
ds2.coords['latitude'] = np.linspace(90., -90., 5)
ds2.coords['longitude'] = np.linspace(0., 360., 6, endpoint=False)
ds2.coords['time'] = pd.to_datetime('2017-10-01')            ### <----- Added this line
ds2.skin_temperature.attrs['GRIB_shortName'] = 'skt'
cfgrib.to_grib(ds2, 'out.grib2', grib_keys={'edition': 2})
KeyError: 1506816000000000000

I realize this is an Alpha capability, but can you please provide any additional information on passing time data when creating a simple grib file. Thanks!

Note: I have cfgrib 0.9.8.3

@alexamici
Copy link
Contributor

alexamici commented Sep 15, 2020

In general you should use arrays of np.datetime64 not pandas time object. But in this case I think you are hitting an incompatibility between newer versions of pandas and xarray, see: pydata/xarray#4283

@honnorat
Copy link

Any news regarding this issue ? I'm stuck too !

Since datetime-like coordinates are always converted to datetime64[ns] in xarray, I can see currently no workaround for this problem.
Except of course pinning pandas version to 1.0.x.

@alexamici
Copy link
Contributor

alexamici commented Nov 30, 2020

This has been fixed in current xarray master, but it is not release yet.

The best fix ATM is still to downgrade pandas to version 1.0.5

@honnorat
Copy link

Here is my test case:

import cfgrib
import numpy as np
import xarray as xr

ds = xr.DataArray(
    np.zeros((1, 6, 6)), name="foo",
    dims=("time", "latitude", "longitude"),
    coords={
        "time": np.asarray(["2000-01-01"], dtype="datetime64[s]"),
        "latitude": np.linspace(0., 1., 6),
        "longitude": np.linspace(0., 1., 6)
    },
).to_dataset()
cfgrib.to_grib(ds.squeeze(), "test.grib2")

Works like a charm with pandas=1.0.5. With pandas>=1.1.0, I get:

KeyError                                  Traceback (most recent call last)
<ipython-input-3-1ce018954770> in <module>
----> 1 cfgrib.to_grib(ds.squeeze(), "test.grib2")

~/usr/conda/envs/dp_11/lib/python3.8/site-packages/cfgrib/xarray_to_grib.py in to_grib(*args, **kwargs)
    267 
    268 def to_grib(*args, **kwargs):
--> 269     return canonical_dataset_to_grib(*args, **kwargs)

~/usr/conda/envs/dp_11/lib/python3.8/site-packages/cfgrib/xarray_to_grib.py in canonical_dataset_to_grib(dataset, path, mode, no_warn, grib_keys, **kwargs)
    263     with open(path, mode=mode) as file:
    264         for data_var in dataset.data_vars.values():
--> 265             canonical_dataarray_to_grib(data_var, file, grib_keys=real_grib_keys, **kwargs)
    266 
    267 

~/usr/conda/envs/dp_11/lib/python3.8/site-packages/cfgrib/xarray_to_grib.py in canonical_dataarray_to_grib(data_var, file, grib_keys, default_grib_keys, **kwargs)
    222     for items in itertools.product(*header_coords_values):
    223         select = {n: v for n, v in zip(coords_names, items)}
--> 224         field_values = data_var.sel(**select).values.flat[:]
    225 
    226         # Missing values handling

~/usr/conda/envs/dp_11/lib/python3.8/site-packages/xarray/core/dataarray.py in sel(self, indexers, method, tolerance, drop, **indexers_kwargs)
   1141 
   1142         """
-> 1143         ds = self._to_temp_dataset().sel(
   1144             indexers=indexers,
   1145             drop=drop,

~/usr/conda/envs/dp_11/lib/python3.8/site-packages/xarray/core/dataset.py in sel(self, indexers, method, tolerance, drop, **indexers_kwargs)
   2103         """
   2104         indexers = either_dict_or_kwargs(indexers, indexers_kwargs, "sel")
-> 2105         pos_indexers, new_indexes = remap_label_indexers(
   2106             self, indexers=indexers, method=method, tolerance=tolerance
   2107         )

~/usr/conda/envs/dp_11/lib/python3.8/site-packages/xarray/core/coordinates.py in remap_label_indexers(obj, indexers, method, tolerance, **indexers_kwargs)
    395     }
    396 
--> 397     pos_indexers, new_indexes = indexing.remap_label_indexers(
    398         obj, v_indexers, method=method, tolerance=tolerance
    399     )

~/usr/conda/envs/dp_11/lib/python3.8/site-packages/xarray/core/indexing.py in remap_label_indexers(data_obj, indexers, method, tolerance)
    273             coords_dtype = data_obj.coords[dim].dtype
    274             label = maybe_cast_to_coords_dtype(label, coords_dtype)
--> 275             idxr, new_idx = convert_label_indexer(index, label, dim, method, tolerance)
    276             pos_indexers[dim] = idxr
    277             if new_idx is not None:

~/usr/conda/envs/dp_11/lib/python3.8/site-packages/xarray/core/indexing.py in convert_label_indexer(index, label, index_name, method, tolerance)
    194                 indexer = index.get_loc(label_value)
    195             else:
--> 196                 indexer = index.get_loc(label_value, method=method, tolerance=tolerance)
    197         elif label.dtype.kind == "b":
    198             indexer = label

~/usr/conda/envs/dp_11/lib/python3.8/site-packages/pandas/core/indexes/datetimes.py in get_loc(self, key, method, tolerance)
    620         else:
    621             # unrecognized type
--> 622             raise KeyError(key)
    623 
    624         try:

KeyError: 946684800000000000

@honnorat
Copy link

honnorat commented Dec 2, 2020

This has been fixed in current xarray master, but it is not release yet.

With xarray 0.16.2, I still have the same error.

The best fix ATM is still to downgrade pandas to version 1.0.5

Yes indeed.

@alexamici
Copy link
Contributor

This has been fixed in current xarray master, but it is not release yet.

With xarray 0.16.2, I still have the same error.

I just tested myself and the fix is not there yet. I'm surprised :(

@honnorat
Copy link

Still the same problem with pandas 1.3:

$ conda list | grep -E 'python |pandas|cfgrib|xarray'
cfgrib                    0.9.9.0            pyhd8ed1ab_1    conda-forge
pandas                    1.3.0            py39hde0f152_0    conda-forge
python                    3.9.6           h49503c6_1_cpython    conda-forge
xarray                    0.18.2             pyhd8ed1ab_0    conda-forge
$ python ./test_cfgrib.py
/home/honnorat/usr/conda/envs/dp_13/lib/python3.9/site-packages/cfgrib/xarray_to_grib.py:255: FutureWarning: GRIB write support is experimental, DO NOT RELY ON IT!
  warnings.warn("GRIB write support is experimental, DO NOT RELY ON IT!", FutureWarning)
Traceback (most recent call last):
  File "/home/honnorat/test/./test_cfgrib.py", line 14, in <module>
    to_grib(ds.squeeze(), "test.grib2")
  File "/home/honnorat/usr/conda/envs/dp_13/lib/python3.9/site-packages/cfgrib/xarray_to_grib.py", line 266, in canonical_dataset_to_grib
    canonical_dataarray_to_grib(data_var, file, grib_keys=real_grib_keys, **kwargs)
  File "/home/honnorat/usr/conda/envs/dp_13/lib/python3.9/site-packages/cfgrib/xarray_to_grib.py", line 225, in canonical_dataarray_to_grib
    field_values = data_var.sel(**select).values.flat[:]
  File "/home/honnorat/usr/conda/envs/dp_13/lib/python3.9/site-packages/xarray/core/dataarray.py", line 1271, in sel
    ds = self._to_temp_dataset().sel(
  File "/home/honnorat/usr/conda/envs/dp_13/lib/python3.9/site-packages/xarray/core/dataset.py", line 2365, in sel
    pos_indexers, new_indexes = remap_label_indexers(
  File "/home/honnorat/usr/conda/envs/dp_13/lib/python3.9/site-packages/xarray/core/coordinates.py", line 421, in remap_label_indexers
    pos_indexers, new_indexes = indexing.remap_label_indexers(
  File "/home/honnorat/usr/conda/envs/dp_13/lib/python3.9/site-packages/xarray/core/indexing.py", line 274, in remap_label_indexers
    idxr, new_idx = convert_label_indexer(index, label, dim, method, tolerance)
  File "/home/honnorat/usr/conda/envs/dp_13/lib/python3.9/site-packages/xarray/core/indexing.py", line 191, in convert_label_indexer
    indexer = index.get_loc(label_value, method=method, tolerance=tolerance)
  File "/home/honnorat/usr/conda/envs/dp_13/lib/python3.9/site-packages/pandas/core/indexes/datetimes.py", line 699, in get_loc
    raise KeyError(key)
KeyError: 946684800000000000

The source of the problem seems to be:

header_coords_values = [data_var.coords[name].values.tolist() for name in coords_names]

The call to tolist() forces numpy to convert the array of np.datetime64 into a list of scalars, which causes the final KeyError. If I remove the call to .tolist() in line 222 above, the problem for my test case above (#156 (comment)) is solved.

However, I have no idea what further impacts this fix would have on other use cases.

@honnorat
Copy link

Well, reading grib2 files (NCEP GFS) with xr.open_dataset and writing them back to grib using my fix above works fine now !

@iainrussell
Copy link
Member

Thank you for your suggestion @honnorat . I have tried it with my local installation, but the tests fail, so we cannot use your change as is. I tried a few further modifications, but with no success so far.

@alexamici
Copy link
Contributor

alexamici commented Jan 28, 2022

This issue should be fixed by #272.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

4 participants