# Multi-temporal Sentinel-2 super-resolution by optimization

> Source: <https://github.com/calebrob6/s2-superres>
> Published: 2026-06-24 09:31:40+00:00

**Jump to: Quick start | How it works | Scripts | Repo layout | Citation**

Reconstruct a high-resolution Sentinel-2 reflectance scene from ~30 low-resolution observations of the same area by exploiting the sub-pixel jitter between repeat passes. The scene is parameterized as a continuous field of 2D Gaussian splats; each S2 observation is the analytic integral of that field over shifted pixel footprints convolved with the sensor PSF. We jointly optimize the splat weights and the per-observation sub-pixel shifts in PyTorch. Pure optimization, with no training data and no neural network.

This repo provides:

- A
**PyTorch SR solver**() that fits a Gaussian splat field to a stack of Sentinel-2 observations using one of four optimization strategies (`sr.py`

`adam`

,`lbfgs`

,`adam_lbfgs`

,`lbfgs_adam`

), with PSF / TV / scale / coarse-to-fine knobs exposed on the command line. - A
**parameterized STAC downloader**() that pulls a cloud-free Sentinel-2 L2A time series over an arbitrary AOI / TOI from the Microsoft Planetary Computer and saves the cropped scenes as Cloud-Optimized GeoTIFFs.`download_data.py`

**Figure 1.** A scene near Redmond, WA reconstructed at 8× from 8 Sentinel-2 observations. (**Left**) One raw S2 observation at the original 10 m grid, nearest-neighbor displayed at the SR output size. (** Middle**) Bicubic upsample of the temporal mean across the 8 observations — clean but soft. (** Right**) The Gaussian-splat fit at 1.25 m, recovering edge structure that no single observation contains. Reproduce with `python experiments/psf_sweep.py`

; the per-σ outputs are committed under [ experiments/psf_sweep/](/calebrob6/s2-superres/blob/main/experiments/psf_sweep).

The achievable resolution gain is bounded by the ratio of the per-observation jitter to the sensor PSF width. For Sentinel-2 the jitter is roughly 1.5 m std and the PSF is roughly 3–5 m wide, so the ~30 shifted observations mostly contain redundant spatial information. Multi-temporal optimization produces clean denoised renders at 2–5 m, but it cannot recover genuine 1 m structure from the data alone. A sharper sensor, larger jitter, or an external structural prior (such as an aerial basemap) would be needed to push further.

```
git clone https://github.com/calebrob6/s2-superres.git
cd s2-superres
pip install -r requirements.txt

# Download the default Redmond, WA dataset (~32 cloud-free scenes, 1024x1024 at 10m)
python download_data.py

# Run super-resolution with the default lbfgs_adam strategy on a 256x256 crop
python sr.py --crop 256 --save-png

# Look at the result
ls output/sr_*/
```

The default `sr.py`

invocation runs LBFGS coarse-to-fine to recover the splat weights and per-observation shifts, then Adam sharpens the result with a TV anneal followed by an unregularized phase. On a single GPU a 256 × 256 LR crop takes a few minutes; full ~1000 × 1000 scenes take longer.

To run on a different area, pass an AOI bbox and a date range to `download_data.py`

:

```
python download_data.py \
    --bbox -122.450,37.700,-122.350,37.800 \
    --start 2024-04-01 --end 2024-10-31 \
    --output-dir data_sf/

python sr.py --data-dir data_sf/ --crop 256 --save-png
```

Each LR observation is modeled as

```
y[t,c,i,j] = sum_k w_k[c] * erf_box(mu_k, sigma_eff, pixel(i,j) - delta_t)
```

where `w_k[c]`

is the per-channel weight of splat `k`

, `erf_box(...)`

is the analytic Gaussian integral over a square LR pixel footprint, and `delta_t`

is a learnable per-observation 2D shift. Both the splat basis and the PSF are Gaussian, so their convolution collapses to `sigma_eff = sqrt(sigma_splat^2 + sigma_psf^2)`

evaluated at the splat-pixel offset. There is no HR pixel grid in the forward model and no numerical integration. The gradient flows analytically through `erf`

to both the weights and the shifts.

The optimizer picks one of four strategies:

