diff --git a/docs/examples/grib_to_netcdf.ipynb b/docs/examples/grib_to_netcdf.ipynb index 681bf814f..41c339948 100644 --- a/docs/examples/grib_to_netcdf.ipynb +++ b/docs/examples/grib_to_netcdf.ipynb @@ -23,7 +23,8 @@ "source": [ "import earthkit.data as ekd\n", "ekd.download_example_file(\"tuv_pl.grib\")\n", - "ds = ekd.from_source(\"file\", \"tuv_pl.grib\")" + "# we only select the temperature fields\n", + "ds = ekd.from_source(\"file\", \"tuv_pl.grib\").sel(param=\"t\")" ] }, { @@ -89,25 +90,9 @@ "variables:\n", "\tdouble t(level, latitude, longitude) ;\n", "\t\tt:_FillValue = NaN ;\n", - "\t\tt:param = \"t\" ;\n", "\t\tt:standard_name = \"air_temperature\" ;\n", "\t\tt:long_name = \"Temperature\" ;\n", - "\t\tt:paramId = 130LL ;\n", "\t\tt:units = \"K\" ;\n", - "\tdouble u(level, latitude, longitude) ;\n", - "\t\tu:_FillValue = NaN ;\n", - "\t\tu:param = \"u\" ;\n", - "\t\tu:standard_name = \"eastward_wind\" ;\n", - "\t\tu:long_name = \"U component of wind\" ;\n", - "\t\tu:paramId = 131LL ;\n", - "\t\tu:units = \"m s**-1\" ;\n", - "\tdouble v(level, latitude, longitude) ;\n", - "\t\tv:_FillValue = NaN ;\n", - "\t\tv:param = \"v\" ;\n", - "\t\tv:standard_name = \"northward_wind\" ;\n", - "\t\tv:long_name = \"V component of wind\" ;\n", - "\t\tv:paramId = 132LL ;\n", - "\t\tv:units = \"m s**-1\" ;\n", "\tint64 level(level) ;\n", "\t\tlevel:units = \"hPa\" ;\n", "\t\tlevel:positive = \"down\" ;\n", @@ -126,6 +111,8 @@ "\t\tlongitude:long_name = \"longitude\" ;\n", "\n", "// global attributes:\n", + "\t\t:param = \"t\" ;\n", + "\t\t:paramId = 130LL ;\n", "\t\t:class = \"od\" ;\n", "\t\t:stream = \"oper\" ;\n", "\t\t:levtype = \"pl\" ;\n", @@ -173,7 +160,7 @@ "tags": [] }, "source": [ - "Alternatively, we can use the ``add_earthkit_attrs=False`` option in :py:meth:`~data.readers.grib.index.GribFieldList.to_xarray`. In this case we can call :py:meth:`xarray.Dataset.to_netcdf` directly on the Xarray dataset without using the ``earthkit`` accessor." + "Alternatively, we can use the ``add_earthkit_attrs=False`` option in :py:meth:`~data.readers.grib.index.GribFieldList.to_xarray`. With this the earthkit attributes are not added to the generated dataset and it is safe to call :py:meth:`to_netcdf ` directly on it." ] }, { @@ -200,25 +187,9 @@ "variables:\n", "\tdouble t(level, latitude, longitude) ;\n", "\t\tt:_FillValue = NaN ;\n", - "\t\tt:param = \"t\" ;\n", "\t\tt:standard_name = \"air_temperature\" ;\n", "\t\tt:long_name = \"Temperature\" ;\n", - "\t\tt:paramId = 130LL ;\n", "\t\tt:units = \"K\" ;\n", - "\tdouble u(level, latitude, longitude) ;\n", - "\t\tu:_FillValue = NaN ;\n", - "\t\tu:param = \"u\" ;\n", - "\t\tu:standard_name = \"eastward_wind\" ;\n", - "\t\tu:long_name = \"U component of wind\" ;\n", - "\t\tu:paramId = 131LL ;\n", - "\t\tu:units = \"m s**-1\" ;\n", - "\tdouble v(level, latitude, longitude) ;\n", - "\t\tv:_FillValue = NaN ;\n", - "\t\tv:param = \"v\" ;\n", - "\t\tv:standard_name = \"northward_wind\" ;\n", - "\t\tv:long_name = \"V component of wind\" ;\n", - "\t\tv:paramId = 132LL ;\n", - "\t\tv:units = \"m s**-1\" ;\n", "\tint64 level(level) ;\n", "\t\tlevel:units = \"hPa\" ;\n", "\t\tlevel:positive = \"down\" ;\n", @@ -237,6 +208,8 @@ "\t\tlongitude:long_name = \"longitude\" ;\n", "\n", "// global attributes:\n", + "\t\t:param = \"t\" ;\n", + "\t\t:paramId = 130LL ;\n", "\t\t:class = \"od\" ;\n", "\t\t:stream = \"oper\" ;\n", "\t\t:levtype = \"pl\" ;\n", @@ -257,6 +230,169 @@ "ds.to_xarray(add_earthkit_attrs=False).to_netcdf(\"_tuv_pl_1.nc\")\n", "ncdump(\"_tuv_pl_1.nc\")" ] + }, + { + "cell_type": "markdown", + "id": "23a28542-078b-488c-8b0b-58c7c1add2ba", + "metadata": {}, + "source": [ + "#### Using to_target" + ] + }, + { + "cell_type": "raw", + "id": "7a1e9600-5b49-4fc8-aef8-ac0b51a42f72", + "metadata": { + "editable": true, + "raw_mimetype": "text/restructuredtext", + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "We can call :func:`to_target` on the fieldlist and the Xarray conversion and writing to NetCDF will happen automatically under the hood using the default options. In this case ``add_earthkit_attrs=False`` is always enforced." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "1bdebe02-5f4f-4e3f-b434-ebab25959577", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "netcdf _tuv_pl_2 {\n", + "dimensions:\n", + "\tlevel = 6 ;\n", + "\tlatitude = 7 ;\n", + "\tlongitude = 12 ;\n", + "variables:\n", + "\tdouble longitude(longitude) ;\n", + "\t\tlongitude:units = \"degrees_east\" ;\n", + "\t\tlongitude:standard_name = \"longitude\" ;\n", + "\t\tlongitude:long_name = \"longitude\" ;\n", + "\t\tlongitude:_FillValue = NaN ;\n", + "\tdouble latitude(latitude) ;\n", + "\t\tlatitude:units = \"degrees_north\" ;\n", + "\t\tlatitude:standard_name = \"latitude\" ;\n", + "\t\tlatitude:long_name = \"latitude\" ;\n", + "\t\tlatitude:_FillValue = NaN ;\n", + "\tdouble t(level, latitude, longitude) ;\n", + "\t\tt:standard_name = \"air_temperature\" ;\n", + "\t\tt:long_name = \"Temperature\" ;\n", + "\t\tt:units = \"K\" ;\n", + "\t\tt:_FillValue = NaN ;\n", + "\tint level(level) ;\n", + "\t\tlevel:units = \"hPa\" ;\n", + "\t\tlevel:positive = \"down\" ;\n", + "\t\tlevel:stored_direction = \"decreasing\" ;\n", + "\t\tlevel:standard_name = \"air_pressure\" ;\n", + "\t\tlevel:long_name = \"pressure\" ;\n", + "\n", + "// global attributes:\n", + "\t\t:param = \"t\" ;\n", + "\t\t:paramId = 130 ;\n", + "\t\t:class = \"od\" ;\n", + "\t\t:stream = \"oper\" ;\n", + "\t\t:levtype = \"pl\" ;\n", + "\t\t:type = \"an\" ;\n", + "\t\t:expver = \"0001\" ;\n", + "\t\t:date = 20180801 ;\n", + "\t\t:time = 1200 ;\n", + "\t\t:domain = \"g\" ;\n", + "\t\t:number = 0 ;\n", + "\t\t:Conventions = \"CF-1.8\" ;\n", + "\t\t:institution = \"ECMWF\" ;\n", + "}\n", + "\n" + ] + } + ], + "source": [ + "ds.to_target(\"file\", \"_tuv_pl_2.nc\")\n", + "ncdump(\"_tuv_pl_2.nc\")" + ] + }, + { + "cell_type": "markdown", + "id": "0ebd29d0-cce4-4d3e-bbd7-2fcc1729c563", + "metadata": {}, + "source": [ + "To control the Xarray conversion we can pass options to the earthkit Xarray engine via ``earthkit_to_xarray_kwargs``." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "669cfc69-f64b-490b-83cc-2415ec42ab7a", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "netcdf _tuv_pl_3 {\n", + "dimensions:\n", + "\tlevel = 6 ;\n", + "\tvalues = 84 ;\n", + "variables:\n", + "\tdouble latitude(values) ;\n", + "\t\tlatitude:units = \"degrees_north\" ;\n", + "\t\tlatitude:standard_name = \"latitude\" ;\n", + "\t\tlatitude:long_name = \"latitude\" ;\n", + "\t\tlatitude:_FillValue = NaN ;\n", + "\tdouble longitude(values) ;\n", + "\t\tlongitude:units = \"degrees_east\" ;\n", + "\t\tlongitude:standard_name = \"longitude\" ;\n", + "\t\tlongitude:long_name = \"longitude\" ;\n", + "\t\tlongitude:_FillValue = NaN ;\n", + "\tdouble t(level, values) ;\n", + "\t\tt:standard_name = \"air_temperature\" ;\n", + "\t\tt:long_name = \"Temperature\" ;\n", + "\t\tt:units = \"K\" ;\n", + "\t\tt:coordinates = \"latitude longitude\" ;\n", + "\t\tt:_FillValue = NaN ;\n", + "\tint level(level) ;\n", + "\t\tlevel:units = \"hPa\" ;\n", + "\t\tlevel:positive = \"down\" ;\n", + "\t\tlevel:stored_direction = \"decreasing\" ;\n", + "\t\tlevel:standard_name = \"air_pressure\" ;\n", + "\t\tlevel:long_name = \"pressure\" ;\n", + "\n", + "// global attributes:\n", + "\t\t:param = \"t\" ;\n", + "\t\t:paramId = 130 ;\n", + "\t\t:class = \"od\" ;\n", + "\t\t:stream = \"oper\" ;\n", + "\t\t:levtype = \"pl\" ;\n", + "\t\t:type = \"an\" ;\n", + "\t\t:expver = \"0001\" ;\n", + "\t\t:date = 20180801 ;\n", + "\t\t:time = 1200 ;\n", + "\t\t:domain = \"g\" ;\n", + "\t\t:number = 0 ;\n", + "\t\t:Conventions = \"CF-1.8\" ;\n", + "\t\t:institution = \"ECMWF\" ;\n", + "}\n", + "\n" + ] + } + ], + "source": [ + "ds.to_target(\"file\", \"_tuv_pl_3.nc\", earthkit_to_xarray_kwargs={\"flatten_values\": True})\n", + "ncdump(\"_tuv_pl_3.nc\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "80f49214-43c6-43be-af56-9d55a16f0f8c", + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { diff --git a/docs/guide/xarray/overview.rst b/docs/guide/xarray/overview.rst index 1be2d050c..28b9ab79d 100644 --- a/docs/guide/xarray/overview.rst +++ b/docs/guide/xarray/overview.rst @@ -113,7 +113,12 @@ For further details see the following notebook: Converting GRIB to NetCDF ---------------------------- -To convert GRIB data to NetCDF first we need to convert GRIB to Xarray with :py:meth:`~data.readers.grib.index.GribFieldList.to_xarray` then generate NetCDF from it with :py:meth:`xarray.Dataset.to_netcdf`. Earthkit-data attaches some special attributes to the generated Xarray dataset that cannot be written to NetCDF. In order to make ``to_netcdf()`` work we need to invoke it on the ``earthkit`` accessor and not directly on the Xarray dataset. +To convert GRIB data to NetCDF first we need to convert GRIB to Xarray with :py:meth:`~data.readers.grib.index.GribFieldList.to_xarray` then generate NetCDF from it with :py:meth:`xarray.Dataset.to_netcdf`. We have 3 options to do this: + +Using the earthkit accessor +++++++++++++++++++++++++++++ + +By default, the earthkit Xarray engine attaches some special attributes to the generated Xarray dataset that cannot be written to NetCDF. In order to make ``to_netcdf()`` work we need to invoke it on the ``earthkit`` accessor and not directly on the Xarray dataset. .. code-block:: python @@ -123,7 +128,10 @@ To convert GRIB data to NetCDF first we need to convert GRIB to Xarray with :py: ds_xr = ds_fl.to_xarray() ds_xr.earthkit.to_netcdf("_from_grib.nc") -Alternatively, we can use the ``add_earthkit_attrs=False`` option in :py:meth:`~data.readers.grib.index.GribFieldList.to_xarray`. In this case we can call :py:meth:`xarray.Dataset.to_netcdf` directly on the Xarray dataset without using the ``earthkit`` accessor. +Using the ``add_earthkit_attrs=False`` option +++++++++++++++++++++++++++++++++++++++++++++++++++ + +Alternatively, we can use the ``add_earthkit_attrs=False`` option in :py:meth:`~data.readers.grib.index.GribFieldList.to_xarray`. With this the earthkit attributes are not added to the generated dataset and it is safe to call :py:meth:`to_netcdf ` directly on it. .. code-block:: python @@ -133,6 +141,31 @@ Alternatively, we can use the ``add_earthkit_attrs=False`` option in :py:meth:`~ ds_xr = ds_fl.to_xarray(add_earthkit_attrs=False) ds_xr.to_netcdf("_from_grib.nc") +Using to_target +++++++++++++++++ + +The third option is to use the :func:`to_target` method to convert GRIB directly to NetCDF. This method will generate an Xarray dataset and write it to a NetCDF file in one step. + +.. code-block:: python + + import earthkit.data as ekd + + ds_fl = ekd.from_source("sample", "pl.grib") + ds.fl.to_target("file", "_from_grib.nc") + + +To control the Xarray conversion we can pass options to the earthkit Xarray engine with ``earthkit_to_xarray_kwargs``. In this case ``add_earthkit_attrs=False`` is always enforced. + +.. code-block:: python + + import earthkit.data as ekd + + ds_fl = ekd.from_source("sample", "pl.grib") + ds.fl.to_target( + "file", "_from_grib.nc", earthkit_to_xarray_kwargs={"flatten_values": True} + ) + + For further details see the following notebook: - :ref:`/examples/grib_to_netcdf.ipynb` diff --git a/docs/release_notes/version_0.15_updates.rst b/docs/release_notes/version_0.15_updates.rst index 957a3678c..2d650c393 100644 --- a/docs/release_notes/version_0.15_updates.rst +++ b/docs/release_notes/version_0.15_updates.rst @@ -94,6 +94,14 @@ New features - Added the :ref:`targets-zarr` target (:pr:`716`). See the :ref:`/examples/grib_to_zarr_target.ipynb` notebook example. - Added new config option ``grib-file-serialisation-policy`` to control how GRIB data on disk is pickled. The options are "path" and "memory". The default is "path". Previously, only "memory" was implemented (:pr:`700`). - Added serialisation to GRIB fields (both on disk and in-memory) (:pr:`700`) +- Enabled specifying earthkit Xarray engine options via the ``earthkit_to_xarray_kwargs`` kwarg in :func:`to_target` when converting GRIB to NetCDF. See the :ref:`/examples/grib_to_netcdf.ipynb` notebook example. E.g. + + .. code-block:: python + + ds.to_target( + "netcdf", "pl.nc", earthkit_to_xarray_kwargs={"flatten_values": True} + ) + Fixes diff --git a/src/earthkit/data/encoders/netcdf.py b/src/earthkit/data/encoders/netcdf.py index 5c064a98a..d503373b9 100644 --- a/src/earthkit/data/encoders/netcdf.py +++ b/src/earthkit/data/encoders/netcdf.py @@ -62,7 +62,17 @@ def _encode( missing_value=9999, **kwargs, ): - return NetCDFEncodedData(data.to_xarray()) + _kwargs = kwargs.copy() + if data is not None: + # TODO: find better way to check if the earthkit engine is used + if hasattr(data, "to_xarray_earthkit"): + earthkit_to_xarray_kwargs = _kwargs.pop("earthkit_to_xarray_kwargs", {}) + earthkit_to_xarray_kwargs["add_earthkit_attrs"] = False + _kwargs = earthkit_to_xarray_kwargs + else: + _kwargs = {} + + return NetCDFEncodedData(data.to_xarray(**_kwargs)) def _encode_field(self, field, **kwargs): return self._encode(field, **kwargs) diff --git a/tests/targets/test_target_file.py b/tests/targets/test_target_file.py index 55445d19d..b3fc8a554 100644 --- a/tests/targets/test_target_file.py +++ b/tests/targets/test_target_file.py @@ -230,7 +230,7 @@ def test_target_file_odb(): assert len(df) == 717 -def test_target_file_grib_to_netcdf(): +def test_target_file_grib_to_netcdf_1(): ds = from_source("file", earthkit_examples_file("test.grib")) # vals_ref = ds.values[:, :4] @@ -240,7 +240,28 @@ def test_target_file_grib_to_netcdf(): ds1 = from_source("file", path) assert len(ds1) == len(ds) assert ds1.metadata("param") == ["2t", "msl"] - # assert np.allclose(ds1.values[:, :4], vals_ref) + + ds2 = ds1.to_xarray() + assert "values" not in ds2.sizes + + +def test_target_file_grib_to_netcdf_2(): + ds = from_source("file", earthkit_examples_file("test.grib")) + # vals_ref = ds.values[:, :4] + + with temp_file() as path: + ds.to_target("file", path, encoder="netcdf", earthkit_to_xarray_kwargs={"flatten_values": True}) + + ds1 = from_source("file", path) + ds2 = ds1.to_xarray() + + for name in ["2t", "msl"]: + assert name in ds2.data_vars + + for name in ["latitude", "longitude"]: + assert name in ds2.coords + + assert "values" in ds2.sizes @pytest.mark.skipif(NO_RIOXARRAY, reason="rioxarray not available")