version for test data V2 but without internal loop over scans

geotiepoints/viiinterpolator.py

@@ -19,39 +19,36 @@
 """Interpolation of geographical tiepoints for the VII products.
 It follows the description provided in document "EPS-SG VII Level 1B Product Format Specification".
 Tiepoints are typically subsampled by a factor 8 with respect to the pixels, along and across the satellite track.
-Due to the bowtie effect, the pixel sampling pattern is discontinuous at the swath edge. As a result, each scan has
-the interpolation is more accurate when performed intra-scan and not inter-scan with discontinuous sampling regions.
-This version uses all tie points for individual scans which are then subsequently re-assembled into granules
-before return as per PFS (although at visual inspection level there appears to be little difference in the resulting
-images). It is also modified to work 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) """
+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.
+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
 import dask.array as da
 import numpy as np
-from satpy.dataset import combine_metadata
 
 # MEAN EARTH RADIUS AS DEFINED BY IUGG
 MEAN_EARTH_RADIUS = 6371008.7714  # [m]
 
 
 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_alt, n_tie_act = data_on_tie_points[0].shape
     dim_alt, dim_act = data_on_tie_points[0].dims
 
@@ -65,57 +62,41 @@ def tie_points_interpolation(data_on_tie_points, scan_alt_tie_points, tie_points
 
     # Compute the dimensions of the pixel points array across and along track
     n_pixel_act = (n_tie_act - 1) * tie_points_factor
-    n_pixel_alt = (scan_alt_tie_points - 1) * tie_points_factor * n_scans
+    n_pixel_alt = (n_tie_alt - 1) * tie_points_factor
 
     # 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 for the current scan
-    n_pixel_alt_scan = (scan_alt_tie_points - 1) * tie_points_factor
-    tie_grid_alt_scan = da.arange(0, n_pixel_alt_scan + 1, tie_points_factor)
-    pixels_grid_alt_scan = da.arange(0, n_pixel_alt_scan)
+    # 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 = []
 
-    data_on_pixel_points = []
+    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))
+    pixel_grid_alt = da.concatenate(pixel_grid_alt)
 
-    # loop over granules
+    # Loop on all arrays
+    data_on_pixel_points = []
     for data in data_on_tie_points:
 
-        # create an xarray to hold the data on pixel grid (note the coords get renamed later so keep the old names)
-        rads = da.zeros((n_pixel_alt, n_pixel_act))
-        pix_act = da.zeros(n_pixel_act)
-        pix_alt = da.zeros(n_pixel_alt)
-
-        data_on_pixel_points_granule = \
-            xr.DataArray(rads, dims=['num_tie_points_alt', 'num_tie_points_act'],
-                         coords={'num_tie_points_alt': pix_alt, 'num_tie_points_act': pix_act})
-        data_on_pixel_points_granule.attrs = combine_metadata(data)
-
-        # loop over scans
-        for j_scan in range(n_scans):
-            index_tie_points_start = j_scan * scan_alt_tie_points
-            index_tie_points_end = (j_scan * scan_alt_tie_points) + scan_alt_tie_points
-            index_pixel_start = j_scan * (scan_alt_tie_points - 1) * tie_points_factor
-            index_pixel_end = (j_scan * (scan_alt_tie_points - 1) * tie_points_factor) + \
-                              (scan_alt_tie_points - 1) * tie_points_factor
+        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")
 
-            data_on_tie_points_scan = data[index_tie_points_start:index_tie_points_end, :]
+        # 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_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_tie_points_scan = data_on_tie_points_scan.assign_coords(
-                {dim_alt: tie_grid_alt_scan, dim_act: tie_grid_act})
-            data_pixel = data_on_tie_points_scan.interp({dim_alt: pixels_grid_alt_scan}, assume_sorted=True) \
-                .interp({dim_act: pixels_grid_act}, assume_sorted=True).drop_vars([dim_alt, dim_act])
-
-            # put the interpolated data in the pixel resolution granule
-            data_on_pixel_points_granule[index_pixel_start:index_pixel_end, :] = data_pixel
-
-        data_on_pixel_points.append(data_on_pixel_points_granule)
+        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.,