# Temporal Attention Multi-Resolution Fusion of Satellite Image Time-Series, applied to Landsat-8/9 and Sentinel-2: all bands, any time, at best spatial resolution

Julien Michel, Jordi Inglada

► To cite this version:

Julien Michel, Jordi Inglada. Temporal Attention Multi-Resolution Fusion of Satellite Image Time-Series, applied to Landsat-8/9 and Sentinel-2: all bands, any time, at best spatial resolution. 2025. [https://hal.science/hal-05101526v3](https://hal.science/hal-05101526v3)

HAL Id: hal-05101526

https://hal.science/hal-05101526v3

Preprint submitted on 13 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Highlights

Temporal Attention Multi-Resolution Fusion of Satellite Image Time-Series, applied to Landsat-8/9 and Sentinel-2: all bands, any time, at best spatial resolution

Julien Michel, Jordi Inglada

* A general formulation of the fusion of Satellite Image Time Series (SITS) is proposed

* A new training strategy and Neural Network architecture solves this general problem

* This solution avoids several unrealistic assumptions found in the literature

* The pre-trained model is evaluated on 4 usual tasks using Landsat and Sentinel-2 SITS

* This single model is on-par or better than ad-hoc existing models, with more benefits
Temporal Attention Multi-Resolution Fusion of Satellite
Image Time-Series, applied to Landsat-8/9 and Sentinel-2:
all bands, any time, at best spatial resolution

Julien Michel<sup>a,*</sup>, Jordi Inglada<sup>a</sup>

<sup>a</sup>Univ Toulouse, CNES, CNRS, INRAE, IRD, CESBIO, Toulouse, France, 18 avenue Edouard Belin, BPI
2801, TOULOUSE Cedex 9, 31401, France

**Abstract**

This paper introduces a general formulation for the fusion of Satellite Image Time Se-
ries (SITS) of variable length from several sensors at different spatial resolutions and
acquisition times over the same geographical area. In this formulation, all the spectral
bands from all the input sensors are predicted at the best input spatial resolution, and at
any observed or non-observed acquisition time requested. To address this general prob-
lem, an advanced Masked Auto-Encoder training strategy is proposed, utilising two
new loss functions: a Linear-Regression Learned Perceptual Image Similarity term to
favor high spatial frequency details, and a mask-contrastive term to ignore clouds and
other non-informative areas in the input data. This strategy is applied to the training of
Temporal Attention Multi-Resolution Fusion of Satellite Image Time-Series (TAMRF-
SITS), a novel Deep Learning architecture designed to implement the proposed gen-
eral formulation. Experiments with joint Landsat-8/9 and Sentinel-2 time-series were
conducted on four different tasks from the literature and demonstrate that a single pre-
trained TAMRF is on par or better than existing ad-hoc methods. For instance, TAMRF
provides a gain of 0.01 surface reflectance Root Mean Square Error on a Spatio-
Temporal Fusion task when compared to competing algorithms, while showing the
highest spatial frequency content. Moreover, the proposed method relaxes unrealistic
assumptions routinely found in the literature, including: same or similar spectral bands
in different sensors, same-day acquisitions, and scale-invariance of the relationship be-

<footnote>*Corresponding author</footnote>
tween high and low resolution images. To the best of our knowledge, our method
is the first to achieve this range of capabilities with a single model, without making
any of these assumptions. The complete source code for training and experiments is
available here: https://github.com/Evoland-Land-Monitoring-Evolution/
tamrsits.

Keywords: Super-Resolution, Spatio-Temporal Fusion, Sharpening, Gap-Filling, Sentinel-2, Landsat, Transformer

1. Introduction

Global coverage of the Earth surface with Satellite Image Time Series (SITS) is a critical component of the continuous monitoring of our planet Essential Climate Variables Bojinski et al. (2014) and Essential Biodiversity Variables Jetz et al. (2019). The Landsat series alone provides 50 years of such constant monitoring Wulder et al. (2022), and for ten years now Sentinel-2 has complemented the Landsat archive with higher revisit, higher spatial resolution images. SITS are used to monitor vegetation Misra et al. (2020), derive Land Cover and Land Use maps Phiri et al. (2020), control the European Common Agricultural Policy Vajsová et al. (2020) or monitor Water Bodies Bhangale et al. (2020).

Spatial and temporal resolutions of SITS are important parameters that define the
reachable scale of analysis both in space and time. The revisit time is even more critical
for passive optical sensors, where cloud occurrences obliterate a vast proportion of the
observations Wilson and Jetz (2016). On the other hand, there is an increasing num-
ber of orbiting sensors that could complement each other to increase the revisit time
and the spatial resolution. In this paper, we consider the case of two or more contem-
porary sensors with global coverage and regular revisit, complementary revisit cycles
and spectral bands, and different spatial resolutions. There are a few combinations
of such sensors already flying, the most popular being Landsat-8/9 with Sentinel-2,
because of their complementary revisit cycles Li and Chen (2020). Another already
available combination is PlanetScope and Sentinel-2 Latte and Lejeune (2020); Sadeh
et al. (2021). In a near future, possible combinations will include the Trishna mission
Lagouarde et al. (2018) and the Land Surface Temperature Monitoring (LSTM) mission Bernard et al. (2023) with Sentinel-2 or Landsat-8/9, as well as Sentinel-2 NG or Landsat-Next. There is therefore a growing need for methods that can jointly leverage those continuous sources of monitoring of the Earth. Ideally one would like to transform all observations from all sensors into measurements as seen by an ubiquitous sensor: all wavelengths, any time, at best spatial resolution.

## 1.1. Existing works

In the literature, this general problem is partly addressed by several families of sub-problems, with different hypotheses and constraints. Those sub-problems and the methods that have been proposed to solve them are detailed in section 1.1.1. Another corpus of work focuses on producing harmonized datasets for identified sets of sensors. To that end, they propose end-to-end processing frameworks including pre-processing and sensor-specific features, in addition to leveraging methods from the above literature. They are summarized in section 1.1.2.

### 1.1.1. Taxonomy of methods

The relevant methodological subdomains are summarized in Fig. 1. They include Single and Multi-Image Super-Resolution, Band-Sharpening, Spatio-Temporal Fusion and Temporal Modeling, all of which are detailed in the following sections. Those works usually focus on machine learning methodological developments within the constraints of existing datasets or subdomain definitions. This proposed taxonomy is complementary to existing reviews such as Wang et al. (2023), where spatial, spectral and temporal reconstructions are analysed independently.

**Single and Multi Image Super-Resolution.** In Single or Multi-Images Super-Resolution (SISR or MISR) Anwar et al. (2020); Liu et al. (2022), images from a higher resolution sensor are used to learn a model improving the spatial resolution of one or several images acquired by a lower resolution sensor. SISR and MISR are inherently ill-posed, since there can be many High-Resolution (HR) images corresponding to the same observed Low-Resolution (LR) image(s). For this reason, successful SISR and MISR usually rely on adversarial training Wang et al. (2018); Ledig et al. (2017) to capture
<figure><img src="image_2.png" /><figcaption>Figure 1: Overview of capabilities of each relevant methodological subdomain, as compared to our proposed method. Input SITS are represented along a time-line showing the acquisition times. The Low Resolution (LR) SITS are represented by small green data-cubes, while the High Resolution (HR) SITS are represented by large blue data-cubes. Single Image Super-Resolution predicts the LR SITS bands at the same acquisition time, but at better spatial resolution. Multi-Image Super-Resolution predicts the LR SITS bands at a single acquisition time but at better spatial resolution, using multiple acquisitions in the input LR SITS. Temporal modelling provides prediction of the LR SITS bands at unobserved acquisition times, without changing the spatial resolution. Band Sharpening provides LR SITS bands at better spatial resolution and at same acquisition time, using simultaneous HR SITS as auxiliary input. Spatio-Temporal Fusion (STF) leverages simultaneous observations of LR SITS and HR SITS in order to predict LR bands at better spatial resolution, at acquisition dates when only LR bands are observed (shaded inputs represent variants of STF). Finally, our proposed method processes joint HR and LR SITS and predicts all HR bands (blue output data-cubes) and LR bands (green output data-cubes) at the HR resolution for any acquisition time.</figcaption></figure>
the target data manifold structure during training. Both SISR Pouliot et al. (2018); Galar et al. (2019); Romero et al. (2020); Nguyen et al. (2023) and MISR Märtens et al. (2019); Molini et al. (2020); Okabayashi et al. (2024); Ibrahim et al. (2025) have been applied to remote-sensing imagery. The models are usually trained by using either simulated or cross-sensor datasets, each having their own strengths and weaknesses, as analyzed in a previous work Michel et al. (2025). Though vastly investigated in the literature, SISR and MISR are not well suited for leveraging the synergies between two global coverage sensors of different spatial resolutions, since they do not integrate the prior HR information given by HR observations at other dates, which always exists in the case of global coverage missions. SISR can not overcome missing data in input images because it lacks temporal context. Finally, both SISR and MISR require HR reference data with similar spectral bands and even same-date acquisitions in case of cross-sensor datasets. This explains why joint Super-Resolution and sensor translation is not represented in Fig. 1: if a different HR sensor is often used to train the super-resolution model, authors usually try to avoid learning sensor-to-sensor differences Michel et al. (2025), and consider that the output of the super-resolution model should be the same sensor as the input with a better spatial resolution. An exception to this is the work of Sambandham et al. Sambandham et al. (2024), which is covered in section 1.1.2.

**Band-Sharpening.** Band-Sharpening (BS) Firozjaei et al. (2022), also known as pan-sharpening Vivone et al. (2015); Ciotola et al. (2025), or reference-based Super-Resolution Lutio et al. (2019); Su et al. (2025) depending on the scientific community, aims at using an auxiliary near-simultaneous High Resolution image (often from the same satellite) to improve the spatial resolution of an image acquired by a lower resolution sensor. As opposed to SISR, where HR images are used as a target in the model supervised training, BS uses HR images as auxiliary inputs to the model. While adding HR inputs certainly helps the model to improve the spatial resolution, those HR inputs can no longer serve as references to train the model against. To overcome this problem, researchers in various remote sensing fields have been using the scale invariance hypothesis, enabling the use of Wald's protocol Wald et al. (1997), which is well il-
lustrated in Fig. 2 of Palsson et al. (2018). It consists in down-scaling by a factor *k* both input and reference data, training the model to perform *r/k* → *r* spatial resolution improvement, and then applying the trained model to perform *r* → *kr* resolution improvement by assuming scale invariance of the model. This practice is commonplace in thermal sharpening, a branch of BS addressing the Thermal Infra Red or Land Surface Temperature (LST) images Gao et al. (2012); Granero-Belinchon et al. (2019), but it is also used in other band-sharpening applications when sensor is equipped with different spatial resolutions across bands Palsson et al. (2018); Salgueiro et al. (2021).

Several studies show that the scale invariance hypothesis is unrealistic and can lead to wrong estimation of local textures Merlin et al. (2010); Nguyen et al. (2022). In thermal sharpening, most methods only learn the mapping from degraded HR images to LR images, and then use the learnt model at High Resolution. During inference, the model only sees the HR image and can not use local LR observations. This imposes to learn a new model for each pair of HR and LR mapping. In order to make sure that predicted dynamic is faithful to the LR observation, residuals between down-sampled predictions and LR observations are often injected to form the final predictions. Nevertheless, the residual signal is inherently an LR signal and may hinder the resolution improvement achieved by the model.

BS seems more equipped that SISR and MISR to solve the problem at stake: it supports LR bands without HR reference, and in some cases might produce gap-less outputs, if several HR acquisitions are used to circumvent clouds occurrences Shao et al. (2019). Yet it relies on the unrealistic scale invariance hypothesis, and is essentially mono-temporal: it can not extrapolate unseen dates or learn the temporal dynamics.

*Spatio-Temporal Fusion.* Spatio-Temporal Fusion (STF) Belgiu and Stein (2019); Xiao et al. (2023) is a set of methods derived from the initial Spatial and Temporal Adaptive Fusion Model (STARFM) Gao et al. (2006). STARFM models the relationship between concomitant Landsat and MODIS surface reflectances at the same location, and extrapolates it to a nearby date where only MODIS is observed. It requires at least one pair of LR and HR simultaneous images close to this target date. Model-based improvement have followed Belgiu and Stein (2019); Zhu et al. (2018), most of which
extend STARFM or try to overcome its limitations by using more complex modelling,
without changing or extending this paired configuration. More recently, the rule-based
method has been traded for Deep-Learning (DL) models, including convolutional lay-
ers Tan et al. (2018, 2019); Liu et al. (2019); Tan et al. (2022), adversarial training
Zhang et al. (024b); Xie et al. (2024), transformers Yang et al. (2022); Li et al. (022b);
Wu and Huang (2024) and vision transformers Chen et al. (2022). It must be stressed
that though the DL models proposed in the literature became more complex, very few
of them break free from the limitation of requiring one or two pairs of simultaneous
LR and HR acquisitions. Notable exceptions are Goyena et al. (2023) and Zhang et al.
(024c) where unpaired STF is investigated. They also do not extend the model tempo-
ral capability beyond two dates. Additionally, several of these methods actually require
the scale invariance hypothesis for training, and some of them also require retraining
for each target date. Finally, it should be stressed that the model is usually not in-
formed of the acquisition dates, which prevents the learning of temporal modelling and
the generalization to other time intervals that were not seen during training. Note that
STF is sometimes called Spatio-Temporal Spectral Fusion when it is combined with
BS in order to process HR and LR images with different spectral bands Shen et al.
(2016).

*Temporal Modeling.* Temporal modeling refers to the interpolation or forecasting of SITS Moskolaï et al. (2021). One strong driver of such research is to mitigate the occurrence of clouds in passive optical data in order to generate gap-less SITS Li et al. (2021), a process also known as gap-filling. Most methods handle SITS from a single sensor, with some notable exceptions Roy et al. (2008). U-TILISE Stucker et al. (2023) combines a UNet architecture with a temporal attention mechanism to perform gap-filling and temporal interpolation. Liu et al. (2024) propose to use a temporal transformer architecture and a masked SITS training scheme to provide gap-free series from the Harmonized Landsat and Sentinel-2 dataset (HLS) dataset. SatFlow Iri-gireddy and Bandaru (2025) uses diffusion models conditioned by acquisition dates in order to reconstruct gap free Landsat SITS. RESTORE-DiT Shu et al. (2025) also uses diffusion models, combined with attention blocks in order to condition the temporal
reconstruction of Sentinel-2 to SAR observations from Sentinel-1. The obvious limita-
tion of temporal modeling is that it does not address the spatial and spectral part of the
problem at stake.

1.1.2. *Producing harmonized datasets*

In addition to those methodological subdomains, a corpus of works aims at lever-
aging those methods to address the broader problem of producing harmonized datasets
from SITS of different sensors. This includes data pre- and post-processing as well as
production strategies.

Some of them are in operational stage, such as the HLS Claverie et al. (2018), where differences between sensors and acquisition conditions are carefully processed into a consistent dataset. Yet HLS does not improve the spatial resolution of Landsat observation, and instead provides matching 30-meter Sentinel-2 observations. A competing operational method aiming at the same goals is Sen2Like Saunier et al. (2022) where an additional fusion step using a simple residual compensation yields 10-meter Landsat predictions.

Aside from HLS and Sen2Like, other more prospective works include Wang et al. (2017), in which the authors propose to use Area-To-Point Regression Kriging (ATPRK), a standard STF method based on geo-statistics, to produce 10-meter Landsat images of corresponding Sentinel-2 bands at times when only Landsat is observed. In Shao et al. (2019), the authors employ the same idea but replace ATPRK with the Super Resolution Convolutional Neural Network Dong et al. (2016), relying on the scale invariance hypothesis and Wald's protocol for training and evaluation. This idea is also used in Latte and Lejeune (2020) to fuse Sentinel-2 and PlanetScope imagery into a 2.5-meter Sentinel-2 product, relying again on the scale invariance hypothesis. Additionally, their methodology requires retraining the model for each Sentinel-2 image because of the varying number of input PlanetScope images. The combination of PlanetScope with Sentinel-2 is also explored in Sadeh et al. (2021), using simple rule-based fusion and linear regression in the context of producing Leaf Area Index high resolution time-series. In Luo et al. (2018) the authors proposed STAIR, later refined into STAIR 2.0 Luo et al. (2020), which is a complete framework for the fu-
sion of SITS from multiple sensors sharing similar spectral bands. The methodology leverages same-day observations to form a series of differences between HR and spatially up-sampled LR. Those differences are then interpolated at acquisition times when only LR is observed, and added to the up-sampled LR observations. STAIR includes many pre-processing steps, including spatial registration, spectral adjustment as well as correction of missing data in Landsat due to faulty detector arrays. A similar idea to STAIR is explored in Wang et al. (2024) where high resolution details are obtained by means of the Smoothing-Sharpening Image Filter. In Mukherjee and Liu (2023), GAN-based SISR is employed to directly predict High Resolution Sentinel-2 bands from low resolution Landsat bands. The same idea is applied in Chen et al. (2021) for the super-resolution of historical Landsat data, using the overlapping period with Sentinel-2 for training, and in Sambandham et al. (2024), albeit without adversarial training. In Chang et al. (2025), the authors leverage the Enhanced Deep Convolutional Spatio-Temporal Fusion Network Tan et al. (2019) in order to fuse Landsat and MODIS data. Most of those works require a substantial amount of cross-calibration, assuming that spectral bands are similar and that same-day acquisitions exist. Most of them have limited temporal modelling, and are unable to produce gap-less output or extrapolate unseen dates for instance.

## 1.2. Contributions

The previous section shows that the different existing methods only address a sub-part of the following more general problem: predict all bands from all-sensors at best observed spatial resolution and for any acquisition time, given the observed SITS from two or more sensors over the same area, without any further assumption. We hypothesize that failing to address this big picture not only restricts the capabilities of those methods, but also requires additional unrealistic assumptions and limitations to solve the sub-part of the problem at stake. Capabilities, limitations and unrealistic assumptions of existing methods are summarized in Table 1.

In this paper, we first introduce a mathematical formulation of the general problem in section 2.1. We then propose a deep-learning architecture (section 2.2) and a training procedure (section 2.3) that allow to solve it by means of Self Supervised Learning
(SSL). Our proposed model, called Time Attention Multi-Resolution Fusion of Satellite Image Time Series (TAMRF-SITS or TAMRF for short) has full capabilities and none of the limitations and unrealistic assumptions of the existing families of methods. Fig. 2 illustrates some of those capabilities on our testing dataset which will be introduced in section 3.1. Though the proposed model is not limited to a specific pair of sensors and can even handle more than two sensors, we demonstrate its properties through the fusion of Sentinel-2 and Landsat-8/9 time-series in section 3. Finally, section 4 proposes several ablation and sensitivity studies allowing to better understand the drivers for TAMRF performances.

<table><caption>Table 1: Capabilities and limitations of the different categories of methods. Temporal refers to Temporal Modeling (see section 1.1.1), SR refers to Super-Resolution (see section 1.1.1), BS refers to Band- Sharpening (see section 3.3.2), and STF refers to Spatio-Temporal Fusion (see section 3.3.3). Parents (✓) indicate that some methods in the literature partially achieve the capability or solve the limitation or unrealistic assumption.</caption><thead><tr><th>Category of methods</th><th>Temporal</th><th>SR</th><th>BS</th><th>STF</th><th>Ours</th></tr></thead><tbody><tr><td><strong>Capability</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Improves spatial resolution</td><td></td><td>✓</td><td>✓</td><td>✓</td><td>✓</td></tr><tr><td>Supports bands without HR reference</td><td></td><td></td><td>✓</td><td>✓</td><td>✓</td></tr><tr><td>Produces gap-less outputs</td><td>✓</td><td></td><td>(✓)</td><td>(✓)</td><td>✓</td></tr><tr><td>Extrapolates unseen dates</td><td>✓</td><td></td><td></td><td></td><td>✓</td></tr><tr><td>Learns temporal dynamic</td><td>✓</td><td></td><td></td><td></td><td>✓</td></tr><tr><td><strong>Limitation solved</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Does not require similar bands in sensors</td><td></td><td></td><td>✓</td><td></td><td>✓</td></tr><tr><td>Does not require same-day acquisitions</td><td>✓</td><td></td><td>(✓)</td><td>(✓)</td><td>✓</td></tr><tr><td>Does not use scale invariance hypothesis</td><td>✓</td><td>✓</td><td></td><td></td><td>✓</td></tr><tr><td>Is not limited to pairs or triplets</td><td>✓</td><td>✓</td><td></td><td></td><td>✓</td></tr><tr><td>Does not require no-data masks as inputs</td><td>(✓)</td><td></td><td></td><td></td><td>✓</td></tr><tr><td>Does not require training for each target date</td><td>✓</td><td>✓</td><td>(✓)</td><td>(✓)</td><td>✓</td></tr></tbody></table>
<figure><img src="image_3.png" /><figcaption>Figure 2: Illustration of the TAMRF model capabilities over Area Of Interest 31TFJ for year 2022 (from the test set, see section 3.1). The model receives 25 Landsat-8/9 observations and 40 Sentinel-2 observations as inputs (highlighted in green), and is asked to predict Landsat-8/9 and Sentinel-2 at different query dates (highlighted in red). For some query dates, additional Sentinel-2 or Landsat-8/9 images (e.g. not part of the model inputs) serve as reference data to compare the prediction with (highlighted in blue). Depending on the availability of clear or cloudy Sentinel-2 or Landsat-8/9 image in the input, several cases are presented. From top to bottom: row 1 shows Landsat-8/9 prediction for a clear Landsat-8/9 input date, row 2 shows Sentinel-2 predictions for a clear Sentinel-2 input date, row 3 shows Sentinel-2 prediction compared to true Sentinel-2 image for a clear Landsat-8/9 input date, row 4 shows Sentinel-2 and Landsat-8/9 prediction for a cloudy Sentinel-2 input date, and row 5 shows Sentinel-2 prediction compared to true Sentinel-2 image for a date when neither Landsat-8/9 nor Sentinel-2 are seen by the model. Note that the model does not use cloud masks. Color compositions are: RGB (B04, B03, B02 for Sentinel-2, B4, B3, B2 for Landsat-8/9), Red Edge (B07, B06, B05 for Sentinel-2), SWIR (B8A, B11, B12 for Sentinel-2), and LST (B10 for Landsat-8/9).</figcaption></figure>
## 2. Proposed method

### 2.1. Problem Formulation

For the sake of simplicity, the problem formulation is presented for two SITS, one with high spatial resolution R called $S_{HR}$ and one with low spatial resolution r called $S_{LR}$. Note that this formulation can easily be extended to any number of SITS. k represents the resolution ratio between the two SITS:

$$r = kR, k \in \mathbb{R}. \qquad (1)$$

The spatial sampling grid of the LR SITS can be defined as $R_{LR} = R(x_0, y_0, h, w, r)$, with:

$$R(x_0, y_0, h, w, r) = \{x_0 + (i + 0.5)r, i \in 0, \dots, w - 1\} \times \{y_0 + (j + 0.5)r, j \in 0, \dots, h - 1\}, (2)$$

where $x_0$ and $y_0$ denote the coordinates of the upper-left corner of the grid, and $h \in \mathbb{N}$ and $w \in \mathbb{N}$ denote the integer height and width of the image. Conversely, the spatial sampling grid of the HR SITS can be written $R_{HR} = R(x_0, y_0, H, W, R)$, with $H = [kh]$ and $W = [kw]$, where $\lceil \cdot \rceil$ represents the ceiling integer. The acquisition times of the LR SITS can be defined as:

$$T_{LR} = \{t_n, n \in 0, \dots, N_{LR} - 1\}, \qquad (3)$$

where $N_{LR} \in \mathbb{N}$ is the number of acquisitions in LR SITS. Acquisition times for HR SITS are likewise noted $T_{HR}$. Finally, $C_{LR}$ denotes the set of spectral bands in LR SITS, and $C_{HR}$ denotes the set of spectral bands in HR SITS, without any further assumption on possible matching bands between $C_{LR}$ and $C_{HR}$. With those notations, both the HR and LR SITS can be seen as a collection of measurements sampled from a universal observation function $F(x, y, r, t, c)$ of spatial location $(x, y)$, spatial resolution $r$, acquisition time $t$ and spectral band $c$:

$$S_{LR} = \{F(x, y, r, t, c) \text{ with } (x, y), t, c \in R_{LR} \times T_{LR} \times C_{LR}\}, \qquad (4)$$

$$S_{HR} = \{F(x, y, r, t, c) \text{ with } (x, y), t, c \in R_{HR} \times T_{HR} \times C_{HR}\}. \qquad (5)$$
For the sake of computation as well as for enforcing spatial and temporal inductive biases, collections of measurements $S_{LR}$ and $S_{HR}$ are usually organized into tensors of shape $[w, h, N_{LR}, \#C_{LR}]$ (resp. $[W, H, N_{HR}, \#C_{HR}]$), where $\#C_{LR}$ (resp. $\#C_{HR}$) denotes the number of elements in $C_{LR}$ (resp. $C_{HR}$). Let $T_{query}$ denote any other set of acquisition times. In this paper, we aim at building a parametric function $\Phi(S_{LR}, S_{HR}, T_{query} | \Theta)$ as well as a training procedure to derive its optimal parameters $\Theta^*$ so that $\Phi$ is able to infer measurements for all bands in $C_{all} = C_{HR} \cup C_{LR}$ at every query acquisition time in $T_{query}$ and sampled on high spatial resolution sampling grid defined by $R_{HR}$:

$$ \Phi(S_{LR}, S_{HR}, T_{query} | \Theta^*) \approx \{F(x, y, r, t, c) \text{ with } (x, y), t, c \in R_{HR} \times T_{query} \times C_{all}\}. \quad (6) $$

Note that the output of $\Phi$ can be further separated into two SITS with $C_{HR}$ and $C_{LR}$ channels:

$$ \Phi(S_{LR}, S_{HR}, T_{query} | \Theta^*) = (S_{LR}^*, S_{HR}^*), \quad (7) $$

with

$$ S_{LR}^* \approx \{F(x, y, r, t, c) \text{ with } (x, y), t, c \in R_{HR} \times T_{query} \times C_{LR}\}, \quad (8) $$

$$ S_{HR}^* \approx \{F(x, y, r, t, c) \text{ with } (x, y), t, c \in R_{HR} \times T_{query} \times C_{HR}\}. \quad (9) $$

## 2.2. Architecture

In this work, function $\Phi$ is implemented through a deep-learning architecture called Temporal Attention Multi-Resolution Fusion of Satellite Image Time Series (TAMRF-SITS or TAMRF for short), leveraging residual Convolutional Neural Networks (CNN) in the spatial dimensions and a Transformer in the temporal dimension. The overall architecture follows a classical encoder - decoder scheme, as presented in Fig. 3. Both $S_{LR}$ and $S_{HR}$ are fed into the TAMRF encoder which outputs a latent representation $S_{latent}$ comprising $D$ features, with shape $[W, H, \#T_{LR} + \#T_{HR}, D]$. It is important to note that unlike many works in auto-encoders, the time dimension is not projected to a fixed grid, and retains the initial number of observations $\#T_{LR} + \#T_{HR}$. To obtain predictions for $T_{query}$ acquisition times, both $T_{query}$ and latent representation $S_{latent}$
are fed into the TAMRF decoder, which outputs $S_{LR}^*$ and $S_{HR}^*$. In the following sections, the architecture of the encoder (section 2.2.1) and the decoder (section 2.2.2) will be detailed. Intuitively, TAMRF can be seen as the combination of a Single Image Super-Resolution network in the spatial dimension and a Transformer in the temporal dimension. Note that simpler alternatives such as Temporal Convolutional Neural Networks (TempCNN) or Long Term Short Term Memory (LSTM) could be considered to process the temporal dimension, with the addition of positional encoding to handle irregular sampling. However, additional tweaking of those architectures would be required to handle a fully generic $T_{query}$ in the decoder (see section 2.2.2), which on the other hand is natively supported by the Transformer architecture. Besides, TempCNN and LSTM are known for their limitations in learning long range relationships, as well as for difficulties in training for LSTM due to its recurrent nature. Nevertheless, the TAMRF architecture is modular and it would be possible to replace the temporal or spatial encoders and decoders with equivalent functional blocks while retaining the benefits of the global architecture presented in this section, as well as the training strategy presented in section 2.3, in solving the ubiquitous task of Eq. 6.

### 2.2.1. Encoder

The workflow of the TAMRF encoder is detailed in Fig. 4 and is composed of two sequential stages. First, each acquisition date of $S_{LR}$ is processed independently by a dedicated spatial LR encoder. The same applies to the acquisition dates of $S_{HR}$, which are processed independently by a dedicated spatial HR encoder. The resulting latent SITS are then concatenated along their temporal dimension and processed by a transformer encoder in the same temporal dimension. Thus, the TAMRF encoder follows the spatial then temporal encoding paradigm also found in SITS-Former Yuan et al. (2022) or UBARN Dumeur et al. (2024).

**Spatial encoders.** The date-wise spatial encoders follow the SRResNet architecture introduced in Ledig et al. (2017) and notably used in ESRGAN Wang et al. (2018). It consists in chaining several Residual-In-Residual-Dense-Block (RRDB), followed by an up-sampling operation achieved by bicubic interpolation. It can be noted that
<figure><img src="image_4.png" /><figcaption>Figure 3: Overview of the proposed method. $S_{LR}$ (in green) and $S_{HR}$ (in blue) are fed into the TAMRF encoder, which outputs latent representation $S^{latent}$ (in red). The latent representation is fed into the TAMRF decoder, along with $T_{query}$, which outputs $S_{LR}^{*}$ (in green) and $S_{HR}^{*}$ (in red). The large downward arrow on the left indicates the processing flow, from top to bottom.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 4: Details of TAMRF encoder. The input S<sub>LR</sub> (in green) and S<sub>HR</sub> (in blue) are first passed through separate date-wise spatial encoders based on the Residual-In-Residual-Dense-Block (RRDB) architecture, including spatial up-sampling of S<sub>LR</sub>. The resulting series S<sub>LR</sub><sup>latent</sup> and S<sub>HR</sub><sup>latent</sup> have matching dimensions, except for the temporal dimension. The F latent features are then extended with temporal positional encoding (L features) and a learnable sensor token (M Features). Therefore, each pixel is converted to a sequence of tokens with a total of D = F + L + M features. This sequence is processed by a pixel-wise transformer encoder along the temporal dimension, forming the latent SITS S<sub> latent </sub>. The large downward arrow on the left indicates the processing flow, from top to bottom.</figcaption></figure>
compared to other up-sampling layers such as Pixel-Shuffle for instance, bicubic interpolation supports noninteger resolution ratios. The up-sampling is not applied in the HR encoder, since its input images already have the target resolution. With respect to the original architecture of SRResNet, the convolution layers after the up-sampling operation are removed. These final layers will form the spatial part of the TAMRF decoder, presented in section 2.2.2.

The LR spatial encoder outputs $F$ features on a $W \times H$ spatial support from the $\#C_{LR}$ input channels on a $w \times h$ spatial support. On the other hand, the HR spatial encoder outputs $F$ features on a $W \times H$ spatial support from the $\#C_{HR}$ input channels on a $W \times H$ spatial support. Each encoder is therefore specific to the set of spectral bands ($C_{LR}$ or $C_{HR}$), and to the spatial resolution it processes. In the outputs of both encoders, all dimensions match except for the temporal dimension.

**Temporal and modality encoder.** For a pixel corresponding to a spatial location, the $F$ features extracted by the date-wise spatial encoders for each date of $S_{HR}$ and each date of $S_{LR}$ are concatenated into a sequence of measurements. To allow the downstream temporal encoder to reason with sensors and acquisition dates, the $F$ features are extended with contextual temporal and sensor information. The temporal context is obtained through the traditional positional encoding used by transformers in Natural Language Processing Vaswani et al. (2017). This positional encoding has been extended to temporal positional encoding in many works on SITS, including Yuan et al. (2022); Dumeur et al. (2024); Guo et al. (2024). It alternates sine and cosine functions with varying frequencies:

$$ \Psi(t, l) = \begin{cases} \sin(t/\lambda^{l/L}) & \text{if } l \text{ is even,} \\ \cos(t/\lambda^{l/L}) & \text{if } l \text{ is odd,} \end{cases} \quad l \in 0, \dots, L-1, \qquad (10) $$

where $L$ is the number of temporal positional features and $\lambda$ is a normalization factor. The $L$ features of the temporal positional encoding are concatenated in the feature dimension with the $F$ features produced by the date-wise spatial encoder. The sensor context is obtained from two dedicated learnable sensor tokens of size $M$, one for $S_{HR}$ and one for $S_{LR}$. The corresponding token is also concatenated to the feature
dimension. A pixel (i.e. spatial location) is therefore described by a sequence of #T<sub>LR</sub> + #T<sub>HR</sub> tokens of D = F + L + M features. Since the D feature space includes temporal positional encoding, the order of dates in the latent SITS does not matter.

This pixel-wise sequence of tokens is then processed by a transformer encoder
module as described in Vaswani et al. (2017): the sequence passes through several
blocks of multi-head self-attention followed by a feed-forward network. The resulting
pixel-wise sequences are rearranged into spatial dimensions to form the S<sup>latent</sup> SITS.

2.2.2. Decoder

The TAMRF decoder operates on the latent SITS S<sup>latent</sup> produced by the encoder as well as on the acquisition time queries T<sub>query</sub>. As shown in Fig. 5, it consists of a pixel-wise temporal decoder followed by a date-wise spatial decoder.

***Temporal decoder.*** The temporal decoder mechanism is similar to the decoder used in the cross-reconstruction task of ALISE Dumeur et al. (2025). The first step is to transform the acquisition time queries *T<sub>query</sub>* into a sequence of #T<sub>query</sub> tokens of *D* features. A temporal positional encoding of size *L* is computed using eq. 10, and concatenated to a learnable placeholder of size *F + M*.

The query sequence and the pixel-wise sequences of $S_{\text{latent}}$ are then fed to a mod-
ified transformer decoder instance: the first self attention block is removed and the
$S_{\text{latent}}$ sequence serves as keys and values for the multi-head cross-attention block,
while the acquisition time query sequence serves as queries. The resulting output pixel-
wise sequence can be seen as $S_{\text{latent}}$ reconstructed at the acquisition times in $T_{\text{query}}$.

***Spatial decoder.*** Following this pixel-wise temporal decoding, the resulting sequences are rearranged into spatial dimensions and fed to the date-wise spatial decoder. This decoder is fairly simple and corresponds to the final two blocks of the SRResNet architecture Wang et al. (2018): two convolutional layers, with activation in-between. The spatial decoder outputs #C<sub>LR</sub> + #C<sub>HR</sub> channels, which can then be split into the two predicted SITS S<sub>LR</sub>* and S<sub>HR</sub>*.
<figure><img src="image_6.png" /><figcaption>Figure 5: Details of TAMRF decoder. The query acquisition times $T_{query}$ are transformed into a sequence of tokens by means of temporal positional encoding and a learnable placeholder. The query sequence and the latent SITS $S^{latent}$ are processed by a pixel-wise temporal decoder, starting with a cross-attention block, which reconstructs $S^{latent}$ at query acquisition times. The resulting SITS is then passed through a simple spatial decoder in order to predict the resulting $S_{LR}^{*}$ and $S_{HR}^{*}$ SITS. The large downward arrow on the left indicates the processing flow, from top to bottom.</figcaption></figure>
## 2.3. Training strategy

The problem formulation proposed in section 2.1 can be functionally implemented by the architecture proposed in section 2.2, but finding the optimal parameters $\Theta^*$ requires a training strategy that pushes toward the desired properties of $S_{LR}^*$ and $S_{HR}^*$.

First, predictions should closely match the actual observations that could have been made on that date. This can be achieved through an SSL strategy, and more specifically a Masked Auto-Encoder (MAE) strategy He et al. (2022); Reed et al. (2023); Liu et al. (2025), described in 2.3.1. The reconstruction loss term of the MAE needs to be adapted to the difference in spatial resolutions of the input SITS as described in section 2.3.2.

Second, despite the probable presence of clouds in both LR and HR input SITS, no cloud masks are required as input to the model. On the other hand predictions should correspond to ground measurement and be free of clouds or missing data. This is enforced by the use of a novel loss term inspired from contrastive learning, which is presented in section 2.3.3.

Third, all predicted bands should have the best spatial resolution among the input sensors. This is enforced by a dedicated loss term that aims at favoring HR details in all predictions, which is detailed in section 2.3.4.

### 2.3.1. Masking Strategy

As opposed to Liu et al. (2024) which adopts a purely random masking strategy, the masking strategy for TAMRF is randomly drawn from a pool of strategies described in Table 2. This allows to favor configurations that are unlikely to happen regularly with random masking, and which are of interest for downstream applications. They include completely missing the HR or LR SITS, or long term forecasting. For each input SITS, a different strategy is drawn, with random parameters if the strategy has parameters, in order to maximize the input variability.

The network is trained to reconstruct both masked and unmasked dates in $S_{HR}$ and $S_{LR}$, and thus during training we use as query dates $T_{query}^{training} = T_{HR} \cup T_{LR}$. Reconstructing unmasked dates during training is necessary to ensure that the network will behave correctly on those dates during inference. It should be reminded that inferring
unmasked dates is of interest with TAMRF, since the model will provide higher resolution versions of unmasked $S_{LR}$ dates and cloud-free predictions of all unmasked inputs. In the following sections, the SITS containing the masked dates will be called $S_{LR}^m$ and $S_{HR}^m$ respectively. The total loss used in order to optimize TAMRF is given by:

$$L = \underbrace{\alpha_m(w_{LR}^m L_{LR}^m + w_{HR}^m L_{HR}^m)}_{\text{masked reconstruction term}} + \underbrace{\alpha_c(w_{LR}^c L_{LR}^c + w_{HR}^c L_{HR}^c)}_{\text{clear reconstruction term}} + \underbrace{\alpha_m(w_{LR}^m \kappa_{LR}^m + w_{HR}^m \kappa_{HR}^m)}_{\text{masked contrastive term}} + \underbrace{\alpha_c(w_{LR}^c \kappa_{LR}^c + w_{HR}^c \kappa_{HR}^c)}_{\text{clear contrastive term}} + \underbrace{\alpha_s w_{HR}^c L_s}_{\text{spatial term}} \quad (11)$$

where superscript c (resp. m) denote clear dates (resp. masked dates), $L_{HR}^c$, $L_{HR}^m$, $L_{LR}^c$ and $L_{LR}^m$ are reconstruction terms detailed in section 2.3.2, $\kappa_{HR}^c$, $\kappa_{HR}^m$, $\kappa_{LR}^c$ and $\kappa_{LR}^m$ are invalid data contrastive terms detailed in section 2.3.3, $L_s$ is the spatial reconstruction loss term detailed in section 2.3.4, $\alpha_m$, $\alpha_c$ and $\alpha_s$ are constant weights designed to balance the different loss terms, and $w_{LR}^m$ (resp. $w_{HR}^m$, $w_{LR}^c$, $w_{HR}^c$) is the ratio of masked LR dates with respect to the total number of LR and HR dates in the SITS. The next sections will detail each of those terms.

<table><caption>Table 2: Pool of strategies for MAE training.</caption><thead><tr><td>Strategy</td><td>Description</td><td>Parameter</td></tr></thead><tbody><tr><td>Random</td><td>Randomly discard dates from S<sub>HR</sub> and S<sub>LR</sub></td><td>Discard rate</td></tr><tr><td>Gaps</td><td>Generate periods where no S<sub>HR</sub> or S<sub>LR</sub> are available</td><td>Length of gaps</td></tr><tr><td>No HR</td><td>Discard S<sub>HR</sub> entirely</td><td>N/A</td></tr><tr><td>No LR</td><td>Discard S<sub>LR</sub> entirely</td><td>N/A</td></tr><tr><td>Forecast</td><td>Discard all dates in S<sub>HR</sub> and S<sub>LR</sub> after a given date</td><td>Date</td></tr><tr><td>Backcast</td><td>Discard all dates in S<sub>HR</sub> and S<sub>LR</sub> before a given date</td><td>Date</td></tr></tbody></table>
### 2.3.2. Reconstruction term

The reconstruction loss terms make use of the Huber loss, introduced in Girshick (2015) and given by:

$$L_{1}^{\text{smooth}}(x) = \begin{cases} 0.5x^2, & \text{if } |x| \le \epsilon \\ \epsilon(|x| - 0.5\epsilon), & \text{otherwise.} \end{cases} \qquad (12)$$

The $L_1^{\text{smooth}}$ loss behaves like the $L_1$ loss when differences between predicted and target values are large, limiting the impact of outliers, and like the $L_2$ when differences are small, which is better for optimisation. In the experiments, it has proven to be beneficial over the standard MSE ($L_2$), especially in early stages of training. Since $L_1^{\text{smooth}}$ will operate on standardized values, we choose $\epsilon = 0.1$.

Two loss terms are defined for each of $S_{HR}$ and $S_{LR}$. Loss terms for $S_{HR}$ are given by:

$$L_{HR}^{c} = \frac{1}{\sum_{R_{HR} \times T_{HR}^{c} \times C_{HR}} V_{HR}^{c}} \sum_{R_{HR} \times T_{HR}^{c} \times C_{HR}} L_{1}^{\text{smooth}}(S_{HR}^{\star} - S_{HR}^{c})V_{HR}^{c}, \qquad (13)$$

$$L_{HR}^{m} = \frac{1}{\sum_{R_{HR} \times T_{HR}^{m} \times C_{HR}} V_{HR}^{m}} \sum_{R_{HR} \times T_{HR}^{m} \times C_{HR}} L_{1}^{\text{smooth}}(S_{HR}^{\star} - S_{HR}^{m})V_{HR}^{m}, \qquad (14)$$

where the sum is performed over all $(x, y)$, $t$, $c \in R_{HR} \times T_{HR}^{c} \times C_{HR}$ (resp. $R_{HR} \times T_{HR}^{m} \times C_{HR}$) and $V_{HR}^{c}$ (resp. $V_{HR}^{m}$) is a binary validity mask that is derived from the metadata of the products and accounts for clouds, shadows, out-of-swath areas, etc. Note that this mask is only used for the loss computation, and therefore only during training, not as input to the network. Using such masks is key to learning quality predictions by avoiding setting cloudy areas as targets for the $L_1^{\text{smooth}}$ loss.

Because $S_{LR}^{\star}$ is of spatial resolution $R$ while $S_{LR}^{c}$ and $S_{LR}^{m}$ are of coarser spatial resolution $r$, the formulation of the $S_{LR}$ loss terms is adapted as follows:

$$L_{LR}^{c} = \frac{1}{\sum_{R_{LR} \times T_{LR}^{c} \times C_{LR}} V_{LR}^{c}} \sum_{R_{LR} \times T_{LR}^{c} \times C_{LR}} L_{1}^{\text{smooth}}((S_{LR}^{\star} * \omega_{\sigma}) \downarrow_{k} - S_{LR}^{c})V_{LR}^{c}, \qquad (15)$$

$$L_{LR}^{m} = \frac{1}{\sum_{R_{LR} \times T_{LR}^{m} \times C_{LR}} V_{LR}^{m}} \sum_{R_{LR} \times T_{LR}^{m} \times C_{LR}} L_{1}^{\text{smooth}}((S_{LR}^{\star} * \omega_{\sigma}) \downarrow_{k} - S_{LR}^{m})V_{LR}^{m}, \qquad (16)$$
where $\omega_\sigma$ is a gaussian spatial kernel designed for smoothing the images in $S_{LR}^*$ and $\downarrow_k$ is the nearest-neighbour downsampling operator of factor $k$. Combining those operators allow bringing high resolution prediction $S_{LR}^*$ back to its original low resolution $r$.

### 2.3.3. Invalid data contrastive term

The use of validity masks $V_{HR}^m$, $V_{HR}^c$, $V_{LR}^m$ and $V_{LR}^c$ to mask the reconstruction loss terms should theoretically help the temporal transformer to learn that input invalid pixels such as cloud, cloud shadow or missing data are uninformative for the temporal reconstruction and should not be used. However, those masks from Level 2 products have commission and omission errors which might actually prevent to learn to ignore those pixels. For instance, if a substantial amount of cloudy pixels are not identified as such in the Level 2 cloud mask, those pixels will serve as active reference data for the reconstruction term of section 2.3.2, and drive the model toward reconstructing cloudy pixels. For this reason, ignoring mask pixels in the reconstruction term is not sufficient to prevent invalid pixels leakage into the predictions. To overcome this issue, we introduce an invalid data contrastive term enforcing that the predictions corresponding to masked reference pixels should be closer to the closest valid pixels in time than to the invalid reference pixel of the predicted date. This is formulated as a triplet margin loss Balntas et al. (2016):

$$ \kappa_{LR} = \max \{ L_2(S_{LR}^*, S_{LR}^{\text{valid}}) - L_2(S_{LR}^*, S_{LR}^{\text{invalid}}) + \text{margin}, 0 \}, \quad (17) $$

where $L_2$ is the MSE, $S_{LR}^{\text{invalid}}$ represents invalid pixels of $S_{LR}$ according to validity mask $V_{LR}$ derived from the image Level 2 product quality masks, $S_{LR}^{\text{valid}}$ represents the closest valid pixel in time according to validity mask $V_{LR}$ derived from the Level 2 product quality masks, and margin is a user-defined threshold. $\kappa_{HR}$ is defined likewise. The behavior of the invalid data contrastive loss is illustrated in Fig. 6. In practice, a threshold on the time distance to the closest valid pixel is used, so that valid pixels that are too far away in time are not used in the computation of the contrastive loss term.
<figure><img src="image_7.png" /><figcaption>Figure 6: Illustration of the invalid data contrastive term. For prediction $S^*$ At time $t$ (noted with $★$ marker), the reference data $S^{invalid}$ (noted with ◦ marker). The closest in time valid reference data $S^{valid}$ (noted with ◆ marker) occurs at time $t_{n-1}$. If $S^*$ is closest to $S^{valid}$ than to $S^{invalid}$ by margin, then the contrastive term is null (Fig. 6a). If $S^*$ is closest to $S^{invalid}$ than to $S^{valid}$ by margin, then the contrastive term is positive, penalising $S^*$ to draw it further away from the invalid data and closer to the valid data (Fig. 6b).</figcaption></figure>

### 2.3.4. Spatial term

The $L_{LR}^c$ and $L_{LR}^m$ reconstruction loss terms operate at initial LR resolution. As such, they can not drive $S_{LR}^*$ toward high spatial resolution details, even if $S_{LR}^*$ has the spatial sampling of $R_{HR}$. The usual solution to this problem is to assume scale invariance and apply Wald's protocol. In our work, this hypothesis is avoided by leveraging the correlation between $S_{HR}$ and $S_{LR}$ SITS in order to favor high spatial frequencies in $S_{LR}^*$.

Indeed, during training $T_{query} = T_{HR} \cup T_{LR}$, and both $C_{LR}$ and $C_{HR}$ bands are predicted for each time in $T_{query}$ (see Eq. 7). Therefore, for each acquisition in $T_{HR}$ (either masked or clear in the MAE strategy), predicted bands in $S_{LR}^*$ can be compared to reference HR images in $S_{HR}$. However, no assumptions have been made regarding corresponding bands in $C_{LR}$ and $C_{HR}$. Therefore, $S_{HR}$ can not be directly used to supervise the spatial reconstruction of $S_{LR}^*$.

To overcome this issue, we project the $C_{HR}$ channels from $S_{HR}$ onto the $C_{LR}$ channels of $S_{LR}^*$, by means of a linear regression. Ridge regularization is employed to make the linear regression more robust to degenerated cases. This linear projection is performed separately for each date in $T_{HR}$, to compensate for the relatively weak estimate provided by linear regression. This is especially true for bands in $C_{LR}$ that do not have a close match in $C_{HR}$ or depend on external factors that are not measured by bands
in C<sub>HR</sub>. The latter case is encountered with the thermal band of Landsat-8/9, whose primary driver is meteorological. This linear regression is denoted Λ(S<sub>LR</sub><sup>*</sup>, S<sub>HR</sub>).

Despite the date-wise estimation of the linear regression, Λ(S<sub>LR</sub><sup>*</sup>, S<sub>HR</sub>) is still a weak estimate of S<sub>LR</sub>. In particular, some local trends might not be correctly accounted for, and using Λ(S<sub>LR</sub><sup>*</sup>, S<sub>HR</sub>) for the direct supervision of S<sub>LR</sub><sup>*</sup> might lead to radiometric distortion of the latter. Since radiometric accuracy is enforced by the reconstruction loss term of Eq. 15, only the high spatial frequency content of Λ(S<sub>LR</sub><sup>*</sup>, S<sub>HR</sub>) should supervise the high spatial frequency content of S<sub>LR</sub><sup>*</sup>. This is achieved through the use of the Learned Perceptual Image Patch Similarity (LPIPS) Zhang et al. (2018). Indeed, in a previous work Michel et al. (2025), we demonstrated that LPIPS is especially well suited for the measurement of difference in spatial frequency content. Furthermore, it has interesting properties such as being relatively blind to radiometric distortion and spatial misalignment. The complete spatial reconstruction loss term is given by:

$$L_{spat} = \text{LPIPS}(S_{LR}^{*}, \Lambda(S_{LR}^{*}, S_{HR})) \quad (18)$$

# 3. Experiments

## 3.1. Dataset

All the experiments have been conducted using a new Landsat-8/9 and Sentinel-2 dataset called Landsat-8/9 to Sentinel-2 (LS2S2), which is publicly available<fnref n="1" />, <fnref n="2" />. This dataset comprises two subdatasets of joint Sentinel-2 and Landsat-8/9 SITS. The first, called eu22 is composed of 64 Areas Of Interest (AOIs) for the training set and 41 AOIs for the testing set, each covering an area of 9.9x9.9 km<fnref n="2" /> (990x990 pixels for Sentinel-2 and 330x330 pixels for Landsat-8/9). It should be noted that testing and training AOIs do not overlap. AOIs of the eu22 split are mostly located in Europe, with a few additional sites in Equatorial Africa and French Guyana. All SITS span the full 2022 year. This yields a total of 2 984 Sentinel-2 images and 1 609 Landsat-8/9 images in the training set, and 1 581 Sentinel-2 images and 736 Landsat-8/9 images in the testing

<footnote n="1">¹https://zenodo.org/records/15471890</footnote>

<footnote n="2">²https://doi.org/10.6096/1029</footnote>
set. The second split, called **wwmy**, is composed of 138 AOIs for the training set and 69 AOIs for the testing set, with the same area size. AOIs were chosen over the entire world within the most populated joint classes between CCI Land-Cover CCI (2017) and the Köppen classification of climates Rubel and Kottek (2010), in order to get a fair sampling of bio-climatic regions and land-cover. They cover 69 unique combinations of Köppen and CCI classes (composed of 20 different CCI classes and 18 different Köppen classes). They still comprise one year of data, but the starting date has been randomly chosen for each AOI between 2022.01.01 (beginning of Landsat-9 availability) and 2024.06.01. This yields a training set of 6 477 Sentinel-2 images and 4 403 Landsat-8/9 images, and a testing set of 2 981 Sentinel-2 images and 2 042 Landsat-8/9 images. Unless stated otherwise, experiments are conducted on the union of the two splits, noted **all**, which comprises 202 training AOIs with 9 461 Sentinel-2 images and 4 590 Landsat-8/9 images, and 110 testing AOIs with 5 984 Sentinel-2 images and 2 778 Landsat-8/9 images. Fig. 7 shows a map of the different AOIs for the training and testing sets of the **eu22** and **wwmy** splits, while figure 8 shows the variability of the temporal sampling of the **wwmy** testing set. All dates with more than 25% of clear pixels over the AOI and the acquisition time range are included in the dataset. Additional statistics about the number of dates per AOI for each sensor are presented in Table 3. Top-of-Canopy surface reflectances from level 2 products are used. The spectral bands included in the dataset are presented in Table 4. It can be observed that the Landsat-8/9 sensor does not have Red Edge or wide Infra-Red bands, and conversely Sentinel-2 sensor does not retrieve Land Surface Temperature (LST). In addition to the spectral bands, corresponding quality masks have been used to derive a validity mask for each date of each sensor. This dataset has been gathered through the OpenEO API Schramm et al. (2021). Fig. 9 gives an example of SITS extracted from the LS2S2 **eu22** dataset. It should be stressed that the LS2S2 **wwmy** split is significantly more challenging than the LS2S2 **eu22** split, given the wider variety of landscapes and cloud occurrences. In particular, due to high latitude sites, it contains significantly more snow: around 18% of the LS2S2 **eu22** test set images contain snow, compared to nearly 26% of the LS2S2 **wwmy** test set. A similar ratio holds for the training sets. Finally, it should be noted that despite the variable spatial resolution across bands for both Landsat-8/9 and
Sentinel-2, for each sensor all bands are up-sampled to the maximum resolution with
the bicubic interpolator (10 meters for Sentinel-2 and 30 meters for Landsat-8/9).

<figure><img src="image_8.png" /><figcaption>Figure 7: Geographical coverage of the LS2S2 dataset, divided into the eu22 split spanning year 2022 and the wwmy split spanning random periods of one year between 2022.01.01 and 2025.06.01.</figcaption></figure>

3.2. *Experimental setup*

3.2.1. *Model and training hyperparameters*

The hyperparameters for the TAMRF model have been set experimentally. Table 6
shows the main hyperparameters for each component, as well as their size in number
of parameters to optimise. With a total of 2.3M parameters, TAMRF is still a relatively
small model compared to the pretrained VGG instance used in LPIPS for instance,
which has 16M parameters.

The LS2S2 **all** training split is subdivided into 6 673 training patches and 239 validation patches, each of 165x165 pixels for Sentinel-2 SITS and 55x55 pixel for Landsat-8/9 SITS. Each patch contains the full time-series. TAMRF is trained with the Adam optimizer Kingma and Ba (2014), with one pair of Sentinel-2 and Landsat-8/9 SITS passed at each step, which can be seen as using a batch size of one. Because
<figure><img src="image_9.png" /><figcaption>Figure 8: Temporal variability of the testing set of the **wwmy** split.</figcaption></figure>

<table><caption>Table 3: Number of 9.9x9.9 km² images in LS2S2 dataset for each sensor and each split. Average, minimum and maximum columns show statistics on the number of acquisition dates per AOI, for each sensor.</caption><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th colspan="4">Sentinel-2</th><th colspan="4">Landsat-8/9</th></tr><tr><th>Total</th><th>Average</th><th>Min</th><th>Max</th><th>Total</th><th>Average</th><th>Min</th><th>Max</th></tr></thead><tbody><tr><td rowspan="2">eu22</td><td>train</td><td>2 984</td><td>44</td><td>11</td><td>94</td><td>1 581</td><td>26</td><td>0</td><td>61</td></tr><tr><td>test</td><td>1 609</td><td>39</td><td>6</td><td>98</td><td>736</td><td>18</td><td>1</td><td>56</td></tr><tr><td rowspan="2">wwmy</td><td>train</td><td>6 477</td><td>47</td><td>7</td><td>112</td><td>4 403</td><td>32</td><td>0</td><td>88</td></tr><tr><td>test</td><td>2 981</td><td>43</td><td>6</td><td>124</td><td>2 042</td><td>32</td><td>0</td><td>80</td></tr></tbody></table>
<table><caption>Table 4: Spectral bands included in the LS2S2 dataset for each sensor. The table layout allows highlighting which bands are similar between sensors and which bands are only available on one sensor. Despite the variable spatial resolution across bands for both Landsat-8/9 and Sentinel-2, for each sensor all bands are up-sampled to the maximum resolution with the bicubic interpolator (10 meters for Sentinel-2 and 30 meters for Landsat).</caption><thead><tr><th colspan="2">Sentinel-2</th><th colspan="2">Landsat-8/9</th><th rowspan="2">Description</th></tr><tr><th>Band</th><th>Resolution</th><th>Band</th><th>Resolution</th></tr></thead><tbody><tr><td>B02</td><td>10.</td><td>B1</td><td>30.</td><td>Deep blue</td></tr><tr><td>B03</td><td>10.</td><td>B2</td><td>30.</td><td>Blue</td></tr><tr><td>B04</td><td>10.</td><td>B3</td><td>30.</td><td>Green</td></tr><tr><td>B05, B06, B07</td><td>20.</td><td>B4</td><td>30.</td><td>Red</td></tr><tr><td>B08</td><td>20.</td><td>B5</td><td>30.</td><td>Red Edge</td></tr><tr><td>B8a</td><td>10.</td><td>B6, B7</td><td></td><td>Near Infra-Red</td></tr><tr><td>B11, B12</td><td>20.</td><td>B10</td><td>100.</td><td>Wide Near Infra-Red<br>Short Wavelength Infra-Red<br>Land Surface Temperature</td></tr></tbody></table>
<figure><img src="image_10.png" /><figcaption>Figure 9: Example of SITS from LS2S2 dataset eu22 training split, for UTM tile 31UFR. Color composition is RGB (B04, B03, B02 for Sentinel-2 and B4, B3, B2 for Landsat). Invalid pixels according to quality masks are highlighted in red.</figcaption></figure>
of memory limitations, if the total number of images for a given step is higher than 50, exceeding images of both sensors are randomly dropped. The masking strategy parameters are summarized in Table 5.

A single TAMRF model is trained with the self-supervised training strategy presented in section 2.3, as TAMRF requires no fine-tuning or re-training when transposed to different AOIs.

TAMRF is trained for 150 epochs, for a total of $10^6$ steps, with validation steps every 2 018 steps. The best model parameters are selected according to the validation loss. The initial learning rate is set to $10^{-4}$. After a linear warm-up of 2 018 steps, the cosine annealing with warm restarts Loshchilov and Hutter (2016) is used to modulate the learning rate between the initial learning rate and a minimum learning rate of $2 \times 10^{-5}$, with an initial period of 40 360 steps and a multiplicative factor of 2. Gradient clipping is set to 0.1, since it proved to speed up the convergence in the experiment. Weights for the different terms of the loss in Eq. 11 were set empirically to $\alpha_c=1$, $\alpha_m=0.5$ and $\alpha_s = 1$, except for the LST (B10) band of Landsat-8/9, for which $\alpha_s = 0.1$. This setting allows balancing the relative importance of the terms. Intuitively, the reconstruction of masked dates is more difficult than the reconstruction of clear dates, hence, masked reconstruction loss terms tend to have higher values. The same applies to the reconstruction of high frequencies of Landsat-8/9 LST band, which is more challenging than the others due to the lower correlation with the HR bands. All codes use PyTorch Paszke et al. (2019). The total training time is around 5 days on a NVIDIA H100 GPU.

### 3.2.2. Metrics

This section describes the metrics used to assess the performances of the trained TAMRF model, and compare it to other models from the literature. In a previous work on cross-sensor SISR Michel et al. (2025), we identified a set of recommended metrics that are able to correctly rank sharpness in a cross-sensor setup, and focus on different key aspects of the SISR performances. Here we retain three of those qualified metrics, each with a distinctive purpose:

* RMSE measures the faithfulness to the radiometry of the target image at ini-
<table><caption>Table 5: Parameters of the MAE strategy described in section 2.3.1. First, a strategy is drawn according to the probability column. Depending on the selected strategy, its parameters are then drawn within the parameter range. U Stands for the uniform distribution. Gaps length and tipping dates are expressed in number of days since the first day of the SITS.</caption><thead><tr><th>Strategy</th><th>Probability</th><th>Parameter range</th></tr></thead><tbody><tr><td>Random</td><td>0.5</td><td>Discard rate ~ U(0.2,0.7)</td></tr><tr><td>Gaps</td><td>0.25</td><td>Length of gaps ~ U(30,90)</td></tr><tr><td>No HR</td><td>0.0625</td><td>N/A</td></tr><tr><td>No LR</td><td>0.0625</td><td>N/A</td></tr><tr><td>Forecast</td><td>0.0625</td><td>Tipping date ~ U(65,300)</td></tr><tr><td>Backcast</td><td>0.0625</td><td>Tipping date ~ U(65,300)</td></tr></tbody></table>

tial resolution,

* BRISQUE is a No-Reference Image Quality (IQ) metric that grades the global IQ,

* Frequency Restoration (FR) compares the spatial frequency content of the predicted and target images at initial resolution.

All metrics are measured separately for each band and for each date, either clear (used as model input) or masked (not used as model input). This forms a self-contained protocol that allows to evaluate several key aspects of the spatial and temporal resolution improvement without relying on additional higher resolution or in-situ reference data.

RMSE. RMSE is computed with respect to the initial resolution of the target data. Since $S_{HR}^*$ has the same spatial resolution R as $S_{HR}$, RMSE$_{HR}$ has a simple formulation:

$$RMSE_{HR} = L_2(S_{HR}^* - S_{HR}). \qquad (19)$$

On the other hand, $S_{LR}^*$ has a spatial resolution of R while $S_{LR}$ has a spatial resolution of r = kR. To be able to compute RMSE$_{LR}$, $S_{LR}^*$ is first smoothed by a Gaussian
<table><caption>Table 6: Summary of model hyperparameters for each component of the TAMRF model. The last column shows the size in number of parameters to optimise for each component.</caption><thead><tr><th>Component</th><th>Key parameters</th><th>#Parameters</th></tr></thead><tbody><tr><td>LR Spatial Encoder</td><td>in channels (C<sub>LR</sub>)<br>out channels (F)<br>nb. blocks<br>up. factor (k)</td><td>8<br>64<br>1<br>3<br>797 K</td></tr><tr><td>HR Spatial Encoder</td><td>in channels (C<sub>HR</sub>)<br>out channels (F)<br>nb. blocks<br>up. factor</td><td>10<br>64<br>1<br>1<br>762 K</td></tr><tr><td>Temporal Encoder</td><td>positional encoding (L)<br>sensor token (M)<br>token size (D = F + M + L)<br>feed-forward<br>nb. layers<br>nb. heads</td><td>64<br>8<br>136<br>256<br>3<br>4<br>435 K</td></tr><tr><td>Temporal Decoder</td><td>token size (D = F + M + L)<br>nb. layers<br>nb. heads</td><td>136<br>1<br>4<br>145 K</td></tr><tr><td>Spatial Decoder</td><td>in channels (D = F + M + L)<br>out channels (C<sub>LR</sub> + C<sub>HR</sub>)</td><td>136<br>18<br>188 K</td></tr><tr><td>Total</td><td></td><td>2.3 M</td></tr></tbody></table>
kernel and down-sampled to a spatial resolution of *r*, similarly to loss term presented in Eq. 13:

$$ RMS E_{LR} = L_2((S_{LR}^* * \omega_\sigma) \downarrow_k -S_{LR}). \qquad (20) $$

It should be noted that both RMS E<sub>HR</sub> and RMS E<sub>LR</sub> are computed according to the validity reference masks already used in losses computation during training (see section 2.3.2).

Since RMS E<sub>LR</sub> is computed at lower resolution *r*, it does not measure the spatial resolution improvement of $S_{LR}^*$. The remaining two metrics are designed for this purpose.

**BRISQUE Score.** BRISQUE Mittal et al. (2012) is a No Reference IQ which builds local features mapped to IQ score by a Support Vector Regressor. It is trained on human annotated scores, yielding a score between 0 and 100, 0 being the best IQ. According to Michel et al. (2025), BRISQUE provides a solid criterion to measure the sharpness of the image as well as its overall quality, even if it tends to slightly favor noise. Using a No Reference IQ metric allows providing an insight on the performances that is not biased by the quality of the reference data, such as unmasked clouds, which are present in the LS2S2 dataset.

**Frequency Restoration.** Frequency Domain Analysis is proposed in Michel et al. (2025) and has notably been used in Ait-Bachir et al. (2025) to measure the improvement of the spatial frequency content of super-resolved images. It consists in analysing the Fourier domain Frequency Attenuation Profile ($F_{AP}$) for bandwidth [$f_m, f_M$]. Let

$$ U_{f_m, f_M} = \{(u, v) : f_m \le \sqrt{u^2 + v^2} < f_M\} \qquad (21) $$

denote the set of discrete spatial frequencies $(u, v)$ that lies within a ring defined by $f_m$ and $f_M$ in Fourier plane, $F_{AP}$ is given by:

$$ F_{AP}[P](f_m, f_M) = \frac{1}{\#U_{f_m, f_M}} \sum_{(u,v) \in U_{f_m, f_M}} |F[P](u, v)|, \qquad (22) $$
where $\#U_{f_m, f_M}$ is the number of elements in $U_{f_m, f_M}$, and $\mathcal{F}[P](u, v)$ is the Discrete Fourier Transform (DFT) of image $P$. $\mathcal{F}_{\mathcal{AP}}[P](f_m, f_M)$ is successively computed over a set of non-overlapping bandwidth intervals as given by the DFT quantization. The resulting set of values is denoted $\mathcal{F}_{\mathcal{AP}}[P](f_n)$, $n \in [0, N]$ with $f_n$ the central frequency of each frequency intervals. $\mathcal{F}_{\mathcal{AP}}[P](f_n)$ is averaged across batches and dataset. Finally, the normalized logarithmic $\mathcal{F}_{\mathcal{AP}}$ is computed, which gives spatial frequency attenuation in decibels:

$$ \mathcal{F}_{\mathcal{AP}}^{\log}[P](f_n) = 10 \cdot (\log_{10}(\mathcal{F}_{\mathcal{AP}}[P](f_n)) - \log_{10}(\mathcal{F}_{\mathcal{AP}}[P](f_0))). \quad (23) $$

Since there are no HR reference images for $S_{LR}$ and conversely no LR input images for $S_{HR}$, we adopt the No-Reference Actual Frequency Restoration metric (Michel et al. (2025), eq. 24), which will be abbreviated Frequency Restoration (FR) in the remaining of the paper:

$$ FR_{LR} = \sum \mathcal{F}_{\mathcal{AP}}^{\log}(S_{LR}^{\star}) - \mathcal{F}_{\mathcal{AP}}^{\log}((S_{LR})^{\uparrow_k}), \quad (24) $$

$$ FR_{HR} = \sum \mathcal{F}_{\mathcal{AP}}^{\log}(S_{HR}^{\star}) - \mathcal{F}_{\mathcal{AP}}^{\log}(S_{HR}), \quad (25) $$

where $(S_{LR})^{\uparrow_k}$ is the bicubic up-sampling of factor $k$ of $S_{LR}$. FR is also expressed in decibels.

### 3.2.3. Dataset particularities

There are three particularities of the LS2S2 dataset which need to be addressed in the model. First, all Landsat-8/9 bands are sampled at 30 meter resolution, and all Sentinel-2 bands are sampled at 10-meter resolution. Therefore, the up-sampling factor $k = 3$. This requires a small adaptation of the up-sampling operation in SRResNet, which usually performs progressive dyadic up-sampling. Second, despite this constant spatial sampling resolution for each sensor, some Sentinel-2 spectral bands have a native resolution of 20 meters, and the LST band of Landsat-8/9 has a native resolution of 100 meters, as stated in Table 4. This is accounted for in the reconstruction loss term (section 2.3.2) and in the computation of the RMSE metric (section 3.2.2) by using
larger Gaussian kernels for blurring those bands, independently of the down-sampling
factor. Finally, the main driver for the Land Surface Temperature band is meteorologi-
cal. Colder conditions or stronger wind will yield lower LST, while warmer conditions
and weaker wind will yield higher LST. Since the TAMRF model is not informed by
meteorological conditions, trying to predict LST at unobserved Landsat-8/9 date is
meaningless. For this reason, the reconstruction loss term (section 2.3.2) and RMSE
metric (section 3.2.2) are only computed for clear Landsat-8/9 dates for the LST band.
Finally, the first date of the SITS is used as the reference date to compute the temporal
position.

3.3. Results

To our best knowledge, there are no other methods in the literature that provide
the same capabilities as TAMRF (see Table 1). In order to provide an overview of
the performances of the proposed method with respect to existing ones, several tasks
have been identified and assessed on the LS2S2 all testing set. They include recon-
struction of Landsat-8/9 and Sentinel-2 dates in gaps of one month in section 3.3.1,
band-sharpening of Sentinel-2 20-meter bands in section 3.3.2, Spatio-Temporal Fu-
sion in section 3.3.3 and Thermal Sharpening in section 3.3.4. It is important to note
that for all those tasks, a single pre-trained instance of the TAMRF model has been
used, without any retraining or fine-tuning for the specific task.

For each task, baseline methods from the literature have been selected for compari-
son. Competing methods have been selected from existing works that address Sentinel-
2 and Landsat-8/9 fusion or spatial resolution enhancement. Only methods for which
the implementation was straightforward or for which the authors provided source code
and pre-trained weights have been considered. For this reason, recent methods such
as STF-Trans Huang et al. (2024), STF-Diff Huang et al. (2024) or RealFusion Guo
et al. (2025) could not be included in the benchmark. Even if those methods advertise
State-of-the-Art performances, they do not break free from the STF formulation and
its limitations. On the other hand, TAMRF proposes a unified and versatile framework,
whose interest is beyond the sole ranking on the STF task.

In the following results, all the images predicted by TAMRF have a spatial reso-
lution of 10 meters, for all bands of Sentinel-2 and Landsat-8/9. No filtering of input
images according to cloud or snow masks is performed, and the complete **all** testing set
is used, including both the **eu22** and the **wwmy** splits. The RMSE is computed follow-
ing Eq. 19 for Sentinel-2 10-meter bands and Eq. 20 for all the remaining Sentinel-2
and Landsat-8/9 bands. In all cases, only valid pixels according to the Level 2 products
mask are used for computing RMSE, thus excluding clouds and other invalid reference
pixels.

As explained in section 3.1, there is almost a quarter of the **wwmy** testing set dates that contains snow. Because the dynamic of snow reflectance is very large with respect to other land cover types, and because snow is a difficult land cover type to predict for unobserved dates, using the snowy dates as target dates in the evaluation of methods results in higher RMSE and lower IQ for all methods and all cases. Moreover, in most of the works published in the literature, snow is absent from the evaluation datasets. In order to provide numbers that can be compared to other works, we therefore choose to exclude the snowy target dates from the results presented in this section. A dedicated analysis on how well TAMRF predicts snowy dates is presented in section 4.5. The reader should keep in mind that this only concerns reference images in the evaluation process, as the snowy dates are used during training and as input to the model during evaluation.

3.3.1. *Gap-Filling*

The first task consists in generating one-month gaps in both Landsat-8/9 and Sentinel-2 SITS for all AOIs in the testing set, following the scheme of Fig. 10. All acquisitions in generated gaps are removed from the model input and kept aside for validation. The model is asked to reconstruct both the missing dates within gaps and the clear input dates, with all bands at 10-meter spatial resolution. TAMRF reconstructions are compared to a naive algorithm corresponding to what is usually used as pre-processing for downstream applications requiring temporally and spatially aligned data cubes. This naive method consists in linearly interpolating SITS from each sensor along the temporal dimension, taking into account the cloud and missing data mask. Then, all low resolution bands are spatially interpolated to 10-meter resolution with bicubic up-sampling.
TAMRF reconstructions are also compared to the reconstructions of U-TILISE Stucker et al. (2023), a gap-filling DL approach that combines hierarchical convolutional spatial encoders and a temporal transformer. For this experiment, we use the weights published by the authors, corresponding to the training of U-TILISE on the EarthNet dataset Requena-Mesa et al. (2021). It should be noted that this instance of U-TILISE is only able to process Sentinel-2 bands B2, B3, B4 and B8, thus only the results for these bands are presented for U-TILISE.

<figure><img src="image_11.png" /><figcaption>Figure 10: Gaps generated for the gap-filling task evaluation. When masked is True, the acquisition is removed from the model input SITS and used for performance evaluation. Note that only AOIs in the eu22 testing set are presented in the figure.</figcaption></figure>

Table 7 shows the mean and standard-deviation of the three evaluation metrics, for a subset of Sentinel-2 and Landsat-8/9 bands (full results are available in appendix 6), separately for masked dates, which are not in the model inputs, and clear dates, which are in the model inputs. Since the RMSE is computed at initial resolution, it is expected that the naive method gives RMSE values that are close to zero for clear dates. When compared to the naive method, TAMRF consistently provides better RMSE on masked dates, with a narrower standard deviation, while being on par with the naive model on clear dates. It also consistently outperforms U-TILISE on those masked dates, while being on par with U-TILISE predictions on clear dates. Moreover, TAMRF provides
better IQ for almost all bands, except for Sentinel-2 B8 and B4 on masked dates, where all the methods perform very similarly. On Landsat-8/9 dates, it lowers the BRISQUE score by a vast amount. Finally, spectral bands with an initial resolution coarser than 10 meters have positive FR for TAMRF, while the naive method provides FR close to zero. This highlights that our model effectively improves the high spatial frequency content for all bands, and explains why it yields RMSE on clear dates that are slightly worse than the naive method: those dates have more high resolution details. Nevertheless, RMSE values are close to the accuracy of Level 2 surface reflectances produced by the MACCS-ATCOR Joint Algorithm (MAJA) Colin et al. (2023) processor.

Fig. 11 shows the annual Normalized Difference Vegetation Index (NDVI) profile derived from the Gap-Filling experiment for an agricultural pixel of AOI 31TCJ. Images of predictions for unobserved dates at start, peak and end of season are provided for each method, along with the reference unobserved Sentinel-2 date. Compared to U-TILISE, TAMRF yields similar IQ on Sentinel-2 visible bands, while being slightly more accurate in terms of radiometry, which is consistent with the trends in Table 7. As expected, profiles from the Naive method fail to correctly capture the peak of vegetation growth, while both U-TILISE and TAMRF are able to retrieve it accurately, despite being in a gap period. TAMRF is slightly closer to the true unobserved maximum. All methods miss the regrowth of the end of season, which is due to the lack of observations beyond 2022. U-TILISE and TAMRF also predict a more gentle drop in NDVI than what the unobserved dates suggest between 2022.08 and 2022.09. Looking at the 10-meter Sentinel-2 predictions for the three unobserved example dates, it can be noted that U-TILISE and TAMRF have similar performances in predicting the land-cover changes. They fail however to retrieve the correct reflectance for bare soils on 2022.06.29 which appear darker in the reference image, probably due to a recent rain event that increased soil moisture. This kind of changes can not be predicted accurately for any of the methods because they lack meteorological conditioning (see section 5.1.2). Finally, it should be stressed that TAMRF also provides 10-meter resolution predictions of the Landsat-8/9 SITS, and that the NDVI profile generated from the Landsat-8/9 is perfectly consistent with the predicted Sentinel-2 profile.

Fig. 12 shows another illustration of predictions of both TAMRF for a target date
<table><caption>Table 7: Comparison between TAMRF model, U-TILISE and Naive interpolation on the gap-filling task, where regular gaps of 30 days are masked from the input SITS and kept appart for validation. Snow target dates are excluded from metrics computation. Only a selection of spectral bands is presented. Full results are available in appendix 6. * if usually higher FR means higher spatial resolution details, in the case of bands that are 10-meter native resolution, the best FR is the closest to zero in absolute value. ↓ (resp. ↑) indicates that the metric should be minimized (resp. maximized). Best mean values for each metric and each band are highlighted in **bold**. Clear designate dates for which the bands were observed by the model, and **masked** designates the dates that were removed from the input SITS as described in Fig. 10.</caption><thead><tr><th rowspan="2">Band</th><th rowspan="2">Method</th><th colspan="4">RMSE↓</th><th colspan="4">BRISQUE↓</th><th colspan="4">FR↑</th></tr><tr><th>Clear mean</th><th>std</th><th>Masked mean</th><th>std</th><th>Clear mean</th><th>std</th><th>Masked mean</th><th>std</th><th>Clear mean</th><th>std</th><th>Masked mean</th><th>std</th></tr></thead><tbody><tr><td rowspan="2">LS B1</td><td>Naive</td><td><strong>0.003</strong></td><td>0.006</td><td>0.079</td><td>0.137</td><td>87.7</td><td>8.7</td><td>87.7</td><td>8.8</td><td>0.6</td><td>2.0</td><td>-0.2</td><td>2.8</td></tr><tr><td>TARMF</td><td>0.008</td><td>0.034</td><td><strong>0.015</strong></td><td>0.048</td><td><strong>44.0</strong></td><td>18.7</td><td><strong>49.7</strong></td><td>21.1</td><td><strong>7.3</strong></td><td>2.5</td><td><strong>7.7</strong></td><td>2.8</td></tr><tr><td rowspan="2">LS B4</td><td>Naive</td><td><strong>0.004</strong></td><td>0.005</td><td>0.079</td><td>0.120</td><td>72.8</td><td>6.5</td><td>73.3</td><td>6.6</td><td>0.7</td><td>1.7</td><td>0.3</td><td>2.7</td></tr><tr><td>TARMF</td><td>0.009</td><td>0.030</td><td><strong>0.022</strong></td><td>0.047</td><td><strong>20.8</strong></td><td>11.3</td><td><strong>21.2</strong></td><td>12.7</td><td><strong>8.0</strong></td><td>2.4</td><td><strong>8.2</strong></td><td>2.7</td></tr><tr><td rowspan="2">LS B5</td><td>Naive</td><td><strong>0.006</strong></td><td>0.004</td><td>0.074</td><td>0.086</td><td>66.9</td><td>7.6</td><td>66.6</td><td>6.7</td><td>0.7</td><td>1.8</td><td>0.3</td><td>2.3</td></tr><tr><td>TARMF</td><td>0.015</td><td>0.025</td><td><strong>0.037</strong></td><td>0.038</td><td><strong>15.2</strong></td><td>7.0</td><td><strong>11.5</strong></td><td>6.1</td><td><strong>7.7</strong></td><td>2.4</td><td><strong>8.1</strong></td><td>2.5</td></tr><tr><td rowspan="2">LS B7</td><td>Naive</td><td><strong>0.005</strong></td><td>0.002</td><td>0.041</td><td>0.042</td><td>68.7</td><td>5.4</td><td>69.1</td><td>5.0</td><td>0.6</td><td>1.6</td><td>0.1</td><td>2.3</td></tr><tr><td>TARMF</td><td>0.009</td><td>0.018</td><td><strong>0.025</strong></td><td>0.027</td><td><strong>18.9</strong></td><td>10.6</td><td><strong>9.8</strong></td><td>4.9</td><td><strong>8.0</strong></td><td>2.4</td><td><strong>8.5</strong></td><td>2.7</td></tr><tr><td rowspan="3">S2 B4</td><td>Naive</td><td><strong>0.000</strong></td><td>0.000</td><td>0.023</td><td>0.033</td><td>19.3</td><td>10.6</td><td><strong>20.4</strong></td><td>10.3</td><td><strong>1.2</strong></td><td>2.6</td><td><strong>0.8</strong></td><td>2.9</td></tr><tr><td>U-TILISE</td><td>0.007</td><td>0.009</td><td>0.021</td><td>0.010</td><td>21.5</td><td>7.8</td><td>25.7</td><td>7.2</td><td>0.2</td><td>1.3</td><td>0.8</td><td>3.2</td></tr><tr><td>TARMF</td><td>0.008</td><td>0.006</td><td><strong>0.019</strong></td><td>0.023</td><td><strong>19.2</strong></td><td>12.7</td><td>24.1</td><td>14.0</td><td>1.1</td><td>2.6</td><td>0.6</td><td>3.0</td></tr><tr><td rowspan="2">S2 B6</td><td>Naive</td><td><strong>0.005</strong></td><td>0.002</td><td>0.031</td><td>0.030</td><td>47.4</td><td>6.3</td><td>48.0</td><td>6.1</td><td>1.2</td><td>2.1</td><td>1.0</td><td>2.3</td></tr><tr><td>TARMF</td><td>0.011</td><td>0.007</td><td><strong>0.027</strong></td><td>0.023</td><td><strong>14.2</strong></td><td>9.8</td><td><strong>16.2</strong></td><td>10.5</td><td><strong>5.5</strong></td><td>1.8</td><td><strong>5.1</strong></td><td>2.3</td></tr><tr><td rowspan="3">S2 B8</td><td>Naive</td><td><strong>0.000</strong></td><td>0.000</td><td>0.036</td><td>0.030</td><td><strong>13.8</strong></td><td>8.2</td><td><strong>14.7</strong></td><td>7.8</td><td>0.7</td><td>1.9</td><td>0.4</td><td>2.3</td></tr><tr><td>U-TILISE</td><td>0.013</td><td>0.017</td><td>0.044</td><td>0.020</td><td><strong>12.4</strong></td><td>6.8</td><td>16.5</td><td>6.4</td><td><strong>0.1*</strong></td><td>1.2</td><td><strong>0.2*</strong></td><td>2.4</td></tr><tr><td>TARMF</td><td>0.012</td><td>0.009</td><td><strong>0.032</strong></td><td>0.024</td><td>14.4</td><td>10.2</td><td>17.7</td><td>10.6</td><td>0.6</td><td>1.9</td><td>0.3</td><td>2.3</td></tr><tr><td rowspan="2">S2 B8a</td><td>Naive</td><td><strong>0.006</strong></td><td>0.003</td><td>0.035</td><td>0.028</td><td>47.3</td><td>6.4</td><td>47.9</td><td>6.3</td><td>1.1</td><td>1.9</td><td>0.9</td><td>2.2</td></tr><tr><td>TARMF</td><td>0.012</td><td>0.009</td><td><strong>0.031</strong></td><td>0.024</td><td><strong>14.7</strong></td><td>9.7</td><td><strong>16.0</strong></td><td>10.4</td><td><strong>5.4</strong></td><td>1.6</td><td><strong>5.1</strong></td><td>2.1</td></tr><tr><td rowspan="2">S2 B12</td><td>Naive</td><td><strong>0.004</strong></td><td>0.002</td><td>0.026</td><td>0.019</td><td>53.0</td><td>7.2</td><td>53.6</td><td>6.8</td><td>1.3</td><td>2.1</td><td>1.1</td><td>2.3</td></tr><tr><td>TARMF</td><td>0.009</td><td>0.008</td><td><strong>0.023</strong></td><td>0.017</td><td><strong>11.8</strong></td><td>9.2</td><td><strong>14.3</strong></td><td>10.0</td><td><strong>6.5</strong></td><td>1.9</td><td><strong>5.9</strong></td><td>2.4</td></tr></tbody></table>
<figure><img src="image_12.png" /><figcaption>Figure 11: NDVI Temporal profile from the gap-filling experiment, for an agricultural pixel of AOI 31TCJ, with predictions from Naive, U-TILISE and TAMRF. Images show reference and predictions for unobserved dates at start, peak and end of season, where color composition is natural colors (B4, B3, B2). Blue (resp. orange) dots denote clear (observed) NDVI for Sentinel-2 (resp. Landsat-8/9). Stars denote masked (unobserved) NDVI.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 12: TAMRF compared to U-TILISE and Naive baseline in interpolating a target date in generated gaps, for tile 31TCJ, for Sentinel-2, 2022-04-10 (top row) and Landsat, 2022-02-10 (bottom row). Sentinel-2 color composition is (B07, B06, B05), while Landsat-8/9 color composition is natural colors (B4, B3, B2). Interpolation artifacts caused by mask transition in the naive method are especially visible on the RGB output. Moreover, the 2022-02-10 prediction of the bottom row is a difficult case, since the closest image (either Landsat-8/9 or Sentinel-2) is 16 days away in time. This explains errors in predicting the state of some agricultural areas, for both TAMRF and naive method.</figcaption></figure>
located in gaps for tile 31TCJ. Note that the prediction task is very difficult for the bottom row example, since the closest clear date is 16 days away in time. This explains the failure for both Naive and TAMRF to predict the state of some agricultural fields. The figure highlights an important benefit of TAMRF that is not reflected in RMSE: naive temporal interpolation relies on Level 2 masks. Depending on those masks, neighboring pixels will be interpolated with different dates, yielding artifacts following the masks delineation, which are clearly visible on the bottom row. TAMRF, on the other hand, does not use the masks as input, and therefore yields images that are free from those artifacts. Unfortunately, the published U-TILISE model can not process the Sentinel-2 20-meter bands, nor the Landsat-8/9 bands.

In addition to predicting consistent surface reflectance values with respect to the reference image, TAMRF improves the spatial resolution of the spectral bands for both Sentinel-2 and Landsat-8/9. In contrast, the naive model only performs spatial bicubic up-sampling, yielding blurry 10-meter resolution data. This highlights the benefits of using TAMRF as both a temporal interpolation method and a spatial resolution enhancement method.

### 3.3.2. Band-Sharpening

The second evaluation task consists in sharpening the 20-meter Sentinel-2 bands (see Table 4) to 10-meter resolution. TAMRF systematically outputs sharpened predictions for 20-meter Sentinel-2 bands, regardless of whether the query date has been observed by the model or not. In this experiment, it is compared to the DSen2 model Lanaras et al. (2018), a network that performs similar Sentinel-2 20-meter bands sharpening for a given Sentinel-2 image. Though this is the only model with publicly available source code and model weights found in the literature, it seems to be representative of the performances of more recent DL-based Sentinel-2 sharpening approaches Wu et al. (2023); Ciotola et al. (2025).

DSen2 is applied to each Sentinel-2 acquisition throughout the testing set. On the other hand TAMRF is fed with the full Sentinel-2 SITS for a given AOI, without the Landsat-8/9 corresponding SITS, and is asked to predict all the input Sentinel-2 dates. In addition to sharpening the 20-meter bands, it will therefore also remove clouds and
no-data areas. The three evaluation metrics are computed with respect to the input
20-meter spectral bands for both methods, taking into account their validity mask to
compute RMSE.

Table 8 shows the comparative performances obtained for a subset of the spectral bands (full results are available in appendix 6), namely one Red Edge band, the Narrow Near Infra-Red band and one of the SWIR bands. It shows that TAMRF is globally on par with DSen2. It tends to yield slightly higher RMSE with respect to the input image, but better IQ according to the BRISQUE score, and slightly better sharpening according to the FR score. This is confirmed by the visual inspection of Fig. 13: Dsen2 predictions look less sharp that those from TAMRF, even if both manage to improve on the bicubic up-sampling of the input Sentinel-2 image. The very good RMSE performances of Dsen2 come from its training on simulated data, using the Wald protocol and scale invariance hypothesis. This is known to yield very consistent radiometric values, since the training data have no distortion, at the expense of lesser IQ and spatial resolution improvement caused by the domain gap between simulated and real data Zhao et al. (2024). Nevertheless, the slightly higher RMSE of the TAMRF prediction is barely noticeable visually. An additional benefit of using TAMRF over DSen2 is of course its ability to seamlessly interpolate cloudy pixels, as demonstrated in Fig. 14.

3.3.3. Spatio-Temporal Fusion

The next task aims at predicting Sentinel-2 at times when only Landsat-8/9 has been observed, by optionally leveraging Sentinel-2 acquisitions from other time steps. Using the LS2S2 test set, this use case is simulated by removing all Sentinel-2 acquisitions occurring on same date as Landsat-8/9 from the model input and keeping them as references for validation. The performances are compared to STAIR Luo et al. (2018) and to our own implementation of the Sen2Like Saunier et al. (2022) fusion rule, which are rule-based fusion methods, as we as to Deep-Harmonization Sambandham et al. (2024) and to DSTFN Wu et al. (2022), the latter being DL based methods. All those methods share the same assumption of matching bands between Landsat-8/9 and Sentinel-2. As such, they do not predict the Red Edge bands, nor the wide Near Infra Red band (see Table 4).
<table><caption>Table 8: Comparison between TAMRF and DSen2 on the sharpening of Sentinel-2 20-meter bands, for a selection of bands. Snow target dates are excluded from metrics computation. ↓ (resp. ↑) indicates that the metric should be minimized (resp. maximized). Best mean values for each band and each metric are highlighted in **bold**. Full results are available in appendix 6.</caption><thead><tr><th rowspan="2">Band</th><th rowspan="2">Method</th><th colspan="2">RMSE↓</th><th colspan="2">BRISQUE↓</th><th colspan="2">FR↑</th></tr><tr><th>mean</th><th>std</th><th>mean</th><th>std</th><th>mean</th><th>std</th></tr></thead><tbody><tr><td rowspan="2">B6</td><td>DSEN2</td><td><strong>0.002</strong></td><td>0.001</td><td>25.907</td><td>6.494</td><td>4.4</td><td>0.6</td></tr><tr><td>TAMRF</td><td>0.011</td><td>0.007</td><td><strong>14.18</strong></td><td>9.401</td><td><strong>5.5</strong></td><td>1.7</td></tr><tr><td rowspan="2">B8a</td><td>DSEN2</td><td><strong>0.003</strong></td><td>0.001</td><td>25.390</td><td>6.360</td><td>4.5</td><td>0.5</td></tr><tr><td>TAMRF</td><td>0.012</td><td>0.009</td><td><strong>14.775</strong></td><td>9.365</td><td><strong>5.3</strong></td><td>1.6</td></tr><tr><td rowspan="2">B12</td><td>DSEN2</td><td><strong>0.002</strong></td><td>0.001</td><td>33.449</td><td>6.225</td><td>3.9</td><td>0.7</td></tr><tr><td>TAMRF</td><td>0.009</td><td>0.008</td><td><strong>11.858</strong></td><td>9.011</td><td><strong>6.5</strong></td><td>2.0</td></tr></tbody></table>

<figure><img src="image_14.png" /><figcaption>Figure 13: Examples of predictions from DSEN2 and TAMRF over tile 31TFJ, 2022-08-07, for Red-Edge color composition (B07, B06, B05) and SWIR color composition (B8a, B11, B12). Inputs are captioned in green, predictions are captioned in red.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure 14: Examples of predictions from DSEN2 and TAMRF over tile 31TFJ, 2022-09-29, where the input Sentinel-2 image is cloudy, for Red-Edge color composition (B07, B06, B05) and SWIR color composition (B8a, B11, B12). Inputs are captioned in green, predictions are captioned in red.</figcaption></figure>

Sen2Like Saunier et al. (2022) is a processing chain providing harmonized surface reflectances by means of spatial registration and careful radiative transfer modelling. Though the standard output is produced at 30-meter resolution, an optional fusion rule is proposed in order to up-sample the prediction to 10-meter resolution. For a given Landsat-8/9 target date, it consists in finding the two closest Sentinel-2 anterior dates with low cloud cover, and build a linear model to interpolate them at the target date. The predicted image is then High Pass filtered to extract the high resolution details, which are added to the target Landsat-8/9 image. Note that though in Sen2Like the 2 dates are searched among the previous predictions in an autoregressive fashion, in our work they are searched among the available Sentinel-2 dates. Another key difference with the original Sen2Like implementation is that we up-sample the 20-meter Sentinel-2 band to 10-meter spatial resolution by means of bicubic operator, for a fair comparison with the other methods that effectively improve the spatial resolution of the 20-meter bands. As such, we do not compare with the full Sen2Like process but only with the temporal extrapolation and spatial fusion methods implemented in Sen2Like. STAIR Luo et al.
(2018) is another processing chain aiming at merging Sentinel-2 and Landsat images. It first performs a linear gap-filling of the Sentinel-2 SITS and Landsat-8/9 SITS independently, as well as bicubic up-sampling of Landsat-8/9 SITS to 10-meter spatial resolution. It then uses same-day pairs in order to build a time-series of differences between Landsat-8/9 and Sentinel-2 bands. This time-series is linearly interpolated at prediction dates where only Landsat-8/9 is available, and added to the Landsat-8/9 image. Deep-Harmonization Sambandham et al. (2024) is an ensemble of five Single Image Super-Resolution Unets trained to map individual Landsat-8/9 acquisition, including the panchromatic channel, to Sentinel-2 resolution and spectral bands. The ensemble includes instances of the same Unet architecture with two different depths and two different up-sampling layers. Being a SISR approach, Deep-Harmonization does not use Sentinel-2 as input during inference and process the Landsat-8/9 image independently for each date. DSTFN Wu et al. (2022) is a Spatio-Temporal Fusion model that does not require 2 pairs of same-day acquisitions framing the target date. Instead, DSTFN only requires the Landsat-8/9 image at target date and an auxiliary Sentinel-2 image at a close date. To compete with these 4 methods, the TAMRF model can be used in different ways. First, it can process the Landsat-8/9 images only, in a similar setup to Deep-Harmonization. Second, it can work on Sentinel-2 images only, resorting to the learnt temporal interpolation in order to perform the task. Finally, it can receive all possible Landsat-8/9 images and all non-simultaneous Sentinel-2 images.

Table 9 shows the performances of the 4 competing algorithms as well as the 3 TAMRF runs with the different input configurations: only LR SITS, only HR SITS and both LR and HR SITS. First, it can be observed that TAMRF systematically yields better RMSE values than the other algorithms. STAIR, Sen2Like fusion method and DSTFN have the worst RMSE performances, while Deep-Harmonization is closer to TAMRF. The Sen2Like fusion algorithm expects carefully harmonized surface reflectances between Landsat-8/9 and Sentinel-2. In this experiment, the fusion process of Sen2Like is instead applied to raw Level 2 surface reflectances, which can explain its poor performance. It can be observed that TAMRF provides good estimates of Sentinel-2 surface reflectances even in the case where only Landsat-8/9 images are seen by the model (TAMRF no hr). With respect to using only the HR images for the
prediction of TAMRF (TAMRF no lr), adding the Landsat-8/9 acquisition of the tar-
get date to the inputs allows to systematically decrease both the mean RMSE and its
standard-deviation, which shows that TAMRF makes use of the Landsat-8/9 low reso-
lution information in the prediction of Sentinel-2 bands. In terms of general image IQ,
TAMRF also always provides the lowest BRISQUE score, surpassing all methods by a
margin of more than 10 BRISQUE points, except for B4 where STAIR is on par with
TAMRF. Among the three variants of TAMRF inputs, using only Landsat-8/9 image
yields worse IQ, on par with the other methods, and using the Landsat-8/9 image of the
target date in addition to the Sentinel-2 SITS yields slightly worse IQ than using only
the Sentinel-2 SITS, albeit with lower RMSE. Finally, in terms of spatial resolution
enhancement, it can be observed that for the 10-meter bands, all methods manage to
retrieve the target resolution (|FR|≈ 0.), except for Deep-Harmonization, which fails to
restore the spatial frequency content of a 10-meter image despite being a SISR method
at its heart. For 20-meter bands however, only TAMRF and to a lower extent DSTFN
manage to provide a spatial frequency content improvement.

Fig. 15 shows sample predictions for all algorithms and provides a visual confir-
mation of those conclusions. The natural color compositions appear sharp except for
Deep-Harmonization, which is both blurry and radiometrically inaccurate. As mea-
sured in Table 9, IQ is poor for Sen2Like and DSTFN. Deep-Harmonization IQ is on
par with TAMRF using only Landsat-8/9 images. Looking a SWIR color compositions,
TAMRF is the only algorithm that provide sharpness and good IQ. DSTFN is sharper
than the other competing algorithms, but bears artifacts similar to aliasing. On STAIR
prediction, an artifact caused by temporal interpolation based on cloud masks can be
observed around the largest of the two lakes. In conclusion, TAMRF provides the best
performances in this spatio-temporal task. It can use only the Landsat-8/9 images and
yield already consistent results, but the best performances are obtained by combining
the Landsat-8/9 image at target date with the SITS of Sentinel-2 images observed at
other instants. Nevertheless, the same pre-trained TAMRF model can do both, and
even only use Sentinel-2 SITS for prediction. Another benefit of TAMRF is that it han-
dles seamlessly the case where the Landsat-8/9 image at target date is cloudy, without
using any additional mask. It is worth noting that none of the competing methods can
<table><caption>Table 9: Comparison of the different methods on the spatio-temporal fusion task. Only a subset of bands is presented. Full results are available in appendix 6. * if usually higher FR means higher spatial resolution details, in the case of bands that are 10-meter native resolution, the best FR is the closest to zero in absolute value. Otherwise, ↓ (resp. ↑) indicates that the metric should be minimized (resp. maximized). Best mean values for each band and metric are highlighted in **bold**, and second best values are **underlined**. † SEN2LIKE does not correspond to the full Sen2Like processing but only to the fusion part of the pipeline, which has been reimplemented following the paper and source code of the original pipeline.</caption><thead><tr><th rowspan="2">Band</th><th rowspan="2">Method</th><th colspan="2">RMSE↓</th><th colspan="2">BRISQUE↓</th><th colspan="2">FR↑</th></tr><tr><th>mean</th><th>std</th><th>mean</th><th>std</th><th>mean</th><th>std</th></tr></thead><tbody><tr><td rowspan="7">B4</td><td>SEN2LIKE†</td><td>0.032</td><td>0.030</td><td>30.1</td><td>17.8</td><td>-0.5*</td><td>2.8</td></tr><tr><td>STAIR</td><td>0.038</td><td>0.040</td><td><u>22.1</u></td><td>11.4</td><td>-0.7</td><td>2.6</td></tr><tr><td>DEEPHARMO</td><td>0.027</td><td>0.010</td><td>44.2</td><td>5.4</td><td>-5.3</td><td>1.0</td></tr><tr><td>DSTFN</td><td>0.033</td><td>0.035</td><td>28.4</td><td>8.7</td><td><u>-0.0</u>*</td><td>2.3</td></tr><tr><td>TAMRF (no hr)</td><td>0.019</td><td>0.009</td><td>49.8</td><td>9.2</td><td>-3.5</td><td>2.5</td></tr><tr><td>TAMRF (no lr)</td><td><u>0.016</u></td><td>0.019</td><td>22.6</td><td>12.6</td><td>0.7</td><td>2.5</td></tr><tr><td>TAMRF (full)</td><td><b>0.014</b></td><td>0.008</td><td><b>21.4</b></td><td>12.8</td><td>0.6</td><td>2.5</td></tr><tr><td rowspan="7">B8a</td><td>SEN2LIKE†</td><td>0.048</td><td>0.044</td><td>58.4</td><td>12.4</td><td>-0.5</td><td>2.3</td></tr><tr><td>STAIR</td><td>0.056</td><td>0.050</td><td>48.4</td><td>5.9</td><td>-0.5</td><td>2.4</td></tr><tr><td>DEEPHARMO</td><td><u>0.026</u></td><td>0.003</td><td>62.4</td><td>4.7</td><td>-3.0</td><td>0.4</td></tr><tr><td>DSTFN</td><td>0.049</td><td>0.048</td><td>25.3</td><td>6.8</td><td>2.9</td><td>2.4</td></tr><tr><td>TAMRF (no hr)</td><td>0.031</td><td>0.019</td><td>47.4</td><td>5.5</td><td>1.1</td><td>1.6</td></tr><tr><td>TAMRF (no lr)</td><td><u>0.026</u></td><td>0.020</td><td><b>15.6</b></td><td>8.4</td><td><b>5.1</b></td><td>1.8</td></tr><tr><td>TAMRF (full)</td><td><b>0.024</b></td><td>0.013</td><td><u>17.4</u></td><td>9.5</td><td><u>4.8</u></td><td>1.8</td></tr><tr><td rowspan="7">B12</td><td>SEN2LIKE†</td><td>0.041</td><td>0.029</td><td>66.8</td><td>13.2</td><td>-0.0</td><td>1.9</td></tr><tr><td>STAIR</td><td>0.042</td><td>0.034</td><td>53.1</td><td>6.4</td><td>-0.1</td><td>1.9</td></tr><tr><td>DEEPHARMO</td><td><u>0.021</u></td><td>0.008</td><td>65.3</td><td>7.3</td><td>-2.3</td><td>0.6</td></tr><tr><td>DSTFN</td><td>0.044</td><td>0.032</td><td>33.2</td><td>9.2</td><td>3.5</td><td>2.1</td></tr><tr><td>TAMRF (no hr)</td><td>0.022</td><td>0.012</td><td>46.6</td><td>6.5</td><td>1.9</td><td>1.7</td></tr><tr><td>TAMRF (no lr)</td><td><u>0.021</u></td><td>0.016</td><td><b>13.6</b></td><td>8.4</td><td><b>5.9</b></td><td>1.8</td></tr><tr><td>TAMRF (full)</td><td><b>0.019</b></td><td>0.012</td><td><u>16.7</u></td><td>9.4</td><td><u>5.6</u></td><td>1.8</td></tr></tbody></table>
make a prediction in this case, since they heavily rely on the LR observation at target
date.

<figure><img src="image_16.png" /><figcaption>Figure 15: Sentinel-2 predictions of the Spatio-Temporal Fusion task for all methods, on tile 31TCJ, 2022-05-10. First two rows shows the visible (B04, B03, B02) color composition, and last two rows shows the SWIR (B8a, B11, B12) color composition.</figcaption></figure>

3.3.4. *Thermal Sharpening*

This last task focuses specifically on sharpening the LST band of Landsat-8/9, for
observed dates. In this task, all available Sentinel-2 and Landsat-8/9 dates are fed to
TAMRF, and the model is asked to predict the clear Landsat-8/9 images. As already
explained, predicting LST for unobserved dates is not feasible without meteorological
conditioning, since meteorological factors are the main drivers for LST and are not
observed in SITS (see section 5.1.2).

We measure specifically the performances of the LST band. As a baseline for comparison, the Data Mining Sharpener (DMS) Gao et al. (2012) is used. DMS is widely used in the Thermal Infra Red community, especially for operational productions. In these experiments, we used the pyDMS<fnref n="3" /> implementation. In order to provide the High Resolution input to DMS, the closest Sentinel-2 image that meets a maximum cloud coverage criterion is selected. Only the 10-meter bands are given to DMS, and the default parameters of the method are used.

It must be stressed that like most Thermal Sharpening methods, DMS uses the residual compensation mechanism: lower resolution differences between the machine learning based prediction and the input LST image are up-sampled to the target resolution and injected in the prediction, following:

$$S_{LR}^{\star rc} = S^{\star} + (S^{\star} \downarrow - S_{LR}) \uparrow, \qquad (26)$$

where $\downarrow$ is a down-sampling operator, $\uparrow$ is an up-sampling operator, and $S_{LR}^{\star rc}$ is the prediction with residual compensation. Injecting back LR residual in the prediction ensures that the prediction dynamic aligns with the LR input image, yielding lower RMSE$_{LR}$ values. This process however is not without drawbacks, as it injects back an LR signal into the prediction, diminishing its high spatial frequency content and general IQ. It also can not be applied for cloudy pixels for instance. This residual compensation can be seen as a post-processing, and can be applied to any Thermal (or non-thermal) sharpening method, here we compare results of DMS and TAMRF both with and without residual compensation.

Table 10 shows the performances of TAMRF and DMS on the thermal sharpening task, with and without residual compensation, evaluated on the **all** testing set, with snowy target dates excluded. The first two rows present the results of both models without the residual compensation. In this case, DMS provides a better RMSE value, more than halving the error made by TAMRF. However, the latter provides better IQ

<footnote n="3">³https://github.com/radosuav/pyDMS</footnote>
and spatial frequency content. The last two rows present the results of both algorithms, this time using residual compensation. As expected it allows to significantly lower the RMSE of both methods, and in this case TAMRF performs better than DMS. Using this compensation slightly diminishes the IQ for TAMRF, while still being better than both runs of DMS. It also has a small impact on FR, which is expected because of the injection of LR signal in the prediction. Nevertheless, TAMRF remains the best method for this metric.

<table><caption>Table 10: Comparison between DMS and TAMRF, with and without residual compensation, on the Thermal Sharpening task of the LST (expressed in Kelvins). Best results for each metric are highlighted in **bold**. Second-best results are *underlined*. ↓ (resp. ↑) indicates that the metric should be minimized (resp. maximized).</caption><thead><tr><th rowspan="2">Method</th><th colspan="2">RMSE↓</th><th colspan="2">BRISQUE↓</th><th colspan="2">FR↑</th></tr><tr><th>mean</th><th>std</th><th>mean</th><th>std</th><th>mean</th><th>std</th></tr></thead><tbody><tr><td>DMS</td><td>0.772</td><td>0.404</td><td>28.769</td><td>25.008</td><td>0.6</td><td>5.6</td></tr><tr><td>TAMRF</td><td>1.772</td><td>2.075</td><td><strong>20.941</strong></td><td>23.156</td><td><strong>2.8</strong></td><td>2.9</td></tr><tr><td>DMS + res. comp.</td><td><u>0.477</u></td><td>0.921</td><td>27.720</td><td>25.237</td><td>0.4</td><td>5.4</td></tr><tr><td>TAMRF + res. comp.</td><td><strong>0.381</strong></td><td>1.441</td><td><u>22.181</u></td><td>21.638</td><td><u>2.5</u></td><td>1.0</td></tr></tbody></table>

Fig. 16 shows sample predictions for both methods. Without residual compensation, the two methods provide consistent sharpening of the spatial resolution, but TAMRF overestimates the LST dynamics while DMS tends to underestimate them. When residual compensation is added, the IQ is noticeably worse, but the dynamics of the LST is preserved. The LR residuals are especially visible in the **DMS + res. comp.** example (lower right image). Finally, several aliasing artifacts can be observed in DMS predictions with and without the compensation, which explains the lower BRISQUE score.

It should be noted that it is always possible to add residual compensation to predictions of other bands, for pixels that are observed and are not cloudy. However, given the performances of TAMRF on the other spectral bands and the degradation of IQ and FRR it incurs, it seems unwise to do so.
<figure><img src="image_17.png" /><figcaption>Figure 16: Predictions of the Thermal Sharpening task for DMS and TAMRF, on tile 31TCJ, 2022-06-03.
From left to right and top to bottom: Landsat-8/9 RGB color composition (B4, B3, B2), Landsat-8/9 LST sharpened by TAMRF without residual compensation (B10), Landsat-8/9 LST sharpened by DMS without residual compensation (B10), Landsat-8/9 LST (B10), Landsat-8/9 LST sharpened by TAMRF with residual compensation (B10), Landsat-8/9 LST sharpened by DMS with residual compensation (B10).</figcaption></figure>
### 3.4. Summary of experiments

In this section, a single pre-trained TAMRF model has been used to solve 4 different downstream tasks. For each of these tasks, the pre-trained model has been compared to algorithms from the literature which have been specifically designed to solve the task at stake. The experiments show that TAMRF is on par or better than those algorithms in all the tasks. Those experiments also demonstrate the benefits of using a single versatile model: cloud removal when performing Sentinel-2 20-meter bands sharpening, or the ability to input Landsat-8/9 SITS, Sentinel-2 SITS or both and always get consistent predictions without modifying and retraining the model. It is important to keep in mind that results presented in this section are measured in the frame of limitations imposed by the dataset as well as the competing algorithms. For instance, RMSE is always measured on pixels of acquired reference dates that are not masked, but TAMRF also produce consistent outputs even when and where there is no reference value to compare to. In the spatio-temporal fusion task, only 6 Sentinel-2 bands are evaluated because the competing algorithms only process those 6 bands, but TAMRF predicts the 10 Sentinel-2 bands and the 8 Landsat-8/9 bands.

## 4. Ablation and sensitivity studies

This section presents various ablation and sensitivity studies in order to assess the role of each element of the proposed method, as well as to better understand the model capabilities. All the variants of the model are trained with the full training procedure presented in section 3.2.1.

Unless stated otherwise, the different variants of the model are evaluated on the same random masking task: for each SITS of the **all** testing set, 50% of the Sentinel-2 dates and 50% of the Landsat-8/9 are removed from the model input and used as reference data to measure performances. Contrary to section 3.3, all target dates of the testing set are used, including snowy target dates.

### 4.1. Model size

In order to evaluate the sensitivity of the model performances to the model hyper-parameters selected in Table 6, in this section two important macroscopic parameters
are varied: the number of RRDB blocks in the spatial encoders on one hand, and the number of Transformer encoder layers in the temporal encoder on the other. Tested configurations are summarized in Table 11.

<table><caption>Table 11: Key parameter changes with respect to Table 6 in order to investigate the influence of model size on the performances.</caption><thead><tr><th>Name</th><th>Key parameters change</th><th>Value</th><th>#Parameters</th></tr></thead><tbody><tr><td>baseline</td><td>N.A.</td><td></td><td>2.3M</td></tr><tr><td>s2t3</td><td>LR Spatial encoder nb. blocks</td><td>2</td><td>3.8M</td></tr><tr><td></td><td>HR Spatial encoder nb. blocks</td><td>2</td><td></td></tr><tr><td>s3t3</td><td>LR Spatial encoder blocks</td><td>3</td><td>5.2M</td></tr><tr><td></td><td>HR Spatial encoder blocks</td><td>3</td><td></td></tr><tr><td>s1t1</td><td>Temporal encoder nb. layers</td><td>1</td><td>2M</td></tr><tr><td>s1t2</td><td>Temporal encoder nb. layers</td><td>2</td><td>2.2M</td></tr></tbody></table>

### 4.1.1. Temporal encoder depth

Figure 17 shows the impact of reducing the number of layers in the temporal transformer encoder on the performances of reconstruction of masked dates. As expected, reducing the depth of the temporal encoder yields higher RMSE values on masked dates, for all bands. Interestingly, the performance degradation is higher for Landsat-8/9 bands, while for many Sentinel-2 bands, it is barely noticeable. If someone was to build a TAMRF model that specifically targets Sentinel-2 reconstructions from both Sentinel-2 and Landsat-8/9, reducing the number of layers in the temporal encoder would yield a lighter model with similar performances. For the all-purposes model proposed in this paper however, the **baseline** value of 3 layers in the temporal encoder seems appropriate. Increasing the number of layers further may yield even better performances, at the expense of an heavier model.

### 4.1.2. Spatial encoders depth

Figure 18 shows the impact of increasing the number of RRDB blocks in the HR and LR spatial encoders on the reconstruction of masked dates. Deeper spatial encoders
<figure><img src="image_18.png" /><figcaption>Figure 17: Impact of the reduction of the number of layers in temporal encoder on the quality of reconstruction of masked dates. **baseline** has 3 cross-attention layers, **s1t2** has two cross-attention layers and **s1t1** has only one cross-attention layer.</figcaption></figure>

seem to have no significant effect on the RMSE, though the **s2t3** run with 2 blocks seems to perform slightly better than the others. However, according to BRISQUE score this configuration has slightly worse IQ, and less frequency restoration according to FR. Increasing the number of blocks to 3 does not provide any significant performance improvement and yields slightly higher BRISQUE and lower FR scores with respect to the **baseline** configuration, which seems to be the appropriate trade-off for this parameter.

## 4.2. Loss Terms

This experiment intends to highlight the role of the spatial term and contrastive term in the total loss of Eq. 11. Each term is ablated separately, forming the **w/o spatial** (no spatial term) and **w/o contrastive** (no contrastive term) experiments, compared to the **baseline** run which uses the full expression of Eq. 11. All other parameters remain the same as presented in section 3.2.

### 4.2.1. Spatial loss term

Figure 19 shows the impact on reconstruction performances of training the TAMRF model without the spatial loss term of Eq. 11. It can be observed that the **w/o spatial** run has lower performances for almost all bands from both Sentinel-2 and Landsat-8/9 and all the metrics. RMSE is slightly better with the spatial term for all bands, the greatest impact of the spatial term is on the BRISQUE score and FR of bands
<figure><img src="image_19.png" /><figcaption>Figure 18: Impact of the increase of the number of blocks in spatial encoders on the quality of reconstruction of masked dates. **baseline** has one RRDB block per encoder, **s2t3** has two RRDB per encoder, and **s3t3** has three RRDB per encoder.</figcaption></figure>
which have a native spatial resolution coarser than 10 meters (all Landsat-8/9 bands and bands B5, B6, B7, B8a, B11 and B12 of Sentinel-2). For those bands there is a dramatic drop in IQ when the model is trained without the spatial term, of more than 10 BRISQUE points for Sentinel-2 20-meter bands and more than 30 BRISQUE points for some Landsat-8/9 bands. The analysis of FR allows to further explain this drop in IQ, since it shows that for those bands the high spatial frequency content of the predicted images is far lower when the spatial term is not used in training: in this case, the model simply fails to improve the effective spatial resolution and merely perform a zooming operation. This is illustrated in Fig. 20, which shows an example of prediction of a Landsat-8/9 date, with both the **baseline** and the **w/o spatial** models. The prediction of the model trained without the spatial term is blurry and exhibits Gibbs artifacts on the edges of some of the agricultural fields. On the other hand, the prediction of the **baseline** model has a better image quality, with a high spatial frequency content similar to a Sentinel-2 images, while being radiometrically faithful to the Landsat-8/9 input image.

### 4.2.2. Contrastive loss term

Figure 21 shows the impact on the reconstruction performances of ablating the contrastive loss term of Eq. 11. It compares the RMSE of clear dates, that are used as input to the model. Contrary to other results presented in this article, this time the RMSE is computed only on invalid pixels, that are masked according to the cloud and quality masks of the Level 2 products. The figure shows that when the contrastive loss term is not used during training, the predictions of TAMRF on observed dates have lower RMSE on invalid pixels, which suggests that TAMRF reconstructs some clouds, cloud shadows and other invalid areas of the input images instead of making clear predictions of the ground beneath them. This is confirmed by the examples shown in Fig. 22: the **w/o contrastive** model leaks clouds and cloud shadows in the prediction, where the **baseline** model learned to ignore clouds and cloud shadows and predicts the surface reflectance instead.
<figure><img src="image_20.png" /><figcaption>Figure 19: Comparison of performances between the **baseline** model and the **w/o spatial** model, where the spatial term has been ablated from Eq. 11. Metrics are computed on all target dates, either clear or masked.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Figure 20: Comparison between the **baseline** model and the **w/o spatial** model where the spatial loss term is ablated, on the prediction of an observed Landsat-8/9 image, acquired on 2022-08-05. Color composition is natural colors (B4, B3, B2).</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>Figure 21: Comparison of the RMSE of **baseline** and **w/o contrastive** models on invalid pixels of clear dates (used as input to the model). Invalid pixels correspond to pixels masked in the cloud masks and quality masks of the Level 2 products.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Figure 22: Examples of predictions of clear Sentinel-2 dates (used as input to the model) that are partly cloudy, on tile 31TCJ, for 2022.03.31 (top row) and 2022.04.15 (bottom row), with model **baseline** and **w/o contrastive**. Color composition is natural colors (B4, B3, B2).</figcaption></figure>
### 4.3. MAE Strategy

This experiment investigates the benefit of diversifying the masking strategy with the Gaps, No HR, No LR, Forecast and Backcast configurations, as presented in Table 5. The **baseline** configuration is compared to a run called **plain mae** which solely uses the classical random masking strategy with a discard rate following the uniform distribution $U(0.2, 0.7)$.

In order to understand the benefit of the diversification of the masking strategies, in this section we evaluate both models on two different tasks of interest for end-users. The first task is a forecasting task: since all SITS in the testing set have a temporal length of one year, the models receive as input the first half year of each SITS and are asked to reconstruct the dates in the second half of the year. The second task aims at reconstructing the Sentinel-2 dates using only the Landsat-8/9 dates as input. Note that with the plain random masking strategy applied for training the **plain mae** model, neither of those scenarios will be seen during training: the range of discard rates forbids masking the Sentinel-2 dates entirely, and even if the range was extended, the odds of drawing this configuration are very low. This also applies to the forecasting scenario, which is very unlikely to occur in plain random masking. Looking at the performances of both models on these two tasks presented in Fig. 23, it can be observed that the **baseline** model consistently performs better on those two tasks than the **plain mae** model. We can conclude that favoring the occurrence of interesting configurations for downstream applications during MAE training allows increasing the model performances on those tasks.

### 4.4. Learnable sensor tokens

The last ablation study in this section investigates the role of the learnable sensor tokens presented in section 2.2.1. The size of sensor tokens $M$ is set to zero, forming the **w/o sensor tokens** run. The results of this ablation are presented in figure 24. It can be observed that the use of learnable sensor tokens is mainly beneficial to the reconstruction of Landsat-8/9 dates, allowing for marginal improvement of the RMSE for Landsat-8/9, at the expense of a slight decrease in FR. On the other hand, it marginally degrades the RMSE of Sentinel-2 bands, albeit to a lesser extent than the gain in RMSE
<figure><img src="image_24.png" /><figcaption>Figure 23: Performances of the **baseline** and **plain mae** models on the forecasting task and missing Sentinel-2 sensor tasks.</figcaption></figure>
for Landsat-8/9. We can infer that sensor tokens allow the model to better distinguish between Sentinel-2 and Landsat-8/9 input dates, favoring the reconstruction of the latter which is less represented in the training set due to its revisit of 16 days for each satellite. Since this improvement is only marginal, it is probable that Sentinel-2 and Landsat-8/9 activate different parts of the latent space, with or without sensor tokens. Nevertheless, the cost of including the learnable sensor tokens is very small with respect to their benefits, as they only add 24 extra learnable parameters and extend the length of the latent representation from 128 to 136 features.

<figure><img src="image_25.png" /><figcaption>Figure 24: Comparison of the performances of **baseline** and **w/o sensor tokens** models on all dates (either clear or masked), on the random masking experiment.</figcaption></figure>

### 4.5. *Impact of snow*

As mentioned in the introduction of section 3.3, target dates with snow were excluded from the metrics computation when comparing TAMRF capabilities in order to provide results that can be compared to other studies which have no snow in their
validation sets. In this section we discuss the quality of the reconstructions of target dates containing snow with TAMRF. For this purpose, we split the target dates into two categories: target dates containing no snow pixels according to Level 2 quality masks on one hand, and target dates containing at least one snow pixel according to those masks on the other. In addition to computing the standard RMSE, we also compute the Normalized RMSE (NRMSE), where RMSE is divided by the standard deviation of the reference data.

Figure 25 shows the RMSE and NRMSE for both categories, separately for clear (available in model inputs) and masked (not available in model inputs) dates. As expected, masked dates with snow are very difficult to predict, and yield large RMSE. Snow is driven by external meteorological factors: if the presence or absence of snow in nearby dates can be a clue for the model, isolated snowfalls or meltdown can not be forecasted by using solely the SITS. This is another incentive to add meteorological conditioning in TAMRF, which is discussed in section 5.1.2. When it comes to observed snowy dates, the RMSE degradation is still large, especially in the visible and red-edge domains. This can be partly explained by the larger dynamic of the surface reflectance of snow in this part of the spectrum, which is often close to 1., while the typical surface reflectance is 0.15. On the other hand, SWIR bands such as B11 and B12 from Sentinel-2 exhibit lower drop in RMSE performances. This is confirmed by the NRMSE trends, which show that a large part of the RMSE drop in visible and red-edge bands is explained by the higher standard deviation of the target reflectance.

Therefore, it can be said that TAMRF has difficulties in predicting unobserved dates with snow, and can predict observed dates with snow with a larger RMSE than for snow-free dates. However, this larger RMSE is in part explained by the higher standard deviation of snow pixel surface reflectances, especially in the visible and red-edge range of the spectrum. Figure 26 illustrates those conclusions. On the first column, it shows the prediction with respect to the reference on an observed Sentinel-2 date, in B2, B3, and B4 (RGB colors). Snowy areas are well handled, but radiometric differences can be observed. On the second column, it shows the prediction with respect to the reference on a date that has not been observed by the model. In this example, TAMRF completely misses the snow and predicts bare soil and vegetation surface
reflectance instead.

Finally, we would like to stress again that if TAMRF may fail to predict unobserved
snow dates in some cases, this is probably the case for most of the methods from the
literature. Moreover, TAMRF has been trained with all dates including snow dates,
and in all cases presented in this paper the snow dates are included in the model inputs.
Snow dates do not prevent the proper training of the model, and the trained model is
able to handle snow input dates correctly and distinguish them from clouds, without
using any of the Level 2 quality masks at inference time.

<figure><img src="image_26.png" /><figcaption>Figure 25: RMSE and NRMSE of TAMRF on target dates with no snow and with snow, for clear (used as input to the model) and masked (not seen in the model input) dates.</figcaption></figure>

4.6. *Number of input dates*

Finally, this section investigates the influence of the number of input images in
the reconstruction performances of TAMRF for unobserved dates. First, a prediction
period of two months is chosen at the center of each SITS for each AOI. Because AOIs
from the **wwmy** testing set have variable time span (see Fig. 8), these 2 months actually
span different period of the year between 2022.01 and 2025.06. The TAMRF model is
asked to reconstruct the masked dates, using series of neighboring dates of increasing
sizes, from 4 neighboring dates per sensor (2 on each side of the prediction period) to
64 neighboring dates per sensor (32 on each side of the prediction period). Only the
<figure><img src="image_27.png" /><figcaption>Figure 26: Examples of prediction of snowy target dates for AOI 51VWH_42 located in eastern Russia. Date 2023.04.20 is part of the model input. Date 203.05.10 is not part of the model input and the closest input date is 10 days away.</figcaption></figure>
AOIs of the **all** testing set that have at least 48 dates to fill the 64 neighboring dates per sensor case are retained for this experiment. This process is illustrated in Fig. 27.

The RMSE of the reconstruction of Landsat-8/9 and Sentinel-2 target dates in the prediction period is presented in Fig. 28. Two phases can be observed for Sentinel-2 reconstructions. The RMSE first improves rapidly with the number of dates, especially for Sentinel-2. It then reaches a plateau with a slight increasing slope after 32 input dates per sensor. For Landsat-8/9 reconstructions the impact of the number of input dates is not as dramatic as for Sentinel-2, though there is also a marginal minimum followed by slowly increasing plateau for most of the bands. Interestingly this minimum seems to happen around 16 dates per sensor, which is less than for Sentinel-2. However, the scale of the increasing slope is very small for both sensors, and therefore there is no reason to restrict the number of input dates in TAMRF, apart for gains in processing time as presented in section 5.1.1.

<figure><img src="image_28.png" /><figcaption>Figure 27: Configurations of the experiment varying the number of input images for TAMRF. The prediction period is highlighted in red: those dates are masked from TAMRF input and the model is asked to reconstruct them, using series of neighboring dates of increasing sizes, from 4 neighboring dates per sensor (2 on each side of the prediction period) to 64 neighboring dates per sensor (32 on each side of the prediction period). Note that only AOIs from the **all** testing set that can provide at least 48 dates per sensor in the 64 neighboring dates per sensor case are used.</figcaption></figure>

# 5. Discussions and Conclusion

## 5.1. Discussion

### 5.1.1. *Model complexity*

Figure 29 gives the processing time and memory consumption to produce an area of 165 × 165 pixels, for both CPU and GPU inference, depending on the number of input
<figure><img src="image_29.png" /><figcaption>Figure 28: RMSE of TAMRF predictions of Sentinel-2 (left) and Landsat-8/9 (right) for masked dates (in red in Fig. 27), depending on the number of input neighboring dates per sensor.</figcaption></figure>

and target dates. As expected, both processing time and memory consumption are mostly linear with respect to the number of input dates. Higher number of target dates only incur marginal overhead, as most of the complexity occurs in the encoder. Since the size of the model is still relatively small, CPU inference, though approximately 35 times slower than GPU inference, is still in a usable range. Because the temporal transformer works at target resolution, memory consumption is high: 20 Gb of memory are required for 128 input images. However, the TAMRF model is embarrassingly parallel in the spatial dimension, and tiled inference can be used to produce larger areas while keeping memory consumption under control. Note that the receptive fields of the spatial encoders should be accounted for in tiling. Memory consumption might be reduced by using optimised transformer formulations in a future work.

### 5.1.2. Generalization to other sensors

TAMRF is a very versatile model that solve many tasks without retraining. It supports variable number of input Sentinel-2 and Landsat-8/9 dates and can predict a variable number of target dates. There is also nothing in the problem formulation or the model architecture that restrict TARMF to only a pair of sensors. However, using different sensors or extending to a greater number of sensors will require adapting the model architecture, by including additional spatial encoders and extending the spatial decoder to predict more spectral bands. This new instance of the model will need to
<figure><img src="image_30.png" /><figcaption>Figure 29: Total processing time and memory consumption to produce an area of 165 × 165 pixels with TAMRF with respect to the total number of input images, for CPU and GPU execution. Measures are averaged over five runs. The ratio between LR and HR dates follows the theoretical ratio between Landsat-8/9 (8 days revisit) and Sentinel-2 (5 days revisit).</figcaption></figure>

be trained for those new sensors. With the advent of foundation models, architectures such as PerceiverIO Jaegle et al. (2021) could alleviate the need for those retraining by encoding the spectral, spatial and temporal context of each band for each sensor, yielding a generic model that could integrate any optical sensor, even if it has not been seen during training. It could also allow generalizing the output spectral bands, removing the need for a fixed set of output bands. Finally, such an architecture would allow introducing conditioning with respect to meteorological data using data sources such as AGERA5, which would enable continuous prediction of LST through time, even for unobserved dates. Meteorological conditioning may also improve the forecasting of snowy dates, and the reflectance of bare soil under different soil moisture conditions. This would require to revisit the architecture of TARMF, but the training process proposed in this paper could still be used to train such a model.

### 5.1.3. Toward new Level 3 products

TAMRF is not trained for a specific task and can be used for various tasks, as demonstrated in section 3.3. One interesting way of using TAMRF is generating joint Landsat-8/9 and Sentinel-2 Level 3 monthly synthesis Hagolle et al. (2021). Those products are generated by weighted averaging of images comprised in each month, taking into account cloud masks. This level of product is very interesting for users, because it simplifies the use of the data by offering a regular, global time grid. How-
ever, because they rely on cloud masks to determine valid pixels, current methods can
yield artifacts such as those shown in Fig. 12. The dates averaged for a given pixel
also depend on the cloud masks, resulting in a different mix of dates for each pixel.
Therefore, current Level 3 products do not correspond to a specific date.

TAMRF can produce regular cloud-free Level 2 predictions with any time-step by
using all observed Landsat-8/9 and Sentinel-2 dates. Fig. 30 shows an example of
such Level 2 products generated with TAMRF by requesting a prediction every 15th
day of the month. Note that all bands, including B8A, B11 and B12 are true 10-meter
resolution, and that there are no cloud mask interpolation artifacts. This could be used
instead of Level 3 products. If one requires real monthly synthesis, the time-step could
be reduced to 10 our 5 days and the predictions could be statistically aggregated. In
both cases, the products for the user would be free of spurious clouds and interpolation
artifacts.

<figure><img src="image_31.png" /><figcaption>Figure 30: Monthly predictions (every 15th of the month) produced by TAMRF from all observed Sentinel-2 and Landsat-8/9 data over the 31TCJ AOIs. The two top rows displays the natural color composition (B4,B3, B2), while the two bottom rows display the SWIR color composition (B8A, B11, B12).</figcaption></figure>
## 5.2. Conclusion

In this paper, we generalise the problem of temporal interpolation, spatial resolution enhancement and spatio-temporal fusion into a single original problem formulation. This problem formulation allows overcoming limitations that have been restricting the capabilities of methods proposed in the literature for decades. We then proposed TAMRF, a DL architecture and training framework that allows to solve the generic problem formulation. This original architecture leverages a combination of Residual CNN for the spatial encoding with a Transformer for temporal encoding. Observations from both sensors are informed by Temporal Positional Encoding and a learnable sensor token, enabling the Transformer to exploit this information in the reconstruction. We demonstrated TAMRF capabilities on a new worldwide, multi-year dataset comprising joint Sentinel-2 and Landsat-8/9 SITS. The TAMRF model has unmatched versatility: it can process any number of Sentinel-2 dates and any number of Landsat-8/9 dates, and it predicts all bands from both sensors at 10-meter spatial resolution, for any target date. All these capabilities are achieved by a single instance of the TAMRF model, trained on a MAE pretext task. We compared TAMRF to existing models in the literature on four tasks: gap-filling, band-sharpening of Sentinel-2 20-meter bands, Spatio-Temporal Fusion, and Thermal sharpening. For all these tasks, TAMRF is on par or better than existing models. TAMRF can also solve those tasks with capabilities that no other method can offer, such as cloud removal when performing Sentinel-2 20-meter bands sharpening, or the ability to input Landsat-8/9 SITS, Sentinel-2 SITS or both and always get consistent predictions without modifying and retraining the model. Finally, we demonstrated the potential of TAMRF in redefining Level 3 processing into a multi-sensor, spatial-resolution enhanced, temporally accurate and artifact free product. An interesting aspect of TAMRF is its modular architecture, which enables swapping the spatial and temporal encoders and decoders with equivalent functional blocks in the future, while retaining the benefits of the global architecture and training strategy. While our work demonstrates good performances of the proposed TAMRF model at a worldwide scale, its performances under different environmental conditions such as biomes or climatic zones remain to be studied. On the methodological side, future work includes generalising TAMRF to any optical sensor and any spatial resolution,
without retraining, which will result in a foundation model for SITS fusion. We also envision to condition TAMRF predictions with exogenous data such as meteorological time-series. The complete source code for training and experiments is available here: https://github.com/Evoland-Land-Monitoring-Evolution/tamrfsits.

## Acknowledgments

* This work was partly performed using HPC resources from GENCI-IDRIS (Grant 2023-AD010114835)

* This work was partly performed using HPC resources from CNES Computing Center.

* The authors acknowledge funding from the EvoLand project (Evolution of the Copernicus Land Service portfolio), grant agreement No 101082130) funded from the European Union's Horizon Europe research and innovation program.

## CRediT

**Julien Michel:** Conceptualization, Formal Analysis, Methodology, Data Curation, Investigation, Software, Writing - Original Draft. **Jordi Inglada:** Conceptualization, Supervision, Writing - Review & Editing.

## References

S. Bojinski, M. Verstraete, T. C. Peterson, C. Richter, A. Simmons, M. Zemp, The Concept of Essential Climate Variables in Support of Climate Research, Applications, and Policy, Bulletin of the American Meteorological Society 95 (2014) 1431 – 1443. URL: https://journals.ametsoc.org/view/journals/bams/95/9/bams-d-13-00047.1.xml. doi:10.1175/BAMS-D-13-00047.1.

W. Jetz, M. A. McGeoch, R. Guralnick, S. Ferrier, J. Beck, M. J. Costello, M. Fernandez, G. N. Geller, P. Keil, C. Merow, C. Meyer, F. E. Muller-Karger, H. M. Pereira, E. C. Regan, D. S. Schmeller, E. Turak, Essential biodiversity variables
for mapping and monitoring species populations, Nature Ecology & Evolution 3 (2019) 539–551. URL: http://dx.doi.org/10.1038/s41559-019-0826-1. doi:10.1038/s41559-019-0826-1.

M. A. Wulder, D. P. Roy, V. C. Radeloff, T. R. Loveland, M. C. Anderson, D. M. Johnson, S. Healey, Z. Zhu, T. A. Scambos, N. Pahlevan, M. Hansen, N. Gorelick, C. J. Crawford, J. G. Masek, T. Hermosilla, J. C. White, A. S. Belward, C. Schaaf, C. E. Woodcock, J. L. Huntington, L. Lymburner, P. Hostert, F. Gao, A. Lyapustin, J.-F. Pekel, P. Strobl, B. D. Cook, Fifty years of landsat science and impacts, Remote Sensing of Environment 280 (2022) 113195. URL: http://dx.doi.org/10.1016/j.rse.2022.113195. doi:10.1016/j.rse.2022.113195.

G. Misra, F. Cawkwell, A. Wingler, Status of Phenological Research Using Sentinel-2 Data: a Review, Remote Sensing 12 (2020) 2760. URL: http://dx.doi.org/10.3390/rs12172760. doi:10.3390/rs12172760.

D. Phiri, M. Simwanda, S. Salekin, V. Nyirenda, Y. Murayama, M. Ranagalage, Sentinel-2 Data for Land Cover/use Mapping: a Review, Remote Sensing 12 (2020) 2291. URL: http://dx.doi.org/10.3390/rs12142291. doi:10.3390/rs12142291.

B. Vajsová, D. Fasbender, C. Wirnhardt, S. Lemajic, W. Devos, Assessing spatial limits of Sentinel-2 data on arable crops in the context of checks by monitoring, Remote Sensing 12 (2020) 2195. URL: http://dx.doi.org/10.3390/rs12142195. doi:10.3390/rs12142195.

U. Bhangale, S. More, T. Shaikh, S. Patil, N. More, Analysis of surface water resources using Sentinel-2 imagery, Procedia Computer Science 171 (2020) 2645–2654. URL: http://dx.doi.org/10.1016/j.procs.2020.04.287. doi:10.1016/j.procs.2020.04.287.

A. M. Wilson, W. Jetz, Remotely sensed high-resolution global cloud dynamics for predicting ecosystem and biodiversity distributions, PLOS Biology 14 (2016) e1002415. URL: http://dx.doi.org/10.1371/journal.pbio.1002415. doi:10.1371/journal.pbio.1002415.
J. Li, B. Chen, Global revisit interval analysis of landsat-8-9 and Sentinel-2a-2b data for terrestrial monitoring, Sensors 20 (2020) 6631. URL: http://dx.doi.org/10.3390/s20226631. doi:10.3390/s20226631.

N. Latte, P. Lejeune, Planetscope Radiometric Normalization and Sentinel-2 Super-Resolution (2.5 m): a Straightforward Spectral-Spatial Fusion of Multi-Satellite Multi-Sensor Images Using Residual Convolutional Neural Networks, Remote Sensing 12 (2020) 2366. URL: http://dx.doi.org/10.3390/rs12152366. doi:10.3390/rs12152366.

Y. Sadeh, X. Zhu, D. Dunkerley, J. P. Walker, Y. Zhang, O. Rozenstein, V. Manivasagam, K. Chenu, Fusion of Sentinel-2 and Planetscope Time-Series Data Into Daily 3 M Surface Reflectance and Wheat Lai Monitoring, International Journal of Applied Earth Observation and Geoinformation 96 (2021) 102260. URL: http://dx.doi.org/10.1016/j.jag.2020.102260. doi:10.1016/j.jag.2020.102260.

J.-P. Lagouarde, B. Bhattacharya, P. Crébassol, P. Gamet, S. S. Babu, G. Boulet, X. Briottet, K. Buddhiraju, S. Cherchali, I. Dadou, G. Dedieu, M. Gouhier, O. Hagolle, M. Irvine, F. Jacob, A. Kumar, K. K. Kumar, B. Laignel, K. Mallick, C. Murthy, A. Olioso, C. Ottlé, M. R. Pandya, P. V. Raju, J.-L. Roujean, M. Sekhar, M. V. Shukla, S. K. Singh, J. Sobrino, R. Ramakrishnan, The Indian-French Trishna Mission: Earth Observation in the Thermal Infrared with High Spatio-Temporal Resolution, in: IGARSS 2018 - 2018 IEEE International Geoscience and Remote Sensing Symposium, 2018, pp. 4078-4081. URL: http://dx.doi.org/10.1109/IGARSS.2018.8518720. doi:10.1109/igarss.2018.8518720.

F. Bernard, I. Manolis, I. Barat, A. B. Alamañac, M. S. Taboada, P. Mingorance, A. Ciapponi, T. Cardone, I. F. Nunez, A. Garcia, P. Hallibert, A. Hammar, P. Henriot, D. M. Codinachs, S. Patti, P. P. Montes, P. Skrzypek, D. Steenari, S. Weixler, S. Deslous, C. Coatantiec, O. A. Trotta, I. C. Vega, D. G. Holgueras, The Copernicus land surface temperature monitoring (LSTM) mission: design, technology and sta-
tus, in: Sensors, Systems, and Next-Generation Satellites XXVII, 2023, p. 8. URL: http://dx.doi.org/10.1117/12.2679705. doi:10.1117/12.2679705.

Q. Wang, Y. Tang, Y. Ge, H. Xie, X. Tong, P. M. Atkinson, A comprehensive review of spatial-temporal-spectral information reconstruction techniques, Science of Remote Sensing 8 (2023) 100102. URL: http://dx.doi.org/10.1016/j.srs.2023.100102. doi:10.1016/j.srs.2023.100102.

S. Anwar, S. Khan, N. Barnes, A Deep Journey Into Super-Resolution, ACM Computing Surveys 53 (2020) 1–34. URL: http://dx.doi.org/10.1145/3390462. doi:10.1145/3390462.

A. Liu, Y. Liu, J. Gu, Y. Qiao, C. Dong, Blind image super-resolution: a survey and beyond, IEEE Transactions on Pattern Analysis and Machine Intelligence (2022) 1–19. URL: http://dx.doi.org/10.1109/TPAMI.2022.3203009. doi:10.1109/tpami.2022.3203009.

X. Wang, K. Yu, S. Wu, J. Gu, Y. Liu, C. Dong, Y. Qiao, C. Change Loy, ES-RGAN: Enhanced Super-Resolution Generative Adversarial Networks, in: Proceedings of the European Conference on Computer Vision (ECCV) Workshops, 2018. URL: https://openaccess.thecvf.com/content_eccv_2018_workshops/w25/html/Wang_ESRGAN_Enhanced_Super-Resolution_Generative_Adversarial_Networks_ECCVW_2018_paper.html.

C. Ledig, L. Theis, F. Huszar, J. Caballero, A. Cunningham, A. Acosta, A. Aitken, A. Tejani, J. Totz, Z. Wang, W. Shi, Photo-Realistic Single Image Super-Resolution Using a Generative Adversarial Network, in: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2017. URL: https://openaccess.thecvf.com/content_cvpr_2017/html/LedigPhoto-Realistic_Single_Image_CVPR_2017_paper.html.

D. Pouliot, R. Latifovic, J. Pasher, J. Duffe, Landsat Super-Resolution Enhancement Using Convolution Neural Networks and Sentinel-2 for Training, Remote Sensing 10 (2018) 394. URL: http://dx.doi.org/10.3390/rs10030394. doi:10.3390/rs10030394.
M. Galar, R. Sesma, C. Ayala, C. Aranda, Super-resolution for sentinel-2 images, The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences XLII-2/W16 (2019) 95-102. URL: http://dx.doi.org/10.5194/isprs-archives-XLII-2-W16-95-2019. doi:10.5194/isprs-archives-xlii-2-w16-95-2019.

L. S. Romero, J. Marcello, V. Vilaplana, Super-resolution of sentinel-2 imagery using generative adversarial networks, Remote Sensing 12 (2020) 2424. URL: http://dx.doi.org/10.3390/rs12152424. doi:10.3390/rs12152424.

N. L. Nguyen, J. Anger, L. Raad, B. Galerne, G. Facciolo, On The Role of Alias and Band-Shift for Sentinel-2 Super-Resolution, in: IGARSS 2023 - 2023 IEEE International Geoscience and Remote Sensing Symposium, 2023, pp. 4294-4297. URL: https://ieeexplore.ieee.org/document/10282805. doi:10.1109/IGARSS52108.2023.10282805.

M. Märtens, D. Izzo, A. Krzic, D. Cox, Super-Resolution of Proba-V Images Using Convolutional Neural Networks, Astrodynamics 3 (2019) 387-402. URL: http://dx.doi.org/10.1007/s42064-019-0059-8. doi:10.1007/s42064-019-0059-8.

A. B. Molini, D. Valsesia, G. Fracastoro, E. Magli, Deepsum: Deep Neural Network for Super-Resolution of Unregistered Multitemporal Images, IEEE Transactions on Geoscience and Remote Sensing 58 (2020) 3644-3656. URL: http://dx.doi.org/10.1109/TGRS.2019.2959248. doi:10.1109/tgrs.2019.2959248.

A. Okabayashi, N. Audebert, S. Donike, C. Pelletier, Cross-sensor super-resolution of irregularly sampled Sentinel-2 time series, in: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR) Workshops, 2024, pp. 502-511. URL: https://openaccess.thecvf.com/content/CVPR2024W/EarthVision/html/Okabayashi_Cross-sensorSuper-resolution_of_ irregularly_sampled_Sentinel-2_time_series_CVPRW_2024_paper.html.
M. R. Ibrahim, R. Benavente, D. Ponsa, F. Lumbreras, Hyda-Net: a Hybrid Dense Attention Network for Remote Sensing Multi-Image Super-Resolution, IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing (2025) 1–24. URL: http://dx.doi.org/10.1109/JSTARS.2025.3547785. doi:10.1109/jstars.2025.3547785.

J. Michel, E. Kalinicheva, J. Inglada, Revisiting Remote Sensing Cross-Sensor Single Image Super-Resolution: the Overlooked Impact of Geometric and Radiometric Distortion, IEEE Transactions on Geoscience and Remote Sensing 63 (2025) 1–22. URL: http://dx.doi.org/10.1109/TGRS.2025.3572548. doi:10.1109/tgrs.2025.3572548.

V. T. Sambandham, K. Kirchheim, F. Ortmeier, S. Mukhopadhaya, Deep Learning-Based Harmonization and Super-Resolution of Landsat-8 and Sentinel-2 Images, ISPRS Journal of Photogrammetry and Remote Sensing 212 (2024) 274–288. URL: http://dx.doi.org/10.1016/j.isprsjprs.2024.04.026. doi:10.1016/j.isprsjprs.2024.04.026.

M. K. Firozjaei, M. Kiavarz, S. K. Alavipanah, Satellite-Derived Land Surface Temperature Spatial Sharpening: a Comprehensive Review on Current Status and Perspectives, European Journal of Remote Sensing 55 (2022) 644–664. URL: http://dx.doi.org/10.1080/22797254.2022.2144764. doi:10.1080/22797254.2022.2144764.

G. Vivone, L. Alparone, J. Chanussot, M. D. Mura, A. Garzelli, G. A. Licciardi, R. Restaino, L. Wald, A Critical Comparison Among Pansharpening Algorithms, IEEE Transactions on Geoscience and Remote Sensing 53 (2015) 2565–2586. URL: http://dx.doi.org/10.1109/TGRS.2014.2361734. doi:10.1109/tgrs.2014.2361734.

M. Ciotola, G. Guarino, G. Vivone, G. Poggi, J. Chanussot, A. Plaza, G. Scarpa, Hyperspectral Pansharpening: Critical Review, Tools, and Future Perspectives, IEEE Geoscience and Remote Sensing Magazine (2025) 2–29. URL: http://dx.doi.org/10.1109/MGRS.2024.3509139. doi:10.1109/mgrs.2024.3509139.
R. d. Lutio, S. D'aronco, J. D. Wegner, K. Schindler, Guided super-resolution as pixel-to-pixel transformation, in: Proceedings of the IEEE/CVF International Conference on Computer Vision, 2019, pp. 8829–8837. URL: https://openaccess.thecvf.com/content_ICCV_2019/html/de_Lutio_Guided_Super-Resolution_As_Pixel-to-Pixel_Transformation_ICCV_2019_paper.html.

H. Su, Y. Li, Y. Xu, X. Fu, S. Liu, A Review of Deep-Learning-Based Super-Resolution: From Methods To Applications, Pattern Recognition 157 (2025) 110935. URL: http://dx.doi.org/10.1016/j.patcog.2024.110935. doi:10.1016/j.patcog.2024.110935.

L. Wald, T. Ranchin, M. Mangolini, Fusion of satellite images of different spatial resolutions: Assessing the quality of resulting images, Photogrammetric engineering and remote sensing 63 (1997) 691–699. URL: https://hal.science/hal-00365304/.

F. Palsson, J. R. Sveinsson, M. O. Ulfarsson, Sentinel-2 image fusion using a deep residual network, Remote Sensing 10 (2018). URL: https://www.mdpi.com/2072-4292/10/8/1290. doi:10.3390/rs10081290.

F. Gao, W. P. Kustas, M. C. Anderson, A data mining approach for sharpening thermal satellite imagery over land, Remote Sensing 4 (2012) 3287–3319. URL: https://www.mdpi.com/2072-4292/4/11/3287. doi:10.3390/rs4113287.

C. Granero-Belinchon, A. Michel, J.-P. Lagouarde, J. A. Sobrino, X. Briottet, Multi-resolution study of thermal unmixing techniques over madrid urban area: Case study of trishna mission, Remote Sensing 11 (2019). URL: https://www.mdpi.com/2072-4292/11/10/1251. doi:10.3390/rs11101251.

L. Salgueiro, J. Marcello, V. Vilaplana, Single-Image Super-Resolution of Sentinel-2 Low Resolution Bands With Residual Dense Convolutional Neural Networks, Remote Sensing 13 (2021) 5007. URL: http://dx.doi.org/10.3390/rs13245007. doi:10.3390/rs13245007.
O. Merlin, B. Duchemin, O. Hagolle, F. Jacob, B. Coudert, G. Chehbouni, G. Dedieu, J. Garatuza, Y. Kerr, Disaggregation of Modis Surface Temperature Over an Agricultural Area Using a Time Series of Formosat-2 Images, Remote Sensing of Environment 114 (2010) 2500–2512. URL: http://dx.doi.org/10.1016/j.rse.2010.05.025. doi:10.1016/j.rse.2010.05.025.

B. M. Nguyen, G. Tian, M.-T. Vo, A. Michel, T. Corpetti, C. Granero-Belinchon, Convolutional Neural Network Modelling for MODIS Land Surface Temperature Super-Resolution, in: 2022 30th European Signal Processing Conference (EUSIPCO), 2022, pp. 1806–1810. URL: http://dx.doi.org/10.23919/EUSIPCO55093.2022.9909569. doi:10.23919/eusipco55093.2022.9909569.

Z. Shao, J. Cai, P. Fu, L. Hu, T. Liu, Deep Learning-Based Fusion of Landsat-8 and Sentinel-2 Images for a Harmonized Surface Reflectance Product, Remote Sensing of Environment 235 (2019) 111425. URL: http://dx.doi.org/10.1016/j.rse.2019.111425. doi:10.1016/j.rse.2019.111425.

M. Belgiu, A. Stein, Spatiotemporal Image Fusion in Remote Sensing, Remote Sensing 11 (2019) 818. URL: http://dx.doi.org/10.3390/rs11070818. doi:10.3390/rs11070818.

J. Xiao, A. K. Aggarwal, N. H. Duc, A. Arya, U. K. Rage, R. Avtar, A Review of Remote Sensing Image Spatiotemporal Fusion: Challenges, Applications and Recent Trends, Remote Sensing Applications: Society and Environment 32 (2023) 101005. URL: http://dx.doi.org/10.1016/j.rsase.2023.101005. doi:10.1016/j.rsase.2023.101005.

F. Gao, J. Masek, M. Schwaller, F. Hall, On the Blending of the Landsat and Modis Surface Reflectance: Predicting Daily Landsat Surface Reflectance, IEEE Transactions on Geoscience and Remote Sensing 44 (2006) 2207–2218. URL: http://dx.doi.org/10.1109/TGRS.2006.872081. doi:10.1109/tgrs.2006.872081.

X. Zhu, F. Cai, J. Tian, T. Williams, Spatiotemporal Fusion of Multisource Remote Sensing Data: Literature Survey, Taxonomy, Principles, Applications, and Future
Directions, Remote Sensing 10 (2018) 527. URL: http://dx.doi.org/10.3390/rs10040527. doi:10.3390/rs10040527.

Z. Tan, P. Yue, L. Di, J. Tang, Deriving High Spatiotemporal Remote Sensing Images Using Deep Convolutional Network, Remote Sensing 10 (2018) 1066. URL: http://dx.doi.org/10.3390/rs10071066. doi:10.3390/rs10071066.

Z. Tan, L. Di, M. Zhang, L. Guo, M. Gao, An enhanced deep convolutional model for spatiotemporal image fusion, Remote Sensing 11 (2019) 2898. URL: http://dx.doi.org/10.3390/rs11242898. doi:10.3390/rs11242898.

X. Liu, C. Deng, J. Chanussot, D. Hong, B. Zhao, StfNet: a two-stream convolutional neural network for spatiotemporal image fusion, IEEE Transactions on Geoscience and Remote Sensing 57 (2019) 6552–6564. URL: http://dx.doi.org/10.1109/TGRS.2019.2907310. doi:10.1109/tgrs.2019.2907310.

Z. Tan, M. Gao, X. Li, L. Jiang, A Flexible Reference-Insensitive Spatiotemporal Fusion Model for Remote Sensing Images Using Conditional Generative Adversarial Network, IEEE Transactions on Geoscience and Remote Sensing 60 (2022) 1–13. URL: http://dx.doi.org/10.1109/TGRS.2021.3050551. doi:10.1109/tgrs.2021.3050551.

Y. Zhang, R. Fan, P. Duan, J. Dong, Z. Lei, Dcdgan-Stf: a Multiscale Deformable Convolution Distillation Gan for Remote Sensing Image Spatiotemporal Fusion (2024b) 1–15. URL: http://dx.doi.org/10.1109/JSTARS.2024.3476153. doi:10.1109/jstars.2024.3476153.

Y. Xie, J. Hu, K. He, L. Cao, K. Yang, L. Chen, The Gan Spatiotemporal Fusion Model Based on Multi-Scale Convolution and Attention Mechanism for Remote Sensing Images, IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing (2024) 1–13. URL: http://dx.doi.org/10.1109/JSTARS.2024.3507815. doi:10.1109/jstars.2024.3507815.

G. Yang, Y. Qian, H. Liu, B. Tang, R. Qi, Y. Lu, J. Geng, Msfusion: Multistage for Remote Sensing Image Spatiotemporal Fusion Based on Texture Transformer and
Convolutional Neural Network, IEEE Journal of Selected Topics in Applied Earth
Observations and Remote Sensing 15 (2022) 4653-4666. URL: http://dx.doi.
org/10.1109/JSTARS.2022.3179415. doi:10.1109/jstars.2022.3179415.

W. Li, D. Cao, M. Xiang, Enhanced multi-stream remote sensing spatiotemporal fusion network based on transformer and dilated convolution, Remote Sensing 14 (2022b) 4544. URL: http://dx.doi.org/10.3390/rs14184544. doi:10.3390/rs14184544.

Y. Wu, M. Huang, A unified generative adversarial network with convolution and transformer for remote sensing image fusion, IEEE Transactions on Geoscience and Remote Sensing (2024) 1-1. URL: http://dx.doi.org/10.1109/TGRS.2024.3441719. doi:10.1109/tgrs.2024.3441719.

G. Chen, P. Jiao, Q. Hu, L. Xiao, Z. Ye, Swinstfm: Remote Sensing Spatiotemporal Fusion Using Swin Transformer, IEEE Transactions on Geoscience and Remote Sensing 60 (2022) 1-18. URL: http://dx.doi.org/10.1109/TGRS.2022.3182809. doi:10.1109/tgrs.2022.3182809.

H. Goyena, U. Pérez-Goya, M. Montesino-SanMartin, A. F. Militino, Q. Wang, P. M. Atkinson, M. D. Ugarte, Unpaired Spatio-Temporal Fusion of Image Patches (USTFIP) From Cloud Covered Images, Remote Sensing of Environment 295 (2023) 113709. URL: http://dx.doi.org/10.1016/j.rse.2023.113709. doi:10.1016/j.rse.2023.113709.

X. Zhang, S. Li, Z. Tan, X. Li, Enhanced wavelet based spatiotemporal fusion networks using cross-paired remote sensing images, ISPRS Journal of Photogrammetry and Remote Sensing 211 (2024c) 281-297. URL: http://dx.doi.org/10.1016/j.isprsjprs.2024.04.016. doi:10.1016/j.isprsjprs.2024.04.016.

H. Shen, X. Meng, L. Zhang, An integrated framework for the spatio-temporal-spectral fusion of remote sensing images, IEEE Transactions on Geoscience and Remote Sensing 54 (2016) 7135-7148. URL: http://dx.doi.org/10.1109/TGRS.2016.2596290. doi:10.1109/tgrs.2016.2596290.
W. R. Moskolaï, W. Abdou, A. Dipanda, Kolyang, Application of deep learning architectures for satellite image time series prediction: a review, Remote Sensing 13 (2021) 4822. URL: http://dx.doi.org/10.3390/rs13234822. doi:10.3390/rs13234822.

S. Li, L. Xu, Y. Jing, H. Yin, X. Li, X. Guan, High-Quality Vegetation Index Product Generation: a Review of Ndvi Time Series Reconstruction Techniques, International Journal of Applied Earth Observation and Geoinformation 105 (2021) 102640. URL: http://dx.doi.org/10.1016/j.jag.2021.102640. doi:10.1016/j.jag.2021.102640.

D. P. Roy, J. Ju, P. Lewis, C. Schaaf, F. Gao, M. Hansen, E. Lindquist, Multi-Temporal Modis-Landsat Data Fusion for Relative Radiometric Normalization, Gap Filling, and Prediction of Landsat Data, Remote Sensing of Environment 112 (2008) 3112-3130. URL: http://dx.doi.org/10.1016/j.rse.2008.03.009. doi:10.1016/j.rse.2008.03.009.

C. Stucker, V. S. F. Garnot, K. Schindler, U-tilise: a sequence-to-sequence model for cloud removal in optical satellite time series, IEEE Transactions on Geoscience and Remote Sensing 61 (2023) 1-16. URL: http://dx.doi.org/10.1109/TGRS.2023.3333391. doi:10.1109/tgrs.2023.3333391.

H. Liu, H. K. Zhang, B. Huang, L. Yan, K. K. Tran, Y. Qiu, X. Zhang, D. P. Roy, Reconstruction of Seamless Harmonized Landsat Sentinel-2 (HLS) Time Series Via Self-Supervised Learning, Remote Sensing of Environment 308 (2024) 114191. URL: http://dx.doi.org/10.1016/j.rse.2024.114191. doi:10.1016/j.rse.2024.114191.

B. Irigireddy, V. Bandaru, SatFlow: Generative model based framework for producing High Resolution Gap Free Remote Sensing Imagery, 2025. URL: https://arxiv.org/abs/2502.01098. arXiv:2502.01098.

Q. Shu, X. Zhu, S. Xu, Y. Wang, D. Liu, Restore-dit: Reliable satellite image time series reconstruction by multimodal sequential diffusion transformer, Remote Sens-
ing of Environment 328 (2025) 114872. URL: http://dx.doi.org/10.1016/j.rse.2025.114872. doi:10.1016/j.rse.2025.114872.

M. Claverie, J. Ju, J. G. Masek, J. L. Dungan, E. F. Vermote, J.-C. Roger, S. V. Skakun, C. Justice, The Harmonized Landsat and Sentinel-2 Surface Reflectance Data Set, Remote Sensing of Environment 219 (2018) 145–161. URL: http://dx.doi.org/10.1016/j.rse.2018.09.002. doi:10.1016/j.rse.2018.09.002.

S. Saunier, B. Pflug, I. Lobos, B. Franch, J. Louis, R. D. L. Reyes, V. Debaecker, E. Cadau, V. Boccia, F. Gascon, S. Kocaman, Sen2like: Paving the Way Towards Harmonization and Fusion of Optical Data, Remote Sensing 14 (2022) 3855. URL: http://dx.doi.org/10.3390/rs14163855. doi:10.3390/rs14163855.

Q. Wang, G. A. Blackburn, A. O. Onojeghuo, J. Dash, L. Zhou, Y. Zhang, P. M. Atkinson, Fusion of Landsat 8 OLI and Sentinel-2 Msi Data, IEEE Transactions on Geoscience and Remote Sensing 55 (2017) 3885–3899. URL: http://dx.doi.org/10.1109/TGRS.2017.2683444. doi:10.1109/tgrs.2017.2683444.

C. Dong, C. C. Loy, K. He, X. Tang, Image Super-Resolution Using Deep Convolutional Networks, IEEE Transactions on Pattern Analysis and Machine Intelligence 38 (2016) 295–307. URL: http://dx.doi.org/10.1109/TPAMI.2015.2439281. doi:10.1109/tpami.2015.2439281.

Y. Luo, K. Guan, J. Peng, Stair: a Generic and Fully-Automated Method To Fuse Multiple Sources of Optical Satellite Data To Generate a High-Resolution, Daily and Cloud-/gap-Free Surface Reflectance Product, Remote Sensing of Environment 214 (2018) 87–99. URL: http://dx.doi.org/10.1016/j.rse.2018.04.042. doi:10.1016/j.rse.2018.04.042.

Y. Luo, K. Guan, J. Peng, S. Wang, Y. Huang, Stair 2.0: a Generic and Automatic Algorithm To Fuse Modis, Landsat, and Sentinel-2 To Generate 10 M, Daily, and Cloud-/gap-Free Surface Reflectance Product, Remote Sensing 12 (2020) 3209. URL: http://dx.doi.org/10.3390/rs12193209. doi:10.3390/rs12193209.
P. Wang, M. Huang, S. Shi, B. Huang, B. Zhou, G. Xu, L. Wang, H. Leung, Landsat-8 and Sentinel-2 image fusion based on multi-scale smoothing-sharpening filter, IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing (2024) 1–14. URL: http://dx.doi.org/10.1109/JSTARS.2024.3469974. doi:10.1109/jstars.2024.3469974.

R. Mukherjee, D. Liu, Spatial and Spectral Translation of Landsat 8 To Sentinel-2 Using Conditional Generative Adversarial Networks, Remote Sensing 15 (2023) 5502. URL: http://dx.doi.org/10.3390/rs15235502. doi:10.3390/rs15235502.

B. Chen, J. Li, Y. Jin, Deep learning for feature-level data fusion: Higher resolution reconstruction of historical landsat archive, Remote Sensing 13 (2021) 167. URL: http://dx.doi.org/10.3390/rs13020167. doi:10.3390/rs13020167.

J. Chang, W. Du, B. Zhang, S. Guo, Y. Yin, Z. Wang, T. Xu, Z. Feng, Based on the Improved Edcstfn Model, Modis, Landsat 8 and Sentinel-2 Data Were Fused To Obtain 10m Dense Time Series Images, IEEE Access (2025) 1–1. URL: http://dx.doi.org/10.1109/ACCESS.2025.3564968. doi:10.1109/access.2025.3564968.

Y. Yuan, L. Lin, Q. Liu, R. Hang, Z.-G. Zhou, Sits-Former: a Pre-Trained Spatio-Spectral-Temporal Representation Model for Sentinel-2 Time Series Classification, International Journal of Applied Earth Observation and Geoinformation 106 (2022) 102651. URL: http://dx.doi.org/10.1016/j.jag.2021.102651. doi:10.1016/j.jag.2021.102651.

I. Dumeur, S. Valero, J. Inglada, Self-supervised spatio-temporal representation learning of satellite image time series, IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing 17 (2024) 4350–4367. URL: http://dx.doi.org/10.1109/JSTARS.2024.3358066. doi:10.1109/jstars.2024.3358066.

A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, L. u. Kaiser, I. Polosukhin, Attention is all you need, in: I. Guyon, U. V. Luxburg, S. Ben-gio, H. Wallach, R. Fergus, S. Vishwanathan, R. Garnett (Eds.), Advances in Neural Information Processing Systems, volume 30, Curran Associates, Inc., 2017.
URL: https://proceedings.neurips.cc/paper_files/paper/2017/file/3f5ee243547dee91fbd053c1c4a845aa-Paper.pdf.

X. Guo, J. Lao, B. Dang, Y. Zhang, L. Yu, L. Ru, L. Zhong, Z. Huang, K. Wu, D. Hu, H. He, J. Wang, J. Chen, M. Yang, Y. Zhang, Y. Li, SkySense: A Multi-Modal Remote Sensing Foundation Model Towards Universal Interpretation for Earth Observation Imagery, in: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), 2024, pp. 27672–27683.
URL: https://openaccess.thecvf.com/content/CVPR2024/html/Guo_SkySense_A_Multi-Modal_Remote_Sensing_Foundation_Model_Towards_Universal_Interpretation_CVPR_2024_paper.html.

I. Dumeur, S. Valero, J. Inglada, Paving the way toward foundation models for irregular and unaligned satellite image time series, IEEE Transactions on Geoscience and Remote Sensing 63 (2025) 1–21. URL: http://dx.doi.org/10.1109/TGRS.2025.3589013. doi:10.1109/tgrs.2025.3589013.

K. He, X. Chen, S. Xie, Y. Li, P. Dollár, R. Girshick, Masked Autoencoders Are Scalable Vision Learners, in: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), 2022, pp. 16000–16009. URL: https://openaccess.thecvf.com/content/CVPR2022/html/He_Masked_Autoencoders_Are_Scalable_Vision_Learners_CVPR_2022_paper.

