# Avant-propos : un panorama de la recherche actuelle sur le quaternaire : le colloque Q10 “ paléoclimats et environnements quaternaires, quoi de neuf sous le soleil? ”

Amaëlle Landais, Maria-fernanda Sanchez-goñi, Nicole Limondin-lozouet

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Amaëlle Landais, Maria-fernanda Sanchez-goñi, Nicole Limondin-lozouet. Avant-propos : un panorama de la recherche actuelle sur le quaternaire : le colloque Q10 “ paléoclimats et environnements quaternaires, quoi de neuf sous le soleil? ”. Quaternaire, 2017, 28 (2), pp.135-136. 10.4000/quaternaire.8080. hal-03226714

HAL Id: hal-03226714

https://hal.science/hal-03226714v1

Submitted on 17 Jul 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Automatic sub-pixel co-registration of Landsat-8 OLI and Sentinel-2A MSI images using phase correlation and machine learning based mapping

Sergii Skakun<fnref n="1" />,<fnref n="2" />, Jean-Claude Roger<fnref n="1" />,<fnref n="2" />, Eric F. Vermote<fnref n="2" />, Jeffrey G. Masek<fnref n="3" />, Christopher O. Justice<fnref n="1" />

## Abstract

This study investigates misregistration issues between Landsat-8/OLI and Sentinel-2A/MSI at 30 m resolution, and between multi-temporal Sentinel-2A images at 10 m resolution using a phase correlation approach and multiple transformation functions. Co-registration of 45 Landsat-8 to Sentinel-2A pairs and 37 Sentinel-2A to Sentinel-2A pairs were analyzed. Phase correlation proved to be a robust approach that allowed us to identify hundreds and thousands of control points on images acquired more than 100 days apart. Overall, misregistration of up to 1.6 pixels at 30 m resolution between Landsat-8 and Sentinel-2A images, and 1.2 pixels and 2.8 pixels at 10 m resolution between multi-temporal Sentinel-2A images from the same and different orbits, respectively, were observed. The non-linear Random Forest regression used for constructing the mapping function showed best results in terms of root mean square error (RMSE), yielding an average RMSE error of 0.07±0.02 pixels at 30 m resolution, and 0.09±0.05 and 0.15±0.06 pixels at 10 m resolution for the same and adjacent Sentinel-2A orbits, respectively, for multiple tiles and multiple conditions. A simpler 1st order polynomial function (affine transformation) yielded RMSE of 0.08±0.02 pixels at 30 m resolution and 0.12±0.06 (same Sentinel-2A orbits) and 0.20±0.09 (adjacent orbits) pixels at 10 m resolution.

## Keywords

sub-pixel co-registration; phase correlation; misregistration; Landsat-8; Sentinel-2; machine learning; random forest

## 1. Introduction

Many applications in climate change and environmental and agricultural monitoring rely heavily on the exploitation of multi-temporal satellite imagery. Multi-temporal satellite images can help to identify and analyze changes in land cover land use (LCLUC) (Justice et

<footnote n="1">¹Department of Geographical Sciences, University of Maryland, College Park, MD, USA</footnote>

<footnote n="2">²NASA Goddard Space Flight Center Code 619, Greenbelt, MD, USA</footnote>

<footnote n="3">³NASA Goddard Space Flight Center Code 618, Greenbelt, MD, USA</footnote>
al. 2015), to capture significant trends in land surface properties, e.g. greenness (Ju and Masek 2016), or to discriminate specific crop types (Shelestov et al. 2017, Skakun et al. 2016), that cannot be identified with a single date image. In order to solve these problems more efficiently at high spatial resolutions (30 m), combined use of freely available Landsat-8 and Sentinel-2 images can offer high temporal frequency of about 1 image every 3–5 days globally.

Landsat-8 was launched in 2013 within the Landsat Program, a joint effort between the U.S. Geological Survey (USGS) and NASA (Roy et al. 2014). The Operational Land Imager (OLI) and Thermal Infrared Sensor (TIRS) instruments onboard the Landsat-8 satellite capture images of the Earth's surface in eleven spectral bands of the electromagnetic spectrum at 30 m spatial resolution (15 m for panchromatic band and 100 m for thermal infrared). The swath of Landsat-8 scene is approximately 185 km, allowing global coverage of the Earth's surface every 16 days. To refine image geolocation, the Landsat-8 processing uses ground control points (GCPs) automatically derived from the Global Land Survey (GLS) Landsat images (Gutman et al. 2013). The Sentinel-2A satellite was launched in 2015 within the European Copernicus program (Drusch et al. 2012). Sentinel-2A has a Multi-Spectral Instrument (MSI) that acquires images of the Earth's surface in thirteen spectral bands at 10 m, 20 m and 60 m spatial resolution. The swath of a Sentinel-2A scene is approximately 290 km, allowing global coverage of the Earth's surface every 10 days. The launch of an identical Sentinel-2B satellite will further reduce revisit time to 5 days globally. Both Sentinel-2A/B satellites will use a Global Reference Image (GRI) derived from orthorectified Sentinel-2 cloud-free images to improve geolocation accuracy and repeatability to meet the requirements of multi-temporal registration of 3m for 10 m bands (Déchoz et al. 2015). The GRI dataset is currently under development and is expected to be completed in 2018 (Storey et al. 2016).

Recent studies are focusing on the combined use of Landsat-8 and Sentinel-2A images to increase temporal coverage; however misregistration issues between Landsat-8 and Sentinel-2A have already been identified (Storey et al. 2016). It has been found that the OLI and MSI misregistration can exceed one 30-meter pixel and, therefore, it is recommended to exploit image registration approaches to further improve alignment between Landsat-8 and Sentinel-2A images (Storey et al. 2016). These approaches should be automatic and computationally efficient in order to perform alignment on a global basis, have sub-pixel accuracy and effectively deal with temporal changes, so these approaches can be further applied with GLS and GRI.

There have been many studies carried out to develop automatic satellite image registration methods (e.g. Gao et al. 2009, Le Moigne et al. 2011, Zitova and Flusser, 2003). The general image-to-image registration workflow consists of automatic generation of control points (CPs) between the *reference* (or *master*) and *sensed* (or *slave*) images, building and evaluating a spatial transformation (mapping function) that aligns the reference and sensed images, and warping the sensed image with radiometric transformation. *Area*-based and *feature*-based approaches are used to automatically derive CPs. Area-based methods, also referred to as correlation-like or template matching, find correspondence between reference and sensed images through a similarity measure, for example cross-correlation (in spatial or
frequency domain) or mutual information. These measures are usually applied on a sliding window basis to derive a dense set of CPs. Feature-based methods aim to find distinctive features on images, for example edges, contours, line intersections, closed boundary regions, and then match the derived features to find correspondences between reference and sensed images. The derived CPs are used to construct and evaluate a mapping function that maps points from the reference image to points in the sensed image. Examples of the mapping function include translation, affine transformation, high-order polynomials, radial basis functions (RBFs), and elastic registration. And finally, radiometric transformation (nearest neighbor, bilinear, splines) should be specified to warp the sensed image to the reference one.

The problem of Landsat-8/OLI and Sentinel-2A/MSI misregistration has already been addressed in several previous studies (Barazzetti et al. 2016, Yan et al. 2016). Yan et al. (2016) used a hierarchical feature-based matching approach to find CPs through construction of image pyramids at various spatial resolutions and an area-based matching approach to further refine and reject non-reliable CPs. Translation, affine transformation and second order polynomial functions were evaluated in the study for three pairs of Landsat-8 and Sentinel-2A images with affine transformation giving the best results in terms of root mean square error (RMSE) of 0.3 pixels at 10 m resolution. Barazzetti et al. (2016) utilized standard software packages to study misregistration between Landsat-8 and Sentinel-2A, and achieved RMSE of up to 1.2 pixels at 15 m spatial resolution.

