Merge pull request #837 from djhoese/bugfix-pyresample-crs

Fix Satpy tests to work with new versions of pyresample

satpy/tests/reader_tests/test_ahi_hsd.py

@@ -71,10 +71,13 @@ def test_region(self, fromfile, np2str):
                             'spare': ''}
 
             area_def = fh.get_area_def(None)
-            self.assertEqual(area_def.proj_dict, {'a': 6378137.0, 'b': 6356752.3,
-                                                  'h': 35785863.0, 'lon_0': 140.7,
-                                                  'proj': 'geos', 'units': 'm'})
-
+            proj_dict = area_def.proj_dict
+            self.assertEqual(proj_dict['a'], 6378137.0)
+            self.assertEqual(proj_dict['b'], 6356752.3)
+            self.assertEqual(proj_dict['h'], 35785863.0)
+            self.assertEqual(proj_dict['lon_0'], 140.7)
+            self.assertEqual(proj_dict['proj'], 'geos')
+            self.assertEqual(proj_dict['units'], 'm')
             self.assertEqual(area_def.area_extent, (592000.0038256244, 4132000.026701824,
                                                     1592000.0102878278, 5132000.033164027))
 
@@ -113,10 +116,13 @@ def test_segment(self, fromfile, np2str):
                             'spare': ''}
 
             area_def = fh.get_area_def(None)
-            self.assertEqual(area_def.proj_dict, {'a': 6378137.0, 'b': 6356752.3,
-                                                  'h': 35785863.0, 'lon_0': 140.7,
-                                                  'proj': 'geos', 'units': 'm'})
-
+            proj_dict = area_def.proj_dict
+            self.assertEqual(proj_dict['a'], 6378137.0)
+            self.assertEqual(proj_dict['b'], 6356752.3)
+            self.assertEqual(proj_dict['h'], 35785863.0)
+            self.assertEqual(proj_dict['lon_0'], 140.7)
+            self.assertEqual(proj_dict['proj'], 'geos')
+            self.assertEqual(proj_dict['units'], 'm')
             self.assertEqual(area_def.area_extent, (-5500000.035542117, -3300000.021325271,
                                                     5500000.035542117, -2200000.0142168473))
 

satpy/tests/reader_tests/test_hrit_base.py

@@ -142,12 +142,13 @@ def test_get_area_extent(self):
 
     def test_get_area_def(self):
         area = self.reader.get_area_def('VIS06')
-        self.assertEqual(area.proj_dict, {'a': 6378169.0,
-                                          'b': 6356583.8,
-                                          'h': 35785831.0,
-                                          'lon_0': 44.0,
-                                          'proj': 'geos',
-                                          'units': 'm'})
+        proj_dict = area.proj_dict
+        self.assertEqual(proj_dict['a'], 6378169.0)
+        self.assertEqual(proj_dict['b'], 6356583.8)
+        self.assertEqual(proj_dict['h'], 35785831.0)
+        self.assertEqual(proj_dict['lon_0'], 44.0)
+        self.assertEqual(proj_dict['proj'], 'geos')
+        self.assertEqual(proj_dict['units'], 'm')
         self.assertEqual(area.area_extent,
                          (-77771774058.38356, -77771774058.38356,
                           30310525626438.438, 3720765401003.719))

satpy/tests/reader_tests/test_seviri_l1b_hrit.py

@@ -349,12 +349,13 @@ def test_get_area_extent(self):
     def test_get_area_def(self):
         """Test getting the area def."""
         area = self.reader.get_area_def(DatasetID('VIS006'))
-        self.assertEqual(area.proj_dict, {'a': 6378169.0,
-                                          'b': 6356583.8,
-                                          'h': 35785831.0,
-                                          'lon_0': 44.0,
-                                          'proj': 'geos',
-                                          'units': 'm'})
+        proj_dict = area.proj_dict
+        self.assertEqual(proj_dict['a'], 6378169.0)
+        self.assertEqual(proj_dict['b'], 6356583.8)
+        self.assertEqual(proj_dict['h'], 35785831.0)
+        self.assertEqual(proj_dict['lon_0'], 44.0)
+        self.assertEqual(proj_dict['proj'], 'geos')
+        self.assertEqual(proj_dict['units'], 'm')
         self.assertEqual(area.area_extent,
                          (-77771774058.38356, -3720765401003.719,
                           30310525626438.438, 77771774058.38356))

satpy/tests/test_config.py

@@ -59,13 +59,25 @@ def test_areas_pyproj(self):
         """Test all areas have valid projections with pyproj."""
         import pyproj
         from pyresample import parse_area_file
+        from pyresample.geometry import SwathDefinition
         from satpy.resample import get_area_file
+        import numpy as np
+        import xarray as xr
 
+        lons = np.array([[0, 0.1, 0.2], [0.05, 0.15, 0.25]])
+        lats = np.array([[0, 0.1, 0.2], [0.05, 0.15, 0.25]])
+        lons = xr.DataArray(lons)
+        lats = xr.DataArray(lats)
+        swath_def = SwathDefinition(lons, lats)
         all_areas = parse_area_file(get_area_file())
         for area_obj in all_areas:
-            if getattr(area_obj, 'optimize_projection', False):
-                # the PROJ.4 is known to not be valid on this DynamicAreaDef
-                continue
+            if hasattr(area_obj, 'freeze'):
+                try:
+                    area_obj = area_obj.freeze(lonslats=swath_def)
+                except RuntimeError:
+                    # we didn't provide enough info to freeze, hard to guess
+                    # in a generic test so just skip this area
+                    continue
             proj_dict = area_obj.proj_dict
             _ = pyproj.Proj(proj_dict)
 
@@ -79,13 +91,32 @@ def test_areas_rasterio(self):
             return unittest.skip("RasterIO 1.0+ required")
 
         from pyresample import parse_area_file
+        from pyresample.geometry import SwathDefinition
         from satpy.resample import get_area_file
+        import numpy as np
+        import xarray as xr
+
+        lons = np.array([[0, 0.1, 0.2], [0.05, 0.15, 0.25]])
+        lats = np.array([[0, 0.1, 0.2], [0.05, 0.15, 0.25]])
+        lons = xr.DataArray(lons)
+        lats = xr.DataArray(lats)
+        swath_def = SwathDefinition(lons, lats)
         all_areas = parse_area_file(get_area_file())
         for area_obj in all_areas:
-            if getattr(area_obj, 'optimize_projection', False):
-                # the PROJ.4 is known to not be valid on this DynamicAreaDef
-                continue
+            if hasattr(area_obj, 'freeze'):
+                try:
+                    area_obj = area_obj.freeze(lonslats=swath_def)
+                except RuntimeError:
+                    # we didn't provide enough info to freeze, hard to guess
+                    # in a generic test so just skip this area
+                    continue
             proj_dict = area_obj.proj_dict
+            if proj_dict['proj'] in ('ob_tran', 'nsper') and \
+                    'wktext' not in proj_dict:
+                # FIXME: rasterio doesn't understand ob_tran unless +wktext
+                # See: https://github.com/pyproj4/pyproj/issues/357
+                # pyproj 2.0+ seems to drop wktext from PROJ dict
+                continue
             _ = CRS.from_dict(proj_dict)
 
 

satpy/tests/writer_tests/test_cf.py

@@ -733,8 +733,14 @@ def test_area2gridmapping(self):
         with mock.patch('satpy.writers.cf_writer.warnings.warn') as warn:
             res, grid_mapping = area2gridmapping(ds)
             warn.assert_called()
-            self.assertDictEqual(dict(pyresample.geometry.proj4_str_to_dict(res.attrs['grid_proj4'])),
-                                 dict(pyresample.geometry.proj4_str_to_dict(proj_str)))
+            proj_dict = pyresample.geometry.proj4_str_to_dict(res.attrs['grid_proj4'])
+            self.assertEqual(proj_dict['lon_0'], 4.535)
+            self.assertEqual(proj_dict['lat_0'], 46.0)
+            self.assertEqual(proj_dict['o_lon_p'], -5.465)
+            self.assertEqual(proj_dict['o_lat_p'], 90.0)
+            self.assertEqual(proj_dict['proj'], 'ob_tran')
+            self.assertEqual(proj_dict['o_proj'], 'stere')
+            self.assertEqual(proj_dict['ellps'], 'WGS84')
             self.assertEqual(grid_mapping, cosmo_expected)
 
     def test_area2lonlat(self):

satpy/writers/mitiff.py

@@ -203,12 +203,19 @@ def _add_proj4_string(self, datasets, first_dataset):
         proj4_string = " Proj string: "
 
         if isinstance(datasets, list):
-            proj4_string += first_dataset.attrs['area'].proj4_string
+            area = first_dataset.attrs['area']
         else:
-            proj4_string += datasets.attrs['area'].proj4_string
+            area = datasets.attrs['area']
+        # Use pyproj's CRS object to get a valid EPSG code if possible
+        # only in newer pyresample versions with pyproj 2.0+ installed
+        if hasattr(area, 'crs') and area.crs.to_epsg() is not None:
+            proj4_string += "+init=EPSG:{}".format(area.crs.to_epsg())
+        else:
+            proj4_string += area.proj_str
 
         x_0 = 0
         y_0 = 0
+        # FUTURE: Use pyproj 2.0+ to convert EPSG to PROJ4 if possible
         if 'EPSG:32631' in proj4_string:
             proj4_string = proj4_string.replace("+init=EPSG:32631",
                                                 "+proj=etmerc +lat_0=0 +lon_0=3 +k=0.9996 +ellps=WGS84 +datum=WGS84")
@@ -246,14 +253,14 @@ def _add_proj4_string(self, datasets, first_dataset):
         if 'units' not in proj4_string:
             proj4_string += ' +units=km'
 
-        if isinstance(datasets, list):
+        if 'x_0' not in proj4_string and isinstance(datasets, list):
             proj4_string += ' +x_0=%.6f' % (
                 (-first_dataset.attrs['area'].area_extent[0] +
                  first_dataset.attrs['area'].pixel_size_x) + x_0)
             proj4_string += ' +y_0=%.6f' % (
                 (-first_dataset.attrs['area'].area_extent[1] +
                  first_dataset.attrs['area'].pixel_size_y) + y_0)
-        else:
+        elif 'x_0' not in proj4_string:
             proj4_string += ' +x_0=%.6f' % (
                 (-datasets.attrs['area'].area_extent[0] +
                  datasets.attrs['area'].pixel_size_x) + x_0)