Merge pull request #53 from sfinkens/feature-rotate

Use satpy to rotate SEVIRI images

level1c4pps/seviri2pps_lib.py

@@ -89,17 +89,16 @@ class UnexpectedSatpyVersion(Exception):
 # MSG4-SEVI-MSG15-1234-NA-20190409121243.927000000Z
 
 
-def load_and_calibrate(filenames, apply_sun_earth_distance_correction):
+def load_and_calibrate(filenames, apply_sun_earth_distance_correction, rotate):
     """Load and calibrate data.
 
     Uses inter-calibration coefficients from Meirink et al.
 
     Args:
         filenames: List of data files
-        filename_info: Corresponding information from the filenames
-        file_format: Specifies the file format (HRIT/Native)
         apply_sun_earth_distance_correction: If True, apply sun-earth-distance-
             correction to visible channels.
+        rotate: Rotate image so that pixel (0, 0) is N-W.
 
     Returns:
         Satpy scene holding calibrated channels
@@ -108,24 +107,48 @@ def load_and_calibrate(filenames, apply_sun_earth_distance_correction):
     parser = SEVIRIFilenameParser()
     file_format, info = parser.parse(os.path.basename(filenames[0]))
 
-    # Get calibration coefficients (time-dependent)
-    coefs = get_calibration_for_time(platform=info['platform_shortname'],
-                                     time=info['start_time'])
+    calib_coefs = get_calibration_for_time(
+        platform=info['platform_shortname'],
+        time=info['start_time']
+    )
+    scn_ = _create_scene(file_format, filenames, calib_coefs)
+    _check_is_seviri_data(scn_)
+    _load_bands(scn_, rotate)
+    _update_scene_attrs(scn_, {'image_rotated': rotate})
 
-    # Load and calibrate data
-    scn_ = Scene(reader=file_format,
-                 filenames=filenames,
-                 reader_kwargs={'calib_mode': CALIB_MODE,
-                                'ext_calib_coefs': coefs})
-    if not scn_.attrs['sensor'] == {'seviri'}:
-        raise ValueError('Not SEVIRI data')
-    scn_.load(BANDNAMES)
     if not apply_sun_earth_distance_correction:
         remove_sun_earth_distance_correction(scn_)
 
     return scn_
 
 
+def _create_scene(file_format, filenames, calib_coefs):
+    return Scene(reader=file_format,
+                 filenames=filenames,
+                 reader_kwargs={'calib_mode': CALIB_MODE,
+                                'ext_calib_coefs': calib_coefs})
+
+
+def _check_is_seviri_data(scene):
+    if not scene.attrs['sensor'] == {'seviri'}:
+        raise ValueError('Not SEVIRI data')
+
+
+def _load_bands(scene, rotate):
+    scene.load(
+        BANDNAMES,
+        upper_right_corner=_get_upper_right_corner(rotate)
+    )
+
+
+def _get_upper_right_corner(rotation_flag):
+    return 'NE' if rotation_flag else 'native'
+
+
+def _update_scene_attrs(scene, attrs):
+    scene.attrs.update(attrs)
+
+
 def remove_sun_earth_distance_correction(scene):
     """Remove sun-earth-distance correction from visible channels.
 
@@ -142,15 +165,6 @@ def remove_sun_earth_distance_correction(scene):
             scene[band] = remove_earthsun_distance_correction(scene[band])
 
 
-def rotate_band(scene, band):
-    """Rotate band by 180 degrees."""
-    scene[band] = scene[band].reindex(x=scene[band].x[::-1],
-                                      y=scene[band].y[::-1])
-    llx, lly, urx, ury = scene[band].attrs['area'].area_extent
-    scene[band].attrs['area'] = scene[band].attrs['area'].copy(
-        area_extent=[urx, ury, llx, lly])
-
-
 def get_lonlats(dataset):
     """Get lat/lon coordinates."""
     lons, lats = dataset.attrs['area'].get_lonlats()