C. J. Reed, R. Gupta, S. Li, S. Brockman, C. Funk, B. Clipp, K. Keutzer, S. Candido, M. Uyttendaele, T. Darrell, Scale-mae: A scale-aware masked autoencoder for multiscale geospatial representation learning, in: Proceedings of the IEEE/CVF International Conference on Computer Vision (ICCV), 2023, pp. 4088–4099.
URL: https://openaccess.thecvf.com/content/ICCV2023/html/Reed_Scale-MAE_A_Scale-Aware_Masked_Autoencoder_for_Multiscale_Geospatial_Representation_Learning_ICCV_2023_paper.html.

H. Liu, J. Gan, X. Fan, Y. Zhang, C. Luo, J. Zhang, G. Jiang, Y. Qian, C. Zhao, H. Ma, Z. Guo, PT-Tuning: Bridging the Gap between Time Series Masked Reconstruction and Forecasting via Prompt Token Tuning, Lecture Notes in Computer Science,
Springer Nature Singapore, 2025, pp. 147–162. URL: http://dx.doi.org/10.1007/978-981-97-5779-4_10. doi:10.1007/978-981-97-5779-4_10.

R. Girshick, Fast R-CNN, in: Proceedings of the IEEE International Conference on Computer Vision (ICCV), 2015. URL: https://openaccess.thecvf.com/content_iccv_2015/html/Girshick_Fast_R-CNN_ICCV_2015_paper.html.

