Merge pull request #1918 from ninahakansson/fix_latlon_interp_aapp

Fix geo interpolation for aapp data
pytroll · Dec 14, 2021 · a2b1e82 · a2b1e82
2 parents 5ec9606 + ab4b750
commit a2b1e82
Show file tree

Hide file tree

Showing 2 changed files with 135 additions and 5 deletions.
diff --git a/satpy/readers/aapp_l1b.py b/satpy/readers/aapp_l1b.py
@@ -233,10 +233,10 @@ def _get_tiepoint_angles_in_degrees(self):
         azidiff40km = self._data["ang"][:, :, 2] * 1e-2
         return sunz40km, satz40km, azidiff40km
 
-    def _interpolate_arrays(self, *input_arrays):
+    def _interpolate_arrays(self, *input_arrays, geolocation=False):
         lines = input_arrays[0].shape[0]
         try:
-            interpolator = self._create_40km_interpolator(lines, *input_arrays)
+            interpolator = self._create_40km_interpolator(lines, *input_arrays, geolocation=geolocation)
         except ImportError:
             logger.warning("Could not interpolate, python-geotiepoints missing.")
             output_arrays = input_arrays
@@ -247,8 +247,12 @@ def _interpolate_arrays(self, *input_arrays):
         return output_arrays
 
     @staticmethod
-    def _create_40km_interpolator(lines, *arrays_40km):
-        from geotiepoints.interpolator import Interpolator
+    def _create_40km_interpolator(lines, *arrays_40km, geolocation=False):
+        if geolocation:
+            # Slower but accurate at datum line
+            from geotiepoints.geointerpolator import GeoInterpolator as Interpolator
+        else:
+            from geotiepoints.interpolator import Interpolator
         cols40km = np.arange(24, 2048, 40)
         cols1km = np.arange(2048)
         rows40km = np.arange(lines)
@@ -273,7 +277,7 @@ def navigate(self, coordinate_id):
     @functools.lru_cache(maxsize=10)
     def _get_all_interpolated_coordinates(self):
         lons40km, lats40km = self._get_coordinates_in_degrees()
-        return self._interpolate_arrays(lons40km, lats40km)
+        return self._interpolate_arrays(lons40km, lats40km, geolocation=True)
 
     def _get_coordinates_in_degrees(self):
         lons40km = self._data["pos"][:, :, 1] * 1e-4

diff --git a/satpy/tests/reader_tests/test_aapp_l1b.py b/satpy/tests/reader_tests/test_aapp_l1b.py
@@ -22,6 +22,7 @@
 import tempfile
 import unittest
 from contextlib import suppress
+from unittest import mock
 
 import numpy as np
 
@@ -153,6 +154,131 @@ def test_navigation(self):
             res = fh.get_dataset(key, info)
             assert(np.all(res == 0))
 
