Merge pull request #1296 from djhoese/bugfix-grib-prime-meridian

Fix grib reader handling for data on 0-360 longitude

satpy/readers/grib.py

@@ -44,10 +44,10 @@
 
 
 class GRIBFileHandler(BaseFileHandler):
-    """File handler for grib files."""
+    """Generic GRIB file handler."""
 
     def __init__(self, filename, filename_info, filetype_info):
-        """Init the file handler."""
+        """Open grib file and do initial message parsing."""
         super(GRIBFileHandler, self).__init__(filename, filename_info, filetype_info)
 
         self._msg_datasets = {}
@@ -143,62 +143,100 @@ def _get_message(self, ds_info):
                 msg = self._idx(**{k: ds_info[k] for k in msg_keys})[0]
             return msg
 
-    def _area_def_from_msg(self, msg):
-        proj_params = msg.projparams.copy()
+    @staticmethod
+    def _correct_cyl_minmax_xy(proj_params, min_lon, min_lat, max_lon, max_lat):
+        proj = Proj(**proj_params)
+        min_x, min_y = proj(min_lon, min_lat)
+        max_x, max_y = proj(max_lon, max_lat)
+        if max_x <= min_x:
+            # wrap around
+            # make 180 longitude the prime meridian
+            # assuming we are going from 0 to 360 longitude
+            proj_params['pm'] = 180
+            proj = Proj(**proj_params)
+            # recompute x/y extents with this new projection
+            min_x, min_y = proj(min_lon, min_lat)
+            max_x, max_y = proj(max_lon, max_lat)
+        return proj_params, (min_x, min_y, max_x, max_y)
+
+    @staticmethod
+    def _get_cyl_minmax_lonlat(lons, lats):
+        min_lon = lons[0]
+        max_lon = lons[-1]
+        min_lat = lats[0]
+        max_lat = lats[-1]
+        if min_lat > max_lat:
+            # lats aren't in the order we thought they were, flip them
+            min_lat, max_lat = max_lat, min_lat
+        return min_lon, min_lat, max_lon, max_lat
+
+    def _get_cyl_area_info(self, msg, proj_params):
+        proj_params['proj'] = 'eqc'
+        lons = msg['distinctLongitudes']
+        lats = msg['distinctLatitudes']
+        shape = (lats.shape[0], lons.shape[0])
+        minmax_lonlat = self._get_cyl_minmax_lonlat(lons, lats)
+        proj_params, minmax_xy = self._correct_cyl_minmax_xy(proj_params, *minmax_lonlat)
+        extents = self._get_extents(*minmax_xy, shape)
+        return proj_params, shape, extents
+
+    @staticmethod
+    def _get_extents(min_x, min_y, max_x, max_y, shape):
+        half_x = abs((max_x - min_x) / (shape[1] - 1)) / 2.
+        half_y = abs((max_y - min_y) / (shape[0] - 1)) / 2.
+        return min_x - half_x, min_y - half_y, max_x + half_x, max_y + half_y
+
+    @staticmethod
+    def _get_corner_xy(proj_params, lons, lats, scans_positively):
+        proj = Proj(**proj_params)
+        x, y = proj(lons, lats)
+        if scans_positively:
+            min_x, min_y = x[0], y[0]
+            max_x, max_y = x[3], y[3]
+        else:
+            min_x, min_y = x[2], y[2]
+            max_x, max_y = x[1], y[1]
+        return min_x, min_y, max_x, max_y
+
+    @staticmethod
+    def _get_corner_lonlat(proj_params, lons, lats):
+        # 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])]
+        # if we have longitudes over 180, assume 0-360
+        if (lons > 180).any():
+            # make 180 longitude the prime meridian
+            proj_params['pm'] = 180
+        return proj_params, lons, lats
+
+    def _get_area_info(self, msg, proj_params):
+        lats, lons = msg.latlons()
+        shape = lats.shape
+        scans_positively = (msg.valid_key('jScansPositively') and
+                            msg['jScansPositively'] == 1)
+        proj_params, lons, lats = self._get_corner_lonlat(
+            proj_params, lons, lats)
+        minmax_xy = self._get_corner_xy(proj_params, lons, lats, scans_positively)
+        extents = self._get_extents(*minmax_xy, shape)
+        return proj_params, shape, extents
+
+    @staticmethod
+    def _correct_proj_params_over_prime_meridian(proj_params):
         # 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
+        return proj_params
 
