Merge pull request #34 from pepephillips/tiepoints_dev

Updated interpolator for vii tie points for test data version V2
pytroll · Feb 18, 2022 · 2f569cc · 2f569cc
2 parents 7c14ae9 + bfe48f4
commit 2f569cc
Show file tree

Hide file tree

Showing 3 changed files with 84 additions and 83 deletions.
diff --git a/AUTHORS.md b/AUTHORS.md
@@ -14,4 +14,5 @@ The following people have made contributions to this project:
 - [Mikhail Itkin (mitkin)](https://github.com/mitkin)
 - [Sauli Joro (sjoro)](https://github.com/sjoro)
 - [Panu Lahtinen (pnuu)](https://github.com/pnuu)
+- [Pepe Phillips (pepephillips)](https://github.com/pepephillips)
 - [Martin Raspaud (mraspaud)](https://github.com/mraspaud)
diff --git a/geotiepoints/tests/test_viiinterpolator.py b/geotiepoints/tests/test_viiinterpolator.py
@@ -17,7 +17,9 @@
 # python-geotiepoints.  If not, see <http://www.gnu.org/licenses/>.
 """Test of the interpolation of geographical tiepoints for the VII products.
 
-It follows the description provided in document "EPS-SG VII Level 1B Product Format Specification".
+It follows the description provided in document "EPS-SG VII Level 1B Product Format Specification V4A".
+This version is compatible for vii (METimage) test data version V2 (Jan 2022). It is not back compatible
+with V1.
 
 """
 
@@ -37,58 +39,70 @@
 
 # Results of latitude/longitude interpolation with simple interpolation on coordinates
 TEST_LON_1 = np.array(
-    [[-12., -11.5, -11., -10.5, -9., -8.5, -8., -7.5],
-     [-9., -8.5, -8., -7.5, -6., -5.5, -5., -4.5],
-     [-6., -5.5, -5., -4.5, -3., -2.5, -2., -1.5],
-     [-3., -2.5, -2., -1.5, 0., 0.5, 1., 1.5],
-     [0., 0.5, 1., 1.5, 3., 3.5, 4., 4.5],
-     [3., 3.5, 4., 4.5, 6., 6.5, 7., 7.5]]
+    [[-12., -11.5, -11., -10.5, -10.0, -9.5],
+     [-10., -9.5, -9., -8.5, -8.0, -7.5],
+     [-8., -7.5, -7., -6.5, -6.0, -5.5],
+     [-6., -5.5, -5., -4.5, -4.0, -3.5],
+     [0., 0.5, 1., 1.5, 2.0, 2.5],
+     [2., 2.5, 3., 3.5, 4.0, 4.5],
+     [4., 4.5, 5., 5.5, 6.0, 6.5],
+     [6., 6.5, 7., 7.5, 8.0, 8.5]]
 )
 TEST_LAT_1 = np.array(
-    [[0., 0.5, 1., 1.5, 3., 3.5, 4., 4.5],
-     [3., 3.5, 4., 4.5, 6., 6.5, 7., 7.5],
-     [6., 6.5, 7., 7.5, 9., 9.5, 10., 10.5],
-     [9., 9.5, 10., 10.5, 12., 12.5, 13., 13.5],
-     [12., 12.5, 13., 13.5, 15., 15.5, 16., 16.5],
-     [15., 15.5, 16., 16.5, 18., 18.5, 19., 19.5]]
+    [[0., 0.5, 1., 1.5, 2., 2.5],
+     [2., 2.5, 3., 3.5, 4., 4.5],
+     [4., 4.5, 5., 5.5, 6., 6.5],
+     [6., 6.5, 7., 7.5, 8., 8.5],
+     [12., 12.5, 13., 13.5, 14.0, 14.5],
+     [14., 14.5, 15., 15.5, 16., 16.5],
+     [16., 16.5, 17., 17.5, 18., 18.5],
+     [18., 18.5, 19., 19.5, 20., 20.5]]
 )
 
 # Results of latitude/longitude interpolation on cartesian coordinates (latitude above 60 degrees)
 TEST_LON_2 = np.array(
-    [[-12., -11.50003808, -11., -10.50011426, -9., -8.50026689, -8., -7.50034342],
-     [-9.00824726, -8.50989187, -8.01100418, -7.51272848, -6.01653996, -5.51842688, -5.01932226, -4.52129225],
-     [-6., -5.50049716, -5., -4.50057447, -3., -2.50073021, -2., -1.50080874],
-     [-3.02492451, -2.52706443, -2.02774808, -1.52997501, -0.03344942, 0.46414517, 0.96366893, 1.4611719],
-     [0., 0.49903263, 1., 1.49895241, 3., 3.49878988, 4., 4.49870746],
-     [2.9578336, 3.45514812, 3.9548757, 4.4520932, 5.94886832, 6.44588569, 6.94581415, 7.44272818]]
+    [[-12., -11.50003808, -11., -10.50011426, -10., -9.50019052],
+     [-10.00243991, -9.5032411, -9.00366173, -8.50454031, -8.00488578, -7.5058423],
+     [-8., -7.50034342, -7., -6.50042016, -6., -5.50049716],
+     [-6.00734362, -5.50845783, -5.00857895, -4.50977302, -4.00981958, -3.51109426],
+     [0., 0.49903263, 1., 1.49895241, 2., 2.49887151],
+     [1.98257947, 2.48080192, 2.98127841, 3.47941324, 3.97996512, 4.47801105],
+     [4., 4.49870746, 5., 5.49862418, 6., 6.49853998],
+     [5.97729789, 6.47516183, 6.97594186, 7.47371256, 7.97456943, 8.47224525]]
 )
 
 TEST_LAT_2 = np.array(
-    [[0., 0.49998096, 1., 1.49994287, 3., 3.49986656, 4., 4.4998283],
-     [2.99588485, 3.49506416, 3.99450923, 4.49364876, 5.99174708, 6.49080542, 6.99035883, 7.4893757],
-     [6., 6.49975143, 7., 7.49971278, 9., 9.49963492, 10., 10.49959566],
-     [8.98756357, 9.48649563, 9.98615477, 10.4850434, 11.98331018, 12.48210974, 12.98187246, 13.4806263],
-     [12., 12.49951634, 13., 13.49947623, 15., 15.49939498, 16., 16.49935377],
-     [14.97896116, 15.47762097, 15.97748548, 16.47609689, 17.97448854, 18.47300011, 18.97296495, 19.47142496]]
+    [[0., 0.49998096, 1., 1.49994287, 2., 2.49990475],
+     [1.99878116, 2.49838091, 2.99817081, 3.49773189, 3.99755935, 4.49708148],
+     [4., 4.4998283, 5., 5.49978993, 6., 6.49975143],
+     [5.99633155, 6.4957749, 6.99571447, 7.49511791, 7.99509473, 8.49445789],
+     [12., 12.49951634, 13., 13.49947623, 14., 14.499435779],
+     [13.99129786, 14.4904098, 14.99064796, 15.48971613, 15.98999196, 16.48901572],
+     [16., 16.49935377, 17., 17.49931213, 18., 18.49927003],
+     [17.98865968, 18.48759253, 18.98798235, 19.48686863, 19.98729684, 20.48613573]]
 )
 
 # Results of latitude/longitude interpolation on cartesian coordinates (longitude with a 360 degrees step)
 TEST_LON_3 = np.array(
-    [[-12., -11.50444038, -11., -10.50459822, -9., -8.50493209, -8., -7.50510905],
-     [-9.17477341, -8.68280267, -8.18102962, -7.68936161, -6.19433153, -5.70332248, -5.20141997, -4.71077058],
-     [-6., -5.50548573, -5., -4.50568668, -3., -2.50611746, -2., -1.506349],
-     [-3.2165963, -2.72673687, -2.22474246, -1.73531828, -0.24232275, 0.24613534, 0.7481615, 1.23608137],
-     [0., 0.49315061, 1., 1.49287934, 3., 3.49228746, 4., 4.49196335],
-     [2.72743411, 3.21414435, 3.71610443, 4.20213182, 5.69115252, 6.17562189, 6.67735289, 7.16092853]]
+    [[-12., -11.50444038, -11., -10.50459822, -10., -9.50476197],
+     [-10.07492627, -9.58101155, -9.07759836, -8.5839056, - 8.0803761, -7.58691614],
+     [-8., -7.50510905, -7., -6.5052934, -6., -5.50548573],
+     [-6.0862821, -5.59332416, -5.08942935, -4.59674283, -4.09272043, -3.60032066],
+     [0., 0.49315061, 1., 1.49287934, 2., 2.49259217],
+     [1.88371709, 2.3739768, 2.87898193, 3.36879304, 3.87395153, 4.36327924],
+     [4., 4.49196335, 5., 5.49161771, 6., 6.49124808],
+     [5.86287282, 6.35111105, 6.85674573, 7.3443667, 7.85016382, 8.33710998]]
 )
 
 TEST_LAT_3 = np.array(
-    [[45., 45.49777998, 46., 46.49770107, 48., 48.49753416, 49., 49.49744569],
-     [47.91286617, 48.40886652, 48.90975264, 49.40560282, 50.90313463, 51.39865815, 51.8996091, 52.39495445],
-     [51., 51.49725738, 52., 52.49715691, 54., 54.50311196, 55., 55.50299846],
-     [54.11364234, 54.61463583, 55.10950687, 55.61043728, 57.10152529, 57.60232834, 58.09766771, 58.59840655],
-     [57., 57.50277937, 58., 58.50267347, 60., 60.50246826, 61., 61.5023687],
-     [60.09019289, 60.59080228, 61.0865661, 61.58711028, 63.07951189, 63.57992468, 64.07607638, 64.57642295]]
+    [[45., 45.49777998, 46., 46.49770107, 47., 47.4976192],
+     [46.96258417, 47.4595462, 47.9612508, 48.4581022, 48.95986481, 49.45660021],
+     [49., 49.49744569, 50., 50.49735352, 51., 51.49725738],
+     [50.95691833, 51.45340364, 51.95534841, 52.45169855, 52.95370691, 53.55477452],
+     [57., 57.50277937, 58., 58.50267347, 59., 59.50256981],
+     [59.04185095, 59.54357898, 60.04020844, 60.54185139, 61.03859839, 61.54015723],
+     [61., 61.5023687, 62., 62.502271, 63., 63.50217506],
+     [63.03546804, 63.53686131, 64.03394419, 64.53525587, 65.03244567, 65.53367647]]
 )
 
 
@@ -103,16 +117,16 @@ def setUp(self):
             np.arange(
                 TEST_VALID_ALT_TIE_POINTS * TEST_ACT_TIE_POINTS,
                 dtype=np.float64,
-            ).reshape(TEST_ACT_TIE_POINTS, TEST_VALID_ALT_TIE_POINTS),
-            dims=('num_tie_points_act', 'num_tie_points_alt'),
+            ).reshape(TEST_VALID_ALT_TIE_POINTS, TEST_ACT_TIE_POINTS),
+            dims=('num_tie_points_alt', 'num_tie_points_act'),
         )
         # The second has an invalid number of n_tie_alt points (not multiple of SCAN_ALT_TIE_POINTS)
         self.invalid_data_for_interpolation = xr.DataArray(
             np.arange(
                 TEST_INVALID_ALT_TIE_POINTS * TEST_ACT_TIE_POINTS,
                 dtype=np.float64,
-            ).reshape(TEST_ACT_TIE_POINTS, TEST_INVALID_ALT_TIE_POINTS),
-            dims=('num_tie_points_act', 'num_tie_points_alt'),
+            ).reshape(TEST_INVALID_ALT_TIE_POINTS, TEST_ACT_TIE_POINTS),
+            dims=('num_tie_points_alt', 'num_tie_points_act'),
         )
         # Then two arrays containing valid longitude and latitude data
         self.longitude = xr.DataArray(
@@ -121,17 +135,17 @@ def setUp(self):
                 11,
                 num=TEST_VALID_ALT_TIE_POINTS * TEST_ACT_TIE_POINTS,
                 dtype=np.float64,
-            ).reshape(TEST_ACT_TIE_POINTS, TEST_VALID_ALT_TIE_POINTS),
-            dims=('num_tie_points_act', 'num_tie_points_alt'),
+            ).reshape(TEST_VALID_ALT_TIE_POINTS, TEST_ACT_TIE_POINTS),
+            dims=('num_tie_points_alt', 'num_tie_points_act'),
         )
         self.latitude = xr.DataArray(
             np.linspace(
                 0,
                 23,
                 num=TEST_VALID_ALT_TIE_POINTS * TEST_ACT_TIE_POINTS,
                 dtype=np.float64,
-            ).reshape(TEST_ACT_TIE_POINTS, TEST_VALID_ALT_TIE_POINTS),
-            dims=('num_tie_points_act', 'num_tie_points_alt'),
+            ).reshape(TEST_VALID_ALT_TIE_POINTS, TEST_ACT_TIE_POINTS),
+            dims=('num_tie_points_alt', 'num_tie_points_act'),
         )
         # Then one containing latitude data above 60 degrees
         self.latitude_over60 = xr.DataArray(
@@ -140,8 +154,8 @@ def setUp(self):
                 68,
                 num=TEST_VALID_ALT_TIE_POINTS * TEST_ACT_TIE_POINTS,
                 dtype=np.float64,
-            ).reshape(TEST_ACT_TIE_POINTS, TEST_VALID_ALT_TIE_POINTS),
-            dims=('num_tie_points_act', 'num_tie_points_alt'),
+            ).reshape(TEST_VALID_ALT_TIE_POINTS, TEST_ACT_TIE_POINTS),
+            dims=('num_tie_points_alt', 'num_tie_points_act'),
         )
         # Then one containing longitude data with a 360 degrees step
         self.longitude_over360 = xr.DataArray(
@@ -150,8 +164,8 @@ def setUp(self):
                 11,
                 num=TEST_VALID_ALT_TIE_POINTS * TEST_ACT_TIE_POINTS,
                 dtype=np.float64,
-            ).reshape(TEST_ACT_TIE_POINTS, TEST_VALID_ALT_TIE_POINTS) % 360.,
-            dims=('num_tie_points_act', 'num_tie_points_alt'),
+            ).reshape(TEST_VALID_ALT_TIE_POINTS, TEST_ACT_TIE_POINTS) % 360.,
+            dims=('num_tie_points_alt', 'num_tie_points_act'),
         )
 
     def tearDown(self):
