Repository: ecmwf/earthkit-data (python)

Issue:
Do not implicitly set generatingProcessIdentifier when writing GRIB numpy fieldlist
### What happened?

When a GRIB numpy fieldlist is written to disk 

```python
ds.save("my_file.grib")
```

the `generatingProcessIdentifier` ecCodes key is implicitly set to 255 for each message. This should be controlled by an option.

### What are the steps to reproduce the bug?

See above

### Version

4.0

### Platform (OS and architecture)

all

### Relevant log output

_No response_

### Accompanying data

_No response_

### Organisation

ECMWF

Pull request:
Do not implicitly set generatingProcessIdentifier when writing GRIB data
Fixes #274 

Previously, when writing a `NumpyFieldList` the `generatingProcessIdentifier` ecCodes GRIB key was always set to 255 (internally) in the output.

With this PR this key is not set implicitly but must be set explicitly in `Metadata.override()`. E.g.:

```python
# ds is a fieldlist
md = ds[0].metadata().override(shortName="msl", generatingProcessIdentifier=255)
r = FieldList.from_numpy(ds[0].values, md)
r.save("my_data.grib")
```


Patch:
diff --git a/earthkit/data/core/__init__.py b/earthkit/data/core/__init__.py
index 9259032..d9a7ad2 100644
--- a/earthkit/data/core/__init__.py
+++ b/earthkit/data/core/__init__.py
@@ -65,12 +65,12 @@ class Base(metaclass=MetaBase):
 
     # I/O
     @abstractmethod
-    def save(self, path):
+    def save(self, path, **kwargs):
         """Writes data into the specified ``path``."""
         self._not_implemented()
 
     @abstractmethod
-    def write(self, f):
+    def write(self, f, **kwargs):
         """Writes data to the ``f`` file object."""
         self._not_implemented()
 
diff --git a/earthkit/data/core/fieldlist.py b/earthkit/data/core/fieldlist.py
index 03ed322..e55c1d6 100644
--- a/earthkit/data/core/fieldlist.py
+++ b/earthkit/data/core/fieldlist.py
@@ -1164,7 +1164,7 @@ class FieldList(Index):
                 )
         return False
 
-    def save(self, filename, append=False):
+    def save(self, filename, append=False, **kwargs):
         r"""Write all the fields into a file.
 
         Parameters
@@ -1177,9 +1177,9 @@ class FieldList(Index):
         """
         flag = "wb" if not append else "ab"
         with open(filename, flag) as f:
-            self.write(f)
+            self.write(f, **kwargs)
 
-    def write(self, f):
+    def write(self, f, **kwargs):
         r"""Write all the fields to a file object.
 
         Parameters
@@ -1188,7 +1188,7 @@ class FieldList(Index):
             The target file object.
         """
         for s in self:
-            s.write(f)
+            s.write(f, **kwargs)
 
     def to_fieldlist(self, backend, **kwargs):
         r"""Convert to a new :class:`FieldList` based on the ``backend``.
diff --git a/earthkit/data/readers/__init__.py b/earthkit/data/readers/__init__.py
index cae4567..a911ec5 100644
--- a/earthkit/data/readers/__init__.py
+++ b/earthkit/data/readers/__init__.py
@@ -60,12 +60,12 @@ class Reader(Base, os.PathLike, metaclass=ReaderMeta):
     def cache_file(self, *args, **kwargs):
         return self.source.cache_file(*args, **kwargs)
 
-    def save(self, path):
+    def save(self, path, **kwargs):
         mode = "wb" if self.binary else "w"
         with open(path, mode) as f:
-            self.write(f)
+            self.write(f, **kwargs)
 
-    def write(self, f):
+    def write(self, f, **kwargs):
         if not self.appendable:
             assert f.tell() == 0
         mode = "rb" if self.binary else "r"
diff --git a/earthkit/data/readers/directory.py b/earthkit/data/readers/directory.py
index ae71856..984dd4e 100644
--- a/earthkit/data/readers/directory.py
+++ b/earthkit/data/readers/directory.py
@@ -83,10 +83,10 @@ class DirectoryReader(Reader):
             merger=self.merger,
         )
 
-    def save(self, path):
+    def save(self, path, **kwargs):
         shutil.copytree(self.path, path)
 
-    def write(self, f):
+    def write(self, f, **kwargs):
         raise NotImplementedError()
 
 
diff --git a/earthkit/data/readers/grib/codes.py b/earthkit/data/readers/grib/codes.py
index ae6616b..a6fa002 100644
--- a/earthkit/data/readers/grib/codes.py
+++ b/earthkit/data/readers/grib/codes.py
@@ -223,8 +223,6 @@ class GribCodesHandle(CodesHandle):
         try:
             assert self.path is None, "Only cloned handles can have values changed"
             eccodes.codes_set_values(self._handle, values.flatten())
-            # This is writing on the GRIB that something has been modified (255=unknown)
-            eccodes.codes_set_long(self._handle, "generatingProcessIdentifier", 255)
         except Exception as e:
             LOG.error("Error setting values")
             LOG.exception(e)
diff --git a/earthkit/data/readers/grib/output.py b/earthkit/data/readers/grib/output.py
index 11bcec7..bf58236 100644
--- a/earthkit/data/readers/grib/output.py
+++ b/earthkit/data/readers/grib/output.py
@@ -123,6 +123,9 @@ class GribOutput:
             k: v for k, v in sorted(metadata.items(), key=lambda x: order(x[0]))
         }
 
+        if "generatingProcessIdentifier" not in metadata:
+            metadata["generatingProcessIdentifier"] = 255
+
         LOG.debug("GribOutput.metadata %s", metadata)
 
         for k, v in metadata.items():
diff --git a/earthkit/data/sources/file.py b/earthkit/data/sources/file.py
index cc118e4..32b197e 100644
--- a/earthkit/data/sources/file.py
+++ b/earthkit/data/sources/file.py
@@ -114,11 +114,11 @@ class FileSource(Source, os.PathLike, metaclass=FileSourceMeta):
     def values(self):
         return self._reader.values
 
-    def save(self, path):
-        return self._reader.save(path)
+    def save(self, path, **kwargs):
+        return self._reader.save(path, **kwargs)
 
-    def write(self, f):
-        return self._reader.write(f)
+    def write(self, f, **kwargs):
+        return self._reader.write(f, **kwargs)
 
     def scaled(self, *args, **kwargs):
         return self._reader.scaled(*args, **kwargs)
diff --git a/earthkit/data/sources/memory.py b/earthkit/data/sources/memory.py
index cce56df..365a145 100644
--- a/earthkit/data/sources/memory.py
+++ b/earthkit/data/sources/memory.py
@@ -51,11 +51,11 @@ class MemoryBaseSource(Source):
     def to_fieldlist(self, *args, **kwargs):
         return self._reader.to_fieldlist(*args, **kwargs)
 
-    def save(self, path):
-        return self._reader.save(path)
+    def save(self, path, **kwargs):
+        return self._reader.save(path, **kwargs)
 
-    def write(self, f):
-        return self._reader.write(f)
+    def write(self, f, **kwargs):
+        return self._reader.write(f, **kwargs)
 
     def scaled(self, *args, **kwargs):
         return self._reader.scaled(*args, **kwargs)
diff --git a/earthkit/data/sources/multi.py b/earthkit/data/sources/multi.py
index a059971..b3e1d99 100644
--- a/earthkit/data/sources/multi.py
+++ b/earthkit/data/sources/multi.py
@@ -84,10 +84,10 @@ class MultiSource(Source):
         string = ",".join(repr(s) for s in self.sources)
         return f"{self.__class__.__name__}({string})"
 
-    def save(self, path):
+    def save(self, path, **kwargs):
         with open(path, "wb") as f:
             for s in self.sources:
-                s.write(f)
+                s.write(f, **kwargs)
 
     def graph(self, depth=0):
         print(" " * depth, self.__class__.__name__, self.merger)
diff --git a/earthkit/data/sources/numpy_list.py b/earthkit/data/sources/numpy_list.py
index 03dd1ce..055caad 100644
--- a/earthkit/data/sources/numpy_list.py
+++ b/earthkit/data/sources/numpy_list.py
@@ -20,7 +20,7 @@ LOG = logging.getLogger(__name__)
 
 
 class NumpyField(Field):
-    r"""Represents a field consisting of an ndarray and metadata object.
+    r"""Represent a field consisting of an ndarray and metadata object.
 
     Parameters
     ----------
@@ -46,10 +46,19 @@ class NumpyField(Field):
     def __repr__(self):
         return f"{self.__class__.__name__}()"
 
-    def write(self, f):
+    def write(self, f, **kwargs):
+        r"""Write the field to a file object.
+
+        Parameters
+        ----------
+        f: file object
+            The target file object.
+        **kwargs: dict, optional
+            Other keyword arguments passed to :meth:`data.writers.grib.GribWriter.write`.
+        """
         from earthkit.data.writers import write
 
-        write(f, self.values, self._metadata, check_nans=True)
+        write(f, self.values, self._metadata, **kwargs)
 
 
 class NumpyFieldListCore(PandasMixIn, XarrayMixIn, FieldList):
@@ -158,7 +167,7 @@ class ListMerger:
 
 
 class NumpyFieldList(NumpyFieldListCore):
-    r"""Represents a list of :obj:`NumpyField <data.sources.numpy_list.NumpyField>`\ s.
+    r"""Represent a list of :obj:`NumpyField <data.sources.numpy_list.NumpyField>`\ s.
 
     The preferred way to create a NumpyFieldList is to use either the
     static :obj:`from_numpy` method or the :obj:`to_fieldlist` method.
diff --git a/earthkit/data/writers/grib.py b/earthkit/data/writers/grib.py
index 6e0fa39..3209110 100644
--- a/earthkit/data/writers/grib.py
+++ b/earthkit/data/writers/grib.py
@@ -14,6 +14,19 @@ class GribWriter(Writer):
     DATA_FORMAT = "grib"
 
     def write(self, f, values, metadata, check_nans=True):
+        r"""Write a GRIB field to a file object.
+
+        Parameters
+        ----------
+        f: file object
+            The target file object.
+        values: ndarray
+            Values of the GRIB field/message.
+        values: :class:`GribMetadata`
+            Metadata of the GRIB field/message.
+        check_nans: bool
+            Replace nans in ``values`` with GRIB missing values when writing to``f``.
+        """
         handle = metadata._handle
         if check_nans:
             import numpy as np
diff --git a/setup.cfg b/setup.cfg
index cb9b572..3f9fe0c 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -43,6 +43,7 @@ install_requires =
     tqdm
     xarray>=0.19.0
     earthkit-meteo>=0.0.1
+    aws-requests-auth
 include_package_data = True
 
 [options.packages.find]


Test patch:
diff --git a/tests/grib/test_grib_output.py b/tests/grib/test_grib_output.py
index 56b9617..8034cd9 100644
--- a/tests/grib/test_grib_output.py
+++ b/tests/grib/test_grib_output.py
@@ -55,6 +55,7 @@ def test_grib_output_latlon():
         assert ds[0].metadata("param") == "2t"
         assert ds[0].metadata("levtype") == "sfc"
         assert ds[0].metadata("edition") == 2
+        assert ds[0].metadata("generatingProcessIdentifier") == 255
 
         assert np.allclose(ds[0].to_numpy(), data, rtol=EPSILON, atol=EPSILON)
 
@@ -79,6 +80,7 @@ def test_grib_output_o96():
         assert ds[0].metadata("param") == "2t"
         assert ds[0].metadata("levtype") == "sfc"
         assert ds[0].metadata("edition") == 2
+        assert ds[0].metadata("generatingProcessIdentifier") == 255
 
         assert np.allclose(ds[0].to_numpy(), data, rtol=EPSILON, atol=EPSILON)
 
@@ -103,6 +105,7 @@ def test_grib_output_o160():
         assert ds[0].metadata("edition") == 2
         assert ds[0].metadata("levtype") == "sfc"
         assert ds[0].metadata("param") == "2t"
+        assert ds[0].metadata("generatingProcessIdentifier") == 255
 
         assert np.allclose(ds[0].to_numpy(), data, rtol=EPSILON, atol=EPSILON)
 
@@ -130,6 +133,7 @@ def test_grib_output_mars_labeling():
         assert ds[0].metadata("levtype") == "sfc"
         assert ds[0].metadata("param") == "msl"
         assert ds[0].metadata("type") == "fc"
+        assert ds[0].metadata("generatingProcessIdentifier") == 255
 
         assert np.allclose(ds[0].to_numpy(), data, rtol=EPSILON, atol=EPSILON)
 
@@ -158,6 +162,7 @@ def test_grib_output_pl(levtype):
         assert ds[0].metadata("level") == 850
         assert ds[0].metadata("levtype") == "pl"
         assert ds[0].metadata("param") == "t"
+        assert ds[0].metadata("generatingProcessIdentifier") == 255
 
         assert np.allclose(ds[0].to_numpy(), data, rtol=EPSILON, atol=EPSILON)
 
@@ -185,6 +190,7 @@ def test_grib_output_tp():
         assert ds[0].metadata("levtype") == "sfc"
         assert ds[0].metadata("edition") == 1
         assert ds[0].metadata("step") == 48
+        assert ds[0].metadata("generatingProcessIdentifier") == 255
 
         assert np.allclose(ds[0].to_numpy(), data, rtol=EPSILON, atol=EPSILON)
 
diff --git a/tests/numpy_fs/test_numpy_fs_write.py b/tests/numpy_fs/test_numpy_fs_write.py
index acb4d7c..5dbefe9 100644
--- a/tests/numpy_fs/test_numpy_fs_write.py
+++ b/tests/numpy_fs/test_numpy_fs_write.py
@@ -13,6 +13,7 @@ import logging
 import os
 
 import numpy as np
+import pytest
 
 from earthkit.data import from_source
 from earthkit.data.core.fieldlist import FieldList
@@ -22,7 +23,8 @@ from earthkit.data.testing import earthkit_examples_file
 LOG = logging.getLogger(__name__)
 
 
-def test_numpy_fs_grib_write_missing():
+@pytest.mark.parametrize("_kwargs", [{}, {"check_nans": True}])
+def test_numpy_fs_grib_write_missing(_kwargs):
     ds = from_source("file", earthkit_examples_file("test.grib"))
 
     assert ds[0].metadata("shortName") == "2t"
@@ -42,14 +44,40 @@ def test_numpy_fs_grib_write_missing():
     assert np.isnan(r[0].values[0])
     assert not np.isnan(r[0].values[1])
 
-    # save to disk
-    tmp = temp_file()
-    r.save(tmp.path)
-    assert os.path.exists(tmp.path)
-    r_tmp = from_source("file", tmp.path)
-    v_tmp = r_tmp[0].values
-    assert np.isnan(v_tmp[0])
-    assert not np.isnan(v_tmp[1])
+    with temp_file() as tmp:
+        r.save(tmp, **_kwargs)
+        assert os.path.exists(tmp)
+        r_tmp = from_source("file", tmp)
+        v_tmp = r_tmp[0].values
+        assert np.isnan(v_tmp[0])
+        assert not np.isnan(v_tmp[1])
+
+
+def test_numpy_fs_grib_write_check_nans_bad():
+    ds = from_source("file", earthkit_examples_file("test.grib"))
+
+    assert ds[0].metadata("shortName") == "2t"
+
+    v = ds[0].values
+    v1 = v + 1
+    assert not np.isnan(v1[0])
+    assert not np.isnan(v1[1])
+    v1[0] = np.nan
+    assert np.isnan(v1[0])
+    assert not np.isnan(v1[1])
+
+    md = ds[0].metadata()
+    md1 = md.override(shortName="msl")
+    r = FieldList.from_numpy(v1, md1)
+
+    assert np.isnan(r[0].values[0])
+    assert not np.isnan(r[0].values[1])
+
+    with temp_file() as tmp:
+        from eccodes import EncodingError
+
+        with pytest.raises(EncodingError):
+            r.save(tmp, check_nans=False)
 
 
 def test_numpy_fs_grib_write_append():
@@ -85,6 +113,34 @@ def test_numpy_fs_grib_write_append():
     assert r_tmp.metadata("shortName") == ["msl", "2d"]
 
 
+def test_numpy_fs_grib_write_generating_proc_id():
+    ds = from_source("file", earthkit_examples_file("test.grib"))
+
+    assert ds[0].metadata("shortName") == "2t"
+
+    v = ds[0].values
+    v1 = v + 1
+    v2 = v + 2
+
+    md = ds[0].metadata()
+    md1 = md.override(shortName="msl", generatingProcessIdentifier=255)
+    md2 = md.override(shortName="2d")
+
+    r1 = FieldList.from_numpy([v1, v2], [md1, md2])
+
+    # save to disk: using generatingProcessIdentifier=255 (default)
+    with temp_file() as tmp:
+        r1.save(tmp)
+        assert os.path.exists(tmp)
+        r_tmp = from_source("file", tmp)
+        assert len(r_tmp) == 2
+        assert r_tmp.metadata("shortName") == ["msl", "2d"]
+        assert r_tmp.metadata("generatingProcessIdentifier") == [
+            255,
+            150,
+        ]
+
+
 if __name__ == "__main__":
     from earthkit.data.testing import main