V. Balntas, E. Riba, D. Ponsa, K. Mikolajczyk, Learning local feature descriptors with triplets and shallow convolutional neural networks, in: Proceedings of the British Machine Vision Conference 2016, 2016, pp. 119.1–119.11. URL: http://dx.doi.org/10.5244/C.30.119. doi:10.5244/c.30.119.

R. Zhang, P. Isola, A. A. Efros, E. Shechtman, O. Wang, The Unreasonable Effectiveness of Deep Features as a Perceptual Metric, in: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2018. URL: https://openaccess.thecvf.com/content_cvpr_2018/html/Zhang_The_Unreasonable_Effectiveness_CVPR_2018_paper.html.

E. L. C. CCI, Product user guide version 2.0, UCL-Geomatics: London, UK 685 (2017). URL: http://maps.elie.ucl.ac.be/CCI/viewer/download/ESACCI-LC-Ph2-PUGv2_2.0.pdf.

F. Rubel, M. Kottek, Observed and projected climate shifts 1901-2100 depicted by world maps of the köppen-geiger climate classification, Meteorologische Zeitschrift 19 (2010) 135–141. URL: http://dx.doi.org/10.1127/0941-2948/2010/0430. doi:10.1127/0941-2948/2010/0430.

M. Schramm, E. Pebesma, M. Milenković, L. Foresta, J. Dries, A. Jacob, W. Wagner, M. Mohr, M. Neteler, M. Kadunc, T. Miksa, P. Kempeneers, J. Verbesselt, B. Gößwein, C. Navacchi, S. Lippens, J. Reiche, The Openeo Api-Harmonising the Use of Earth Observation Cloud Services Using Virtual Data Cube Function alities, Remote Sensing 13 (2021) 1125. URL: http://dx.doi.org/10.3390/rs13061125. doi:10.3390/rs13061125.
D. P. Kingma, J. Ba, Adam: a Method for Stochastic Optimization (2014). URL: https://arxiv.org/abs/1412.6980. doi:10.48550/ARXIV.1412.6980.