In this paper, we explore a phase correlation approach to automatically generate a dense grid of CPs when registering Landsat-8 to Sentinel-2A images, as well as multiple mapping functions including those based on machine learning approaches. We also address the issue of multi-temporal misregistration between Sentinel-2A images. (For assessment of registration accuracy of multi-temporal Landsat-8 images, we refer readers to Storey et al. (2014)). Our analysis shows misregistration magnitudes of up to 3 pixels at 10 m resolution can be observed. This issue has not been reported in previous studies to our knowledge (except Sentinel-2 Data Quality Reports, see ESA (2016a), and reports delivered at the Land Cover Land Use Change (LCLUC) Multi-Source Land Imaging (MuSLI) Science Team 2016, http://lcluc.umd.edu/meetings/2016-lcluc-spring-science-team-meeting-18-19-april-and-musli-science-team-meeting-20-21), and should be further addressed especially for users dealing with Sentinel-2A time-series at 10 m spatial resolution. As in Yan et al. (2016), we use near-infrared (NIR) bands from Landsat-8 (band 5, 0.85 – 0.88 um) and Sentinel-2A (band 08, 0.842 um) to find CPs on reference and sensed images, since the NIR provides a wide dynamic range of values for multiple land cover types and is less sensitive to atmospheric effects. Co-registration of Landsat-8 to Sentinel-2A was performed at 30 m spatial resolution while co-registration between Sentinel-2A images was undertaken at 10 m spatial resolution.

## 2. Data description

### 2.1. Landsat-8 and Sentinel-2A products description

We used a standard Landsat-8 Level-1 terrain corrected (L1T) product distributed by USGS through the EarthExplorer system (Roy et al. 2014). The product is provided in the World-
wide Reference System (WRS-2) of path and row coordinates. The size of the Landsat-8
scene is approximately 185 km × 180 km. The product is provided with the corresponding
metafile to convert digital numbers (DNs) into the top-of-atmosphere reflectance (TOA)
values. For Sentinel-2, we used a standard Level-1C (L1C) product which is radiometrically
and geometrically corrected with ortho-rectification, and provided with the TOA reflectance
values (ESA 2016b). The product is delivered in tiles, or granules, of approximately 110 km
× 110 km size. Both Landsat-8 L1T and Sentinel-2A L1C products are provided in the
Universal Transverse Mercator (UTM) projection with the World Geodetic System 1984
(WGS84) datum. Each Sentinel-2 tile is assigned a UTM zone, and overlapping tiles
covering the same geographic region might have different UTM zones assigned. The
Sentinel-2 tiling grid is referenced to the U.S. Military Grid Reference System (MGRS). A
tile identifier consists of five signs: two numbers and three letters, e.g. 20HNH. The first two
numbers in the tile identifier correspond to the UTM zone while the remaining three letters
correspond to the tile position. ESA provides a kml file with tile coverage and their
identifiers (https://sentinel.esa.int/documents/247904/1955685/
S2A_OPER_GIP_TILPAR_MPC__20151209T095117_V20150622T000000_21000101T00
0000_B00.kml).

It should be also noted that pixel value is for the center of the pixel for the Landsat-8 L1T product, while it is for the upper left corner of the pixel for the Sentinel-2A L1C product. In this work, we used Sentinel-2 tile system as a reference, i.e. Landsat-8 data were subset for the corresponding Sentinel-2 tile with the nearest neighborhood resampling and co-registered to the reference Sentinel-2 scene using the proposed approach. Depending on the application and user needs other reference systems can be specified, for example Web Enabled Landsat Data (WELD) (Roy et al. 2010) and this proposed approach can be easily adapted to it.

## 2.2. Landsat-8 and Sentinel-2A test data description

This study was carried out for five Sentinel-2 tiles in three countries: Argentina (Sentinel-2 tile grid numbers 20HNH and 20HPH), US (14SKF) and Ukraine (36UUU, 34UFU) (see Figure 1 for examples). The selected tiles cover intensive agriculture regions, where changes are rapid due to seasonal crop development, and a mountain region in the Carpathians (tile 34UFU) where surface elevation varies approximately from 100 m to 1000 m. All selected regions feature considerable difference in Landsat-8 and Sentinel-2A image acquisition dates ranging from July 4, 2015 to July 25, 2016, as well as variable cloud conditions.

Overall co-registration of 45 Landsat-8 to Sentinel-2A pairs and 37 Sentinel-2A to
Sentinel-2A pairs were analyzed. Table 1 gives details on the Landsat-8 and Sentinel-2A
imagery used in the study.

# 3. Methodology

**3.1. General overview**

The proposed approach follows the general concept of automatic image-to-image
registration outlined in (Zitova and Flusser 2003). It has the following steps (Figure 2):
image pre-processing; automatic identification of CPs; filtering of CPs; building and evaluating a transformation, and image warping.

## 3.2. Image pre-processing

Landsat-8 images were converted from DNs to TOA reflectance values using calibration coefficients in the metadata file, so both Landsat-8 and Sentinel-2A images were generated with TOA reflectance values. The steps outlined in this and following subsections will also be valid, should co-registration be applied for atmospherically corrected products (Vermote et al. 2016).

For Landsat-8 to Sentinel-2A co-registration, Sentinel-2A band 08 (NIR) was resampled from 10 m to 30 m using averaging, and Landsat-8 data were subset to the corresponding Sentinel-2A tile using the nearest neighborhood resampling technique. Sentinel-2A to Sentinel-2A co-registration was performed at the original 10 m spatial resolution without further resampling.

## 3.3. Automatic generation of CPs

In this study, we exploited a phase-only correlation image matching method introduced by Guizar-Sicairos et al. (2008). It uses a cross-correlation approach in the frequency domain by means of the Fourier transform and exploits a computationally efficient procedure based on nonlinear optimization and Discrete Fourier Transforms (DFTs) to detect sub-pixel shifts between reference and sensed images. For the detailed description of the algorithm, we refer readers to Guizar-Sicairos et al. (2008).

The phase-correlation algorithm allows detection of translation between reference and sensed images, and therefore is routinely applied using a moving square window. The size of the window and the step are selected empirically. Window size should be large enough to capture similarities on the reference and sensed images, and small enough to have a dense grid of CPs to accurately construct a transform function. In this study, window size and step were selected 100 and 50 pixels, respectively when co-registering Landsat-8 to Sentinel-2A, and 64 and 32 when co-registering multi-temporal Sentinel-2A images.

Compared to other area-based methods (e.g. cross-correlation in the spatial domain), the phase-correlation image usually contains a sharp peak corresponding to the dominant shift between images, and is usually more robust to temporal changes between reference and sensed images (Kravchenko et al. 2014). Compared to feature-based methods, the phase-correlation approach with a moving window allows detection of a dense grid of CPs, especially in cases where features cannot be reliably identified and detected, including at the 30 m spatial resolution imagery.

No matter what method is applied for automatic generation of CPs, filtering is necessary to remove unreliable CPs. First, a peak cross-correlation normalized magnitude is used for initial rejection of CPs. In our study, this value was set to 0.5. After that, a RANdom SAmple Consensus (RANSAC) algorithm (Fischler and Bolles 1981) is run for the linear transformation model to detect inliers and outliers. The RANSAC is a widely used algorithm in computer vision and image processing to detect strong outliers with large deviations
(Brown and Lowe 2007). However, one should be careful to not aggressively remove outliers that can be actually inliers. In our study, we ran a conservative approach for removing outliers, i.e. the RANSAC parameters were set in such a way to remove only outliers with a high confidence level. In our particular case, the number of RANSAC trials was set to 100, and a confidence level for selecting outliers was set to 0.99. The detected inlier CPs are further used to construct and evaluate a transformation function.

### 3.4. Transformation function

The goal of constructing a transformation function $\mathcal{R}$ is to find correspondence between points in the reference image $\mathbf{x}_r = (x_r, y_r)$ and points in the sensed image $\mathbf{x}_s = (x_s, y_s)$:

$$ (x_s, y_s) = F(x_r, y_r). \qquad (1) $$

Function $\mathcal{R}$ is constructed using the set of CPs identified in the previous steps. In this study we compared three different approaches to creating the transformation function: polynomial models, radial basis functions, and random forest regression trees (see below for details). Regardless of transformation approach, all available CPs were randomly split into a training (calibration) set (80%) and a testing set (20%). The training set was used to build the model and identify its parameters, while the testing set was used to evaluate the model on independent data. We will denote CPs from the testing set with $(x_{r,l}, y_{r,l})$ and $(x_{s,l}, y_{s,l})$ where $l = \overline{1, L}$ and $L$ is the number of points with corresponding shifts:

$$ \Delta_{x,l} = x_{r,l} - x_{s,l}, \Delta_{y,l} = y_{r,l} - y_{s,l} \qquad (2) $$

The quality of the transformation is evaluated using a root mean square error (RMSE) between reference values $(x_{s,l}, y_{s,l})$ and estimated values $(\hat{x}_{s,l}, \hat{y}_{s,l})$ by transformation function $\mathcal{R}$ using the testing set:

$$ (\hat{x}_{s,l}, \hat{y}_{s,l}) = F(x_{r,l}, y_{r,l}), \qquad (3) $$

$$ RMSE = \sqrt{\frac{1}{L} \sum_{l=1}^{L} ((\hat{x}_{s,l} - x_{s,l})^2 + (\hat{y}_{s,l} - y_{s,l})^2)}. \qquad (4) $$

#### 3.4.1. Polynomial models—Selection of the type of the transformation function
depends on a-priori knowledge of expected geometric deformations and distortions between reference and sensed images, and required registration accuracy (Zitova and Flusser 2003). Polynomial functions of the n-th degree have the following form:

$$ x_s = P_{x,n}(x_r, y_r) = \sum_{i=0}^{n} \sum_{j+1 \le n} a_{ij} x_r^j y_r^j \qquad (5) $$
$$y_s = P_{y,n}(x_r, y_r) = \sum_{i=0}^{n} \sum_{j \le n} b_{ij} x_r^i y_r^j. \qquad (6)$$

In the case of n=0, the models (5)–(6) are simple translation models where the same values of shift, namely $a_{00}$ and $b_{00}$, are applied in x and y directions, respectively:

$$x_s = P_{x,0}(x_r, y_r) = a_{00}, \qquad (7)$$

$$y_s = P_{y,0}(x_r, y_r) = b_{00}. \qquad (8)$$

A linear model, also referred to as affine, can be reduced to the following form:

$$x_s = P_{x,1}(x_r, y_r) = a_{00} + a_{10}x_r + a_{01}y_r, \qquad (9)$$

$$y_s = P_{y,1}(x_r, y_r) = b_{00} + b_{10}x_r + b_{01}y_r. \qquad (10)$$

Model parameters $a_{ij}$ and $b_{ij}$ are estimated through the ordinary least square (OLS) method, by minimizing the sum of the squares of the differences between predicted values of the model and reference values.

**3.4.2. Radial Basis Functions (RBFs)**—The transformation function based on Radial Basis Functions (RBFs) has the following form (Zitova and Flusser 2003):

$$x_s = a_0 + a_1 x_r + a_2 y_r + \sum_{k=1}^{K} w_k K(\mathbf{x}_r, \mathbf{x}_k), \qquad (11)$$

$$y_s = b_0 + b_1 x_r + b_2 y_r + \sum_{k=1}^{K} \tilde{w}_k K(\mathbf{x}_r, \tilde{\mathbf{x}}_k), \qquad (12)$$

where $K(\cdot, \cdot)$ is the kernel function with parameters (centers) $\mathbf{x}_k$ and $\tilde{\mathbf{x}}_k$ and weights $w_k$ and $\tilde{w}_k$.

In this study, we used two types of kernels, namely Gaussian and thin-plate splines (TPS):

$$\text{Gaussian: } K(\mathbf{x}, \mathbf{x}_k) = \exp(-\|\mathbf{x} - \mathbf{x}_k\|^2), \qquad (13)$$

$$\text{TPS: } K(\mathbf{x}, \mathbf{x}_k) = \|\mathbf{x} - \mathbf{x}_k\|^2 \ln(\|\mathbf{x} - \mathbf{x}_k\|). \qquad (14)$$

There are several ways of selecting a set of centers $\mathbf{x}_k$: randomly, on a regular grid, or adaptively through clustering. We used the k-means clustering approach (Lloyd 1982, Forgy 1965) to adaptively select centers in the models (11)–(12). We varied number of clusters K
and found values from 1 to 10 producing best results. Increasing the number over 10 did not
improve results.

And finally, weights $w_k$ and $\tilde{w}_k$ in models (11)–(12) were estimated from training data using
the RANSAC algorithm (Fischler and Bolles 1981).

As with polynomial models, RBF models are the global mapping functions, however, they
are able to account for local non-linear distortions (Zitova and Flusser 2003).

**3.4.3. Random forest (RF) regression—**Random forest (RF) is a machine learning
algorithm that represents an ensemble of decision trees (DTs) (Breiman 2001). A DT
classifier or regression model is built from a set of data using the concept of information
entropy. At each node of the tree, one attribute of the data, that most effectively splits its set
of samples into subsets enriched in one class or the other, is selected. Its criterion is the
normalized information gain that results from choosing an attribute for splitting the data.
The attribute with the highest normalized information gain is chosen to make the decision.
The algorithm then recurs on the smaller sublists. One disadvantage of the DT classifier is
the considerable sensitivity to the input dataset, so that a small change to the training data
can result in a very different set of subsets (Bishop 2006). In order to overcome
disadvantages of a single DT, an ensemble of DTs is used to form a random forest. Each DT
represents an independent expert (or weak classifier) in the RF that is trained based on
different input datasets that are generated through a bagging procedure (Bishop 2006). RF
can be used for building classification and regression models.

In this study, the RF regression was used to build a transformation function. We used points
from the reference image $\mathbf{x}_r = (x_r, y_r)$ as inputs to the RF regression model with a
polynomial preprocessing. For example, in case of the 2^nd^ degree polynomial function, the
following features were input to the RF model: $x_r, y_r, x_r^2, y_r^2, x_ry_r$. The RF model was further
trained to predict points in the sensed image $\mathbf{x}_s = (x_s, y_s)$. In this study, the number of DTs in
the RF models was kept low (about 5) to avoid overfitting. The optimal number of DTs in
the RF model in terms of RMSE error was identified through the cross-validation procedure.

As with RBFs based mapping functions, RF belongs to the class of global mapping models
that can account for local non-linear distortions.

**4. Results**

The use of phase correlation allowed us to generate hundreds and thousands of CPs when
co-registering Landsat-8 to Sentinel-2A at 30 m spatial resolution (hereafter referred as
LandSen30), and when co-registering multi-temporal Sentinel-2A images at 10 m spatial
resolution (hereafter referred as SenSen10) (Figure 3). The average misregistration between
Landsat-8 and Sentinel-2A calculated on a tile basis, using identified CPs, varied from 0.11
pixels to 1.35 pixels among tiles considered in the study with the maximum misregistration
value varying from 0.25 to 1.59 pixels (per tile). Misregistration between Landsat-8 and
Sentinel-2A was stable in time over the same tile (Table 2, Figure 4) with average standard
deviation of the misregistration through the time varying from 0.03 pixels to 0.16 pixels with average of 0.10 pixels (at 30 m resolution).

In the SenSen10 case, performance depended on whether the reference and sensed Sentinel-2A images were acquired from the same or different (adjacent) orbits. In case of the same orbits, the average misregistration calculated on a tile basis varied from 0.05 pixels to 0.46 pixels among tiles (with the maximum misregistration up to 1.21 pixels) with average 0.23±0.12 pixels at 10 m resolution. In case of different Sentinel-2 orbits over the same tile, the average misregistration calculated on a tile basis varied from 0.14 pixels to 1.45 pixels among tiles (with the maximum misregistration up to 2.83 pixels) with average 0.61±0.42 pixels at 10 m resolution (Table 3, Figure 5). This might be the result of the satellite yaw bias that was corrected within the recent baseline processing version 02.04 (ESA, 2016a). Overall, our estimates of the multi-temporal misregistration in the Sentinel-2A imagery were consistent with the Sentinel-2 Data Quality Reports (ESA, 2016a).

When building a transformation function, all considered approaches showed a similar performance (Table 4 and Table 5) with the complex non-linear RF regression slightly outperforming other methods, namely, a simple translation (Eq. 7–8), 1st order polynomial (Eq. 9–10), Gaussian RBFs (Eq. 11–12, 13) and TPS (Eq. 11–12, 14). For the RF regression, RMSE values varied from 0.02 to 0.12 pixels for the LandSen30 case at 30 m resolution, and from 0.025 to 0.22 pixels for the SenSen10 case at 10 m resolution for the same orbits and from 0.05 to 0.26 pixels for adjacent orbits. The RF model was built using a 1st order polynomial function for input CP coordinates. Increasing the order of the polynomial function did not improve performance of the RF model. It means that the RF was able to capture non-linearity when building a transformation function between CPs on the reference and sensed images.

Our results obtained for the 1st order polynomial function were comparable to previous studies by Barazzetti et al. (2016) and Yan et al. (2016).

Temporal fluctuations of the RMSE error were analysed to explore temporal stability of the constructed transformation functions for multi-temporal Sentinel-2A images co-registration. Figure 7 shows dependence of the RMSE error with time for the 1st degree polynomial function for different Sentinel-2A orbits. For the case of the same orbits, there are not too many variations in time, except three cases with RMSE values over 0.2 pixels. These are due to heavy cloud contamination presented in the sensed imagery. As to the adjacent orbits, the tile T20HNH shows a trend which is due to the difference between the sensed and reference images: the difference is up to 103 days. For the T36UUU case, high RMSE error can be attributed to the 70 day difference between the reference and sensed images acquired over highly intensive agriculture region that features a lot of changes within this time period (Figure 1, D). Also, reduction of the RMSE error for 36UUU tile comparing to the 20HNH tile can be related to the improvements made in the baseline processing version 02.04 (ESA, 2016a).

Results on correcting misregistration between Landsat-8 and Sentinel-2A are shown in Figure 8.
The proposed workflow was implemented in Python programming language using the Geospatial Data Abstraction Library (GDAL) for managing geospatial datasets, the 'skimage package' for phase correlation implementation, and the 'sklearn package' for the RF regression and RANSAC implementations, and was computationally efficient. Typical processing times on the Dell Laptop with processor Inter® Core™ i7-4810MQ CPU @ 2.80GHz with 16 Gb RAM are presented in Table 7.

## 5. Conclusions

Since many applications involving satellite imagery require the use of multi-temporal datasets, misregistration issues can lead to incorrect results. This study investigated misregistration issues between Landsat-8/OLI and Sentinel-2A/MSI at 30 m resolution, and between multi-temporal Sentinel-2A images at 10 m resolution using a phase correlation approach and multiple transformation functions. Phase correlation proved to be a robust approach that allowed us to identify hundreds and thousands of control points on images acquired more than 100 days apart. Overall, misregistration of up to 1.6 pixels at 30 m resolution between Landsat-8 and Sentinel-2A images, and 1.2 pixels (for the same orbit) and 2.8 pixels (for different orbits) at 10 m resolution between multi-temporal Sentinel-2A images were observed. The non-linear RF regression used for constructing the mapping function showed best results in terms of error, yielding the average RMSE error of 0.07±0.02 pixels at 30 m resolution and 0.09±0.05 and 0.15±0.06 pixels at 10 m resolution for the same and adjacent Sentinel-2A orbits, respectively, for multiple tiles and multiple conditions. On the other hand, a simple linear model such as 1st order polynomial function can provide an error of up to 0.08±0.02 pixels at 30 m resolution and 0.12±0.06 (same Sentinel-2A orbits) and 0.20±0.09 (adjacent orbits) pixels at 10 m resolution.

## References

Barazzetti L, Cuca B, and Previtali M 2016 "Evaluation of registration accuracy between Sentinel-2 and Landsat 8." In Fourth International Conference on Remote Sensing and Geoinformation of the Environment, pp. 968809–968809. International Society for Optics and Photonics. doi: 10.1117/12.2241765.

Bishop CM 2006 Pattern Recognition and Machine Learning. New York: Springer.

Breiman L 2001 "Random forests." Machine Learning 45 (1): 5–32. doi:10.1023/A:1010933404324.

Brown M, and Lowe DG 2007 "Automatic panoramic image stitching using invariant features." International Journal of Computer Vision 74 (1): 59–73. doi:10.1007/s11263-006-0002-3.

Dechoz C, Poulain V, Massera S, Languille F, Greslou D, de Lussy F, Gaudel A, L'Helguen C, Picard C, and Trémas T 2015 "Sentinel 2 global reference image." In SPIE Remote Sensing, pp. 96430A–96430A. International Society for Optics and Photonics. doi:10.1117/12.2195046.

Drusch M, Del Bello U, Carlier S, Colin O, Fernandez V, Gascon F, Hoersch B, et al. 2012 "Sentinel-2: ESA's optical high-resolution mission for GMES operational services." Remote Sensing of Environment 120: 25–36. doi:10.1016/j.rse.2011.11.026.

ESA. 2016a "S2 MPC. Data Quality Report." Ref: S2-PDGS-MPC-DQR, Issue: 08, Date: 05/10/2016. https://earth.esa.int/documents/247904/685211/Sentinel-2-Data-Quality-Report.

ESA. 2016b "Sentinel-2 Products Specification Document." Ref: S2-PDGS-TAS-DI-PSD, Issue: 14.0, Date: 15/07/2016. https://earth.esa.int/web/sentinel/user-guides/sentinel-2-msi/document-library.

Fischler MA, and Bolles RC 1981 "Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography." Communications of the ACM 24 (6): 381–395. doi:10.1145/358669.358692.
Forgy EW 1965 "Cluster analysis of multivariate data: efficiency versus interpretability of classifications." Biometrics 21: 768-769.

Gao F, Masek J, and Wolfe RE 2009 "Automated registration and orthorectification package for Landsat and Landsat-like data processing." Journal of Applied Remote Sensing 3 (1): 033515-033515. doi:10.1117/1.3104620.

Guizar-Sicairos M, Thurman ST, and Fienup JR 2008 "Efficient subpixel image registration algorithms." Optics Letters 33 (2): 156-158. doi:10.1364/OL.33.000156. [PubMed: 18197224]

Gutman G, Huang C, Chander G, Noojipady P, and Masek JG 2013 "Assessment of the NASA-USGS global land survey (GLS) datasets." Remote Sensing of Environment 134: 249-265. doi:10.1016/j.rse.2013.02.026.

Ju J, and Masek JG 2016 "The vegetation greenness trend in Canada and US Alaska from 1984-2012 Landsat data." Remote Sensing of Environment 176: 1-16. doi:10.1016/j.rse.2016.01.001.

Justice C, Gutman G, and Vadrevu KP 2015 "NASA land cover and land use change (LCLUC): An interdisciplinary research program." Journal of Environmental Management 148: 4-9. doi: 10.1016/j.jenvman.2014.12.004. [PubMed: 25500156]

Kravchenko O, Lavrenyuk M, and Kussul N 2014 "Orthorectification of Sich-2 satellite images using elastic models." In 2014 IEEE Geoscience and Remote Sensing Symposium, 2281-2284. IEEE. doi:10.1109/IGARSS.2014.6946925.

Le Moigne J, Netanyahu NS, and Eastman RD 2011 Image registration for remote sensing. Cambridge: Cambridge University Press.

Lloyd S 1982 "Least squares quantization in PCM." IEEE Transactions on Information Theory 28 (2): 129-137. doi:10.1109/TIT.1982.1056489.

Roy DP, Wulder MA, Loveland TR, Woodcock CE, Allen RG, Anderson MC, Helder D, et al. 2014 "Landsat-8: Science and product vision for terrestrial global change research." Remote Sensing of Environment 145: 154-172. doi:10.1016/j.rse.2014.02.001.

Roy DP, Ju J., Kline K, Scaramuzza PL, Kovalskyy V, Hansen M, Loveland TR, Vermote E, and Zhang C 2010 "Web-enabled Landsat Data (WELD): Landsat ETM+ composited mosaics of the conterminous United States." Remote Sensing of Environment 114 (1): 35-49. doi:10.1016/j.rse.2009.08.011.

Shelestov A, Lavreniuk M, Kussul N, Novikov A, and Skakun S 2017 "Exploring Google Earth Engine Platform for Big Data Processing: Classification of Multi-temporal Satellite Imagery for Crop Mapping." Frontiers in Earth Science 5(17): 1-7. doi:10.3389/feart.2017.00017.

Skakun S, Kussul N, Shelestov AY, Lavreniuk M, and Kussul O 2016 "Efficiency Assessment of Multi-temporal C-Band Radarsat-2 Intensity and Landsat-8 Surface Reflectance Satellite Imagery for Crop Classification in Ukraine." IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing 9 (8): 3712-3719. doi:10.1109/JSTARS.2015.2454297.

Storey J, Choate M, and Lee K 2014 "Landsat 8 Operational Land Imager on-orbit geometric calibration and performance." Remote Sensing 6 (11): 11127-11152. doi:10.3390/rs61111127.

Storey J, Roy DP, Masek J, Gascon F, Dwyer J, and Choate M 2016 "A note on the temporary misregistration of Landsat-8 Operational Land Imager (OLI) and Sentinel-2 Multi Spectral Instrument (MSI) imagery." Remote Sensing of Environment 186: 121-122. doi:10.1016/j.rse.2016.08.025.

Vermote E, Justice C, Claverie M, and Franch B 2016 "Preliminary analysis of the performance of the Landsat 8/OLI land surface reflectance product." Remote Sensing of Environment 185: 46-56. doi:10.1016/j.rse.2016.04.008.

Zitova B, and Flusser J 2003 "Image registration methods: a survey." Image and Vision Computing 21 (11): 977-1000. doi:10.1016/S0262-8856(03)00137-9.

Yan L, Roy DP, Zhang H, Li J, and Huang H 2016 "An automated approach for sub-pixel registration of Landsat-8 Operational Land Imager (OLI) and Sentinel-2 Multi Spectral Instrument (MSI) imagery." Remote Sensing 8 (6): 520. doi:10.3390/rs8060520.
<figure><img src="image_4.png" /></figure>

(A)
<figure><img src="image_5.png" /></figure>

(B)

<figure><img src="image_6.png" /></figure>

(C)

<figure><img src="image_7.png" /></figure>

(D)
<figure><img src="image_8.png" /><figcaption>**Figure 1.**</figcaption></figure>

Examples of TOA true color images acquired by Sentinel-2A/MSI (**left**) and Landsat-8/OLI (**right**): (A) tile 20HNH over Argentina, dates of Sentinel-2A/MSI and Landsat-8/OLI acquisitions are 2015358 and 2015361 respectively; (B) tile 20HPH over Argentina, acquisitions dates are 2015358 and 2015258; (C) tile 14SKF over Texas, US, dates of acquisitions are 2016012 and 2016104; (D) tile T36UUU over Ukraine, acquisition dates are 2016169 and 2016108; (E) tile T34UFU over the Carpathian Mountains, Ukraine, acquisitions dates are 2016198 and 2016063.
<figure><img src="image_9.png" /><figcaption>**Figure 2.**
General workflow of Landsat-8 and Sentinel-2A image co-registration.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 3.</figcaption></figure>

Location of CPs shown in the form of vectors outlining the direction and magnitude of shifts ($\Delta_x$ and $\Delta_y$ (Eq. 2)) found between Landsat-8 image acquired on 2016021 (21-Jan-2016), and Sentinel-2A image acquired on 2015358 (24-Dec-2015) and used as a reference image, over the study area in Argentina, tile T20HNH. Vector lengths were multiplied by 100 for visual clarity. Overall, 1634 CPs were found using the phase correlation approach in this case. The background is a Landsat-8 TOA NIR (band 5) image with TOA reflectance values scaled from 0.05 to 0.65.
<figure><img src="image_11.png" /><figcaption>**Figure 4.**</figcaption></figure>

Distribution of misregistration values $\Delta_x$ and $\Delta_y$ (Eq. 2) when co-registering Landsat-8 to Sentinel-2A images for different tiles used in the study. Units are shown in pixel values at 30 m spatial resolution.
<figure><img src="image_12.png" /><figcaption>**Figure 5.** Distribution of misregistration values $\Delta_X$ and $\Delta_Y$ when co-registering multi-temporal Sentinel-2A images from the same orbits for different tiles used in the study. Units are shown in pixel values at 10 m spatial resolution.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 6.**</figcaption></figure>

Distribution of misregistration values $\Delta_x$ and $\Delta_y$ when co-registering multi-temporal
Sentinel-2A images from the adjacent orbits for different tiles used in the study. Units are
shown in pixel values at 10 m spatial resolution.
<figure><img src="image_14.png" /><figcaption>Figure 7.</figcaption></figure>

Changes of RMSE error of building a 1st degree polynomial transformation function when registering multi-temporal Sentinel-2A images over the time for different tiles and different Sentinel-2A orbits: (A) same orbits; (B) adjacent orbits.
<figure><img src="image_15.png" /><figcaption>**Figure 8.**</figcaption></figure>

A 30 m “chessboard” composed of alternating Landsat-8 (acquired on 20-Dec-2015) and Sentinel-2A (24-Dec-2015) images before (left panel) and after co-registration (right panel). Near infrared images from band 5 (Landsat-8) and band 8 (Sentinel-2A) were used to produce these “chessboard”. TOA reflectance values were scaled from 0.05 to 0.55. This subset covers the area in the south-east part of the tile 20HNH over Argentina (Figure 1, A). Misregistrations between satellite images can be seen in the irrigated fields (circles, *middle left image*) and in the bridge over the lake (*bottom left image*) with corrections applied and misregistration disappearing in the *right images* (*middle* and *bottom*). Middle and bottom subset images are shown in corresponding boxes on the top images.
<caption>Table 1.</caption>

Description of data used in the study. Acquisition dates are given in the format YYYYDOY (where DOY is the day of the year).

<table><thead><tr><td>Country</td><td>Tile number</td><td>Acquisition date of Sentinel-2A reference image</td><td>Acquisition dates of Landsat-8 co-registered images</td><td>Acquisition dates of Sentinel-2A co-registered images</td></tr></thead><tbody><tr><td>Argentina</td><td>20HNH</td><td>2015358</td><td>2015354, 2015185, 2015201, 2015242, 2015249, 2015258, 2015290, 2015306, 2015329, 2015338, 2015345, 2015361, 2016021, 2016037, 2016053</td><td>2015341, 2016006, 2016013, 2016016, 2016023, 2016026, 2016036, 2016043, 2016046, 2016063, 2016065, 2016073, 2016083, 2016093, 2016096</td></tr><tr><td>Argentina</td><td>20HPH</td><td>2015358</td><td>2015242, 2015258, 2015290, 2015306, 2015338, 2015354, 2016021, 2016037, 2016053</td><td>2016003, 2016013, 2016023, 2016043, 2016063, 2016073, 2016083, 2016093</td></tr><tr><td>US (Texas)</td><td>14SKF</td><td>2016012</td><td>2015245, 2015261, 2015293, 2015309, 2015325, 2015341, 2015357, 2016024, 2016040, 2016056, 2016072, 2016088, 2016104</td><td>2016042, 2016072, 2016132</td></tr><tr><td>Ukraine</td><td>36UUU</td><td>2016169</td><td>2016076, 2016092, 2016108, 2016156, 2016172, 2016188</td><td>2016096, 2016109, 2016119, 2016156, 2016166, 2016179, 2016196, 2016199, 2016206</td></tr><tr><td>Ukraine</td><td>34UFU</td><td>2016198</td><td>2016063, 2016182</td><td>2016048, 2016208</td></tr></tbody></table>
<caption>Table 2.</caption>

Results of identifying CPs on the sensed (Landsat-8) and reference (Sentinel-2A) images using phase correlation approach at 30 m spatial resolution. Values for shifts $\Delta_x$ and $\Delta_y$ (Eq. 2) are shown in pixel units.

<table><thead><tr><th rowspan="2">Tile</th><th rowspan="2">Date of sensed image acquisition</th><th rowspan="2">Number of CPs</th><th colspan="4">Δ<sub>x</sub></th><th colspan="4">Δ<sub>y</sub></th></tr><tr><th>mean</th><th>std</th><th>min</th><th>max</th><th>mean</th><th>std</th><th>min</th><th>max</th></tr></thead><tbody><tr><td rowspan="12">20HNH</td><td>2015185</td><td>694</td><td>1.15</td><td>0.08</td><td>0.96</td><td>1.34</td><td>0.24</td><td>0.06</td><td>0.12</td><td>0.35</td></tr><tr><td>2015201</td><td>652</td><td>1.18</td><td>0.11</td><td>0.85</td><td>1.46</td><td>0.22</td><td>0.08</td><td>0</td><td>0.45</td></tr><tr><td>2015242</td><td>718</td><td>0.91</td><td>0.09</td><td>0.74</td><td>1.07</td><td>0.04</td><td>0.05</td><td>-0.07</td><td>0.13</td></tr><tr><td>2015249</td><td>512</td><td>1.14</td><td>0.11</td><td>0.88</td><td>1.34</td><td>0.13</td><td>0.06</td><td>-0.03</td><td>0.25</td></tr><tr><td>2015258</td><td>679</td><td>0.97</td><td>0.13</td><td>0.63</td><td>1.31</td><td>0.03</td><td>0.05</td><td>-0.1</td><td>0.16</td></tr><tr><td>2015290</td><td>939</td><td>0.98</td><td>0.09</td><td>0.72</td><td>1.23</td><td>-0.01</td><td>0.08</td><td>-0.24</td><td>0.17</td></tr><tr><td>2015306</td><td>810</td><td>0.99</td><td>0.08</td><td>0.8</td><td>1.17</td><td>0.02</td><td>0.05</td><td>-0.1</td><td>0.19</td></tr><tr><td>2015329</td><td>454</td><td>1.34</td><td>0.11</td><td>1.02</td><td>1.59</td><td>0.05</td><td>0.09</td><td>-0.16</td><td>0.28</td></tr><tr><td>2015338</td><td>1854</td><td>1</td><td>0.09</td><td>0.76</td><td>1.25</td><td>-0.03</td><td>0.11</td><td>-0.31</td><td>0.2</td></tr><tr><td>2015345</td><td>866</td><td>1.22</td><td>0.05</td><td>1.1</td><td>1.35</td><td>0.14</td><td>0.02</td><td>0.08</td><td>0.19</td></tr><tr><td>2015354</td><td>2338</td><td>0.96</td><td>0.13</td><td>0.59</td><td>1.32</td><td>-0.03</td><td>0.07</td><td>-0.24</td><td>0.13</td></tr><tr><td>2015361</td><td>1268</td><td>1.22</td><td>0.06</td><td>1.08</td><td>1.36</td><td>0.07</td><td>0.04</td><td>-0.03</td><td>0.16</td></tr><tr><td>2016021</td><td>1419</td><td>0.96</td><td>0.11</td><td>0.69</td><td>1.27</td><td>0.01</td><td>0.06</td><td>-0.17</td><td>0.18</td></tr><tr><td>2016037</td><td>1017</td><td>0.95</td><td>0.12</td><td>0.6</td><td>1.29</td><td>0.12</td><td>0.07</td><td>-0.06</td><td>0.29</td></tr><tr><td>2016053</td><td>876</td><td>0.81</td><td>0.15</td><td>0.41</td><td>1.19</td><td>0.08</td><td>0.06</td><td>-0.07</td><td>0.21</td></tr><tr><td rowspan="12">20HPH</td><td>2015242</td><td>726</td><td>0.88</td><td>0.11</td><td>0.53</td><td>1.14</td><td>-0.06</td><td>0.07</td><td>-0.27</td><td>0.1</td></tr><tr><td>2015258</td><td>652</td><td>0.93</td><td>0.13</td><td>0.54</td><td>1.19</td><td>-0.01</td><td>0.05</td><td>-0.13</td><td>0.09</td></tr><tr><td>2015290</td><td>848</td><td>0.89</td><td>0.1</td><td>0.56</td><td>1.13</td><td>-0.08</td><td>0.07</td><td>-0.27</td><td>0.06</td></tr><tr><td>2015306</td><td>709</td><td>0.98</td><td>0.07</td><td>0.73</td><td>1.16</td><td>-0.03</td><td>0.04</td><td>-0.17</td><td>0.06</td></tr><tr><td>2015338</td><td>2637</td><td>0.79</td><td>0.12</td><td>0.46</td><td>1.08</td><td>-0.09</td><td>0.07</td><td>-0.28</td><td>0.07</td></tr><tr><td>2015354</td><td>2912</td><td>0.82</td><td>0.11</td><td>0.53</td><td>1.08</td><td>-0.13</td><td>0.07</td><td>-0.3</td><td>0.04</td></tr><tr><td>2016021</td><td>1634</td><td>0.77</td><td>0.1</td><td>0.5</td><td>1</td><td>-0.08</td><td>0.07</td><td>-0.28</td><td>0.09</td></tr><tr><td>2016037</td><td>1070</td><td>0.82</td><td>0.1</td><td>0.54</td><td>1.05</td><td>0.02</td><td>0.07</td><td>-0.18</td><td>0.18</td></tr><tr><td>2016053</td><td>574</td><td>0.71</td><td>0.11</td><td>0.41</td><td>0.91</td><td>0</td><td>0.07</td><td>-0.21</td><td>0.16</td></tr><tr><td rowspan="6">36UUU</td><td>2016076</td><td>770</td><td>0.28</td><td>0.08</td><td>0.08</td><td>0.47</td><td>0.6</td><td>0.07</td><td>0.44</td><td>0.76</td></tr><tr><td>2016092</td><td>695</td><td>0.33</td><td>0.1</td><td>0.08</td><td>0.62</td><td>0.58</td><td>0.06</td><td>0.41</td><td>0.73</td></tr><tr><td>2016108</td><td>856</td><td>0.29</td><td>0.09</td><td>0.1</td><td>0.55</td><td>0.6</td><td>0.06</td><td>0.45</td><td>0.79</td></tr><tr><td>2016156</td><td>796</td><td>0.27</td><td>0.09</td><td>0.08</td><td>0.51</td><td>0.57</td><td>0.06</td><td>0.38</td><td>0.72</td></tr><tr><td>2016172</td><td>2379</td><td>0.31</td><td>0.08</td><td>0.14</td><td>0.5</td><td>0.6</td><td>0.05</td><td>0.46</td><td>0.73</td></tr><tr><td>2016188</td><td>818</td><td>0.37</td><td>0.07</td><td>0.2</td><td>0.55</td><td>0.61</td><td>0.06</td><td>0.39</td><td>0.75</td></tr><tr><td rowspan="3">14SKF</td><td>2015245</td><td>721</td><td>-0.45</td><td>0.14</td><td>-0.75</td><td>-0.12</td><td>0.13</td><td>0.1</td><td>-0.09</td><td>0.39</td></tr><tr><td>2015261</td><td>472</td><td>-0.43</td><td>0.18</td><td>-0.89</td><td>0.01</td><td>0.05</td><td>0.11</td><td>-0.22</td><td>0.37</td></tr><tr><td>2015293</td><td>721</td><td>-0.48</td><td>0.14</td><td>-0.76</td><td>-0.06</td><td>-0.24</td><td>0.11</td><td>-0.53</td><td>0.12</td></tr><tr><td>2015309</td><td>1195</td><td>-0.38</td><td>0.13</td><td>-0.65</td><td>-0.05</td><td>0.36</td><td>0.08</td><td>0.18</td><td>0.55</td></tr></tbody></table>
<table><thead><tr><th rowspan="2">Tile</th><th rowspan="2">Date of sensed image acquisition</th><th rowspan="2">Number of CPs</th><th colspan="4">Δ<sub>x</sub></th><th colspan="4">Δ<sub>y</sub></th></tr><tr><th>mean</th><th>std</th><th>min</th><th>max</th><th>mean</th><th>std</th><th>min</th><th>max</th></tr></thead><tbody><tr><td></td><td>2015325</td><td>1095</td><td>-0.38</td><td>0.14</td><td>-0.67</td><td>-0.04</td><td>-0.12</td><td>0.12</td><td>-0.49</td><td>0.19</td></tr><tr><td></td><td>2015341</td><td>703</td><td>-0.03</td><td>0.08</td><td>-0.21</td><td>0.22</td><td>0.06</td><td>0.06</td><td>-0.06</td><td>0.18</td></tr><tr><td></td><td>2015357</td><td>1100</td><td>-0.36</td><td>0.14</td><td>-0.65</td><td>0</td><td>-0.16</td><td>0.1</td><td>-0.42</td><td>0.12</td></tr><tr><td></td><td>2016024</td><td>729</td><td>-0.05</td><td>0.11</td><td>-0.28</td><td>0.26</td><td>0.11</td><td>0.09</td><td>-0.11</td><td>0.32</td></tr><tr><td></td><td>2016040</td><td>1258</td><td>-0.05</td><td>0.09</td><td>-0.23</td><td>0.19</td><td>-0.08</td><td>0.07</td><td>-0.29</td><td>0.11</td></tr><tr><td></td><td>2016056</td><td>1117</td><td>-0.39</td><td>0.12</td><td>-0.66</td><td>-0.08</td><td>0.33</td><td>0.09</td><td>0.12</td><td>0.6</td></tr><tr><td></td><td>2016072</td><td>752</td><td>-0.05</td><td>0.08</td><td>-0.24</td><td>0.18</td><td>0.19</td><td>0.06</td><td>0.02</td><td>0.34</td></tr><tr><td></td><td>2016088</td><td>1085</td><td>-0.37</td><td>0.1</td><td>-0.59</td><td>-0.08</td><td>0.14</td><td>0.08</td><td>-0.04</td><td>0.36</td></tr><tr><td></td><td>2016104</td><td>992</td><td>-0.06</td><td>0.08</td><td>-0.22</td><td>0.16</td><td>0.19</td><td>0.07</td><td>0.01</td><td>0.4</td></tr><tr><td>34UFU</td><td>2016063</td><td>636</td><td>0.61</td><td>0.07</td><td>0.46</td><td>0.78</td><td>0.15</td><td>0.1</td><td>-0.16</td><td>0.46</td></tr><tr><td></td><td>2016182</td><td>1498</td><td>0.76</td><td>0.04</td><td>0.66</td><td>0.87</td><td>0.14</td><td>0.05</td><td>0.05</td><td>0.29</td></tr></tbody></table>
<caption>Table 3.</caption>

Results of identifying CPs on the sensed (Sentinel-2A) and reference (Sentinel-2A) images using phase correlation approach at 10 m spatial resolution. Values for shifts $\Delta_x$ and $\Delta_y$ (Eq. 2) are shown in pixel units.

<table><thead><tr><th rowspan="2">Tile</th><th rowspan="2">Date of sensed image acquisition (baseline processing version, orbit relative to the orbit of reference image)</th><th rowspan="2">Number of CPs</th><th colspan="4">Δ<sub>x</sub></th><th colspan="4">Δ<sub>y</sub></th></tr><tr><th>mean</th><th>std</th><th>min</th><th>max</th><th>mean</th><th>std</th><th>min</th><th>max</th></tr></thead><tbody><tr><td rowspan="12">20HNH</td><td>2016013 (v02.01, same)</td><td>20622</td><td>-0.03</td><td>0.04</td><td>-0.16</td><td>0.08</td><td>0.04</td><td>0.04</td><td>-0.04</td><td>0.15</td></tr><tr><td>2016023 (v02.01, same)</td><td>21934</td><td>0.06</td><td>0.04</td><td>-0.04</td><td>0.16</td><td>0.05</td><td>0.04</td><td>-0.06</td><td>0.15</td></tr><tr><td>2016043 (v02.01, same)</td><td>15453</td><td>-0.03</td><td>0.06</td><td>-0.17</td><td>0.1</td><td>0.02</td><td>0.06</td><td>-0.13</td><td>0.18</td></tr><tr><td>2016063 (v02.01, same)</td><td>6133</td><td>0.12</td><td>0.11</td><td>-0.37</td><td>0.41</td><td>0.19</td><td>0.12</td><td>-0.14</td><td>0.45</td></tr><tr><td>2016073 (v02.01, same)</td><td>12919</td><td>-0.11</td><td>0.07</td><td>-0.23</td><td>0.03</td><td>0.21</td><td>0.1</td><td>-0.03</td><td>0.48</td></tr><tr><td>2016083 (v02.01, same)</td><td>13439</td><td>-0.27</td><td>0.11</td><td>-0.6</td><td>0.03</td><td>0.11</td><td>0.1</td><td>-0.18</td><td>0.42</td></tr><tr><td>2016093 (v02.01, same)</td><td>4160</td><td>-0.16</td><td>0.17</td><td>-0.61</td><td>0.31</td><td>0.35</td><td>0.31</td><td>-0.54</td><td>1.09</td></tr><tr><td>2015341 (v02.00, adjacent)</td><td>15056</td><td>0.22</td><td>0.09</td><td>0</td><td>0.43</td><td>-0.42</td><td>0.19</td><td>-0.93</td><td>0.12</td></tr><tr><td>2016006 (v02.01, adjacent)</td><td>17858</td><td>0.08</td><td>0.04</td><td>-0.01</td><td>0.16</td><td>-0.76</td><td>0.18</td><td>-1.2</td><td>-0.3</td></tr><tr><td>2016016 (v02.01, adjacent)</td><td>18265</td><td>0.24</td><td>0.08</td><td>0.05</td><td>0.44</td><td>-0.6</td><td>0.3</td><td>-1.3</td><td>0.15</td></tr><tr><td>2016026 (v02.01, adjacent)</td><td>16795</td><td>0.41</td><td>0.18</td><td>-0.12</td><td>0.92</td><td>-0.53</td><td>0.24</td><td>-1.15</td><td>0.1</td></tr><tr><td>2016036 (v02.01, adjacent)</td><td>2281</td><td>-0.04</td><td>0.12</td><td>-0.34</td><td>0.28</td><td>-0.05</td><td>0.16</td><td>-0.42</td><td>0.34</td></tr><tr><td>2016046 (v02.01, adjacent)</td><td>8566</td><td>0.3</td><td>0.14</td><td>-0.08</td><td>0.68</td><td>-0.34</td><td>0.15</td><td>-0.69</td><td>0</td></tr><tr><td>2016065 (v02.01, adjacent)</td><td>9997</td><td>0.23</td><td>0.12</td><td>-0.01</td><td>0.49</td><td>-0.36</td><td>0.23</td><td>-1.04</td><td>0.33</td></tr><tr><td>2016096 (v02.01, adjacent)</td><td>1060</td><td>-0.02</td><td>0.17</td><td>-0.47</td><td>0.47</td><td>0.02</td><td>0.16</td><td>-0.35</td><td>0.4</td></tr><tr><td rowspan="12">20HPH</td><td>2016003 (v02.01, same)</td><td>2286</td><td>0.01</td><td>0.15</td><td>-0.35</td><td>0.38</td><td>0</td><td>0.14</td><td>-0.42</td><td>0.47</td></tr><tr><td>2016013 (v02.01, same)</td><td>29887</td><td>0</td><td>0.04</td><td>-0.08</td><td>0.08</td><td>0.01</td><td>0.04</td><td>-0.08</td><td>0.11</td></tr><tr><td>2016023 (v02.01, same)</td><td>20264</td><td>0.07</td><td>0.05</td><td>-0.05</td><td>0.19</td><td>0</td><td>0.06</td><td>-0.18</td><td>0.15</td></tr><tr><td>2016043 (v02.01, same)</td><td>10197</td><td>-0.04</td><td>0.07</td><td>-0.24</td><td>0.19</td><td>-0.01</td><td>0.07</td><td>-0.19</td><td>0.21</td></tr><tr><td>2016063 (v02.01, same)</td><td>7065</td><td>0.24</td><td>0.12</td><td>-0.02</td><td>0.51</td><td>0.15</td><td>0.11</td><td>-0.18</td><td>0.45</td></tr><tr><td>2016073 (v02.01, same)</td><td>7492</td><td>-0.09</td><td>0.07</td><td>-0.27</td><td>0.12</td><td>0.03</td><td>0.07</td><td>-0.12</td><td>0.17</td></tr><tr><td>2016083 (v02.01, same)</td><td>9020</td><td>-0.24</td><td>0.08</td><td>-0.43</td><td>-0.05</td><td>-0.07</td><td>0.11</td><td>-0.35</td><td>0.22</td></tr><tr><td>2016093 (v02.01, same)</td><td>1981</td><td>-0.06</td><td>0.15</td><td>-0.39</td><td>0.25</td><td>0.05</td><td>0.17</td><td>-0.38</td><td>0.47</td></tr><tr><td>36UUU</td><td>2016109 (v02.01, same)</td><td>18989</td><td>0.18</td><td>0.12</td><td>-0.16</td><td>0.47</td><td>0.2</td><td>0.24</td><td>-0.5</td><td>0.93</td></tr><tr><td>2016119 (v02.01, same)</td><td>21363</td><td>0.08</td><td>0.11</td><td>-0.26</td><td>0.45</td><td>-0.17</td><td>0.29</td><td>-1.16</td><td>0.61</td></tr><tr><td>2016179 (v02.04, same)</td><td>19300</td><td>0</td><td>0.04</td><td>-0.08</td><td>0.08</td><td>0.21</td><td>0.05</td><td>0.08</td><td>0.33</td></tr><tr><td>2016199 (v02.04, same)</td><td>16647</td><td>-0.07</td><td>0.06</td><td>-0.24</td><td>0.08</td><td>0.14</td><td>0.06</td><td>-0.03</td><td>0.3</td></tr><tr><td>2016096 (v02.01, adjacent)</td><td>10046</td><td>-0.1</td><td>0.13</td><td>-0.47</td><td>0.29</td><td>1.44</td><td>0.51</td><td>-0.13</td><td>2.82</td></tr><tr><td>2016156 (v02.02, adjacent)</td><td>8145</td><td>-0.32</td><td>0.13</td><td>-0.68</td><td>0.1</td><td>1.19</td><td>0.4</td><td>-0.04</td><td>2.02</td></tr><tr><td>2016166 (v02.02, adjacent)</td><td>10436</td><td>0.01</td><td>0.04</td><td>-0.09</td><td>0.14</td><td>0.13</td><td>0.05</td><td>-0.02</td><td>0.24</td></tr><tr><td>2016196 (v02.04, adjacent)</td><td>8530</td><td>-0.06</td><td>0.07</td><td>-0.22</td><td>0.08</td><td>0.25</td><td>0.08</td><td>0.01</td><td>0.47</td></tr><tr><td>2016206 (v02.04, adjacent)</td><td>4905</td><td>-0.17</td><td>0.09</td><td>-0.38</td><td>0.07</td><td>0</td><td>0.08</td><td>-0.23</td><td>0.2</td></tr><tr><td>14SKF</td><td>2016042 (v02.01, same)</td><td>28459</td><td>-0.17</td><td>0.04</td><td>-0.27</td><td>-0.05</td><td>0.1</td><td>0.06</td><td>-0.02</td><td>0.22</td></tr></tbody></table>
<table><thead><tr><th rowspan="2">Tile</th><th rowspan="2">Date of sensed image acquisition (baseline processing version, orbit relative to the orbit of reference image)</th><th rowspan="2">Number of CPs</th><th colspan="4">Δ<sub>x</sub></th><th colspan="4">Δ<sub>y</sub></th></tr><tr><th>mean</th><th>std</th><th>min</th><th>max</th><th>mean</th><th>std</th><th>min</th><th>max</th></tr></thead><tbody><tr><td></td><td>2016072 (v02.01, same)</td><td>19183</td><td>-0.37</td><td>0.09</td><td>-0.56</td><td>-0.13</td><td>0.21</td><td>0.09</td><td>-0.08</td><td>0.45</td></tr><tr><td></td><td>2016132 (v02.01, same)</td><td>18091</td><td>-0.26</td><td>0.1</td><td>-0.58</td><td>0.05</td><td>0.35</td><td>0.17</td><td>-0.22</td><td>0.94</td></tr><tr><td>34UFU</td><td>2016048 (v02.01, same)</td><td>16528</td><td>0.28</td><td>0.08</td><td>0.06</td><td>0.55</td><td>0.05</td><td>0.18</td><td>-0.37</td><td>0.61</td></tr><tr><td></td><td>2016208 (v02.04, same)</td><td>26453</td><td>0.16</td><td>0.03</td><td>0.06</td><td>0.25</td><td>0.22</td><td>0.04</td><td>0.08</td><td>0.33</td></tr></tbody></table>
<caption>Table 4.</caption>

Average and standard deviation of the RMSE error (Eq. 4) calculated for different transformation functions using CPs from testing set when co-registering Landsat-8 and Sentinel-2A images. RMSE values are shown in pixel units at 30 m spatial resolution.

<table><thead><tr><th rowspan="2">Tile</th><th colspan="2">Translation (Eq. 7–8)</th><th colspan="2">1<sup>st</sup> order polynomial (Eq. 9–10)</th><th colspan="2">Gaussian RBFs (Eq. 11–12, 13)</th><th colspan="2">RF regression</th><th colspan="2">TPS (Eq. 11–12, 14)</th></tr><tr><th>Mean</th><th>Standard deviation</th><th>Mean</th><th>Standard deviation</th><th>Mean</th><th>Standard deviation</th><th>Mean</th><th>Standard deviation</th><th>Mean</th><th>Standard deviation</th></tr></thead><tbody><tr><td>20HNH</td><td>0.119</td><td>0.031</td><td>0.091</td><td>0.026</td><td>0.093</td><td>0.027</td><td>0.084</td><td>0.024</td><td>0.090</td><td>0.025</td></tr><tr><td>20HPH</td><td>0.123</td><td>0.014</td><td>0.078</td><td>0.016</td><td>0.081</td><td>0.017</td><td>0.073</td><td>0.018</td><td>0.079</td><td>0.018</td></tr><tr><td>36UUU</td><td>0.108</td><td>0.011</td><td>0.072</td><td>0.015</td><td>0.074</td><td>0.015</td><td>0.059</td><td>0.014</td><td>0.073</td><td>0.015</td></tr><tr><td>14SKF</td><td>0.145</td><td>0.037</td><td>0.094</td><td>0.018</td><td>0.095</td><td>0.018</td><td>0.074</td><td>0.018</td><td>0.094</td><td>0.017</td></tr><tr><td>34UFU</td><td>0.095</td><td>0.045</td><td>0.056</td><td>0.034</td><td>0.060</td><td>0.038</td><td>0.044</td><td>0.030</td><td>0.054</td><td>0.033</td></tr></tbody></table>
<caption>Table 5.</caption>

Average and standard deviation of the RMSE error (Eq. 4) calculated for different transformation functions using CPs from testing set when co-registering multi-temporal Sentinel-2A images from the same orbit. RMSE values are shown in pixel units at 10 m spatial resolution.

<table><thead><tr><th rowspan="2">Tile</th><th colspan="2">Translation (Eq. 7–8)</th><th colspan="2">1<sup>st</sup> order polynomial (Eq. 9–10)</th><th colspan="2">Gaussian RBFs (Eq. 11–12, 13)</th><th colspan="2">RF regression</th><th colspan="2">TPS (Eq. 11–12, 14)</th></tr><tr><th>Mean</th><th>Standard deviation</th><th>Mean</th><th>Standard deviation</th><th>Mean</th><th>Standard deviation</th><th>Mean</th><th>Standard deviation</th><th>Mean</th><th>Standard deviation</th></tr></thead><tbody><tr><td>20HNH</td><td>0.141</td><td>0.104</td><td>0.125</td><td>0.076</td><td>0.125</td><td>0.076</td><td>0.105</td><td>0.060</td><td>0.126</td><td>0.074</td></tr><tr><td>20HPH</td><td>0.133</td><td>0.064</td><td>0.128</td><td>0.062</td><td>0.127</td><td>0.062</td><td>0.114</td><td>0.059</td><td>0.129</td><td>0.064</td></tr><tr><td>36UUU</td><td>0.181</td><td>0.126</td><td>0.114</td><td>0.046</td><td>0.114</td><td>0.046</td><td>0.088</td><td>0.035</td><td>0.112</td><td>0.048</td></tr><tr><td>14SKF</td><td>0.133</td><td>0.066</td><td>0.123</td><td>0.050</td><td>0.118</td><td>0.046</td><td>0.089</td><td>0.036</td><td>0.123</td><td>0.051</td></tr><tr><td>34UFU</td><td>0.122</td><td>0.101</td><td>0.092</td><td>0.086</td><td>0.091</td><td>0.084</td><td>0.066</td><td>0.059</td><td>0.093</td><td>0.088</td></tr></tbody></table>
<caption>Table 6.</caption>

The same as Table 5, but for adjacent Sentinel-2A orbits.

<table><thead><tr><th rowspan="2">Tile</th><th colspan="2">Translation (Eq. 7–8)</th><th colspan="2">1<sup>st</sup> order polynomial (Eq. 9–10)</th><th colspan="2">Gaussian RBFs (Eq. 11–12, 13)</th><th colspan="2">RF regression</th><th colspan="2">TPS (Eq. 11–12, 14)</th></tr><tr><th>Mean</th><th>Standard deviation</th><th>Mean</th><th>Standard deviation</th><th>Mean</th><th>Standard deviation</th><th>Mean</th><th>Standard deviation</th><th>Mean</th><th>Standard deviation</th></tr></thead><tbody><tr><td>20HNH</td><td>0.239</td><td>0.048</td><td>0.207</td><td>0.031</td><td>0.202</td><td>0.028</td><td>0.164</td><td>0.032</td><td>0.205</td><td>0.031</td></tr><tr><td>36UUU</td><td>0.248</td><td>0.212</td><td>0.191</td><td>0.139</td><td>0.189</td><td>0.138</td><td>0.138</td><td>0.087</td><td>0.193</td><td>0.142</td></tr></tbody></table>
<caption>**Table 7.**</caption>

Typical processing times (in s) for the main steps of the proposed approach.

<table><thead><tr><td>Processing step</td><td>SenSen10</td><td>LandSen30</td></tr></thead><tbody><tr><td>CPs identification</td><td>135</td><td>22</td></tr><tr><td>Precomputing shifts with 1<sup>st</sup> order polynomial function</td><td>5</td><td>0.5</td></tr><tr><td>Precomputing shifts with the RF regression</td><td>220</td><td>8</td></tr><tr><td>Warping (single scene)</td><td>210</td><td>25</td></tr></tbody></table>