Fix generating output for surface GRIB data with reduced Gaussian gri…

…ds (#449)

* Fix writing reduced Gaussian surface GRIB data

docs/release_notes/version_0.10_updates.rst

@@ -2,6 +2,15 @@ Version 0.10 Updates
 /////////////////////////
 
 
+Version 0.10.2
+===============
+
+Fixes
+++++++
+
+- Fixed issue when generating output with reduced Gaussian surface GRIB data did not work
+
+
 Version 0.10.1
 ===============
 

src/earthkit/data/readers/grib/output.py

@@ -343,9 +343,10 @@ def _gg_field(self, values, metadata):
             else:
                 levtype = "sfc"
 
-        if octahedral:
+        if octahedral or levtype == "sfc":
             return f"reduced_gg_{levtype}_grib{edition}"
         else:
+
             return f"reduced_gg_{levtype}_{N}_grib{edition}"
 
 

tests/grib/test_grib_output.py

@@ -79,7 +79,7 @@ def test_grib_output_latlon():
     sys.version_info < (3, 10),
     reason="ignore_cleanup_errors requires Python 3.10 or later",
 )
-def test_grib_output_o96():
+def test_grib_output_o96_sfc():
     data = np.random.random((40320,))
 
     with tempfile.TemporaryDirectory(ignore_cleanup_errors=True) as tmp:
@@ -91,11 +91,19 @@ def test_grib_output_o96():
 
         ds = earthkit.data.from_source("file", path)
 
-        assert ds[0].metadata("date") == 20010101
-        assert ds[0].metadata("param") == "2t"
-        assert ds[0].metadata("levtype") == "sfc"
-        assert ds[0].metadata("edition") == 2
-        assert ds[0].metadata("generatingProcessIdentifier") == 255
+        ref = {
+            "date": 20010101,
+            "param": "2t",
+            "levtype": "sfc",
+            "edition": 2,
+            "generatingProcessIdentifier": 255,
+            "gridType": "reduced_gg",
+            "N": 96,
+            "isOctahedral": 1,
+        }
+
+        for k, v in ref.items():
+            assert ds[0].metadata(k) == v, f"{k}: {ds[0].metadata(k)}!={v}"
 
         assert np.allclose(ds[0].to_numpy(), data, rtol=EPSILON, atol=EPSILON)
 
@@ -104,7 +112,7 @@ def test_grib_output_o96():
     sys.version_info < (3, 10),
     reason="ignore_cleanup_errors requires Python 3.10 or later",
 )
-def test_grib_output_o160():
+def test_grib_output_o160_sfc():
     data = np.random.random((108160,))
 
     with tempfile.TemporaryDirectory(ignore_cleanup_errors=True) as tmp:
@@ -116,11 +124,52 @@ def test_grib_output_o160():
 
         ds = earthkit.data.from_source("file", path)
 
-        assert ds[0].metadata("date") == 20010101
-        assert ds[0].metadata("edition") == 2
-        assert ds[0].metadata("levtype") == "sfc"
-        assert ds[0].metadata("param") == "2t"
-        assert ds[0].metadata("generatingProcessIdentifier") == 255
+        ref = {
+            "date": 20010101,
+            "param": "2t",
+            "levtype": "sfc",
+            "edition": 2,
+            "generatingProcessIdentifier": 255,
+            "gridType": "reduced_gg",
+            "N": 160,
+            "isOctahedral": 1,
+        }
+
+        for k, v in ref.items():
+            assert ds[0].metadata(k) == v, f"{k}: {ds[0].metadata(k)}!={v}"
+
+        assert np.allclose(ds[0].to_numpy(), data, rtol=EPSILON, atol=EPSILON)
+
+
+@pytest.mark.skipif(
+    sys.version_info < (3, 10),
+    reason="ignore_cleanup_errors requires Python 3.10 or later",
+)
+def test_grib_output_n96_sfc():
+    data = np.random.random(50662)
+
+    with tempfile.TemporaryDirectory(ignore_cleanup_errors=True) as tmp:
+        path = os.path.join(tmp, "a.grib")
+
+        f = earthkit.data.new_grib_output(path, date=20010101)
+        f.write(data, param="2t")
+        f.close()
+
+        ds = earthkit.data.from_source("file", path)
+
+        ref = {
+            "date": 20010101,
+            "param": "2t",
+            "levtype": "sfc",
+            "edition": 2,
+            "generatingProcessIdentifier": 255,
+            "gridType": "reduced_gg",
+            "N": 96,
+            "isOctahedral": 0,
+        }
+
+        for k, v in ref.items():
+            assert ds[0].metadata(k) == v, f"{k}: {ds[0].metadata(k)}!={v}"
 
         assert np.allclose(ds[0].to_numpy(), data, rtol=EPSILON, atol=EPSILON)
 
@@ -149,6 +198,9 @@ def test_grib_output_mars_labeling():
         assert ds[0].metadata("param") == "msl"
         assert ds[0].metadata("type") == "fc"
         assert ds[0].metadata("generatingProcessIdentifier") == 255
+        assert ds[0].metadata("gridType") == "reduced_gg"
+        assert ds[0].metadata("N") == 96
+        assert ds[0].metadata("isOctahedral") == 1
 
         assert np.allclose(ds[0].to_numpy(), data, rtol=EPSILON, atol=EPSILON)
 
@@ -158,7 +210,7 @@ def test_grib_output_mars_labeling():
     reason="ignore_cleanup_errors requires Python 3.10 or later",
 )
 @pytest.mark.parametrize("levtype", [{}, {"levtype": "pl"}])
-def test_grib_output_pl(levtype):
+def test_grib_output_o96_pl(levtype):
     data = np.random.random((40320,))
 
     with tempfile.TemporaryDirectory(ignore_cleanup_errors=True) as tmp:
@@ -178,6 +230,9 @@ def test_grib_output_pl(levtype):
         assert ds[0].metadata("levtype") == "pl"
         assert ds[0].metadata("param") == "t"
         assert ds[0].metadata("generatingProcessIdentifier") == 255
+        assert ds[0].metadata("gridType") == "reduced_gg"
+        assert ds[0].metadata("N") == 96
+        assert ds[0].metadata("isOctahedral") == 1
 
         assert np.allclose(ds[0].to_numpy(), data, rtol=EPSILON, atol=EPSILON)
 
@@ -198,14 +253,16 @@ def test_grib_output_tp():
         f.close()
 
         ds = earthkit.data.from_source("file", path)
-        print(ds[0])
 
         assert ds[0].metadata("date") == 20010101
         assert ds[0].metadata("param") == "tp"
         assert ds[0].metadata("levtype") == "sfc"
         assert ds[0].metadata("edition") == 1
         assert ds[0].metadata("step", astype=int) == 48
         assert ds[0].metadata("generatingProcessIdentifier") == 255
+        assert ds[0].metadata("gridType") == "regular_ll"
+        assert ds[0].metadata("Ni") == 360
+        assert ds[0].metadata("Nj") == 181
 
         assert np.allclose(ds[0].to_numpy(), data, rtol=EPSILON, atol=EPSILON)