I. Loshchilov, F. Hutter, SGDR: Stochastic Gradient Descent With Warm Restarts (2016). URL: https://arxiv.org/abs/1608.03983. doi:10.48550/ARXIV.1608.03983.

A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Kopf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, S. Chintala, Pytorch: An imperative style, high-performance deep learning library, in: H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, R. Garnett (Eds.), Advances in Neural Information Processing Systems 32, Curran Associates, Inc., 2019, pp. 8024–8035. URL: http://papers.neurips.cc/paper/9015-pytorch-an-imperative-style-high-performance-deep-learning-library.pdf.

A. Mittal, A. K. Moorthy, A. C. Bovik, No-Reference Image Quality Assessment in the Spatial Domain, IEEE Transactions on Image Processing 21 (2012) 4695–4708. URL: http://dx.doi.org/10.1109/TIP.2012.2214050. doi:10.1109/tip.2012.2214050.

R. Ait-Bachir, C. Granero-Belinchon, A. Michel, J. Michel, X. Briottet, L. Drumetz, Land surface temperature super-resolution with a scale-invariance-free neural approach: Application to modis, IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing 18 (2025) 14480–14494. URL: http://dx.doi.org/10.1109/JSTARS.2025.3573610. doi:10.1109/jstars.2025.3573610.

