Spatial Rectification Algorithm

This chapter describes the algorithm used in the function rectify_dataset() of module xcube.core.resampling. The function geometrically transforms spatial data variables of a given dataset from an irregular source grid mapping into new data variables for a given regular target grid mapping and returns them as a new dataset.

Problem Description

The following figure shows a Sentinel-3 OLCI Level-1b scene in its original satellite perspective. In addition to the measured reflectances, the data product also provides the latitude and longitude for each pixel as two individual coordinate images. To the right, the absolute value of the gradient vectors of the latitude and longitude |(∇ lon)² + (∇ lat)²| are shown:

The given coordinates, latitude and longitude, are terrain corrected with respect to a digital elevation model (DEM) that approximates the Earth’s true fractal surface. Therefore, the gradient vectors are not monotonically varying over the scene, instead they represent the roughness of the DEM. Areas of high surface roughness are indicated by the red circles in the figure above. The figure below should explain the cause of the effect.

A rectification is the transformation of satellite imagery from its original viewing geometry into a target geometry that forms a regular grid in a defined coordinate reference system (CRS) with uniform spatial resolution for each pixel in each dimension. For the Sentinel-3 OLCI Level-1b scene above, the rectified measurement image for the geographic projection (CRS EPSG:4326) is shown here:

Algorithm Description

The input to the rectification algorithm is satellite imagery in satellite viewing perspective. In addition, two images – one for each spatial dimension – provide the terrain corrected spatial coordinates. Thus, for each source pixel we have a given spatial coordinate pair x,y, which is assumed to refer to a source image pixel’s center at i+½,j+½. The expected algorithm input is:

• N source measurement images;

• 2 source coordinate images (one for each spatial dimension) comprising terrain corrected coordinates x,y in units of the CRS for each pixel i,j;

The output produced by the algorithm is

• N target measurement images generated from source measurement images by projection.

• Target image geometry given as pixel size Δx=Δy and the coordinate offset of the upper left target image pixel x0,y0. Pixel size and offsets are given in units of the CRS.

• 2 target lookup images (one for each spatial dimension) comprising fractional source pixel coordinates i+½+u, j+½+v (explained below);

In the following, the CRS of the target and source is assumed to be the same. This is not a limitation, see remarks below.

While the gradient of the coordinate images, Δx and Δy per pixel, is generally not constant in the source image geometry, we demand that the target pixel size Δx = Δy = const for all pixels in the generated target images.

A fast and simple algorithm to perform the rectification is to visit each source pixel i,j, collect the spatial coordinates P = x,y, span two triangles between four adjacent coordinates, (P1, P2, P3) and (P2, P4, P3), and “paint” fractional source pixel coordinates into a new target lookup image. The lookup image can than be used to retrieve the pixel values from a source measurement image for a given pixel of the corresponding target measurement image, either by nearest neighbor lookup or by interpolation.

The true Earth surface is unknown in between any given coordinates points P(i+½, j+½) and its neighborhood, and there is no defined “best guess” for any point P(i+u, j+v) with 0 ≤ u ≤ 1 and 0 ≤ v ≤ 1. Hence, we use triangulation for its simplicity so that any in-between P is found by triangular interpolation.

From the coordinates (P1, P2, P3) of the first source triangle, the bounding box in pixel coordinates in the target image can be exactly determined, because the target grid is regular, P = x0 + i Δx, y0 - j Δy, for each target pixel i,j, target pixel size Δx,Δy, and x0,y0 being the coordinates of the upper left pixel at i=½, j=½. Given P and the plane given by (P1, P2, P3) the parameters u, v can be computed from triangular interpolation form P = P1 + u (P2 – P1) + v (P3 – P1). If 0 ≤ u ≤ 1, 0 ≤ u ≤ 1, u+v ≤ 1, then P is a point within the triangle.

At the same time, u and v are the fractions of source pixel coordinate, i + ½ + u and j + ½ + v, which will be both stored in two target lookup images.

After all source pixels have been processed, the resulting target lookup images containing the fractional source pixel indexes, i + ½ + u and j + ½ + v, can be used to efficiently map a source measurement image V into the target measurement image. The pixel values of the source measurement image are given as V1 = V(i, j), V2 = V(i+1, j), V3 = V(i ,j+1), and V4 = V(i+1, j+1) here and are represented by different colour values, e.g., measurements such as radiances, reflectances, or higher level data variables:

In the simplest case, as shown above, a nearest neighbor lookup is performed to determine the pixel value V for the target measurement image according to:

V = V2 if u > ½; V3 if v > ½; V1 else

The fractions u, v can also be used to perform the triangular interpolation between the three adjacent source measurements pixels. Note that the u, v were found from the three coordinate pairs that formed the original source triangle:

V = VA + u (V2 − V1) + v (V3 − V1)

Using bilinear interpolation between the four adjacent source measurements pixels also takes the fourth source coordinate P4 into account:

V = VA + v (VB − VA)

with

VA = V1 + u (V2 − V1)
VB = V3 + u (V4 − V3)

Remarks

(1) The target pixel size should be less or equal the source pixel size, so that triangle patches in the source refer to multiple pixels in the target image. However, if the target pixel size is distinctly smaller than the source pixel size, and the source has a low spatial resolution, results will be inaccurate, as curved source pixel boundaries need to be taken into account for many projections.

(2) If the target pixel size is greater than the source pixel size, target bounding boxes refer to single source pixels only. However, if the target pixel size is distinctly larger than the source pixel size, a prepended down-sampling or convolution of source images should be taken into account, e.g., using a Gaussian filter. Note, the algorithm can also be easily adapted to aggregate values of source pixels that refer to same target pixels.

(3) The algorithm description assumes that the source CRS and target CRS be the same. Different reference systems can be easily and efficiently supported by initially transforming each coordinate P(i, j) = (x, y) in the source coordinate images into the desired target CRS, i.e., P_new = transform_coord(P, CRS_S, CRS_T).

(4) If x, y are decimal longitude and latitude, and the north or south poles are in the scene, than the algorithm will fail. One can get around this problem by transforming source coordinates into a another suitable CRS first or by transforming longitude values x into complex numbers and optionally normalizing latitudes y to the same range from -1 to +1:

x’ = cos(x) + i sin(x)
y’ = 2y / π

(5) The algorithm can be rewritten to use bilinear surface patches comprising the adjacent four coordinates (P1, P2, P2, P4) instead of two triangles. The bilinear interpolation between the four coordinates is

PA = P2 + u (P2 − P1)
PB = P3 + u (P4 − P3)
P = PA + v (PB − PA)

Given that P is a known point on that surface the above can be rewritten as the two equations in u and v

0 = a + u b + v c + u v d
0 = e + u f + v g + u v h

with the solutions for u and v (from Wolfram|Alpha¹):

and with the differences

a = P1.x − P.x
b = P2.x − P1.x
c = P3.x − P1.x
d = P4.x − P3.x
e = P1.y − P.y
f = P2.y − P1.y
g = P3.y − P1.y
h = P4.y − P3.y

¹ Wolfram Alpha LLC, 2024. Wolfram|Alpha. https://www.wolframalpha.com/input?i2d=true&i=0+%3D+a+%2B+u++b+%2B+v++c+%2B+u++v++d%5C%2844%29+0+%3D+e+%2B+u++f+%2B+v++g+%2B+u++v++h+for+u%5C%2844%29+v (access April 8, 2024).