From b4f9ae5a6db59cea04bd973a818e0ffcc13e157e Mon Sep 17 00:00:00 2001
From: Martin Raspaud <martin.raspaud@smhi.se>
Date: Thu, 30 Apr 2020 14:52:58 +0200
Subject: [PATCH 1/2] Fix cf writing of 3d arrays

---
 satpy/tests/writer_tests/test_cf.py | 26 ++++++++++++++++++++++++++
 satpy/writers/cf_writer.py          |  3 ++-
 2 files changed, 28 insertions(+), 1 deletion(-)

diff --git a/satpy/tests/writer_tests/test_cf.py b/satpy/tests/writer_tests/test_cf.py
index 001d18714f..d6b3147d23 100644
--- a/satpy/tests/writer_tests/test_cf.py
+++ b/satpy/tests/writer_tests/test_cf.py
@@ -907,6 +907,7 @@ def test_area2lonlat(self):
         """Test the conversion from areas to lon/lat."""
         import pyresample.geometry
         import xarray as xr
+        import dask.array as da
         from satpy.writers.cf_writer import area2lonlat
 
         area = pyresample.geometry.AreaDefinition(
@@ -934,6 +935,31 @@ def test_area2lonlat(self):
         self.assertDictContainsSubset({'name': 'longitude', 'standard_name': 'longitude', 'units': 'degrees_east'},
                                       lon.attrs)
 
+        area = pyresample.geometry.AreaDefinition(
+            'seviri',
+            'Native SEVIRI grid',
+            'geos',
+            "+a=6378169.0 +h=35785831.0 +b=6356583.8 +lon_0=0 +proj=geos",
+            10, 10,
+            [-5570248.686685662, -5567248.28340708, 5567248.28340708, 5570248.686685662]
+        )
+        lons_ref, lats_ref = area.get_lonlats()
+        dataarray = xr.DataArray(data=da.from_array(np.arange(3*10*10).reshape(3, 10, 10), chunks=(1, 5, 5)),
+                                 dims=('bands', 'y', 'x'), attrs={'area': area})
+        res = area2lonlat(dataarray)
+
+        # original should be unmodified
+        self.assertNotIn('longitude', dataarray.coords)
+        self.assertEqual(set(res.coords), {'longitude', 'latitude'})
+        lat = res['latitude']
+        lon = res['longitude']
+        self.assertTrue(np.all(lat.data == lats_ref))
+        self.assertTrue(np.all(lon.data == lons_ref))
+        self.assertDictContainsSubset({'name': 'latitude', 'standard_name': 'latitude', 'units': 'degrees_north'},
+                                      lat.attrs)
+        self.assertDictContainsSubset({'name': 'longitude', 'standard_name': 'longitude', 'units': 'degrees_east'},
+                                      lon.attrs)
+
 
 def suite():
     """Test suite for this writer's tests."""
diff --git a/satpy/writers/cf_writer.py b/satpy/writers/cf_writer.py
index 10ae9d8c69..30140db0e6 100644
--- a/satpy/writers/cf_writer.py
+++ b/satpy/writers/cf_writer.py
@@ -167,7 +167,8 @@ def area2lonlat(dataarray):
     """Convert an area to longitudes and latitudes."""
     dataarray = dataarray.copy()
     area = dataarray.attrs['area']
-    chunks = getattr(dataarray.data, 'chunks', None)
+    ignore_dims = {dim: 0 for dim in dataarray.dims if dim not in ['x', 'y']}
+    chunks = getattr(dataarray.isel(**ignore_dims), 'chunks', None)
     lons, lats = area.get_lonlats(chunks=chunks)
     lons = xr.DataArray(lons, dims=['y', 'x'],
                         attrs={'name': "longitude",

From b5bd072f5e5757df533a797b5b1ad9fe9b2a0c45 Mon Sep 17 00:00:00 2001
From: Martin Raspaud <martin.raspaud@smhi.se>
Date: Mon, 4 May 2020 10:58:55 +0200
Subject: [PATCH 2/2] Clean up cf's area2lonlat

---
 satpy/writers/cf_writer.py | 22 ++++++++++------------
 1 file changed, 10 insertions(+), 12 deletions(-)

diff --git a/satpy/writers/cf_writer.py b/satpy/writers/cf_writer.py
index 30140db0e6..80e9bc7491 100644
--- a/satpy/writers/cf_writer.py
+++ b/satpy/writers/cf_writer.py
@@ -170,18 +170,16 @@ def area2lonlat(dataarray):
     ignore_dims = {dim: 0 for dim in dataarray.dims if dim not in ['x', 'y']}
     chunks = getattr(dataarray.isel(**ignore_dims), 'chunks', None)
     lons, lats = area.get_lonlats(chunks=chunks)
-    lons = xr.DataArray(lons, dims=['y', 'x'],
-                        attrs={'name': "longitude",
-                               'standard_name': "longitude",
-                               'units': 'degrees_east'},
-                        name='longitude')
-    lats = xr.DataArray(lats, dims=['y', 'x'],
-                        attrs={'name': "latitude",
-                               'standard_name': "latitude",
-                               'units': 'degrees_north'},
-                        name='latitude')
-    dataarray['longitude'] = lons
-    dataarray['latitude'] = lats
+    dataarray['longitude'] = xr.DataArray(lons, dims=['y', 'x'],
+                                          attrs={'name': "longitude",
+                                                 'standard_name': "longitude",
+                                                 'units': 'degrees_east'},
+                                          name='longitude')
+    dataarray['latitude'] = xr.DataArray(lats, dims=['y', 'x'],
+                                         attrs={'name': "latitude",
+                                                'standard_name': "latitude",
+                                                'units': 'degrees_north'},
+                                         name='latitude')
     return dataarray