H. Huang, W. He, H. Zhang, Y. Xia, L. Zhang, Stfdiff: Remote Sensing Image Spatiotemporal Fusion With Diffusion Models, Information Fusion 111 (2024) 102505. URL: http://dx.doi.org/10.1016/j.inffus.2024.102505. doi:10.1016/j.inffus.2024.102505.

D. Guo, Z. Li, X. Gao, M. Gao, C. Yu, C. Zhang, W. Shi, Realfusion:
a Reliable Deep Learning-Based Spatiotemporal Fusion Framework for Gen-
erating Seamless Fine-Resolution Imagery, Remote Sensing of Environment
321 (2025) 114689. URL: http://dx.doi.org/10.1016/j.rse.2025.114689.
doi:10.1016/j.rse.2025.114689.

C. Requena-Mesa, V. Benson, M. Reichstein, J. Runge, J. Denzler, Earthnet2021: A large-scale dataset and challenge for earth surface forecasting as a guided video prediction task., in: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR) Workshops, 2021, pp. 1132–1142. URL: https://openaccess.thecvf.com/content/CVPR2021W/EarthVision/html/Requena-Mesa_EarthNet2021_A_Large-Scale_Dataset_and_Challenge_for_Earth

J. Colin, O. Hagolle, L. Landier, S. Coustance, P. Kettig, A. Meygret, J. Osman, E. Ver-
mote, Assessment of the Performance of the Atmospheric Correction Algorithm
Maja for Sentinel-2 Surface Reflectance Estimates, Remote Sensing 15 (2023) 2665.
URL: http://dx.doi.org/10.3390/rs15102665. doi:10.3390/rs15102665.