| Strategy | Description |
|---|---|
`adam` |
Joint Adam on weights and shifts at one splat scale. Simple baseline. |
`lbfgs` |
Block-coordinate LBFGS, alternating between weights and shifts, run coarse-to-fine across `--c2f-levels` . |
`adam_lbfgs` |
Short Adam warmup at the coarsest level (escapes bad basins), then LBFGS C2F. |
`lbfgs_adam` |
LBFGS C2F to convergence, then a long Adam phase that anneals TV to zero and then runs unregularized for sharpening. Default; highest measured edge contrast. |

The shared core is in [ src/](/calebrob6/s2-superres/blob/main/src):

—`src/splat_field.py`

`SplatField`

`nn.Module`

with the analytic`erf`

forward model.—`src/strategies.py`

`run_strategy(...)`

dispatches to the four strategies above.—`src/sharpness.py`

`edge_contrast`

,`laplacian_var`

,`gradient_energy`

for evaluating real-data outputs.— phase-correlation shift initialization with parabolic sub-pixel refinement.`src/shift_estimation.py`

— STAC GeoTIFF stack loader with cloud masking and COG writer.`src/data_loader.py`

Searches the Microsoft Planetary Computer STAC catalog and saves cropped Sentinel-2 L2A scenes as COGs.

```
python download_data.py \
    --bbox -122.186,47.629,-122.056,47.719 \
    --start 2025-04-01 --end 2025-10-31 \
    --max-cloud-cover 1.0 \
    --bands B02,B03,B04,B08 \
    --min-scenes 32 \
    --output-dir data/
```

Pass `--basemap`

to also fetch an ESRI World Imagery aerial mosaic over the same AOI for visual comparison. Run `python download_data.py --help`

for the full flag list.

```
python sr.py \
    --data-dir data/ \
    --output-dir output/my_run/ \
    --strategy lbfgs_adam \
    --psf 0.5 \
    --tv 1e-3 \
    --c2f-levels 2,4,8 \
    --n-obs 8 \
    --crop 256 \
    --save-png
```

Useful flags:

| Flag | Description |
|---|---|
`--strategy` |
One of `{adam, lbfgs, adam_lbfgs, lbfgs_adam}` . |
`--psf` |
Sensor PSF Gaussian σ in LR pixels (Sentinel-2 is approximately 0.4–0.6). |
`--tv` |
Huber-TV regularization weight on the splat weights. |
`--c2f-levels` |
Comma-separated splat-scale ladder for the LBFGS strategies. |
`--n-obs` |
Number of lowest-cloud observations to use. |
`--crop` |
Center-crop size in LR pixels (`0` = full scene). |
`--save-cog` |
Write a georeferenced Cloud-Optimized GeoTIFF (only when `--crop 0` ). |

Run `python sr.py --help`

for everything.

Sweeps `sigma_psf ∈ {0.3, 0.4, 0.5, 0.6, 0.7, 0.8}`

over the central 256 × 256 crop with `lbfgs_adam`

and otherwise-default settings. Writes per-σ RGB PNGs, a 2 × 4 comparison grid (input + bicubic + 6 σ values), an `edge_contrast`

vs `sigma_psf`

plot, and a metrics CSV to [ experiments/psf_sweep/](/calebrob6/s2-superres/blob/main/experiments/psf_sweep). The PNGs in this README come from this script.

```
python experiments/psf_sweep.py
```

See [ experiments/psf_sweep/README.md](/calebrob6/s2-superres/blob/main/experiments/psf_sweep/README.md) for the most recent results.

```
sr.py                       SR entry point
download_data.py            Parameterized STAC downloader
src/                        Reusable library (SplatField, strategies, sharpness, ...)
experiments/                Experiment scripts and their committed PNG/CSV outputs
images/                     Figures embedded in this README
```

`data/`

and `output/`

are gitignored. Fill them locally with `download_data.py`

and `sr.py`

runs.

If you use this repo, please cite it:

```
@misc{robinson2026s2superres,
  author       = {Robinson, Caleb},
  title        = {{s2-superres}: multi-temporal {Sentinel-2} super-resolution by optimization},
  year         = {2026},
  howpublished = {\url{https://github.com/calebrob6/s2-superres}}
}
```

MIT. See [ LICENSE](/calebrob6/s2-superres/blob/main/LICENSE).