-        if proj_params['proj'] == 'cyl':
-            proj_params['proj'] = 'eqc'
-            proj = Proj(**proj_params)
-            lons = msg['distinctLongitudes']
-            lats = msg['distinctLatitudes']
-            min_lon = lons[0]
-            max_lon = lons[-1]
-            min_lat = lats[0]
-            max_lat = lats[-1]
-            if min_lat > max_lat:
-                # lats aren't in the order we thought they were, flip them
-                # we also need to flip the data in the data loading section
-                min_lat, max_lat = max_lat, min_lat
-            shape = (lats.shape[0], lons.shape[0])
-            min_x, min_y = proj(min_lon, min_lat)
-            max_x, max_y = proj(max_lon, max_lat)
-            if max_x < min_x and 'over' not in proj_params:
-                # wrap around
-                proj_params['over'] = True
-                proj = Proj(**proj_params)
-                max_x, max_y = proj(max_lon, max_lat)
-            pixel_size_x = (max_x - min_x) / (shape[1] - 1)
-            pixel_size_y = (max_y - min_y) / (shape[0] - 1)
-            extents = (
-                min_x - pixel_size_x / 2.,
-                min_y - pixel_size_y / 2.,
-                max_x + pixel_size_x / 2.,
-                max_y + pixel_size_y / 2.,
-            )
-        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
+    def _area_def_from_msg(self, msg):
+        proj_params = msg.projparams.copy()
+        proj_params = self._correct_proj_params_over_prime_meridian(proj_params)
 
