diff --git a/geotiepoints/viiinterpolator.py b/geotiepoints/viiinterpolator.py index e17b5ce..d49dcfa 100644 --- a/geotiepoints/viiinterpolator.py +++ b/geotiepoints/viiinterpolator.py @@ -19,17 +19,19 @@ """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] @@ -37,21 +39,16 @@ 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.,