C. Lanaras, J. Bioucas-Dias, S. Galliani, E. Baltsavias, K. Schindler, Super-
Resolution of Sentinel-2 Images: Learning a Globally Applicable Deep Neu-
ral Network, ISPRS Journal of Photogrammetry and Remote Sensing 146
(2018) 305–319. URL: http://dx.doi.org/10.1016/j.isprsjprs.2018.09.
018. doi:10.1016/j.isprsjprs.2018.09.018.

J. Wu, L. Lin, C. Zhang, T. Li, X. Cheng, F. Nan, Generating sentinel-2 all-band 10-m data by sharpening 20/60-m bands: a hierarchical fusion network, ISPRS Journal of Photogrammetry and Remote Sensing 196 (2023) 16–31. URL: http://dx.doi.org/10.1016/j.isprsjprs.2022.12.017. doi:10.1016/j.isprsjprs.2022.12.017.

M. Ciotola, G. Guarino, A. Mazza, G. Poggi, G. Scarpa, A comprehensive bench-
marking framework for sentinel-2 sharpening: Methods, dataset, and evaluation
metrics, Remote Sensing 17 (2025) 1983. URL: http://dx.doi.org/10.3390/
rs17121983. doi:10.3390/rs17121983.
N. Zhao, C. Zhang, H. Zhang, Z. Jiang, How down-sampling affects supervised-learning-based image super-resolutions, in: SPIE Future Sensing Technologies 2024, 2024. URL: http://dx.doi.org/10.1117/12.3022989. doi:10.1117/12.3022989.