@@ -172,26 +186,12 @@ def test_tie_points_interpolation(self):
         num_scans = TEST_VALID_ALT_TIE_POINTS // TEST_SCAN_ALT_TIE_POINTS
         scan_alt_points_interp = (TEST_SCAN_ALT_TIE_POINTS - 1) * TEST_TIE_POINTS_FACTOR
 
-        # It is easier to check the delta between interpolated points, which must be 1/8 of the original delta
-        # Across the track, it is possible to check the delta on the entire array
-        delta_axis_0 = 1.0 * TEST_VALID_ALT_TIE_POINTS / TEST_TIE_POINTS_FACTOR
-        self.assertTrue(np.allclose(
-            np.diff(result_valid, axis=0),
-            np.ones((act_points_interp - 1, num_scans * scan_alt_points_interp)) * delta_axis_0
-            ))
-
-        delta_axis_1 = 1.0 / TEST_TIE_POINTS_FACTOR
-        # Along the track, it is necessary to check the delta on each scan separately
-        for i in range(num_scans):
-            first_index = i*(TEST_SCAN_ALT_TIE_POINTS-1)*TEST_TIE_POINTS_FACTOR
-            last_index = (i+1)*(TEST_SCAN_ALT_TIE_POINTS-1)*TEST_TIE_POINTS_FACTOR
-            result_per_scan = result_valid[:, first_index:last_index]
-            self.assertTrue(np.allclose(
-                np.diff(result_per_scan, axis=1),
-                np.ones((act_points_interp, (TEST_SCAN_ALT_TIE_POINTS-1)*TEST_TIE_POINTS_FACTOR - 1)) * delta_axis_1
-                ))
-
-        self.assertEqual(len(result_valid.coords), 0)
+        # Across the track
+        delta_axis_0 = [0., 0.5, 1., 1.5, 2., 2.5]
+        self.assertTrue(np.allclose(result_valid[0, :], delta_axis_0))
+        # Along track
+        delta_axis_1 = [0., 2., 4., 6., 12., 14., 16., 18]
+        self.assertTrue(np.allclose(result_valid[:, 0], delta_axis_1))
 
         # Test the interpolation routine with invalid input
         pytest.raises(ValueError, tie_points_interpolation,

diff --git a/geotiepoints/viiinterpolator.py b/geotiepoints/viiinterpolator.py
@@ -21,10 +21,12 @@
 Tiepoints are typically subsampled by a factor 8 with respect to the pixels, along and across the satellite track.
 Because of the bowtie effect, tiepoints at the scan edges are not continuous between neighbouring scans,
 therefore each scan has its own edge tiepoints in the along-track direction.
-Each scan typically extends on 3 tiepoints in the along-track direction.
-At the edges of a given scan (both along and across track) the tie points lie outside the original data point raster
-and are therefore excluded from the interpolation grid.
+However, for computational efficiency, the edge tie points that lie outside the original data point raster
+are excluded from the interpolation grid which is then carried out per granule, rather than per scan
+at a (very) small geolocation accuracy cost at the swath edge (investigation to quantify this ongoing).
 The interpolation functions are implemented for xarray.DataArrays as input.
+This version works with vii test data V2 to be released Jan 2022 which has the data stored
+in alt, act (row,col) format instead of act,alt (col,row)
 """
 
 import xarray as xr
@@ -37,22 +39,18 @@
 
 def tie_points_interpolation(data_on_tie_points, scan_alt_tie_points, tie_points_factor):
     """Interpolate the data from the tie points to the pixel points.
-
     The data are provided as a list of xarray DataArray objects, allowing to interpolate on several arrays
     at the same time; however the individual arrays must have exactly the same dimensions.
-
     Args:
         data_on_tie_points: list of xarray DataArray objects containing the values defined on the tie points.
         scan_alt_tie_points: number of tie points along the satellite track for each scan
         tie_points_factor: sub-sampling factor of tie points wrt pixel points
-
     Returns:
         list of xarray DataArray objects containing the interpolated values on the pixel points.
-
     """
     # Extract the dimensions of the tie points array across and along track
-    n_tie_act, n_tie_alt = data_on_tie_points[0].shape
-    dim_act, dim_alt = data_on_tie_points[0].dims
+    n_tie_alt, n_tie_act = data_on_tie_points[0].shape
+    dim_alt, dim_act = data_on_tie_points[0].dims
 
     # Check that the number of tie points along track is multiple of the number of tie points per scan
     if n_tie_alt % scan_alt_tie_points != 0:
@@ -68,12 +66,13 @@ def tie_points_interpolation(data_on_tie_points, scan_alt_tie_points, tie_points
 
     # Create the grids used for interpolation across the track
     tie_grid_act = da.arange(0, n_pixel_act + 1, tie_points_factor)
-    pixels_grid_act = da.arange(0, n_pixel_act)
+    pixel_grid_act = da.arange(0, n_pixel_act)
 
     # Create the grids used for the interpolation along the track (must not include the spurious points between scans)
     tie_grid_alt = da.arange(0, n_pixel_alt + 1, tie_points_factor)
     n_pixel_alt_per_scan = (scan_alt_tie_points - 1) * tie_points_factor
     pixel_grid_alt = []
+
     for j_scan in range(n_scans):
         start_index_scan = j_scan * scan_alt_tie_points * tie_points_factor
         pixel_grid_alt.append(da.arange(start_index_scan, start_index_scan + n_pixel_alt_per_scan))
@@ -83,20 +82,21 @@ def tie_points_interpolation(data_on_tie_points, scan_alt_tie_points, tie_points
     data_on_pixel_points = []
     for data in data_on_tie_points:
 
-        if data.shape != (n_tie_act, n_tie_alt) or data.dims != (dim_act, dim_alt):
+        if data.shape != (n_tie_alt, n_tie_act) or data.dims != (dim_alt, dim_act):
             raise ValueError("The dimensions of the arrays are not consistent")
 
         # Interpolate using the xarray interp function twice: first across, then along the scan
         # (much faster than interpolating directly in the two dimensions)
-        data = data.assign_coords({dim_act: tie_grid_act, dim_alt: tie_grid_alt})
-        data_pixel = data.interp({dim_act: pixels_grid_act}, assume_sorted=True) \
-                         .interp({dim_alt: pixel_grid_alt}, assume_sorted=True).drop_vars([dim_act, dim_alt])
+        data = data.assign_coords({dim_alt: tie_grid_alt, dim_act: tie_grid_act})
+        data_pixel = data.interp({dim_alt: pixel_grid_alt}, assume_sorted=True) \
+                         .interp({dim_act: pixel_grid_act}, assume_sorted=True).drop_vars([dim_alt, dim_act])
 
         data_on_pixel_points.append(data_pixel)
 
     return data_on_pixel_points
 
 
+
 def tie_points_geo_interpolation(longitude, latitude,
                                  scan_alt_tie_points, tie_points_factor,
                                  lat_threshold_use_cartesian=60.,