Merge pull request #476 from djhoese/bugfix-grib-area

Fix handling of navigation in a grib file with lons greater than 180

satpy/readers/grib.py

@@ -139,6 +139,11 @@ def _get_message(self, ds_info):
 
     def _area_def_from_msg(self, msg):
         proj_params = msg.projparams.copy()
+        # correct for longitudes over 180
+        for lon_param in ['lon_0', 'lon_1', 'lon_2']:
+            if proj_params.get(lon_param, 0) > 180:
+                proj_params[lon_param] -= 360
+
         if proj_params['proj'] == 'cyl':
             proj_params['proj'] = 'eqc'
             proj = Proj(**proj_params)
@@ -171,9 +176,20 @@ def _area_def_from_msg(self, msg):
         else:
             lats, lons = msg.latlons()
             shape = lats.shape
+            # take the corner points only
+            lons = lons[([0, 0, -1, -1], [0, -1, 0, -1])]
+            lats = lats[([0, 0, -1, -1], [0, -1, 0, -1])]
+            # correct for longitudes over 180
+            lons[lons > 180] -= 360
+
             proj = Proj(**proj_params)
-            min_x, min_y = proj(lons[-1, 0], lats[-1, 0])
-            max_x, max_y = proj(lons[0, -1], lats[0, -1])
+            x, y = proj(lons, lats)
+            x[np.abs(x) >= 1e30] = np.nan
+            y[np.abs(y) >= 1e30] = np.nan
+            min_x = np.nanmin(x)
+            max_x = np.nanmax(x)
+            min_y = np.nanmin(y)
+            max_y = np.nanmax(y)
             extents = (min_x, min_y, max_x, max_y)
 
         return geometry.AreaDefinition(

satpy/tests/reader_tests/test_grib.py

@@ -22,13 +22,17 @@
 class FakeMessage(object):
     """Fake message returned by pygrib.open().message(x)."""
 
-    def __init__(self, values, proj_params=None, **attrs):
+    def __init__(self, values, proj_params=None, latlons=None, **attrs):
         super(FakeMessage, self).__init__()
         self.attrs = attrs
         self.values = values
         if proj_params is None:
             proj_params = {'a': 6371229, 'b': 6371229, 'proj': 'cyl'}
         self.projparams = proj_params
+        self._latlons = latlons
+
+    def latlons(self):
+        return self._latlons
 
     def __getitem__(self, item):
         return self.attrs[item]
@@ -37,7 +41,7 @@ def __getitem__(self, item):
 class FakeGRIB(object):
     """Fake GRIB file returned by pygrib.open."""
 
-    def __init__(self, messages=None):
+    def __init__(self, messages=None, proj_params=None, latlons=None):
         super(FakeGRIB, self).__init__()
         if messages is not None:
             self._messages = messages
@@ -62,6 +66,8 @@ def __init__(self, messages=None):
                     minimum=100.,
                     maximum=200.,
                     typeOfLevel='isobaricInhPa',
+                    proj_params=proj_params,
+                    latlons=latlons,
                 ),
                 FakeMessage(
                     values=np.arange(25.).reshape((5, 5)),
@@ -82,6 +88,8 @@ def __init__(self, messages=None):
                     minimum=100.,
                     maximum=200.,
                     typeOfLevel='isobaricInhPa',
+                    proj_params=proj_params,
+                    latlons=latlons,
                 ),
                 FakeMessage(
                     values=np.arange(25.).reshape((5, 5)),
@@ -102,6 +110,8 @@ def __init__(self, messages=None):
                     minimum=100.,
                     maximum=200.,
                     typeOfLevel='isobaricInhPa',
+                    proj_params=proj_params,
+                    latlons=latlons,
                 ),
             ]
         self.messages = len(self._messages)
@@ -176,10 +186,46 @@ def test_load_all(self, pg):
             self.assertEqual(v.attrs['units'], 'K')
             self.assertIsInstance(v, xr.DataArray)
 
+    @mock.patch('satpy.readers.grib.pygrib')
+    def test_load_all_lcc(self, pg):
+        """Test loading all test datasets with lcc projections"""
+        lons = np.array([
+            [12.19, 0, 0, 0, 14.34208538],
+            [0, 0, 0, 0, 0],
+            [0, 0, 0, 0, 0],
+            [0, 0, 0, 0, 0],
+            [54.56534318, 0, 0, 0, 57.32843565]])
+        lats = np.array([
+            [-133.459, 0, 0, 0, -65.12555139],
+            [0, 0, 0, 0, 0],
+            [0, 0, 0, 0, 0],
+            [0, 0, 0, 0, 0],
+            [-152.8786225, 0, 0, 0, -49.41598659]])
+        pg.open.return_value = FakeGRIB(
+            proj_params={
+                'a': 6371229, 'b': 6371229, 'proj': 'lcc',
+                'lon_0': 265.0, 'lat_0': 25.0,
+                'lat_1': 25.0, 'lat_2': 25.0},
+            latlons=(lats, lons))
+        from satpy.readers import load_reader
+        from satpy import DatasetID
+        r = load_reader(self.reader_configs)
+        loadables = r.select_files_from_pathnames([
+            'gfs.t18z.sfluxgrbf106.grib2',
+        ])
+        r.create_filehandlers(loadables)
+        datasets = r.load([
+            DatasetID(name='t', level=100),
+            DatasetID(name='t', level=200),
+            DatasetID(name='t', level=300)])
+        self.assertEqual(len(datasets), 3)
+        for v in datasets.values():
+            self.assertEqual(v.attrs['units'], 'K')
+            self.assertIsInstance(v, xr.DataArray)
+
 
 def suite():
-    """The test suite for test_viirs_l1b.
-    """
+    """The test suite for test_grib."""
     loader = unittest.TestLoader()
     mysuite = unittest.TestSuite()
     mysuite.addTest(loader.loadTestsFromTestCase(TestGRIBReader))