-            proj = Proj(**proj_params)
-            x, y = proj(lons, lats)
-            if msg.valid_key('jScansPositively') and msg['jScansPositively'] == 1:
-                min_x, min_y = x[0], y[0]
-                max_x, max_y = x[3], y[3]
-            else:
-                min_x, min_y = x[2], y[2]
-                max_x, max_y = x[1], y[1]
-            half_x = abs((max_x - min_x) / (shape[1] - 1)) / 2.
-            half_y = abs((max_y - min_y) / (shape[0] - 1)) / 2.
-            extents = (min_x - half_x, min_y - half_y, max_x + half_x, max_y + half_y)
+        if proj_params['proj'] in ('cyl', 'eqc'):
+            # eqc projection that goes from 0 to 360
+            proj_params, shape, extents = self._get_cyl_area_info(msg, proj_params)
+        else:
+            proj_params, shape, extents = self._get_area_info(msg, proj_params)
 
         return geometry.AreaDefinition(
             'on-the-fly grib area',

satpy/tests/reader_tests/test_grib.py

@@ -19,14 +19,50 @@
 
 import os
 import sys
-import unittest
 from unittest import mock
 
 import numpy as np
 import xarray as xr
+import pytest
 
 from satpy.dataset import DataQuery
 
+# Parameterized cases
+TEST_ARGS = ('proj_params', 'lon_corners', 'lat_corners')
+TEST_PARAMS = (
+    (None, None, None),  # cyl default case
+    (
+        {
+            'a': 6371229, 'b': 6371229, 'proj': 'lcc',
+            'lon_0': 265.0, 'lat_0': 25.0,
+            'lat_1': 25.0, 'lat_2': 25.0
+        },
+        [-133.459, -65.12555139, -152.8786225, -49.41598659],
+        [12.19, 14.34208538, 54.56534318, 57.32843565]
+    ),
+)
+
+
+def _round_trip_projection_lonlat_check(area):
+    """Check that X/Y coordinates can be transformed multiple times.
+
+    Many GRIB files include non-standard projects that work for the
+    initial transformation of X/Y coordinates to longitude/latitude,
+    but may fail in the reverse transformation. For example, an eqc
+    projection that goes from 0 longitude to 360 longitude. The X/Y
+    coordinates may accurately go from the original X/Y metered space
+    to the correct longitude/latitude, but transforming those coordinates
+    back to X/Y space will produce the wrong result.
+
+    """
+    from pyproj import Proj
+    p = Proj(area.crs)
+    x, y = area.get_proj_vectors()
+    lon, lat = p(x, y, inverse=True)
+    x2, y2 = p(lon, lat)
+    np.testing.assert_almost_equal(x, x2)
+    np.testing.assert_almost_equal(y, y2)
+
 
 class FakeMessage(object):
     """Fake message returned by pygrib.open().message(x)."""
@@ -157,12 +193,12 @@ def __exit__(self, exc_type, exc_val, exc_tb):
         pass
 
 
-class TestGRIBReader(unittest.TestCase):
+class TestGRIBReader:
     """Test GRIB Reader."""
 
     yaml_file = "grib.yaml"
 
-    def setUp(self):
+    def setup_method(self):
         """Wrap pygrib to read fake data."""
         from satpy.config import config_search_paths
         self.reader_configs = config_search_paths(os.path.join('readers', self.yaml_file))
@@ -174,23 +210,61 @@ def setUp(self):
         self.orig_pygrib = pygrib
         sys.modules['pygrib'] = mock.MagicMock()
 
-    def tearDown(self):
+    def teardown_method(self):
         """Re-enable pygrib import."""
         sys.modules['pygrib'] = self.orig_pygrib
 
-    @mock.patch('satpy.readers.grib.pygrib')
-    def test_init(self, pg):
+    def _get_test_datasets(self, dataids, fake_pygrib=None):
+        from satpy.readers import load_reader
+        if fake_pygrib is None:
+            fake_pygrib = FakeGRIB()
+
+        with mock.patch('satpy.readers.grib.pygrib') as pg:
+            pg.open.return_value = fake_pygrib
+            r = load_reader(self.reader_configs)
+            loadables = r.select_files_from_pathnames([
+                'gfs.t18z.sfluxgrbf106.grib2',
+            ])
+            r.create_filehandlers(loadables)
+            datasets = r.load(dataids)
+        return datasets
+
+    @staticmethod
+    def _get_fake_pygrib(proj_params, lon_corners, lat_corners):
+        latlons = None
+        if lon_corners is not None:
+            lats = np.array([
+                [lat_corners[0], 0, 0, 0, lat_corners[1]],
+                [0, 0, 0, 0, 0],
+                [0, 0, 0, 0, 0],
+                [0, 0, 0, 0, 0],
+                [lat_corners[2], 0, 0, 0, lat_corners[3]]])
+            lons = np.array([
+                [lon_corners[0], 0, 0, 0, lon_corners[1]],
+                [0, 0, 0, 0, 0],
+                [0, 0, 0, 0, 0],
+                [0, 0, 0, 0, 0],
+                [lon_corners[2], 0, 0, 0, lon_corners[3]]])
+            latlons = (lats, lons)
+
+        fake_pygrib = FakeGRIB(
+            proj_params=proj_params,
+            latlons=latlons)
+        return fake_pygrib
+
+    def test_init(self):
         """Test basic init with no extra parameters."""
-        pg.open.return_value = FakeGRIB()
         from satpy.readers import load_reader
-        r = load_reader(self.reader_configs)
-        loadables = r.select_files_from_pathnames([
-            'gfs.t18z.sfluxgrbf106.grib2',
-        ])
-        self.assertEqual(len(loadables), 1)
-        r.create_filehandlers(loadables)
-        # make sure we have some files
-        self.assertTrue(r.file_handlers)
+        with mock.patch('satpy.readers.grib.pygrib') as pg:
+            pg.open.return_value = FakeGRIB()
+            r = load_reader(self.reader_configs)
+            loadables = r.select_files_from_pathnames([
+                'gfs.t18z.sfluxgrbf106.grib2',
+            ])
+            assert len(loadables) == 1
+            r.create_filehandlers(loadables)
+            # make sure we have some files
+            assert r.file_handlers
 
     def test_file_pattern(self):
         """Test matching of file patterns."""
@@ -205,59 +279,31 @@ def test_file_pattern(self):
 
         r = load_reader(self.reader_configs)
         files = r.select_files_from_pathnames(filenames)
-        self.assertEqual(len(files), 4)
+        assert len(files) == 4
 
-    @mock.patch('satpy.readers.grib.pygrib')
-    def test_load_all(self, pg):
+    @pytest.mark.parametrize(TEST_ARGS, TEST_PARAMS)
+    def test_load_all(self, proj_params, lon_corners, lat_corners):
         """Test loading all test datasets."""
-        pg.open.return_value = FakeGRIB()
-        from satpy.readers import load_reader
-        r = load_reader(self.reader_configs)
-        loadables = r.select_files_from_pathnames([
-            'gfs.t18z.sfluxgrbf106.grib2',
-        ])
-        r.create_filehandlers(loadables)
-        datasets = r.load([
-            DataQuery(name='t', level=100, modifiers=tuple()),
-            DataQuery(name='t', level=200, modifiers=tuple()),
-            DataQuery(name='t', level=300, modifiers=tuple())])
-        self.assertEqual(len(datasets), 3)
-        for v in datasets.values():
-            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
-        r = load_reader(self.reader_configs)
-        loadables = r.select_files_from_pathnames([
-            'gfs.t18z.sfluxgrbf106.grib2',
-        ])
-        r.create_filehandlers(loadables)
-        datasets = r.load([
+        fake_pygrib = self._get_fake_pygrib(proj_params, lon_corners, lat_corners)
+        dataids = [
             DataQuery(name='t', level=100, modifiers=tuple()),
             DataQuery(name='t', level=200, modifiers=tuple()),
-            DataQuery(name='t', level=300, modifiers=tuple())])
-        self.assertEqual(len(datasets), 3)
+            DataQuery(name='t', level=300, modifiers=tuple())
+        ]
+        datasets = self._get_test_datasets(dataids, fake_pygrib)
+
+        assert len(datasets) == 3
         for v in datasets.values():
-            self.assertEqual(v.attrs['units'], 'K')
-            self.assertIsInstance(v, xr.DataArray)
+            assert v.attrs['units'] == 'K'
+            assert isinstance(v, xr.DataArray)
+
+    @pytest.mark.parametrize(TEST_ARGS, TEST_PARAMS)
+    def test_area_def_crs(self, proj_params, lon_corners, lat_corners):
+        """Check that the projection is accurate."""
+        fake_pygrib = self._get_fake_pygrib(proj_params, lon_corners, lat_corners)
+        dataids = [DataQuery(name='t', level=100, modifiers=tuple())]
+        datasets = self._get_test_datasets(dataids, fake_pygrib)
+        area = datasets['t'].attrs['area']
+        if not hasattr(area, 'crs'):
+            pytest.skip("Can't test with pyproj < 2.0")
+        _round_trip_projection_lonlat_check(area)