Merge pull request #446 from pytroll/bugfix-nc_nwcsaf_coordinates

Bugfix nc_nwcsaf_msg reader projection units handling

satpy/readers/nc_nwcsaf.py

@@ -210,18 +210,10 @@ def get_area_def(self, dsid):
         if dsid.name.endswith('_pal'):
             raise NotImplementedError
 
-        try:
-            proj_str = self.nc.attrs['gdal_projection'] + ' +units=km'
-        except TypeError:
-            proj_str = self.nc.attrs['gdal_projection'].decode() + ' +units=km'
+        proj_str, area_extent = self._get_projection()
 
         nlines, ncols = self.nc[dsid.name].shape
 
-        area_extent = (float(self.nc.attrs['gdal_xgeo_up_left']) / 1000,
-                       float(self.nc.attrs['gdal_ygeo_low_right']) / 1000,
-                       float(self.nc.attrs['gdal_xgeo_low_right']) / 1000,
-                       float(self.nc.attrs['gdal_ygeo_up_left']) / 1000)
-
         area = get_area_def('some_area_name',
                             "On-the-fly area",
                             'geosmsg',
@@ -271,6 +263,32 @@ def end_time(self):
             return datetime.strptime(self.nc.attrs['time_coverage_end'],
                                      '%Y%m%dT%H%M%S%fZ')
 
+    def _get_projection(self):
+        """Get projection from the NetCDF4 attributes"""
+        try:
+            proj_str = self.nc.attrs['gdal_projection']
+        except TypeError:
+            proj_str = self.nc.attrs['gdal_projection'].decode()
+
+        # Check the a/b/h units
+        radius_a = proj_str.split('+a=')[-1].split()[0]
+        if float(radius_a) > 10e3:
+            units = 'm'
+            scale = 1.0
+        else:
+            units = 'km'
+            scale = 1e3
+
+        if 'units' not in proj_str:
+            proj_str = proj_str + ' +units=' + units
+
+        area_extent = (float(self.nc.attrs['gdal_xgeo_up_left']) / scale,
+                       float(self.nc.attrs['gdal_ygeo_low_right']) / scale,
+                       float(self.nc.attrs['gdal_xgeo_low_right']) / scale,
+                       float(self.nc.attrs['gdal_ygeo_up_left']) / scale)
+
+        return proj_str, area_extent
+
 
 def remove_empties(variable):
     """Remove empty objects from the *variable*'s attrs."""

satpy/tests/reader_tests/__init__.py

@@ -35,7 +35,8 @@
                                       test_grib, test_hrit_goes, test_ahi_hsd,
                                       test_iasi_l2, test_generic_image,
                                       test_scmi, test_hrit_jma, test_nc_goes,
-                                      test_nc_slstr, test_nc_olci, test_viirs_flood)
+                                      test_nc_slstr, test_nc_olci,
+                                      test_viirs_flood, test_nc_nwcsaf)
 
 
 if sys.version_info < (2, 7):
@@ -75,5 +76,6 @@ def suite():
     mysuite.addTests(test_hrit_jma.suite())
     mysuite.addTests(test_nc_goes.suite())
     mysuite.addTests(test_nc_slstr.suite())
+    mysuite.addTests(test_nc_nwcsaf.suite())
 
     return mysuite

satpy/tests/reader_tests/test_nc_nwcsaf.py

@@ -0,0 +1,95 @@
+#!/usr/bin/python
+# Copyright (c) 2018.
+#
+
+# Author(s):
+#   Panu Lahtinen <[email protected]>
+
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+# General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see <http://www.gnu.org/licenses/>.
+#
+
+"""Unittests for NWC SAF reader.
+"""
+
+import sys
+if sys.version_info < (2, 7):
+    import unittest2 as unittest
+else:
+    import unittest
+
+try:
+    from unittest import mock
+except ImportError:
+    import mock
+
+PROJ_KM = {'gdal_projection': '+proj=geos +a=6378.137000 +b=6356.752300 +lon_0=0.000000 +h=35785.863000',
+           'gdal_xgeo_up_left': -5569500.0,
+           'gdal_ygeo_up_left': 5437500.0,
+           'gdal_xgeo_low_right': 5566500.0,
+           'gdal_ygeo_low_right': 2653500.0}
+PROJ = {'gdal_projection': '+proj=geos +a=6378137.000 +b=6356752.300 +lon_0=0.000000 +h=35785863.000',
+        'gdal_xgeo_up_left': -5569500.0,
+        'gdal_ygeo_up_left': 5437500.0,
+        'gdal_xgeo_low_right': 5566500.0,
+        'gdal_ygeo_low_right': 2653500.0}
+
+
+class TestNcNWCSAF(unittest.TestCase):
+
+    @mock.patch('satpy.readers.nc_nwcsaf.unzip_file')
+    @mock.patch('satpy.readers.nc_nwcsaf.xr')
+    def setUp(self, xr_, unzip):
+        from satpy.readers.nc_nwcsaf import NcNWCSAF
+        xr_.return_value = mock.Mock(attrs={})
+        unzip.return_value = ''
+        self.scn = NcNWCSAF('filename', {}, {})
+
+    def tearDown(self):
+        pass
+
+    def test_get_projection(self):
+        # a, b and h in kilometers
+        self.scn.nc.attrs = PROJ_KM
+        proj_str, area_extent = self.scn._get_projection()
+        self.assertTrue('+units=km' in proj_str)
+        self.assertAlmostEqual(area_extent[0],
+                               PROJ_KM['gdal_xgeo_up_left'] / 1000.)
+        self.assertAlmostEqual(area_extent[1],
+                               PROJ_KM['gdal_ygeo_low_right'] / 1000.)
+        self.assertAlmostEqual(area_extent[2],
+                               PROJ_KM['gdal_xgeo_low_right'] / 1000.)
+        self.assertAlmostEqual(area_extent[3],
+                               PROJ_KM['gdal_ygeo_up_left'] / 1000.)
+
+        # a, b and h in meters
+        self.scn.nc.attrs = PROJ
+        proj_str, area_extent = self.scn._get_projection()
+        self.assertTrue('+units=m' in proj_str)
+        self.assertAlmostEqual(area_extent[0], PROJ['gdal_xgeo_up_left'])
+        self.assertAlmostEqual(area_extent[1], PROJ['gdal_ygeo_low_right'])
+        self.assertAlmostEqual(area_extent[2], PROJ['gdal_xgeo_low_right'])
+        self.assertAlmostEqual(area_extent[3], PROJ['gdal_ygeo_up_left'])
+
+
+def suite():
+    """The test suite for test_writers."""
+    loader = unittest.TestLoader()
+    my_suite = unittest.TestSuite()
+    my_suite.addTest(loader.loadTestsFromTestCase(TestNcNWCSAF))
+
+    return my_suite
+
+
+if __name__ == '__main__':
+    unittest.main()