@@ -548,15 +562,10 @@ def process_one_scan(tslot_files, out_path, rotate=True, engine='h5netcdf',
     tic = time.time()
     scn_ = load_and_calibrate(
         tslot_files,
-        apply_sun_earth_distance_correction=apply_sun_earth_distance_correction
+        apply_sun_earth_distance_correction=apply_sun_earth_distance_correction,
+        rotate=rotate
     )
 
-    # By default pixel (0,0) is S-E. Rotate bands so that (0,0) is N-W.
-    if rotate:
-        for band in BANDNAMES:
-            rotate_band(scn_, band)
-    scn_.attrs['image_rotated'] = rotate
-
     # Find lat/lon data
     lons, lats = get_lonlats(scn_['IR_108'])
 

level1c4pps/tests/test_seviri2pps.py

@@ -24,7 +24,6 @@
 
 import datetime as dt
 import numpy as np
-from pyresample.geometry import AreaDefinition
 import unittest
 try:
     from unittest import mock
@@ -37,40 +36,43 @@
 import level1c4pps.calibration_coefs as calib
 
 
+def get_fake_scene():
+    scene = Scene()
+    start_time = dt.datetime(2020, 1, 1, 12)
+    scene['VIS006'] = xr.DataArray(
+        [[1, 2],
+         [3, 4]],
+        dims=('y', 'x'),
+        attrs={'calibration': 'reflectance',
+               'sun_earth_distance_correction_applied': True,
+               'start_time': start_time}
+    )
+    scene['IR_108'] = xr.DataArray(
+        [[5, 6],
+         [7, 8]],
+        dims=('y', 'x'),
+        attrs={'calibration': 'brightness_temperature',
+               'start_time': start_time}
+    )
+    scene.attrs['sensor'] = {'seviri'}
+    return scene
+
+
 class TestSeviri2PPS(unittest.TestCase):
     """Test for SEVIRI converter."""
-
 
     @mock.patch('level1c4pps.seviri2pps_lib.Scene')
     def test_load_and_calibrate(self, mocked_scene):
         """Test loading and calibrating the data."""
 
-        # Create test scene
-        scene = Scene()
-        start_time = dt.datetime(2020, 1, 1, 12)
-        scene['VIS006'] = xr.DataArray(
-            [[1, 2],
-             [3, 4]],
-            dims=('y', 'x'),
-            attrs={'calibration': 'reflectance',
-                   'sun_earth_distance_correction_applied': True,
-                   'start_time': start_time}
-        )
-        scene['IR_108'] = xr.DataArray(
-            [[5, 6],
-             [7, 8]],
-            dims=('y', 'x'),
-            attrs={'calibration': 'brightness_temperature',
-                   'start_time': start_time}
-        )
-        scene.attrs['sensor'] = {'seviri'}
-        mocked_scene.return_value = scene
+        mocked_scene.return_value = get_fake_scene()
 
         # Load and calibrate
         filenames = ['MSG4-SEVI-MSG15-1234-NA-20190409121243.927000000Z']
         res = seviri2pps.load_and_calibrate(
             filenames,
-            apply_sun_earth_distance_correction=False
+            apply_sun_earth_distance_correction=False,
+            rotate=False
         )
 
         # Compare results and expectations
@@ -90,39 +92,19 @@ def test_load_and_calibrate(self, mocked_scene):
             res['VIS006'].attrs['sun_earth_distance_correction_applied'],
         )
 
-    def test_rotate_band(self):
-        """Test rotation of bands."""
-        area = AreaDefinition(area_id='test',
-                              description='test',
-                              proj_id='test',
-                              projection={'proj': 'geos', 'h': 12345},
-                              width=3,
-                              height=3,
-                              area_extent=[1001, 1002, -1003, -1004])
-        data = xr.DataArray(data=[[1, 2, 3],
-                                  [4, 5, 6],
-                                  [7, 8, 9]],
-                            dims=('y', 'x'),
-                            coords=[('y', [1.1, 0, -1.1]), ('x', [1, 0, -1])],
-                            attrs={'area': area})
-        scene = {'data': data}
-
-        # Rotate
-        seviri2pps.rotate_band(scene, 'data')
-
-        # Check results
-        self.assertTupleEqual(scene['data'].attrs['area'].area_extent,
-                              (-1003, -1004, 1001, 1002))
-        np.testing.assert_array_equal(scene['data']['x'], [-1, 0, 1])
-        np.testing.assert_array_equal(scene['data']['y'], [-1.1, 0, 1.1])
-        np.testing.assert_array_equal(scene['data'], [[9, 8, 7],
-                                                      [6, 5, 4],
-                                                      [3, 2, 1]])
-        lons, lats = scene['data'].attrs['area'].get_lonlats()
-        self.assertTrue(lons[0, 0] < 0)
-        self.assertTrue(lons[0, 2] > 0)
-        self.assertTrue(lats[0, 0] > 0)
-        self.assertTrue(lons[2, 0] < 0)
+    @mock.patch('level1c4pps.seviri2pps_lib.Scene')
+    def test_load_and_calibrate_with_rotation(self, mocked_scene):
+        scene = get_fake_scene()
+        scene.load = mock.MagicMock()
+        mocked_scene.return_value = scene
+        filenames = ['MSG4-SEVI-MSG15-1234-NA-20190409121243.927000000Z']
+        res = seviri2pps.load_and_calibrate(
+            filenames,
+            apply_sun_earth_distance_correction=False,
+            rotate=True
+        )
+        scene.load.assert_called_with(mock.ANY, upper_right_corner='NE')
+        self.assertTrue(res.attrs['image_rotated'])
 
     def test_get_lonlats(self):
         """Test getting lat/lon coordinates."""