J. Wu, L. Lin, T. Li, Q. Cheng, C. Zhang, H. Shen, Fusing Landsat 8 and Sentinel-2 Data for 10-m Dense Time-Series Imagery Using a Degradation-Term Constrained Deep Network, International Journal of Applied Earth Observation and Geoinformation 108 (2022) 102738. URL: http://dx.doi.org/10.1016/j.jag.2022.102738. doi:10.1016/j.jag.2022.102738.

A. Jaegle, S. Borgeaud, J.-B. Alayrac, C. Doersch, C. Ionescu, D. Ding, S. Kop-pula, D. Zoran, A. Brock, E. Shelhamer, O. Hénaff, M. M. Botvinick, A. Zis-serman, O. Vinyals, J. Carreira, Perceiver Io: a General Architecture for Struc-tured Inputs & Outputs (2021). URL: https://arxiv.org/abs/2107.14795. doi:10.48550/ARXIV.2107.14795.

O. Hagolle, J. Colin, S. Coustance, P. Kettig, P. D'Angelo, S. Auer, G. Doxani, C. Desjardins, Sentinel-2 Surface Reflectance Products Generated By CNES and DLR: Methods, Validation and Applications, ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences V-1-2021 (2021) 9-15. URL: http://dx.doi.org/10.5194/isprs-annals-V-1-2021-9-2021. doi:10.5194/isprs-annals-v-1-2021-9-2021.