+    def test_interpolation(self):
+        """Test reading the lon and lats."""
+        with tempfile.TemporaryFile() as tmpfile:
+            self._header.tofile(tmpfile)
+            tmpfile.seek(22016, 0)
+            self._data.tofile(tmpfile)
+            fh = AVHRRAAPPL1BFile(tmpfile, self.filename_info, self.filetype_info)
+            lons40km = np.array([
+                [-115.9773, -122.3054, -127.7482, -132.464, -136.5788, -140.1951,
+                 -143.3961, -146.2497, -148.8112, -151.1259, -153.2309, -155.1568,
+                 -156.9291, -158.5689, -160.0941, -161.5196, -162.8584, -164.1212,
+                 -165.3176, -166.4557, -167.5426, -168.5846, -169.5872, -170.5555,
+                 -171.4937, -172.406, -173.296, -174.1671, -175.0224, -175.865,
+                 -176.6976, -177.523, -178.3439, -179.1628, -179.9825, 179.1944,
+                 178.3651, 177.5267, 176.6761, 175.8098, 174.9242, 174.0149,
+                 173.0773, 172.1057, 171.0935, 170.0326, 168.9128, 167.7211,
+                 166.4397, 165.0436, 163.4946],
+                [-115.9639, -122.2967, -127.7441, -132.4639, -136.5824, -140.2018,
+                 -143.4055, -146.2614, -148.8249, -151.1413, -153.2478, -155.175,
+                 -156.9484, -158.5892, -160.1152, -161.5415, -162.8809, -164.1443,
+                 -165.3412, -166.4797, -167.567, -168.6094, -169.6123, -170.5808,
+                 -171.5192, -172.4317, -173.3219, -174.1931, -175.0486, -175.8913,
+                 -176.724, -177.5494, -178.3703, -179.1893, 179.991, 179.168,
+                 178.3388, 177.5005, 176.6499, 175.7838, 174.8983, 173.9892,
+                 173.0518, 172.0805, 171.0685, 170.0079, 168.8885, 167.6972,
+                 166.4164, 165.0209, 163.4726],
+                [-115.9504, -122.288, -127.7399, -132.4639, -136.5859, -140.2084,
+                 -143.4148, -146.2731, -148.8386, -151.1567, -153.2647, -155.1932,
+                 -156.9677, -158.6095, -160.1363, -161.5634, -162.9034, -164.1674,
+                 -165.3648, -166.5038, -167.5915, -168.6341, -169.6374, -170.6061,
+                 -171.5448, -172.4575, -173.3478, -174.2192, -175.0748, -175.9176,
+                 -176.7503, -177.5758, -178.3968, -179.2157, 179.9646, 179.1416,
+                 178.3124, 177.4742, 176.6238, 175.7577, 174.8724, 173.9635,
+                 173.0263, 172.0552, 171.0436, 169.9833, 168.8643, 167.6734,
+                 166.3931, 164.9982, 163.4507]])
+            lats40km = np.array([
+                [78.6613, 78.9471, 79.0802, 79.1163, 79.0889, 79.019, 78.9202,
+                 78.8016, 78.6695, 78.528, 78.38, 78.2276, 78.0721, 77.9145,
+                 77.7553, 77.5949, 77.4335, 77.2712, 77.1079, 76.9435, 76.7779,
+                 76.6108, 76.4419, 76.2708, 76.0973, 75.921, 75.7412, 75.5576,
+                 75.3696, 75.1764, 74.9776, 74.7721, 74.5592, 74.3379, 74.1069,
+                 73.865, 73.6106, 73.342, 73.057, 72.7531, 72.4273, 72.076,
+                 71.6945, 71.2773, 70.8171, 70.3046, 69.7272, 69.0676, 68.3014,
+                 67.3914, 66.2778],
+                [78.6703, 78.9565, 79.0897, 79.1259, 79.0985, 79.0286, 78.9297,
+                 78.8111, 78.6789, 78.5373, 78.3892, 78.2367, 78.0811, 77.9233,
+                 77.764, 77.6035, 77.442, 77.2796, 77.1162, 76.9518, 76.7861,
+                 76.6188, 76.4498, 76.2787, 76.1051, 75.9287, 75.7488, 75.5651,
+                 75.377, 75.1838, 74.9848, 74.7793, 74.5663, 74.3448, 74.1138,
+                 73.8718, 73.6173, 73.3486, 73.0635, 72.7595, 72.4336, 72.0821,
+                 71.7005, 71.2832, 70.8229, 70.3102, 69.7326, 69.0729, 68.3065,
+                 67.3963, 66.2825],
+                [78.6794, 78.9658, 79.0993, 79.1355, 79.1082, 79.0381, 78.9392,
+                 78.8205, 78.6882, 78.5465, 78.3984, 78.2458, 78.0901, 77.9322,
+                 77.7728, 77.6122, 77.4506, 77.2881, 77.1246, 76.96, 76.7942,
+                 76.6269, 76.4578, 76.2866, 76.1129, 75.9364, 75.7564, 75.5727,
+                 75.3844, 75.1911, 74.9921, 74.7864, 74.5734, 74.3518, 74.1207,
+                 73.8786, 73.624, 73.3552, 73.0699, 72.7658, 72.4398, 72.0882,
+                 71.7065, 71.2891, 70.8286, 70.3158, 69.7381, 69.0782, 68.3116,
+                 67.4012, 66.2872]])
+            fh._get_coordinates_in_degrees = mock.MagicMock()
+            fh._get_coordinates_in_degrees.return_value = (lons40km, lats40km)
+            (lons, lats) = fh._get_all_interpolated_coordinates()
+            lon_data = lons.compute()
+            self.assertTrue(np.max(lon_data) <= 180)
+            # Not longitdes between -110, 110 in indata
+            self.assertTrue(np.all(np.abs(lon_data) > 110))
+
+    def test_interpolation_angles(self):
+        """Test reading the lon and lats."""
+        with tempfile.TemporaryFile() as tmpfile:
+            self._header.tofile(tmpfile)
+            tmpfile.seek(22016, 0)
+            self._data.tofile(tmpfile)
+            fh = AVHRRAAPPL1BFile(tmpfile, self.filename_info, self.filetype_info)
+
+            sunz40km = np.array(
+                [[122.42, 121.72, 121.14, 120.63, 120.19, 119.79, 119.43, 119.1, 118.79, 118.51,
+                  118.24, 117.99, 117.76, 117.53, 117.31, 117.1, 116.9, 116.71, 116.52, 116.33,
+                  116.15, 115.97, 115.79, 115.61, 115.44, 115.26, 115.08, 114.91, 114.73, 114.55,
+                  114.36, 114.18, 113.98, 113.79, 113.58, 113.37, 113.15, 112.92, 112.68, 112.43,
+                  112.15, 111.87, 111.55, 111.22, 110.85, 110.44, 109.99, 109.47, 108.88, 108.18,
+                  107.33],
+                 [122.41, 121.71, 121.13, 120.62, 120.18, 119.78, 119.42, 119.09, 118.78, 118.5,
+                  118.24, 117.99, 117.75, 117.52, 117.31, 117.1, 116.9, 116.7, 116.51, 116.32,
+                  116.14, 115.96, 115.78, 115.6, 115.43, 115.25, 115.08, 114.9, 114.72, 114.54,
+                  114.36, 114.17, 113.98, 113.78, 113.57, 113.36, 113.14, 112.91, 112.67, 112.42,
+                  112.15, 111.86, 111.55, 111.21, 110.84, 110.43, 109.98, 109.46, 108.87, 108.17,
+                  107.32]])
+            satz40km = np.array(
+                [[6.623e+01, 6.281e+01, 5.960e+01, 5.655e+01, 5.360e+01, 5.075e+01, 4.797e+01,
+                  4.524e+01, 4.256e+01, 3.992e+01, 3.731e+01, 3.472e+01, 3.216e+01, 2.962e+01,
+                  2.710e+01, 2.460e+01, 2.210e+01, 1.962e+01, 1.714e+01, 1.467e+01, 1.221e+01,
+                  9.760e+00, 7.310e+00, 4.860e+00, 2.410e+00, 3.000e-02, 2.470e+00, 4.920e+00,
+                  7.370e+00, 9.820e+00, 1.227e+01, 1.474e+01, 1.720e+01, 1.968e+01, 2.216e+01,
+                  2.466e+01, 2.717e+01, 2.969e+01, 3.223e+01, 3.479e+01, 3.737e+01, 3.998e+01,
+                  4.263e+01, 4.531e+01, 4.804e+01, 5.082e+01, 5.368e+01, 5.662e+01, 5.969e+01,
+                  6.290e+01, 6.633e+01],
+                 [6.623e+01, 6.281e+01, 5.960e+01, 5.655e+01, 5.360e+01, 5.075e+01, 4.797e+01,
+                  4.524e+01, 4.256e+01, 3.992e+01, 3.731e+01, 3.472e+01, 3.216e+01, 2.962e+01,
+                  2.710e+01, 2.460e+01, 2.210e+01, 1.962e+01, 1.714e+01, 1.467e+01, 1.221e+01,
+                  9.760e+00, 7.310e+00, 4.860e+00, 2.410e+00, 3.000e-02, 2.470e+00, 4.920e+00,
+                  7.370e+00, 9.820e+00, 1.227e+01, 1.474e+01, 1.720e+01, 1.968e+01, 2.216e+01,
+                  2.466e+01, 2.717e+01, 2.969e+01, 3.223e+01, 3.479e+01, 3.737e+01, 3.998e+01,
+                  4.263e+01, 4.531e+01, 4.804e+01, 5.082e+01, 5.368e+01, 5.662e+01, 5.969e+01,
+                  6.290e+01, 6.633e+01]])
+            azidiff40km = np.array([
+                [56.9, 56.24, 55.71, 55.27, 54.9, 54.57, 54.29, 54.03, 53.8, 53.59,
+                 53.4, 53.22, 53.05, 52.89, 52.74, 52.6, 52.47, 52.34, 52.22, 52.1,
+                 51.98, 51.87, 51.76, 51.65, 51.55, 128.55, 128.65, 128.76, 128.86, 128.96,
+                 129.07, 129.17, 129.27, 129.38, 129.49, 129.6, 129.72, 129.83, 129.95, 130.08,
+                 130.21, 130.35, 130.5, 130.65, 130.81, 130.99, 131.18, 131.39, 131.63, 131.89,
+                 132.19],
+                [56.9, 56.24, 55.72, 55.28, 54.9, 54.58, 54.29, 54.03, 53.8, 53.59,
+                 53.4, 53.22, 53.05, 52.89, 52.75, 52.6, 52.47, 52.34, 52.22, 52.1,
+                 51.98, 51.87, 51.76, 51.65, 51.55, 128.55, 128.65, 128.75, 128.86, 128.96,
+                 129.06, 129.17, 129.27, 129.38, 129.49, 129.6, 129.71, 129.83, 129.95, 130.08,
+                 130.21, 130.35, 130.49, 130.65, 130.81, 130.99, 131.18, 131.39, 131.62, 131.89,
+                 132.19]])
+            fh._get_tiepoint_angles_in_degrees = mock.MagicMock()
+            fh._get_tiepoint_angles_in_degrees.return_value = (sunz40km, satz40km, azidiff40km)
+            (sunz, satz, azidiff) = fh._get_all_interpolated_angles()
+            self.assertTrue(np.max(sunz) <= 123)
+            self.assertTrue(np.max(satz) <= 70)
+
 
 class TestAAPPL1BChannel3AMissing(unittest.TestCase):
     """Test the filehandler when channel 3a is missing."""