Merge pull request #339 from hazbottles/master

fixed meteosat native georeferencing

satpy/readers/native_msg.py

@@ -28,7 +28,8 @@
 References:
     MSG Level 1.5 Native Format File Definition
     https://www.eumetsat.int/website/wcm/idc/idcplg?IdcService=GET_FILE&dDocName=PDF_FG15_MSG-NATIVE-FORMAT-15&RevisionSelectionMethod=LatestReleased&Rendition=Web
-
+    MSG Level 1.5 Image Data Format Description
+    https://www.eumetsat.int/website/wcm/idc/idcplg?IdcService=GET_FILE&dDocName=PDF_TEN_05105_MSG_IMG_DATA&RevisionSelectionMethod=LatestReleased&Rendition=Web
 """
 
 import os
@@ -272,6 +273,17 @@ def get_area_extent(self, dsid):
 
         if dsid.name != 'HRV':
 
+            # following calculations assume grid origin is south-east corner
+            # section 7.2.4 of MSG Level 1.5 Image Data Format Description
+            origins = {0: 'NW', 1: 'SW', 2: 'SE', 3: 'NE'}
+            grid_origin = data15hd['ImageDescription'][
+                "ReferenceGridVIS_IR"]["GridOrigin"]
+            if grid_origin != 2:
+                raise NotImplementedError(
+                    'Grid origin not supported number: {}, {} corner'
+                    .format(grid_origin, origins[grid_origin])
+                )
+
             center_point = 3712/2
 
             north = int(sec15hd["NorthLineSelectedRectangle"]['Value'])
@@ -283,11 +295,25 @@ def get_area_extent(self, dsid):
                 "ReferenceGridVIS_IR"]["ColumnDirGridStep"] * 1000.0
             line_step = data15hd['ImageDescription'][
                 "ReferenceGridVIS_IR"]["LineDirGridStep"] * 1000.0
-
-            ll_c = (center_point - west - 0.5) * column_step
-            ll_l = (south - center_point + 1.5) * line_step
-            ur_c = (center_point - east + 0.5) * column_step
-            ur_l = (north - center_point + 0.5) * line_step
+            # section 3.1.4.2 of MSG Level 1.5 Image Data Format Description
+            earth_model = data15hd['GeometricProcessing']['EarthModel'][
+                'TypeOfEarthModel']
+            if earth_model == 2:
+                ns_offset = 0  # north +ve
+                we_offset = 0  # west +ve
+            elif earth_model == 1:
+                ns_offset = -0.5  # north +ve
+                we_offset = 0.5  # west +ve
+            else:
+                raise NotImplementedError(
+                    'unrecognised earth model: {}'.format(earth_model)
+                )
+
+            # section 3.1.5 of MSG Level 1.5 Image Data Format Description
+            ll_c = (center_point - west - 0.5 + we_offset) * column_step
+            ll_l = (south - center_point - 0.5 + ns_offset) * line_step
+            ur_c = (center_point - east + 0.5 + we_offset) * column_step
+            ur_l = (north - center_point + 0.5 + ns_offset) * line_step
 
             area_extent = (ll_c, ll_l, ur_c, ur_l)
 

satpy/tests/reader_tests/test_native_msg.py

@@ -27,7 +27,10 @@
 import sys
 
 import numpy as np
-from satpy.readers.native_msg import get_available_channels
+from satpy.readers.native_msg import (
+    NativeMSGFileHandler,
+    get_available_channels,
+)
 
 if sys.version_info < (2, 7):
     import unittest2 as unittest
@@ -104,6 +107,75 @@ def tearDown(self):
         pass
 
 
+class TestNativeMSGAreaExtent(unittest.TestCase):
+    """Test NativeMSGFileHandler.get_area_extent
+    The expected results have been verified by manually
+    inspecting the output of geoferenced imagery.
+    """
+    @staticmethod
+    def get_mock_file_handler(earth_model):
+        """
+        Mocked NativeMSGFileHandler with sufficient attributes for
+        NativeMSGFileHandler.get_area_extent to be able to execute.
+        """
+        header = {
+            '15_DATA_HEADER': {
+                'ImageDescription': {
+                    'ReferenceGridVIS_IR': {
+                        'ColumnDirGridStep': 3.0004032,
+                        'LineDirGridStep': 3.0004032,
+                        'GridOrigin': 2,  # south-east corner
+                    }
+                },
+                'GeometricProcessing': {
+                    'EarthModel': {'TypeOfEarthModel': earth_model}
+                }
+            },
+            '15_SECONDARY_PRODUCT_HEADER': {
+                'NorthLineSelectedRectangle': {'Value': 3712},
+                'EastColumnSelectedRectangle': {'Value': 1},
+                'WestColumnSelectedRectangle': {'Value': 3712},
+                'SouthLineSelectedRectangle': {'Value': 1},
+            }
+
+        }
+        return mock.Mock(header=header)
+
+    def setUp(self):
+        pass
+
+    def test_earthmodel1(self):
+        """TypeOfEarthModel = 1, need to offset by 0.5 pixels"""
+        calc_area_extent = NativeMSGFileHandler.get_area_extent(
+            self.get_mock_file_handler(earth_model=1),
+            mock.Mock(name='VIS006')  # mocked dsid (not 'HRV')
+        )
+        expected_area_extent = (
+            -5568748.275756836, -5568748.275756836,
+            5568748.275756836, 5568748.275756836
+        )
+        assertNumpyArraysEqual(
+            np.array(calc_area_extent), np.array(expected_area_extent)
+        )
+
+    def test_earthmodel2(self):
+        """TypeOfEarthModel = 2, do not offset"""
+        calc_area_extent = NativeMSGFileHandler.get_area_extent(
+            self.get_mock_file_handler(earth_model=2),
+            mock.Mock(name='VIS006')  # mocked dsid (not 'HRV')
+        )
+        expected_area_extent = (
+            -5570248.477339745, -5567248.074173927,
+            5567248.074173927, 5570248.477339745
+        )
+        assertNumpyArraysEqual(
+            np.array(calc_area_extent), np.array(expected_area_extent)
+        )
+
+    def tearDown(self):
+        pass
+
+
 def suite():
     """The test suite for test_scene.
     """