## 6. Appendix: Complete results
<table><caption>Table 12: Comparison between TAMRF model, U-TILISE and Naive interpolation on the gap-filling task, where regular gaps of 30 days are masked from the input SITS and kept appart for validation. Snow target dates are excluded from metrics computation. * if usually higher FR means higher spatial resolution details, in the case of bands that are 10-meter native resolution, the best FR is the closest to zero in absolute value. Otherwise, ↓ (resp. ↑) indicates that the metric should be minimized (resp. maximized). Clear designate dates for which the bands were observed by the model, and masked designate the dates that were removed from the input SITS as described in Fig. 10.</caption><thead><tr><th rowspan="2">Band</th><th rowspan="2">Method</th><th colspan="4">RMSE↓</th><th colspan="4">BRISQUE↓</th><th colspan="4">FR↑</th></tr><tr><th>Clear mean</th><th>std</th><th>Masked mean</th><th>std</th><th>Clear mean</th><th>std</th><th>Masked mean</th><th>std</th><th>Clear mean</th><th>std</th><th>Masked mean</th><th>std</th></tr></thead><tbody><tr><td rowspan="2">LS B1</td><td>Naive</td><td><strong>0.003</strong></td><td>0.006</td><td>0.079</td><td>0.137</td><td>87.7</td><td>8.7</td><td>87.7</td><td>8.8</td><td>0.6</td><td>2.0</td><td>-0.2</td><td>2.8</td></tr><tr><td>TARMF</td><td>0.008</td><td>0.034</td><td><strong>0.015</strong></td><td>0.048</td><td>44.0</td><td><strong>18.7</strong></td><td>49.7</td><td><strong>21.1</strong></td><td><strong>7.3</strong></td><td>2.5</td><td><strong>7.7</strong></td><td>2.8</td></tr><tr><td rowspan="2">LS B2</td><td>Naive</td><td><strong>0.003</strong></td><td>0.006</td><td>0.079</td><td>0.135</td><td>84.4</td><td>8.4</td><td>84.5</td><td>8.3</td><td>0.5</td><td>1.8</td><td>0.1</td><td>2.8</td></tr><tr><td>TARMF</td><td>0.008</td><td>0.033</td><td><strong>0.016</strong></td><td>0.048</td><td><strong>37.4</strong></td><td>17.6</td><td><strong>42.2</strong></td><td>20.0</td><td><strong>7.5</strong></td><td>2.5</td><td><strong>7.8</strong></td><td>2.7</td></tr><tr><td rowspan="2">LS B3</td><td>Naive</td><td><strong>0.004</strong></td><td>0.005</td><td>0.076</td><td>0.122</td><td>78.4</td><td>7.0</td><td>78.4</td><td>6.4</td><td>0.7</td><td>1.8</td><td>0.4</td><td>2.8</td></tr><tr><td>TARMF</td><td>0.009</td><td>0.030</td><td><strong>0.018</strong></td><td>0.045</td><td><strong>24.9</strong></td><td>12.4</td><td><strong>27.4</strong></td><td>13.4</td><td><strong>8.0</strong></td><td>2.6</td><td><strong>8.2</strong></td><td>2.9</td></tr><tr><td rowspan="2">LS B4</td><td>Naive</td><td><strong>0.004</strong></td><td>0.005</td><td>0.079</td><td>0.120</td><td>72.8</td><td>6.5</td><td>73.3</td><td>6.6</td><td>0.7</td><td>1.7</td><td>0.3</td><td>2.7</td></tr><tr><td>TARMF</td><td>0.009</td><td>0.030</td><td><strong>0.022</strong></td><td>0.047</td><td><strong>20.8</strong></td><td>11.3</td><td><strong>21.2</strong></td><td>12.7</td><td><strong>8.0</strong></td><td>2.4</td><td><strong>8.2</strong></td><td>2.7</td></tr><tr><td rowspan="2">LS B5</td><td>Naive</td><td><strong>0.006</strong></td><td>0.004</td><td>0.074</td><td>0.086</td><td>66.9</td><td>7.6</td><td>66.6</td><td>6.7</td><td>0.7</td><td>1.8</td><td>0.3</td><td>2.3</td></tr><tr><td>TARMF</td><td>0.015</td><td>0.025</td><td><strong>0.037</strong></td><td>0.038</td><td><strong>15.2</strong></td><td>7.0</td><td><strong>11.5</strong></td><td>6.1</td><td><strong>7.7</strong></td><td>2.4</td><td><strong>8.1</strong></td><td>2.5</td></tr><tr><td rowspan="2">LS B6</td><td>Naive</td><td><strong>0.006</strong></td><td>0.003</td><td>0.049</td><td>0.050</td><td>66.9</td><td>5.6</td><td>66.9</td><td>5.0</td><td>0.6</td><td>1.5</td><td>0.2</td><td>2.4</td></tr><tr><td>TARMF</td><td>0.011</td><td>0.020</td><td><strong>0.030</strong></td><td>0.028</td><td><strong>17.6</strong></td><td>10.1</td><td><strong>9.4</strong></td><td>4.8</td><td><strong>8.2</strong></td><td>2.5</td><td><strong>8.4</strong></td><td>2.7</td></tr><tr><td rowspan="2">LS B7</td><td>Naive</td><td><strong>0.005</strong></td><td>0.002</td><td>0.041</td><td>0.042</td><td>68.7</td><td>5.4</td><td>69.1</td><td>5.0</td><td>0.6</td><td>1.6</td><td>0.1</td><td>2.3</td></tr><tr><td>TARMF</td><td>0.009</td><td>0.018</td><td><strong>0.025</strong></td><td>0.027</td><td><strong>18.9</strong></td><td>10.6</td><td><strong>9.8</strong></td><td>4.9</td><td><strong>8.0</strong></td><td>2.4</td><td><strong>8.5</strong></td><td>2.7</td></tr><tr><td rowspan="2">LS LST</td><td>Naive</td><td><strong>0.425</strong></td><td>0.430</td><td></td><td></td><td>81.7</td><td>8.4</td><td></td><td></td><td>0.8</td><td>3.3</td><td></td><td></td></tr><tr><td>TARMF</td><td>1.943</td><td>3.002</td><td></td><td></td><td><strong>24.4</strong></td><td>23.9</td><td></td><td></td><td><strong>3.1</strong></td><td>3.5</td><td></td><td></td></tr><tr><td rowspan="3">S2 B2</td><td>Naive</td><td><strong>0.000</strong></td><td>0.000</td><td>0.019</td><td>0.034</td><td>24.1</td><td>11.5</td><td>26.1</td><td>12.0</td><td>1.5</td><td>3.0</td><td>1.1</td><td>3.5</td></tr><tr><td>TARMF</td><td>0.008</td><td>0.006</td><td><strong>0.015</strong></td><td>0.023</td><td>31.7</td><td>16.0</td><td>35.2</td><td>17.4</td><td>1.0</td><td>3.0</td><td><strong>0.7*</strong></td><td>3.5</td></tr><tr><td>U-TILISE</td><td>0.006</td><td>0.009</td><td><strong>0.015</strong></td><td>0.007</td><td><strong>24.8</strong></td><td>9.4</td><td><strong>34.1</strong></td><td>10.1</td><td><strong>0.1*</strong></td><td>1.6</td><td>1.2</td><td>3.8</td></tr><tr><td rowspan="3">S2 B3</td><td>Naive</td><td><strong>0.000</strong></td><td>0.000</td><td>0.020</td><td>0.032</td><td>18.7</td><td>9.2</td><td>20.1</td><td>9.3</td><td>1.4</td><td>2.9</td><td>1.0</td><td>3.2</td></tr><tr><td>TARMF</td><td>0.007</td><td>0.006</td><td><strong>0.016</strong></td><td>0.023</td><td>20.0</td><td>11.9</td><td>24.9</td><td>13.5</td><td>1.2</td><td>2.8</td><td><strong>0.7*</strong></td><td>3.3</td></tr><tr><td>U-TILISE</td><td>0.007</td><td>0.009</td><td><strong>0.017</strong></td><td>0.008</td><td><strong>17.8</strong></td><td>6.6</td><td><strong>22.8</strong></td><td>6.2</td><td><strong>0.2*</strong></td><td>1.5</td><td>0.9</td><td>3.4</td></tr><tr><td rowspan="3">S2 B4</td><td>Naive</td><td><strong>0.000</strong></td><td>0.000</td><td>0.023</td><td>0.033</td><td>19.3</td><td>10.6</td><td>20.4</td><td>10.3</td><td>1.2</td><td>2.6</td><td>0.8</td><td>2.9</td></tr><tr><td>TARMF</td><td>0.008</td><td>0.006</td><td><strong>0.019</strong></td><td>0.023</td><td><strong>19.2</strong></td><td>12.7</td><td><strong>24.1</strong></td><td>14.0</td><td>1.1</td><td>2.6</td><td><strong>0.6*</strong></td><td>3.0</td></tr><tr><td>U-TILISE</td><td>0.007</td><td>0.009</td><td><strong>0.021</strong></td><td>0.010</td><td>21.5</td><td>7.8</td><td>25.7</td><td>7.2</td><td><strong>0.2*</strong></td><td>1.3</td><td>0.8</td><td>3.2</td></tr><tr><td rowspan="2">S2 B5</td><td>Naive</td><td><strong>0.004</strong></td><td>0.002</td><td>0.024</td><td>0.032</td><td>48.7</td><td>6.2</td><td>49.2</td><td>6.1</td><td>1.5</td><td>2.5</td><td>1.1</td><td>2.6</td></tr><tr><td>TARMF</td><td>0.008</td><td>0.006</td><td><strong>0.020</strong></td><td>0.023</td><td><strong>14.9</strong></td><td>10.3</td><td><strong>17.5</strong></td><td>11.5</td><td><strong>5.7</strong></td><td>2.1</td><td><strong>5.2</strong></td><td>2.6</td></tr><tr><td rowspan="2">S2 B6</td><td>Naive</td><td><strong>0.005</strong></td><td>0.002</td><td>0.031</td><td>0.030</td><td>47.4</td><td>6.3</td><td>48.0</td><td>6.1</td><td>1.2</td><td>2.1</td><td>1.0</td><td>2.3</td></tr><tr><td>TARMF</td><td>0.011</td><td>0.007</td><td><strong>0.027</strong></td><td>0.023</td><td><strong>14.2</strong></td><td>9.8</td><td><strong>16.2</strong></td><td>10.5</td><td><strong>5.5</strong></td><td>1.8</td><td><strong>5.1</strong></td><td>2.3</td></tr><tr><td rowspan="2">S2 B7</td><td>Naive</td><td><strong>0.006</strong></td><td>0.003</td><td>0.034</td><td>0.030</td><td>46.7</td><td>6.4</td><td>47.4</td><td>6.2</td><td>1.1</td><td>2.0</td><td>0.9</td><td>2.3</td></tr><tr><td>TARMF</td><td>0.012</td><td>0.008</td><td><strong>0.030</strong></td><td>0.024</td><td><strong>14.6</strong></td><td>9.7</td><td><strong>16.0</strong></td><td>10.5</td><td><strong>5.4</strong></td><td>1.7</td><td><strong>5.0</strong></td><td>2.2</td></tr><tr><td rowspan="3">S2 B8</td><td>Naive</td><td><strong>0.000</strong></td><td>0.000</td><td>0.036</td><td>0.030</td><td>13.8</td><td>8.2</td><td>14.7</td><td>7.8</td><td>0.7</td><td>1.9</td><td>0.4</td><td>2.3</td></tr><tr><td>TARMF</td><td>0.012</td><td>0.009</td><td><strong>0.032</strong></td><td>0.024</td><td>14.4</td><td>10.2</td><td>17.7</td><td>10.6</td><td>0.6</td><td>1.9</td><td>0.3</td><td>2.3</td></tr><tr><td>U-TILISE</td><td>0.013</td><td>0.017</td><td><strong>0.044</strong></td><td>0.020</td><td><strong>12.4</strong></td><td>6.8</td><td><strong>16.5</strong></td><td>6.4</td><td><strong>0.1*</strong></td><td>1.2</td><td><strong>0.2*</strong></td><td>2.4</td></tr><tr><td rowspan="2">S2 B8a</td><td>Naive</td><td><strong>0.006</strong></td><td>0.003</td><td>0.035</td><td>0.028</td><td>47.3</td><td>6.4</td><td>47.9</td><td>6.3</td><td>1.1</td><td>1.9</td><td>0.9</td><td>2.2</td></tr><tr><td>TARMF</td><td>0.012</td><td>0.009</td><td><strong>0.031</strong></td><td>0.024</td><td><strong>14.7</strong></td><td>9.7</td><td><strong>16.0</strong></td><td>10.4</td><td><strong>5.4</strong></td><td>1.6</td><td><strong>5.1</strong></td><td>2.1</td></tr><tr><td rowspan="2">S2 B11</td><td>Naive</td><td><strong>0.004</strong></td><td>0.002</td><td>0.030</td><td>0.021</td><td>53.8</td><td>7.4</td><td>54.3</td><td>6.9</td><td>1.3</td><td>2.1</td><td>1.1</td><td>2.4</td></tr><tr><td>TARMF</td><td>0.010</td><td>0.009</td><td><strong>0.027</strong></td><td>0.018</td><td><strong>10.9</strong></td><td>9.1</td><td><strong>13.3</strong></td><td>9.7</td><td><strong>6.6</strong></td><td>1.9</td><td><strong>6.1</strong></td><td>2.4</td></tr><tr><td rowspan="2">S2 B12</td><td>Naive</td><td><strong>0.004</strong></td><td>0.002</td><td>0.026</td><td>0.019</td><td>53.0</td><td>7.2</td><td>53.6</td><td>6.8</td><td>1.3</td><td>2.1</td><td>1.1</td><td>2.3</td></tr><tr><td>TARMF</td><td>0.009</td><td>0.008</td><td><strong>0.023</strong></td><td>0.017</td><td><strong>11.8</strong></td><td>9.2</td><td><strong>14.3</strong></td><td>10.0</td><td><strong>6.5</strong></td><td>1.9</td><td><strong>5.9</strong></td><td>2.4</td></tr></tbody></table>
<table><caption>Table 13: Comparison between TAMRF and DSen2 on the sharpening of Sentinel-2 20m bands, for a selection of bands. Snow target dates are excluded from metrics computation. ↓ (resp. ↑) indicates that the metric should be minimized (resp. maximized). Best mean values for each band and each metric are highlighted in **bold**.</caption><thead><tr><th rowspan="2">Band</th><th rowspan="2">Method</th><th colspan="2">RMSE↓</th><th colspan="2">BRISQUE↓</th><th colspan="2">FR↑</th></tr><tr><th>mean</th><th>std</th><th>mean</th><th>std</th><th>mean</th><th>std</th></tr></thead><tbody><tr><td rowspan="2">B5</td><td>DSEN2</td><td><strong>0.002</strong></td><td>0.001</td><td>29.033</td><td>5.846</td><td>4.1</td><td>0.6</td></tr><tr><td>TAMRF</td><td>0.008</td><td>0.006</td><td><strong>14.763</strong></td><td>10.091</td><td><strong>5.6</strong></td><td>2.1</td></tr><tr><td rowspan="2">B6</td><td>DSEN2</td><td><strong>0.002</strong></td><td>0.001</td><td>25.907</td><td>6.494</td><td>4.4</td><td>0.6</td></tr><tr><td>TAMRF</td><td>0.011</td><td>0.007</td><td><strong>14.180</strong></td><td>9.401</td><td><strong>5.5</strong></td><td>1.7</td></tr><tr><td rowspan="2">B7</td><td>DSEN2</td><td><strong>0.003</strong></td><td>0.001</td><td>24.930</td><td>6.058</td><td>4.4</td><td>0.5</td></tr><tr><td>TAMRF</td><td>0.012</td><td>0.008</td><td><strong>14.617</strong></td><td>9.443</td><td><strong>5.3</strong></td><td>1.7</td></tr><tr><td rowspan="2">B8a</td><td>DSEN2</td><td><strong>0.003</strong></td><td>0.001</td><td>25.390</td><td>6.360</td><td>4.5</td><td>0.5</td></tr><tr><td>TAMRF</td><td>0.012</td><td>0.009</td><td><strong>14.775</strong></td><td>9.365</td><td><strong>5.3</strong></td><td>1.6</td></tr><tr><td rowspan="2">B11</td><td>DSEN2</td><td><strong>0.002</strong></td><td>0.001</td><td>33.575</td><td>6.501</td><td>3.8</td><td>0.7</td></tr><tr><td>TAMRF</td><td>0.010</td><td>0.009</td><td><strong>10.941</strong></td><td>8.930</td><td><strong>6.6</strong></td><td>1.9</td></tr><tr><td rowspan="2">B12</td><td>DSEN2</td><td><strong>0.002</strong></td><td>0.001</td><td>33.449</td><td>6.225</td><td>3.9</td><td>0.7</td></tr><tr><td>TAMRF</td><td>0.009</td><td>0.008</td><td><strong>11.858</strong></td><td>9.011</td><td><strong>6.5</strong></td><td>2.0</td></tr></tbody></table>
<table><caption>Table 14: Full results table of the comparison of the different methods on the spatio-temporal fusion task. Only a subset of bands is presented. Full results are available in appendix 6. * if usually higher FR means higher spatial resolution details, in the case of bands that are 10-meter native resolution, the best FR is the closest to zero in absolute value. Otherwise, ↓ (resp. ↑) indicates that the metric should be minimized (resp. maximized). Best mean values for each band and metric are highlighted in **bold**, and second best values are **underlined**. † SEN2LIKE does not correspond to the full Sen2Like processing but only to the fusion part of the pipeline, which has been reimplemented following the paper and source code of the original pipeline.</caption><thead><tr><th>Band</th><th>Method</th><th>RMSE↓</th><th>BRISQUE↓</th><th>FR↑</th><th></th></tr><tr><th></th><th></th><th>mean</th><th>std</th><th>mean</th><th>std</th></tr></thead><tbody><tr><td rowspan="6">B2</td><td>SEN2LIKE†</td><td>0.028</td><td>0.028</td><td>34.3</td><td>18.9</td><td><strong>0.1*</strong></td><td>3.5</td></tr><tr><td>STAIR</td><td>0.032</td><td>0.038</td><td><strong>27.7</strong></td><td>13.0</td><td><strong>-0.6*</strong></td><td>3.2</td></tr><tr><td>DEEPHARMO</td><td>0.034</td><td>0.014</td><td>50.7</td><td>7.0</td><td>-6.1</td><td>1.4</td></tr><tr><td>DSTFN</td><td>0.029</td><td>0.033</td><td>33.5</td><td>12.9</td><td><strong>-0.1*</strong></td><td>2.8</td></tr><tr><td>TAMRF (no lr)</td><td><u>0.013</u></td><td>0.019</td><td>33.1</td><td>16.3</td><td>0.7</td><td>3.2</td></tr><tr><td>TAMRF (no hr)</td><td>0.016</td><td>0.008</td><td>56.3</td><td>12.2</td><td>-3.7</td><td>3.0</td></tr><tr><td></td><td>TAMRF (full)</td><td><strong>0.012</strong></td><td>0.007</td><td><strong>32.5</strong></td><td>16.8</td><td><strong>0.6*</strong></td><td>3.1</td></tr><tr><td rowspan="6">B3</td><td>SEN2LIKE†</td><td>0.028</td><td>0.028</td><td>30.8</td><td>17.3</td><td><strong>-0.4*</strong></td><td>3.2</td></tr><tr><td>STAIR</td><td>0.034</td><td>0.038</td><td><strong>22.7</strong></td><td>11.0</td><td>-0.8</td><td>3.0</td></tr><tr><td>DEEPHARMO</td><td>0.033</td><td>0.011</td><td>46.9</td><td>6.1</td><td>-5.7</td><td>1.2</td></tr><tr><td>DSTFN</td><td>0.028</td><td>0.033</td><td>28.4</td><td>7.3</td><td><strong>0.0*</strong></td><td>2.5</td></tr><tr><td>TAMRF (no lr)</td><td><u>0.014</u></td><td>0.018</td><td><strong>23.5</strong></td><td>12.0</td><td>0.8</td><td>2.9</td></tr><tr><td>TAMRF (no hr)</td><td>0.016</td><td>0.008</td><td>53.0</td><td>9.9</td><td>-3.7</td><td>2.7</td></tr><tr><td></td><td>TAMRF (full)</td><td><strong>0.012</strong></td><td>0.007</td><td>23.9</td><td>12.8</td><td>0.6</td><td>2.8</td></tr><tr><td rowspan="6">B4</td><td>SEN2LIKE†</td><td>0.032</td><td>0.030</td><td>30.1</td><td>17.8</td><td><strong>-0.5*</strong></td><td>2.8</td></tr><tr><td>STAIR</td><td>0.038</td><td>0.040</td><td><strong>22.1</strong></td><td>11.4</td><td>-0.7</td><td>2.6</td></tr><tr><td>DEEPHARMO</td><td>0.027</td><td>0.010</td><td>44.2</td><td>5.4</td><td>-5.3</td><td>1.0</td></tr><tr><td>DSTFN</td><td>0.033</td><td>0.035</td><td>28.4</td><td>8.7</td><td><strong>-0.0*</strong></td><td>2.3</td></tr><tr><td>TAMRF (no lr)</td><td><u>0.016</u></td><td>0.019</td><td>22.6</td><td>12.6</td><td>0.7</td><td>2.5</td></tr><tr><td>TAMRF (no hr)</td><td>0.019</td><td>0.009</td><td>49.8</td><td>9.2</td><td>-3.5</td><td>2.5</td></tr><tr><td></td><td>TAMRF (full)</td><td><strong>0.014</strong></td><td>0.008</td><td><strong>21.4</strong></td><td>12.8</td><td>0.6</td><td>2.5</td></tr><tr><td rowspan="6">B8a</td><td>SEN2LIKE†</td><td>0.048</td><td>0.044</td><td>58.4</td><td>12.4</td><td>-0.5</td><td>2.3</td></tr><tr><td>STAIR</td><td>0.056</td><td>0.050</td><td>48.4</td><td>5.9</td><td>-0.5</td><td>2.4</td></tr><tr><td>DEEPHARMO</td><td>0.026</td><td>0.003</td><td>62.4</td><td>4.7</td><td>-3.0</td><td>0.4</td></tr><tr><td>DSTFN</td><td>0.049</td><td>0.048</td><td>25.3</td><td>6.8</td><td>2.9</td><td>2.4</td></tr><tr><td>TAMRF (no lr)</td><td><u>0.026</u></td><td>0.020</td><td><strong>15.6</strong></td><td>8.4</td><td><strong>5.1</strong></td><td>1.8</td></tr><tr><td>TAMRF (no hr)</td><td>0.031</td><td>0.019</td><td>47.4</td><td>5.5</td><td>1.1</td><td>1.6</td></tr><tr><td></td><td>TAMRF (full)</td><td><strong>0.024</strong></td><td>0.013</td><td><strong>17.4</strong></td><td>9.5</td><td><strong>4.8</strong></td><td>1.8</td></tr><tr><td rowspan="6">B11</td><td>SEN2LIKE†</td><td>0.046</td><td>0.036</td><td>67.3</td><td>12.3</td><td>-0.0</td><td>2.1</td></tr><tr><td>STAIR</td><td>0.050</td><td>0.043</td><td>54.0</td><td>6.1</td><td>-0.1</td><td>2.1</td></tr><tr><td>DEEPHARMO</td><td>0.024</td><td>0.007</td><td>64.1</td><td>6.3</td><td>-2.2</td><td>0.4</td></tr><tr><td>DSTFN</td><td>0.048</td><td>0.040</td><td>28.4</td><td>9.6</td><td>3.8</td><td>2.3</td></tr><tr><td>TAMRF (no lr)</td><td><u>0.024</u></td><td>0.017</td><td><strong>12.7</strong></td><td>8.3</td><td><strong>6.1</strong></td><td>2.0</td></tr><tr><td>TAMRF (no hr)</td><td>0.027</td><td>0.014</td><td>47.0</td><td>5.9</td><td>2.2</td><td>1.8</td></tr><tr><td></td><td>TAMRF (full)</td><td><strong>0.022</strong></td><td>0.013</td><td><strong>16.0</strong></td><td>9.1</td><td><strong>5.8</strong></td><td>1.8</td></tr><tr><td rowspan="6">B12</td><td>SEN2LIKE†</td><td>0.041</td><td>0.029</td><td>66.8</td><td>13.2</td><td>-0.0</td><td>1.9</td></tr><tr><td>STAIR</td><td>0.042</td><td>0.034</td><td>53.1</td><td>6.4</td><td>-0.1</td><td>1.9</td></tr><tr><td>DEEPHARMO</td><td>0.021</td><td>0.008</td><td>65.3</td><td>7.3</td><td>-2.3</td><td>0.6</td></tr><tr><td>DSTFN</td><td>0.044</td><td>0.032</td><td>33.2</td><td>9.2</td><td>3.5</td><td>2.1</td></tr><tr><td>TAMRF (no lr)</td><td><u>0.021</u></td><td>0.016</td><td><strong>13.6</strong></td><td>8.4</td><td><strong>5.9</strong></td><td>1.8</td></tr><tr><td>TAMRF (no hr)</td><td>0.022</td><td>0.012</td><td>46.6</td><td>6.5</td><td>1.9</td><td>1.7</td></tr><tr><td></td><td>TAMRF (full)</td><td><strong>0.019</strong></td><td>0.012</td><td><strong>16.7</strong></td><td>9.4</td><td><strong>5.6</strong></td><td>1.8</td></tr></tbody></table>