CF-compliant storage for lon/lat case

satpy/tests/utils.py

@@ -325,9 +325,9 @@ def make_fake_scene(content_dict, daskify=False, area=True,
         area (bool or BaseDefinition): Can be ``True``, ``False``, or an
             instance of ``pyresample.geometry.BaseDefinition`` such as
             ``AreaDefinition`` or ``SwathDefinition``.  If ``True``, which is
-            the default, automatically generate areas.  If ``False``, values
-            will not have assigned areas.  If an instance of
-            ``pyresample.geometry.BaseDefinition``, those instances will be
+            the default, automatically generate areas with the name "test-area".
+            If ``False``, values will not have assigned areas.  If an instance
+            of ``pyresample.geometry.BaseDefinition``, those instances will be
             used for all generated fake datasets.  Warning: Passing an area as
             a string (``area="germ"``) is not supported.
         common_attrs (Mapping): optional, additional attributes that will
@@ -341,31 +341,45 @@ def make_fake_scene(content_dict, daskify=False, area=True,
     sc = Scene()
     for (did, arr) in content_dict.items():
         extra_attrs = common_attrs.copy()
-        if isinstance(area, BaseDefinition):
-            extra_attrs["area"] = area
-        elif area:
-            extra_attrs["area"] = create_area_def(
-                    "test-area",
-                    {"proj": "eqc", "lat_ts": 0, "lat_0": 0, "lon_0": 0,
-                     "x_0": 0, "y_0": 0, "ellps": "sphere", "units": "m",
-                     "no_defs": None, "type": "crs"},
-                    units="m",
-                    shape=arr.shape,
-                    resolution=1000,
-                    center=(0, 0))
-        if isinstance(arr, DataArray):
-            sc[did] = arr.copy()  # don't change attributes of input
-            sc[did].attrs.update(extra_attrs)
-        else:
-            if daskify:
-                arr = da.from_array(arr)
-            sc[did] = DataArray(
-                    arr,
-                    dims=("y", "x"),
-                    attrs=extra_attrs)
+        if area:
+            extra_attrs["area"] = _get_fake_scene_area(arr, area)
+        sc[did] = _get_did_for_fake_scene(area, arr, extra_attrs, daskify)
     return sc
 
 
+def _get_fake_scene_area(arr, area):
+    """Get area for fake scene.  Helper for make_fake_scene."""
+    if isinstance(area, BaseDefinition):
+        return area
+    return create_area_def(
+        "test-area",
+        {"proj": "eqc", "lat_ts": 0, "lat_0": 0, "lon_0": 0,
+         "x_0": 0, "y_0": 0, "ellps": "sphere", "units": "m",
+         "no_defs": None, "type": "crs"},
+        units="m",
+        shape=arr.shape,
+        resolution=1000,
+        center=(0, 0))
+
+
+def _get_did_for_fake_scene(area, arr, extra_attrs, daskify):
+    """Add instance to fake scene.  Helper for make_fake_scene."""
+    from satpy.resample import add_crs_xy_coords
+    if isinstance(arr, DataArray):
+        new = arr.copy()  # don't change attributes of input
+        new.attrs.update(extra_attrs)
+    else:
+        if daskify:
+            arr = da.from_array(arr)
+        new = DataArray(
+                arr,
+                dims=("y", "x"),
+                attrs=extra_attrs)
+    if area:
+        new = add_crs_xy_coords(new, extra_attrs["area"])
+    return new
+
+
 def assert_attrs_equal(attrs, attrs_exp, tolerance=0):
     """Test that attributes are equal.
 

satpy/tests/writer_tests/test_cf.py

@@ -17,6 +17,7 @@
 # satpy.  If not, see <http://www.gnu.org/licenses/>.
 """Tests for the CF writer."""
 
+import logging
 import os
 import tempfile
 import unittest
@@ -1132,6 +1133,79 @@ def test_global_attr_history_and_Conventions(self):
                 self.assertIn('Created by pytroll/satpy on', f.attrs['history'])
 
 
+def test_lonlat_storage(tmp_path):
+    """Test correct storage for area with lon/lat units."""
+    import xarray as xr
+    from pyresample import create_area_def
+
+    from ..utils import make_fake_scene
+    scn = make_fake_scene(
+            {"ketolysis": np.arange(25).reshape(5, 5)},
+            daskify=True,
+            area=create_area_def("mavas", 4326, shape=(5, 5),
+                                 center=(0, 0), resolution=(1, 1)))
+
+    filename = os.fspath(tmp_path / "test.nc")
+    scn.save_datasets(filename=filename, writer="cf", include_lonlats=False)
+    with xr.open_dataset(filename) as ds:
+        assert ds["ketolysis"].attrs["grid_mapping"] == "mavas"
+        assert ds["mavas"].attrs["grid_mapping_name"] == "latitude_longitude"
+        assert ds["x"].attrs["units"] == "degrees_east"
+        assert ds["y"].attrs["units"] == "degrees_north"
+        assert ds["mavas"].attrs["longitude_of_prime_meridian"] == 0.0
+        np.testing.assert_allclose(ds["mavas"].attrs["semi_major_axis"], 6378137.0)
+        np.testing.assert_allclose(ds["mavas"].attrs["inverse_flattening"], 298.257223563)
+
+
+def test_da2cf_lonlat():
+    """Test correct da2cf encoding for area with lon/lat units."""
+    import xarray as xr
+    from pyresample import create_area_def
+
+    from satpy.resample import add_crs_xy_coords
+    from satpy.writers.cf_writer import CFWriter
+
+    area = create_area_def("mavas", 4326, shape=(5, 5),
+                           center=(0, 0), resolution=(1, 1))
+    da = xr.DataArray(
+        np.arange(25).reshape(5, 5),
+        dims=("y", "x"),
+        attrs={"area": area})
+    da = add_crs_xy_coords(da, area)
+    new_da = CFWriter.da2cf(da)
+    assert new_da["x"].attrs["units"] == "degrees_east"
+    assert new_da["y"].attrs["units"] == "degrees_north"
+
+
+def test_is_projected(caplog):
+    """Tests for private _is_projected function."""
+    import xarray as xr
+
+    from satpy.writers.cf_writer import CFWriter
+
+    # test case with units but no area
+    da = xr.DataArray(
+        np.arange(25).reshape(5, 5),
+        dims=("y", "x"),
+        coords={"x": xr.DataArray(np.arange(5), dims=("x",), attrs={"units": "m"}),
+                "y": xr.DataArray(np.arange(5), dims=("y",), attrs={"units": "m"})})
+    assert CFWriter._is_projected(da)
+
+    da = xr.DataArray(
+        np.arange(25).reshape(5, 5),
+        dims=("y", "x"),
+        coords={"x": xr.DataArray(np.arange(5), dims=("x",), attrs={"units": "degrees_east"}),
+                "y": xr.DataArray(np.arange(5), dims=("y",), attrs={"units": "degrees_north"})})
+    assert not CFWriter._is_projected(da)
+
+    da = xr.DataArray(
+        np.arange(25).reshape(5, 5),
+        dims=("y", "x"))
+    with caplog.at_level(logging.WARNING):
+        assert CFWriter._is_projected(da)
+    assert "Failed to tell if data are projected." in caplog.text
+
+
 class TestCFWriterData(unittest.TestCase):
     """Test case for CF writer where data arrays are needed."""
 

satpy/writers/cf_writer.py

@@ -34,6 +34,11 @@
 * You can select the netCDF backend using the ``engine`` keyword argument. If `None` if follows
   :meth:`~xarray.Dataset.to_netcdf` engine choices with a preference for 'netcdf4'.
 * For datasets with area definition you can exclude lat/lon coordinates by setting ``include_lonlats=False``.
+  If the area has a projected CRS, units are assumed to be in metre.  If the
+  area has a geographic CRS, units are assumed to be in degrees.  The writer
+  does not verify that the CRS is supported by the CF conventions.  One
+  commonly used projected CRS not supported by the CF conventions is the
+  equirectangular projection, such as EPSG 4087.
 * By default non-dimensional coordinates (such as scanline timestamps) are prefixed with the corresponding
   dataset name. This is because they are likely to be different for each dataset. If a non-dimensional
   coordinate is identical for all datasets, the prefix can be removed by setting ``pretty=True``.
@@ -640,6 +645,9 @@ def da2cf(dataarray, epoch=EPOCH, flatten_attrs=False, exclude_attrs=None, compr
 
         CFWriter._remove_satpy_attributes(new_data)
 
+        new_data = CFWriter._encode_time(new_data, epoch)
+        new_data = CFWriter._encode_coords(new_data)
+
         # Remove area as well as user-defined attributes
         for key in ['area'] + exclude_attrs:
             new_data.attrs.pop(key, None)
@@ -655,9 +663,6 @@ def da2cf(dataarray, epoch=EPOCH, flatten_attrs=False, exclude_attrs=None, compr
         if compression is not None:
             new_data.encoding.update(compression)
 
-        new_data = CFWriter._encode_time(new_data, epoch)
-        new_data = CFWriter._encode_coords(new_data)
-
         if 'long_name' not in new_data.attrs and 'standard_name' not in new_data.attrs:
             new_data.attrs['long_name'] = new_data.name
         if 'prerequisites' in new_data.attrs:
@@ -685,14 +690,73 @@ def _cleanup_attrs(new_data):
 
     @staticmethod
     def _encode_coords(new_data):
+        """Encode coordinates."""
+        if not new_data.coords.keys() & {"x", "y", "crs"}:
+            # there are no coordinates
+            return new_data
+        is_projected = CFWriter._is_projected(new_data)
+        if is_projected:
+            new_data = CFWriter._encode_xy_coords_projected(new_data)
+        else:
+            new_data = CFWriter._encode_xy_coords_geographic(new_data)
+        if 'crs' in new_data.coords:
+            new_data = new_data.drop_vars('crs')
+        return new_data
+
+    @staticmethod
+    def _is_projected(new_data):
+        """Guess whether data are projected or not."""
+        crs = CFWriter._try_to_get_crs(new_data)
+        if crs:
+            return crs.is_projected
+        units = CFWriter._try_get_units_from_coords(new_data)
+        if units:
+            if units.endswith("m"):
+                return True
+            if units.startswith("degrees"):
+                return False
+        logger.warning("Failed to tell if data are projected. Assuming yes.")
+        return True
+
+    @staticmethod
+    def _try_to_get_crs(new_data):
+        """Try to get a CRS from attributes."""
+        if "area" in new_data.attrs:
+            if isinstance(new_data.attrs["area"], AreaDefinition):
+                return new_data.attrs["area"].crs
+            # at least one test case passes an area of type str
+            logger.warning(
+                f"Could not tell CRS from area of type {type(new_data.attrs['area']).__name__:s}. "
+                "Assuming projected CRS.")
+        if "crs" in new_data.coords:
+            return new_data.coords["crs"].item()
+
+    @staticmethod
+    def _try_get_units_from_coords(new_data):
+        for c in "xy":
+            if "units" in new_data.coords[c].attrs:
+                return new_data.coords[c].attrs["units"]
+
+    @staticmethod
+    def _encode_xy_coords_projected(new_data):
+        """Encode coordinates, assuming projected CRS."""
         if 'x' in new_data.coords:
             new_data['x'].attrs['standard_name'] = 'projection_x_coordinate'
             new_data['x'].attrs['units'] = 'm'
         if 'y' in new_data.coords:
             new_data['y'].attrs['standard_name'] = 'projection_y_coordinate'
             new_data['y'].attrs['units'] = 'm'
-        if 'crs' in new_data.coords:
-            new_data = new_data.drop_vars('crs')
+        return new_data
+
+    @staticmethod
+    def _encode_xy_coords_geographic(new_data):
+        """Encode coordinates, assuming geographic CRS."""
+        if 'x' in new_data.coords:
+            new_data['x'].attrs['standard_name'] = 'longitude'
+            new_data['x'].attrs['units'] = 'degrees_east'
+        if 'y' in new_data.coords:
+            new_data['y'].attrs['standard_name'] = 'latitude'
+            new_data['y'].attrs['units'] = 'degrees_north'
         return new_data
 
     @staticmethod