Merge pull request #901 from mraspaud/fix-compact-viirs-pole-interpol…

…ation

Fix compact viirs angle interpolation at the poles

satpy/readers/viirs_compact.py

@@ -184,7 +184,8 @@ def read_geo(self, key, info):
                  ('dnb_lunar_azimuth_angle', 'dnb_lunar_zenith_angle'):
                  ("LunarAzimuthAngle", "LunarZenithAngle"),
                  }
-
+        if self.lons is None or self.lats is None:
+            self.lons, self.lats = self.navigate()
         for pair, fkeys in pairs.items():
             if key.name in pair:
                 if (self.cache.get(pair[0]) is None
@@ -324,15 +325,16 @@ def navigate(self):
     def angles(self, azi_name, zen_name):
         """Compute the angle datasets."""
         all_lat = da.from_array(self.geostuff["Latitude"])
+        all_lon = da.from_array(self.geostuff["Longitude"])
         all_zen = da.from_array(self.geostuff[zen_name])
         all_azi = da.from_array(self.geostuff[azi_name])
 
         res = []
-
         param_start = 0
         for tpz_size, nb_tpz, start in zip(self.tpz_sizes, self.nb_tpzs,
                                            self.group_locations):
             lat = all_lat[:, start:start + nb_tpz + 1]
+            lon = all_lon[:, start:start + nb_tpz + 1]
             zen = all_zen[:, start:start + nb_tpz + 1]
             azi = all_azi[:, start:start + nb_tpz + 1]
 
@@ -344,12 +346,13 @@ def angles(self, azi_name, zen_name):
             if (np.max(azi) - np.min(azi) > 5) or (np.min(zen) < 10) or (
                     np.max(abs(lat)) > 80):
                 expanded = []
-                for data in angle2xyz(azi, zen):
+                cart = convert_from_angles(azi, zen, lon, lat)
+                for data in cart:
                     expanded.append(expand_array(
                         data, self.scans, c_align, c_exp, self.scan_size,
                         tpz_size, nb_tpz, self.track_offset, self.scan_offset))
 
-                azi, zen = xyz2angle(*expanded)
+                azi, zen = convert_to_angles(*expanded, lon=self.lons, lat=self.lats)
                 res.append((azi, zen))
             else:
                 expanded = []
@@ -363,6 +366,32 @@ def angles(self, azi_name, zen_name):
         return da.hstack(azi), da.hstack(zen)
 
 
+def convert_from_angles(azi, zen, lon, lat):
+    """Convert the angles to cartesian coordinates."""
+    lon = np.deg2rad(lon)
+    lat = np.deg2rad(lat)
+    x, y, z, = angle2xyz(azi, zen)
+    # Conversion to ECEF is recommended by the provider, but no significant
+    # difference has been seen.
+    # x, y, z = (-np.sin(lon) * x + np.cos(lon) * y,
+    #            -np.sin(lat) * np.cos(lon) * x - np.sin(lat) * np.sin(lon) * y + np.cos(lat) * z,
+    #            np.cos(lat) * np.cos(lon) * x + np.cos(lat) * np.sin(lon) * y + np.sin(lat) * z)
+    return x, y, z
+
+
+def convert_to_angles(x, y, z, lon, lat):
+    """Convert the cartesian coordinates to angles."""
+    lon = np.deg2rad(lon)
+    lat = np.deg2rad(lat)
+    # Conversion to ECEF is recommended by the provider, but no significant
+    # difference has been seen.
+    # x, y, z = (-np.sin(lon) * x - np.sin(lat) * np.cos(lon) * y + np.cos(lat) * np.cos(lon) * z,
+    #            np.cos(lon) * x - np.sin(lat) * np.sin(lon) * y + np.cos(lat) * np.sin(lon) * z,
+    #            np.cos(lat) * y + np.sin(lat) * z)
+    azi, zen = xyz2angle(x, y, z, acos=True)
+    return azi, zen
+
+
 def expand_array(data,
                  scans,
                  c_align,
@@ -391,7 +420,8 @@ def expand_array(data,
     data_c = data[1:scans * 2:2, np.newaxis, 1:, np.newaxis]
     data_d = data[1:scans * 2:2, np.newaxis, :-1, np.newaxis]
 
-    fdata = ((1 - a_track) *
-             ((1 - a_scan) * data_a + a_scan * data_b) + a_track * (
-                 (1 - a_scan) * data_d + a_scan * data_c))
+    fdata = ((1 - a_track)
+             * ((1 - a_scan) * data_a + a_scan * data_b)
+             + a_track
+             * ((1 - a_scan) * data_d + a_scan * data_c))
     return fdata.reshape(scans * scan_size, nties * tpz_size)

satpy/tests/test_utils.py

@@ -136,6 +136,7 @@ def test_angle2xyz(self):
         self.assertAlmostEqual(z__, sqrt(1) / 2)
 
     def test_xyz2lonlat(self):
+        """Test xyz2lonlat."""
         lon, lat = xyz2lonlat(1, 0, 0)
         self.assertAlmostEqual(lon, 0)
         self.assertAlmostEqual(lat, 0)
@@ -144,6 +145,10 @@ def test_xyz2lonlat(self):
         self.assertAlmostEqual(lon, 90)
         self.assertAlmostEqual(lat, 0)
 
+        lon, lat = xyz2lonlat(0, 0, 1, asin=True)
+        self.assertAlmostEqual(lon, 0)
+        self.assertAlmostEqual(lat, 90)
+
         lon, lat = xyz2lonlat(0, 0, 1)
         self.assertAlmostEqual(lon, 0)
         self.assertAlmostEqual(lat, 90)
@@ -153,6 +158,7 @@ def test_xyz2lonlat(self):
         self.assertAlmostEqual(lat, 0)
 
     def test_xyz2angle(self):
+        """Test xyz2angle."""
         azi, zen = xyz2angle(1, 0, 0)
         self.assertAlmostEqual(azi, 90)
         self.assertAlmostEqual(zen, 90)
@@ -165,6 +171,10 @@ def test_xyz2angle(self):
         self.assertAlmostEqual(azi, 0)
         self.assertAlmostEqual(zen, 0)
 
+        azi, zen = xyz2angle(0, 0, 1, acos=True)
+        self.assertAlmostEqual(azi, 0)
+        self.assertAlmostEqual(zen, 0)
+
         azi, zen = xyz2angle(sqrt(2) / 2, sqrt(2) / 2, 0)
         self.assertAlmostEqual(azi, 45)
         self.assertAlmostEqual(zen, 90)
@@ -178,6 +188,7 @@ def test_xyz2angle(self):
         self.assertAlmostEqual(zen, 90)
 
     def test_proj_units_to_meters(self):
+        """Test proj units to meters conversion."""
         prj = '+asd=123123123123'
         res = proj_units_to_meters(prj)
         self.assertEqual(res, prj)
@@ -196,6 +207,7 @@ def test_proj_units_to_meters(self):
 
     @mock.patch('satpy.utils.warnings.warn')
     def test_get_satpos(self, warn_mock):
+        """Test getting the satellite position."""
         orb_params = {'nadir_longitude': 1,
                       'satellite_actual_longitude': 1.1,
                       'satellite_nominal_longitude': 1.2,
@@ -244,8 +256,7 @@ def test_get_satpos(self, warn_mock):
 
 
 def suite():
-    """The test suite.
-    """
+    """Test suite."""
     loader = unittest.TestLoader()
     mysuite = unittest.TestSuite()
     mysuite.addTest(loader.loadTestsFromTestCase(TestUtils))

satpy/utils.py

@@ -16,8 +16,7 @@
 #
 # You should have received a copy of the GNU General Public License
 # along with satpy.  If not, see <http://www.gnu.org/licenses/>.
-"""Module defining various utilities.
-"""
+"""Module defining various utilities."""
 
 import logging
 import os
@@ -37,21 +36,21 @@
 
 
 class OrderedConfigParser(object):
-
     """Intercepts read and stores ordered section names.
+
     Cannot use inheritance and super as ConfigParser use old style classes.
     """
 
     def __init__(self, *args, **kwargs):
+        """Initialize the instance."""
         self.config_parser = configparser.ConfigParser(*args, **kwargs)
 
     def __getattr__(self, name):
+        """Get the attribute."""
         return getattr(self.config_parser, name)
 
     def read(self, filename):
-        """Reads config file
-        """
-
+        """Read config file."""
         try:
             conf_file = open(filename, 'r')
             config = conf_file.read()
@@ -65,17 +64,15 @@ def read(self, filename):
         return self.config_parser.read(filename)
 
     def sections(self):
-        """Get sections from config file
-        """
-
+        """Get sections from config file."""
         try:
             return self.section_keys
         except:  # noqa: E722
             return self.config_parser.sections()
 
 
 def ensure_dir(filename):
-    """Checks if the dir of f exists, otherwise create it."""
+    """Check if the dir of f exists, otherwise create it."""
     directory = os.path.dirname(filename)
     if directory and not os.path.isdir(directory):
         os.makedirs(directory)
@@ -92,8 +89,7 @@ def trace_on():
 
 
 def logging_on(level=logging.WARNING):
-    """Turn logging on.
-    """
+    """Turn logging on."""
     global _is_logging_on
 
     if not _is_logging_on:
@@ -112,8 +108,7 @@ def logging_on(level=logging.WARNING):
 
 
 def logging_off():
-    """Turn logging off.
-    """
+    """Turn logging off."""
     logging.getLogger('').handlers = [logging.NullHandler()]
 
 
@@ -136,7 +131,7 @@ def trace(self, message, *args, **kwargs):
 
 
 def in_ipynb():
-    """Are we in a jupyter notebook?"""
+    """Check if we are in a jupyter notebook."""
     try:
         return 'ZMQ' in get_ipython().__class__.__name__
     except NameError:
@@ -156,10 +151,13 @@ def lonlat2xyz(lon, lat):
     return x, y, z
 
 
-def xyz2lonlat(x, y, z):
+def xyz2lonlat(x, y, z, asin=False):
     """Convert cartesian to lon lat."""
     lon = np.rad2deg(np.arctan2(y, x))
-    lat = np.rad2deg(np.arctan2(z, np.sqrt(x ** 2 + y ** 2)))
+    if asin:
+        lat = np.rad2deg(np.arcsin(z))
+    else:
+        lat = np.rad2deg(np.arctan2(z, np.sqrt(x ** 2 + y ** 2)))
     return lon, lat
 
 
@@ -173,10 +171,13 @@ def angle2xyz(azi, zen):
     return x, y, z
 
 
-def xyz2angle(x, y, z):
+def xyz2angle(x, y, z, acos=False):
     """Convert cartesian to azimuth and zenith."""
     azi = np.rad2deg(np.arctan2(x, y))
-    zen = 90 - np.rad2deg(np.arctan2(z, np.sqrt(x ** 2 + y ** 2)))
+    if acos:
+        zen = np.rad2deg(np.arccos(z))
+    else:
+        zen = 90 - np.rad2deg(np.arctan2(z, np.sqrt(x ** 2 + y ** 2)))
     return azi, zen
 
 
@@ -217,7 +218,6 @@ def sunzen_corr_cos(data, cos_zen, limit=88., max_sza=95.):
     0. Both ``data`` and ``cos_zen`` should be 2D arrays of the same shape.
 
     """
-
     # Convert the zenith angle limit to cosine of zenith angle
     limit_rad = np.deg2rad(limit)
     limit_cos = np.cos(limit_rad)
@@ -295,6 +295,7 @@ def get_satpos(dataset):
 
     Returns:
         Geodetic longitude, latitude, altitude
+
     """
     try:
         orb_params = dataset.attrs['orbital_parameters']