Super-resolve Sentinel-2 image resolution while ensuring spatial and radiometric consistency¶
Notebook prepared by Ben Maathuis. ITC-University of Twente, Enschede. The Netherlands
To upscale your Sentinel-2 data from 10 meters to 2.5 meters using Python, you can utilize the official, open-source opensr-model Python package managed by the ESA OpenSR ecosystem. This tool applies a Latent Diffusion Super-Resolution model specifically optimized for Sentinel-2's 10m bands (RGB + NIR). See also: https://github.com/ESAOpenSR and https://opensr.eu/
The step-by-step notebook prepared below demonstrates how to prepare your local data, set up the environment, and execute the 4x super-resolution process using a standard PyTorch machine learning workflow.
Here use is made of the Copernicus DataSpace Ecosystem (https://dataspace.copernicus.eu/). Before you continue register and create an account (for free). For online selection of cloud free images check the Copernicus Browser, available at: https://browser.dataspace.copernicus.eu/
Based on the lower left coordinate specified, which can be retrieved using the Copernicus Browser a small area of interest will be retrieved. The lower left coordinate is provided in degrees, and for an approximate window of 128 by 128 lines and columns the upper right coordinate pair is calculated based on the Sentinel-2 MSI VIS spatial resolution (10 metres) which is subsequently converted to degrees to get the upper right corner coordinate.
Ensure that the required packages are installed, especially torch, OmegaConf and opensr_model
Installing the packages, setting the folders and connecting to the Copernicus data space ecosystem¶
In [1]:
#import required libraries - eventually install if not available
import torch
import rasterio
from omegaconf import OmegaConf
import opensr_model
import numpy as np
import os
from io import StringIO
import requests
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
from matplotlib.ticker import FormatStrFormatter
import cartopy
from cartopy import crs as ccrs, feature as cfeature
import ilwis
import openeo
from openeo.rest.auth.config import RefreshTokenStore
import math
#suppress warnings
import warnings
warnings.filterwarnings('ignore')
warnings.simplefilter('ignore')
In [2]:
#within the notebook folder a directory is created to store your processing results
work_dir = os.getcwd()+'/result'
print("current dir is: %s" % (os.getcwd()))
print("current working directory is:",work_dir)
if os.path.isdir(work_dir):
print("Folder exists")
else:
print("Folder doesn't exists")
os.mkdir(work_dir)
current dir is: d:\jupyter\notebook_scripts\Special\Super_Resolution current working directory is: d:\jupyter\notebook_scripts\Special\Super_Resolution/result Folder exists
In [3]:
#set the working directory for ILWISPy
ilwis.setWorkingCatalog(work_dir)
print(work_dir)
d:\jupyter\notebook_scripts\Special\Super_Resolution/result
In [4]:
ilwis.version()
Out[4]:
'1.0 build 20260708'
In [5]:
#uncomment line below to remove previous login credentials
#RefreshTokenStore().remove()
In [6]:
#ensure openeo is installed and you have registered
connection = openeo.connect("openeo.dataspace.copernicus.eu").authenticate_oidc()
Authenticated using refresh token.
Download and visualize selected Sentinel L2A image subset¶
In [7]:
connection.describe_collection("SENTINEL2_L2A")
Out[7]:
We want an AoI size of 128 lines by 128 columns, coordinate system = EPSG:4326 and spatial resolution is 10 meters, we use pixel centre - for 128 pixels, there are 127 intervals between pixel centres. Specify the lower left corner pixel in the code field below
In [8]:
# Specify the lower-left pixel centre, for interactive selection use the Copernicus Browser and note the coordinate in the lower right hand corner - see link above
lon_ll = 6.75479
lat_ll = 52.29420
In [9]:
#Specify you selected date, here the date of a cloud free AoI should be selected - select the date using the Copernicus Browser - see link above
date_start ="2026-08-04"
date_end = date_start
In [10]:
# Sentinel-2 MSI Resolution
resolution = 10 # metres/pixel
# Image dimensions of selected window
ncols = 128
nrows = 128
# Approximate conversion from metres to degrees
meters_per_degree_lat = 111320
meters_per_degree_lon = (
111320 * np.cos(np.deg2rad(lat_ll))
)
dlon = resolution / meters_per_degree_lon
dlat = resolution / meters_per_degree_lat
# Upper-right pixel centre
lon_ur = lon_ll + (ncols - 1) * dlon
lat_ur = lat_ll + (nrows - 1) * dlat
print(f"Pixel size: {dlon:.8f}° lon × {dlat:.8f}° lat")
print(f"Lower-left: ({lon_ll:.8f}, {lat_ll:.8f})")
print(f"Upper-right: ({lon_ur:.8f}, {lat_ur:.8f})")
Pixel size: 0.00014688° lon × 0.00008983° lat Lower-left: (6.75479000, 52.29420000) Upper-right: (6.77344340, 52.30560855)
In [11]:
#note the date
t = [date_start, date_end]
s2_cube = connection.load_collection("SENTINEL2_L2A",
spatial_extent={'west': lon_ll, 'east': lon_ur, 'south': lat_ll, 'north': lat_ur, "crs": "EPSG:4326" },
temporal_extent= t,
bands=['B02', 'B03', 'B04', 'B08'],
max_cloud_cover=100,
)
In [12]:
s2_cube.download(work_dir+"/S2_selected.tif")
In [13]:
#read the S2 image in ilwispy and note the number of spectral bands and image size
S2_in = ilwis.RasterCoverage ('S2_selected.tif')
print(S2_in.size())
Size(132, 132, 4)
In [14]:
#stretch the rasterbands using a loop
multiple_stretch = []
multiple_bands = ilwis.do('selection',S2_in,"rasterbands(0..3)")
ls = ilwis.do('linearstretch',multiple_bands, 1) #using an upper and lower data limit defined by the cumulative 1 and 99 % thresholds
mb_stretch = ilwis.do('setvaluerange', ls, 0, 255, 1) #set ouput to byte range
In [15]:
#store the results as an ILWIS maplist - display the map using the ilwis386 desktop software
mb_stretch.store('S2_selected_stretch.mpl')
In [16]:
#load the individual spectral channels
Blues = ilwis.do('selection',mb_stretch,"rasterbands(0)")
Greens = ilwis.do('selection',mb_stretch,"rasterbands(1)")
Reds = ilwis.do('selection',mb_stretch,"rasterbands(2)")
In [17]:
#transform the spectral channels from ilwis format to a numpy array using the iterator
Blues_2np = np.fromiter(iter(Blues), np.ubyte, Blues.size().linearSize())
Blues_2np = Blues_2np.reshape((Blues.size().ysize, Blues.size().xsize))
Greens_2np = np.fromiter(iter(Greens), np.ubyte, Blues.size().linearSize())
Greens_2np = Greens_2np.reshape((Blues.size().ysize, Blues.size().xsize))
Reds_2np = np.fromiter(iter(Reds), np.ubyte, Blues.size().linearSize())
Reds_2np = Reds_2np.reshape((Blues.size().ysize, Blues.size().xsize))
In [18]:
#create a numpy natural color 3D data stack
ncol = np.dstack((Reds_2np, Greens_2np, Blues_2np))
Check the image obtained
In [19]:
img_extent = (lon_ll, lon_ur, lat_ll, lat_ur)
fig = plt.figure(figsize=(10, 14))
ax = plt.axes(projection=ccrs.PlateCarree())
plt.title('Selected Sentinel-2 composite of Area of Interest of ' + str(date_start))
# Data raster
ax.imshow(
ncol,
extent=img_extent,
transform=ccrs.PlateCarree()
)
# Coordinate grid
gl = ax.gridlines(
crs=ccrs.PlateCarree(),
draw_labels=True,
linewidth=0.6,
color="gray",
alpha=0.7,
linestyle="--",
zorder=5
)
# Grid spacing
gl.xlocator = mticker.MultipleLocator(0.005)
gl.ylocator = mticker.MultipleLocator(0.005)
# Labels
gl.top_labels = False
gl.right_labels = False
gl.xlabel_style = {
"size": 10,
}
gl.ylabel_style = {
"size": 10,
}
plt.show()
Ensure the submap has 128 lines and 128 columns for further processing
In [20]:
# sub map creation to ensure that the final AoI has 128 lines by 128 columns
rcSelect = ilwis.do('selection',S2_in,'boundingbox(1 1, 128 128)')
print(rcSelect.size().xsize)
print(rcSelect.size().ysize)
print(rcSelect.size().zsize)
128 128 4
Save the results, as a tiff file
In [21]:
rcSelect.store("AoI.tif", "GTiff", "gdal")
Start the Spatial Super Resolution Processing procedure¶
In [22]:
# -------------------------------------------------------------
# 1. Environment Setup & Configuration Setup
# -------------------------------------------------------------
# Determine device acceleration (GPU is heavily recommended for Latent Diffusion)
device = "cuda" if torch.cuda.is_available() else "cpu"
print(f"Using processing engine: {device.upper()}")
# Fetch the official neural network configuration mapping from ESA OpenSR
config_url = (
"https://raw.githubusercontent.com/"
"ESAOpenSR/opensr-model/refs/heads/main/"
"opensr_model/configs/config_10m.yaml"
)
response = requests.get(config_url)
config = OmegaConf.load(StringIO(response.text))
# Instantiate the model and load its pretrained 4x upscaling weights
print("Downloading and preparing the Latent Diffusion model...")
model = opensr_model.SRLatentDiffusion(config, device=device)
model.load_pretrained(config.ckpt_version)
model.eval() # Put neural network into evaluation mode
Using processing engine: CPU Downloading and preparing the Latent Diffusion model... LatentDiffusion: Running in eps-prediction mode DiffusionWrapper has 113.63 M params. Keeping EMAs of 308. making attention of type 'vanilla' with 512 in_channels Working with z of shape (1, 4, 128, 128) = 65536 dimensions. making attention of type 'vanilla' with 512 in_channels Normalization disabled. Loaded pretrained weights from: opensr-ldsrs2_v1_0_0.ckpt
Out[22]:
SRLatentDiffusion(
(model): LatentDiffusion(
(model): DiffusionWrapper(
(diffusion_model): UNetModel(
(time_embed): Sequential(
(0): Linear(in_features=160, out_features=640, bias=True)
(1): SiLU()
(2): Linear(in_features=640, out_features=640, bias=True)
)
(input_blocks): ModuleList(
(0): TimestepEmbedSequential(
(0): Conv2d(8, 160, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
(1-2): 2 x TimestepEmbedSequential(
(0): ResBlock(
(in_layers): Sequential(
(0): GroupNorm32(32, 160, eps=1e-05, affine=True, bias=True)
(1): SiLU()
(2): Conv2d(160, 160, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
(h_upd): Identity()
(x_upd): Identity()
(emb_layers): Sequential(
(0): SiLU()
(1): Linear(in_features=640, out_features=160, bias=True)
)
(out_layers): Sequential(
(0): GroupNorm32(32, 160, eps=1e-05, affine=True, bias=True)
(1): SiLU()
(2): Dropout(p=0, inplace=False)
(3): Conv2d(160, 160, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
(skip_connection): Identity()
)
)
(3): TimestepEmbedSequential(
(0): Downsample(
(op): Conv2d(160, 160, kernel_size=(3, 3), stride=(2, 2), padding=(1, 1))
)
)
(4): TimestepEmbedSequential(
(0): ResBlock(
(in_layers): Sequential(
(0): GroupNorm32(32, 160, eps=1e-05, affine=True, bias=True)
(1): SiLU()
(2): Conv2d(160, 320, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
(h_upd): Identity()
(x_upd): Identity()
(emb_layers): Sequential(
(0): SiLU()
(1): Linear(in_features=640, out_features=320, bias=True)
)
(out_layers): Sequential(
(0): GroupNorm32(32, 320, eps=1e-05, affine=True, bias=True)
(1): SiLU()
(2): Dropout(p=0, inplace=False)
(3): Conv2d(320, 320, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
(skip_connection): Conv2d(160, 320, kernel_size=(1, 1), stride=(1, 1))
)
)
(5): TimestepEmbedSequential(
(0): ResBlock(
(in_layers): Sequential(
(0): GroupNorm32(32, 320, eps=1e-05, affine=True, bias=True)
(1): SiLU()
(2): Conv2d(320, 320, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
(h_upd): Identity()
(x_upd): Identity()
(emb_layers): Sequential(
(0): SiLU()
(1): Linear(in_features=640, out_features=320, bias=True)
)
(out_layers): Sequential(
(0): GroupNorm32(32, 320, eps=1e-05, affine=True, bias=True)
(1): SiLU()
(2): Dropout(p=0, inplace=False)
(3): Conv2d(320, 320, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
(skip_connection): Identity()
)
)
(6): TimestepEmbedSequential(
(0): Downsample(
(op): Conv2d(320, 320, kernel_size=(3, 3), stride=(2, 2), padding=(1, 1))
)
)
(7-8): 2 x TimestepEmbedSequential(
(0): ResBlock(
(in_layers): Sequential(
(0): GroupNorm32(32, 320, eps=1e-05, affine=True, bias=True)
(1): SiLU()
(2): Conv2d(320, 320, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
(h_upd): Identity()
(x_upd): Identity()
(emb_layers): Sequential(
(0): SiLU()
(1): Linear(in_features=640, out_features=320, bias=True)
)
(out_layers): Sequential(
(0): GroupNorm32(32, 320, eps=1e-05, affine=True, bias=True)
(1): SiLU()
(2): Dropout(p=0, inplace=False)
(3): Conv2d(320, 320, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
(skip_connection): Identity()
)
)
(9): TimestepEmbedSequential(
(0): Downsample(
(op): Conv2d(320, 320, kernel_size=(3, 3), stride=(2, 2), padding=(1, 1))
)
)
(10): TimestepEmbedSequential(
(0): ResBlock(
(in_layers): Sequential(
(0): GroupNorm32(32, 320, eps=1e-05, affine=True, bias=True)
(1): SiLU()
(2): Conv2d(320, 640, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
(h_upd): Identity()
(x_upd): Identity()
(emb_layers): Sequential(
(0): SiLU()
(1): Linear(in_features=640, out_features=640, bias=True)
)
(out_layers): Sequential(
(0): GroupNorm32(32, 640, eps=1e-05, affine=True, bias=True)
(1): SiLU()
(2): Dropout(p=0, inplace=False)
(3): Conv2d(640, 640, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
(skip_connection): Conv2d(320, 640, kernel_size=(1, 1), stride=(1, 1))
)
(1): AttentionBlock(
(norm): GroupNorm32(32, 640, eps=1e-05, affine=True, bias=True)
(qkv): Conv1d(640, 1920, kernel_size=(1,), stride=(1,))
(attention): QKVAttentionLegacy()
(proj_out): Conv1d(640, 640, kernel_size=(1,), stride=(1,))
)
)
(11): TimestepEmbedSequential(
(0): ResBlock(
(in_layers): Sequential(
(0): GroupNorm32(32, 640, eps=1e-05, affine=True, bias=True)
(1): SiLU()
(2): Conv2d(640, 640, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
(h_upd): Identity()
(x_upd): Identity()
(emb_layers): Sequential(
(0): SiLU()
(1): Linear(in_features=640, out_features=640, bias=True)
)
(out_layers): Sequential(
(0): GroupNorm32(32, 640, eps=1e-05, affine=True, bias=True)
(1): SiLU()
(2): Dropout(p=0, inplace=False)
(3): Conv2d(640, 640, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
(skip_connection): Identity()
)
(1): AttentionBlock(
(norm): GroupNorm32(32, 640, eps=1e-05, affine=True, bias=True)
(qkv): Conv1d(640, 1920, kernel_size=(1,), stride=(1,))
(attention): QKVAttentionLegacy()
(proj_out): Conv1d(640, 640, kernel_size=(1,), stride=(1,))
)
)
)
(middle_block): TimestepEmbedSequential(
(0): ResBlock(
(in_layers): Sequential(
(0): GroupNorm32(32, 640, eps=1e-05, affine=True, bias=True)
(1): SiLU()
(2): Conv2d(640, 640, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
(h_upd): Identity()
(x_upd): Identity()
(emb_layers): Sequential(
(0): SiLU()
(1): Linear(in_features=640, out_features=640, bias=True)
)
(out_layers): Sequential(
(0): GroupNorm32(32, 640, eps=1e-05, affine=True, bias=True)
(1): SiLU()
(2): Dropout(p=0, inplace=False)
(3): Conv2d(640, 640, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
(skip_connection): Identity()
)
(1): AttentionBlock(
(norm): GroupNorm32(32, 640, eps=1e-05, affine=True, bias=True)
(qkv): Conv1d(640, 1920, kernel_size=(1,), stride=(1,))
(attention): QKVAttentionLegacy()
(proj_out): Conv1d(640, 640, kernel_size=(1,), stride=(1,))
)
(2): ResBlock(
(in_layers): Sequential(
(0): GroupNorm32(32, 640, eps=1e-05, affine=True, bias=True)
(1): SiLU()
(2): Conv2d(640, 640, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
(h_upd): Identity()
(x_upd): Identity()
(emb_layers): Sequential(
(0): SiLU()
(1): Linear(in_features=640, out_features=640, bias=True)
)
(out_layers): Sequential(
(0): GroupNorm32(32, 640, eps=1e-05, affine=True, bias=True)
(1): SiLU()
(2): Dropout(p=0, inplace=False)
(3): Conv2d(640, 640, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
(skip_connection): Identity()
)
)
(output_blocks): ModuleList(
(0-1): 2 x TimestepEmbedSequential(
(0): ResBlock(
(in_layers): Sequential(
(0): GroupNorm32(32, 1280, eps=1e-05, affine=True, bias=True)
(1): SiLU()
(2): Conv2d(1280, 640, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
(h_upd): Identity()
(x_upd): Identity()
(emb_layers): Sequential(
(0): SiLU()
(1): Linear(in_features=640, out_features=640, bias=True)
)
(out_layers): Sequential(
(0): GroupNorm32(32, 640, eps=1e-05, affine=True, bias=True)
(1): SiLU()
(2): Dropout(p=0, inplace=False)
(3): Conv2d(640, 640, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
(skip_connection): Conv2d(1280, 640, kernel_size=(1, 1), stride=(1, 1))
)
(1): AttentionBlock(
(norm): GroupNorm32(32, 640, eps=1e-05, affine=True, bias=True)
(qkv): Conv1d(640, 1920, kernel_size=(1,), stride=(1,))
(attention): QKVAttentionLegacy()
(proj_out): Conv1d(640, 640, kernel_size=(1,), stride=(1,))
)
)
(2): TimestepEmbedSequential(
(0): ResBlock(
(in_layers): Sequential(
(0): GroupNorm32(32, 960, eps=1e-05, affine=True, bias=True)
(1): SiLU()
(2): Conv2d(960, 640, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
(h_upd): Identity()
(x_upd): Identity()
(emb_layers): Sequential(
(0): SiLU()
(1): Linear(in_features=640, out_features=640, bias=True)
)
(out_layers): Sequential(
(0): GroupNorm32(32, 640, eps=1e-05, affine=True, bias=True)
(1): SiLU()
(2): Dropout(p=0, inplace=False)
(3): Conv2d(640, 640, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
(skip_connection): Conv2d(960, 640, kernel_size=(1, 1), stride=(1, 1))
)
(1): AttentionBlock(
(norm): GroupNorm32(32, 640, eps=1e-05, affine=True, bias=True)
(qkv): Conv1d(640, 1920, kernel_size=(1,), stride=(1,))
(attention): QKVAttentionLegacy()
(proj_out): Conv1d(640, 640, kernel_size=(1,), stride=(1,))
)
(2): Upsample(
(conv): Conv2d(640, 640, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
)
(3): TimestepEmbedSequential(
(0): ResBlock(
(in_layers): Sequential(
(0): GroupNorm32(32, 960, eps=1e-05, affine=True, bias=True)
(1): SiLU()
(2): Conv2d(960, 320, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
(h_upd): Identity()
(x_upd): Identity()
(emb_layers): Sequential(
(0): SiLU()
(1): Linear(in_features=640, out_features=320, bias=True)
)
(out_layers): Sequential(
(0): GroupNorm32(32, 320, eps=1e-05, affine=True, bias=True)
(1): SiLU()
(2): Dropout(p=0, inplace=False)
(3): Conv2d(320, 320, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
(skip_connection): Conv2d(960, 320, kernel_size=(1, 1), stride=(1, 1))
)
)
(4): TimestepEmbedSequential(
(0): ResBlock(
(in_layers): Sequential(
(0): GroupNorm32(32, 640, eps=1e-05, affine=True, bias=True)
(1): SiLU()
(2): Conv2d(640, 320, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
(h_upd): Identity()
(x_upd): Identity()
(emb_layers): Sequential(
(0): SiLU()
(1): Linear(in_features=640, out_features=320, bias=True)
)
(out_layers): Sequential(
(0): GroupNorm32(32, 320, eps=1e-05, affine=True, bias=True)
(1): SiLU()
(2): Dropout(p=0, inplace=False)
(3): Conv2d(320, 320, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
(skip_connection): Conv2d(640, 320, kernel_size=(1, 1), stride=(1, 1))
)
)
(5): TimestepEmbedSequential(
(0): ResBlock(
(in_layers): Sequential(
(0): GroupNorm32(32, 640, eps=1e-05, affine=True, bias=True)
(1): SiLU()
(2): Conv2d(640, 320, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
(h_upd): Identity()
(x_upd): Identity()
(emb_layers): Sequential(
(0): SiLU()
(1): Linear(in_features=640, out_features=320, bias=True)
)
(out_layers): Sequential(
(0): GroupNorm32(32, 320, eps=1e-05, affine=True, bias=True)
(1): SiLU()
(2): Dropout(p=0, inplace=False)
(3): Conv2d(320, 320, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
(skip_connection): Conv2d(640, 320, kernel_size=(1, 1), stride=(1, 1))
)
(1): Upsample(
(conv): Conv2d(320, 320, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
)
(6-7): 2 x TimestepEmbedSequential(
(0): ResBlock(
(in_layers): Sequential(
(0): GroupNorm32(32, 640, eps=1e-05, affine=True, bias=True)
(1): SiLU()
(2): Conv2d(640, 320, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
(h_upd): Identity()
(x_upd): Identity()
(emb_layers): Sequential(
(0): SiLU()
(1): Linear(in_features=640, out_features=320, bias=True)
)
(out_layers): Sequential(
(0): GroupNorm32(32, 320, eps=1e-05, affine=True, bias=True)
(1): SiLU()
(2): Dropout(p=0, inplace=False)
(3): Conv2d(320, 320, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
(skip_connection): Conv2d(640, 320, kernel_size=(1, 1), stride=(1, 1))
)
)
(8): TimestepEmbedSequential(
(0): ResBlock(
(in_layers): Sequential(
(0): GroupNorm32(32, 480, eps=1e-05, affine=True, bias=True)
(1): SiLU()
(2): Conv2d(480, 320, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
(h_upd): Identity()
(x_upd): Identity()
(emb_layers): Sequential(
(0): SiLU()
(1): Linear(in_features=640, out_features=320, bias=True)
)
(out_layers): Sequential(
(0): GroupNorm32(32, 320, eps=1e-05, affine=True, bias=True)
(1): SiLU()
(2): Dropout(p=0, inplace=False)
(3): Conv2d(320, 320, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
(skip_connection): Conv2d(480, 320, kernel_size=(1, 1), stride=(1, 1))
)
(1): Upsample(
(conv): Conv2d(320, 320, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
)
(9): TimestepEmbedSequential(
(0): ResBlock(
(in_layers): Sequential(
(0): GroupNorm32(32, 480, eps=1e-05, affine=True, bias=True)
(1): SiLU()
(2): Conv2d(480, 160, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
(h_upd): Identity()
(x_upd): Identity()
(emb_layers): Sequential(
(0): SiLU()
(1): Linear(in_features=640, out_features=160, bias=True)
)
(out_layers): Sequential(
(0): GroupNorm32(32, 160, eps=1e-05, affine=True, bias=True)
(1): SiLU()
(2): Dropout(p=0, inplace=False)
(3): Conv2d(160, 160, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
(skip_connection): Conv2d(480, 160, kernel_size=(1, 1), stride=(1, 1))
)
)
(10-11): 2 x TimestepEmbedSequential(
(0): ResBlock(
(in_layers): Sequential(
(0): GroupNorm32(32, 320, eps=1e-05, affine=True, bias=True)
(1): SiLU()
(2): Conv2d(320, 160, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
(h_upd): Identity()
(x_upd): Identity()
(emb_layers): Sequential(
(0): SiLU()
(1): Linear(in_features=640, out_features=160, bias=True)
)
(out_layers): Sequential(
(0): GroupNorm32(32, 160, eps=1e-05, affine=True, bias=True)
(1): SiLU()
(2): Dropout(p=0, inplace=False)
(3): Conv2d(160, 160, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
(skip_connection): Conv2d(320, 160, kernel_size=(1, 1), stride=(1, 1))
)
)
)
(out): Sequential(
(0): GroupNorm32(32, 160, eps=1e-05, affine=True, bias=True)
(1): SiLU()
(2): Conv2d(160, 4, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
)
)
(model_ema): LitEma()
(first_stage_model): AutoencoderKL(
(encoder): Encoder(
(conv_in): Conv2d(4, 128, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
(down): ModuleList(
(0): Module(
(block): ModuleList(
(0-1): 2 x ResnetBlock(
(norm1): GroupNorm(32, 128, eps=1e-06, affine=True, bias=True)
(conv1): Conv2d(128, 128, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
(norm2): GroupNorm(32, 128, eps=1e-06, affine=True, bias=True)
(dropout): Dropout(p=0.0, inplace=False)
(conv2): Conv2d(128, 128, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
)
(attn): ModuleList()
(downsample): Downsample(
(conv): Conv2d(128, 128, kernel_size=(3, 3), stride=(2, 2))
)
)
(1): Module(
(block): ModuleList(
(0): ResnetBlock(
(norm1): GroupNorm(32, 128, eps=1e-06, affine=True, bias=True)
(conv1): Conv2d(128, 256, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
(norm2): GroupNorm(32, 256, eps=1e-06, affine=True, bias=True)
(dropout): Dropout(p=0.0, inplace=False)
(conv2): Conv2d(256, 256, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
(nin_shortcut): Conv2d(128, 256, kernel_size=(1, 1), stride=(1, 1))
)
(1): ResnetBlock(
(norm1): GroupNorm(32, 256, eps=1e-06, affine=True, bias=True)
(conv1): Conv2d(256, 256, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
(norm2): GroupNorm(32, 256, eps=1e-06, affine=True, bias=True)
(dropout): Dropout(p=0.0, inplace=False)
(conv2): Conv2d(256, 256, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
)
(attn): ModuleList()
(downsample): Downsample(
(conv): Conv2d(256, 256, kernel_size=(3, 3), stride=(2, 2))
)
)
(2): Module(
(block): ModuleList(
(0): ResnetBlock(
(norm1): GroupNorm(32, 256, eps=1e-06, affine=True, bias=True)
(conv1): Conv2d(256, 512, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
(norm2): GroupNorm(32, 512, eps=1e-06, affine=True, bias=True)
(dropout): Dropout(p=0.0, inplace=False)
(conv2): Conv2d(512, 512, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
(nin_shortcut): Conv2d(256, 512, kernel_size=(1, 1), stride=(1, 1))
)
(1): ResnetBlock(
(norm1): GroupNorm(32, 512, eps=1e-06, affine=True, bias=True)
(conv1): Conv2d(512, 512, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
(norm2): GroupNorm(32, 512, eps=1e-06, affine=True, bias=True)
(dropout): Dropout(p=0.0, inplace=False)
(conv2): Conv2d(512, 512, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
)
(attn): ModuleList()
)
)
(mid): Module(
(block_1): ResnetBlock(
(norm1): GroupNorm(32, 512, eps=1e-06, affine=True, bias=True)
(conv1): Conv2d(512, 512, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
(norm2): GroupNorm(32, 512, eps=1e-06, affine=True, bias=True)
(dropout): Dropout(p=0.0, inplace=False)
(conv2): Conv2d(512, 512, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
(attn_1): AttnBlock(
(norm): GroupNorm(32, 512, eps=1e-06, affine=True, bias=True)
(q): Conv2d(512, 512, kernel_size=(1, 1), stride=(1, 1))
(k): Conv2d(512, 512, kernel_size=(1, 1), stride=(1, 1))
(v): Conv2d(512, 512, kernel_size=(1, 1), stride=(1, 1))
(proj_out): Conv2d(512, 512, kernel_size=(1, 1), stride=(1, 1))
)
(block_2): ResnetBlock(
(norm1): GroupNorm(32, 512, eps=1e-06, affine=True, bias=True)
(conv1): Conv2d(512, 512, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
(norm2): GroupNorm(32, 512, eps=1e-06, affine=True, bias=True)
(dropout): Dropout(p=0.0, inplace=False)
(conv2): Conv2d(512, 512, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
)
(norm_out): GroupNorm(32, 512, eps=1e-06, affine=True, bias=True)
(conv_out): Conv2d(512, 8, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
(decoder): Decoder(
(conv_in): Conv2d(4, 512, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
(mid): Module(
(block_1): ResnetBlock(
(norm1): GroupNorm(32, 512, eps=1e-06, affine=True, bias=True)
(conv1): Conv2d(512, 512, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
(norm2): GroupNorm(32, 512, eps=1e-06, affine=True, bias=True)
(dropout): Dropout(p=0.0, inplace=False)
(conv2): Conv2d(512, 512, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
(attn_1): AttnBlock(
(norm): GroupNorm(32, 512, eps=1e-06, affine=True, bias=True)
(q): Conv2d(512, 512, kernel_size=(1, 1), stride=(1, 1))
(k): Conv2d(512, 512, kernel_size=(1, 1), stride=(1, 1))
(v): Conv2d(512, 512, kernel_size=(1, 1), stride=(1, 1))
(proj_out): Conv2d(512, 512, kernel_size=(1, 1), stride=(1, 1))
)
(block_2): ResnetBlock(
(norm1): GroupNorm(32, 512, eps=1e-06, affine=True, bias=True)
(conv1): Conv2d(512, 512, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
(norm2): GroupNorm(32, 512, eps=1e-06, affine=True, bias=True)
(dropout): Dropout(p=0.0, inplace=False)
(conv2): Conv2d(512, 512, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
)
(up): ModuleList(
(0): Module(
(block): ModuleList(
(0): ResnetBlock(
(norm1): GroupNorm(32, 256, eps=1e-06, affine=True, bias=True)
(conv1): Conv2d(256, 128, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
(norm2): GroupNorm(32, 128, eps=1e-06, affine=True, bias=True)
(dropout): Dropout(p=0.0, inplace=False)
(conv2): Conv2d(128, 128, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
(nin_shortcut): Conv2d(256, 128, kernel_size=(1, 1), stride=(1, 1))
)
(1-2): 2 x ResnetBlock(
(norm1): GroupNorm(32, 128, eps=1e-06, affine=True, bias=True)
(conv1): Conv2d(128, 128, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
(norm2): GroupNorm(32, 128, eps=1e-06, affine=True, bias=True)
(dropout): Dropout(p=0.0, inplace=False)
(conv2): Conv2d(128, 128, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
)
(attn): ModuleList()
)
(1): Module(
(block): ModuleList(
(0): ResnetBlock(
(norm1): GroupNorm(32, 512, eps=1e-06, affine=True, bias=True)
(conv1): Conv2d(512, 256, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
(norm2): GroupNorm(32, 256, eps=1e-06, affine=True, bias=True)
(dropout): Dropout(p=0.0, inplace=False)
(conv2): Conv2d(256, 256, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
(nin_shortcut): Conv2d(512, 256, kernel_size=(1, 1), stride=(1, 1))
)
(1-2): 2 x ResnetBlock(
(norm1): GroupNorm(32, 256, eps=1e-06, affine=True, bias=True)
(conv1): Conv2d(256, 256, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
(norm2): GroupNorm(32, 256, eps=1e-06, affine=True, bias=True)
(dropout): Dropout(p=0.0, inplace=False)
(conv2): Conv2d(256, 256, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
)
(attn): ModuleList()
(upsample): Upsample(
(conv): Conv2d(256, 256, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
)
(2): Module(
(block): ModuleList(
(0-2): 3 x ResnetBlock(
(norm1): GroupNorm(32, 512, eps=1e-06, affine=True, bias=True)
(conv1): Conv2d(512, 512, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
(norm2): GroupNorm(32, 512, eps=1e-06, affine=True, bias=True)
(dropout): Dropout(p=0.0, inplace=False)
(conv2): Conv2d(512, 512, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
)
(attn): ModuleList()
(upsample): Upsample(
(conv): Conv2d(512, 512, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
)
)
(norm_out): GroupNorm(32, 128, eps=1e-06, affine=True, bias=True)
(conv_out): Conv2d(128, 4, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
(quant_conv): Conv2d(8, 8, kernel_size=(1, 1), stride=(1, 1))
(post_quant_conv): Conv2d(4, 4, kernel_size=(1, 1), stride=(1, 1))
)
(cond_stage_model): Identity()
)
)
In [23]:
# -------------------------------------------------------------
# Loading and Standardizing the 10m Sentinel-2 AoI
# -------------------------------------------------------------
input_path = work_dir+"/AoI.tif" # 4-band file: Blue, Green, Red, NIR
output_path = work_dir+"/AoI_2_5m_upscaled.tif"
# with rasterio.open(input_path) as src:
# meta = src.meta.copy()
# # Read bands 1, 2, 3, 4 corresponding to Blue, Green, Red, NIR
# img_array = src.read([1, 2, 3, 4])
S2_sel_in = ilwis.RasterCoverage (input_path)
print(S2_sel_in.size())
#print(rcSelect.size().zsize)
#load the individual spectral channels 1, 2, 3, 4 corresponding to Blue, Green, Red, NIR
B = ilwis.do('selection',S2_sel_in,"rasterbands(0)")
G = ilwis.do('selection',S2_sel_in,"rasterbands(1)")
R = ilwis.do('selection',S2_sel_in,"rasterbands(2)")
N = ilwis.do('selection',S2_sel_in,"rasterbands(3)")
#transform the spectral channels from ilwis format to a numpy array using the iterator
B_2np = np.fromiter(iter(B), np.uint32, B.size().linearSize())
B_2np = B_2np.reshape((B.size().ysize, B.size().xsize))
G_2np = np.fromiter(iter(G), np.uint32, B.size().linearSize())
G_2np = G_2np.reshape((B.size().ysize, B.size().xsize))
R_2np = np.fromiter(iter(R), np.uint32, B.size().linearSize())
R_2np = R_2np.reshape((B.size().ysize, B.size().xsize))
N_2np = np.fromiter(iter(N), np.uint32, B.size().linearSize())
N_2np = N_2np.reshape((B.size().ysize, B.size().xsize))
img_array = np.dstack((B_2np, G_2np, R_2np, N_2np))
# rasters_np[name] = img_array
# Convert list of arrays to one NumPy array
#img_array = np.stack(img_array, axis=0)
# Move bands from last axis to first
img_array = np.moveaxis(img_array, -1, 0)
print("input img_array shape:", img_array.shape)
# Transform spatial array to a PyTorch tensor scaled between 0.0 and 1.0
# The model expects dimensions shaped as: (1, Bands, Height, Width)
img_tensor = torch.from_numpy(img_array).float() / 10000.0 # Normalized from S2 L2A SR scale
img_tensor = img_tensor.unsqueeze(0).to(device)
print("Model input:", img_tensor.shape)
Size(128, 128, 4) input img_array shape: (4, 128, 128) Model input: torch.Size([1, 4, 128, 128])
In [24]:
#img_array
In [25]:
# -------------------------------------------------------------
# Executing Super-Resolution (Processing)
# -------------------------------------------------------------
print("Processing imagery from 10m down to 2.5m resolution...")
with torch.no_grad():
# Pass tensor through the model pipeline
output_tensor = model.forward(img_tensor)
# Strip batch dimension, move back to CPU, and scale back to reflectance values
output_array = output_tensor.squeeze(0).cpu().numpy()
output_array = (output_array * 10000.0).astype("uint16")
Processing imagery from 10m down to 2.5m resolution...
In [26]:
#note the new size obtained
output_array.shape
Out[26]:
(4, 512, 512)
In [27]:
# -------------------------------------------------------------
# 1. load the selected 128 lines by 128 columns original 10 m image submap
# Input order: Blue, Green, Red, NIR
# -------------------------------------------------------------
print("Original image size = ",rcSelect.size())
#load the individual spectral channels and convert to numpy array
Borg = ilwis.do('selection',rcSelect,"rasterbands(0)")
b2_2np = np.fromiter(iter(Borg), np.uint16, Borg.size().linearSize())
B2org_2np = b2_2np.reshape((Borg.size().ysize, Borg.size().xsize))
Gorg = ilwis.do('selection',rcSelect,"rasterbands(1)")
b3_2np = np.fromiter(iter(Gorg), np.uint16, Borg.size().linearSize())
B3org_2np = b3_2np.reshape((Borg.size().ysize, Borg.size().xsize))
Rorg = ilwis.do('selection',rcSelect,"rasterbands(2)")
b4_2np = np.fromiter(iter(Rorg), np.uint16, Borg.size().linearSize())
B4org_2np = b4_2np.reshape((Borg.size().ysize, Borg.size().xsize))
Norg = ilwis.do('selection',rcSelect,"rasterbands(3)")
b8_2np = np.fromiter(iter(Norg), np.uint16, Borg.size().linearSize())
B8org_2np = b8_2np.reshape((Borg.size().ysize, Borg.size().xsize))
#create a numpy 3D data stack
original_rgb = np.dstack((B4org_2np, B3org_2np, B2org_2np))#note B8org_2np not used, change accordingly for false color composite
scaling = 1500 #change the value is visualization is too light or dark
# Visualization applying the scaling factor
original_rgb_display = np.clip(
original_rgb / scaling,
0,
1
)
# -------------------------------------------------------------
# 2. Create RGB from OpenSR output
# output_array shape: (4, 512, 512) - see above
# Input order: Blue, Green, Red, NIR
# -------------------------------------------------------------
rgb = np.moveaxis(
output_array[[2, 1, 0]],
0,
-1
)
rgb_display = np.clip(
rgb / scaling,
0,
1
)
output_array_display = np.moveaxis(output_array, 0, -1)
print("SR enhanced image size = Size",output_array_display.shape)
# -------------------------------------------------------------
# 3. Plot side by side
# -------------------------------------------------------------
fig, axes = plt.subplots(
1, 2,
figsize=(14, 7)
)
# Original
axes[0].imshow(original_rgb_display)
axes[0].set_title("Original Sentinel-2 (10 m)")
axes[0].axis("off")
# OpenSR
axes[1].imshow(rgb_display)
axes[1].set_title("OpenSR Super-Resolved (2.5 m)")
axes[1].axis("off")
plt.tight_layout()
plt.show()
Original image size = Size(128, 128, 4) SR enhanced image size = Size (512, 512, 4)
In [28]:
# Use cartopy and show with geographical coordinates
# Approximate extent based on pixel centres
img_extent = (lon_ll, lon_ur, lat_ll, lat_ur)
fig, axes = plt.subplots(1, 2, figsize=(14, 7))
axes[0].imshow(
original_rgb_display,
extent=img_extent,
)
axes[1].imshow(
rgb_display,
extent=img_extent,
)
# for ax in axes:
# ax.ticklabel_format(
# axis="y",
# style="plain",
# useOffset=False
# )
# ax.tick_params(
# axis="both",
# labelsize=8
# )
for ax in axes:
ax.xaxis.set_major_formatter(FormatStrFormatter("%.4f"))
ax.yaxis.set_major_formatter(FormatStrFormatter("%.4f"))
ax.tick_params(axis="both", labelsize=8)
axes[0].set_title("Original Sentinel-2 (10 m)")
axes[1].set_title("OpenSR Super-Resolved (2.5 m)")
plt.tight_layout()
plt.show()
Save the resulting super-resolution image as geotif¶
In [29]:
#Calculate the UTM zone from the lower left coordinate provided
def get_utm_epsg(longitude, latitude):
zone = math.floor((longitude + 180) / 6) + 1
if latitude >= 0:
epsg = 32600 + zone
else:
epsg = 32700 + zone
return zone, epsg
zone, epsg = get_utm_epsg(lon_ll, lat_ll)
print(f"UTM zone: {zone}")
print(f"EPSG: {epsg}")
UTM zone: 32 EPSG: 32632
In [30]:
#get info on the original image sizes, image extent and the coordinate system used
print(S2_sel_in.size())
print(S2_sel_in.envelope())
coordSys = S2_sel_in.coordinateSystem()
coordSys.toWKT()
Size(128, 128, 4) 346890.000000 5796120.000000 348170.000000 5797400.000000
Out[30]:
'PROJCS["aoi.tif",GEOCS["aoi.tif",DATUM[" WGS 84",[DWGS84],ELLIPSOID["WGS 84",6378137.000000000000,298.257223563000],PRIMEM["Greenwich",0, AUTHORITY["EPSG",8901"]]],PROJECTION["Transverse_Mercator"],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],PARAMETER["scale",0.9996],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",9],UNIT[meter,1.0]]'
Create empty raster using the epsg zone and the extent from the image evelope information, note the resolution improvement is a factor of 4
In [31]:
size_new = S2_sel_in.size().xsize * 4
print(size_new)
grf_new = ilwis.GeoReference(
f"code=georef:type=corners, "
f"csy=epsg:{epsg}, "
f"envelope={S2_sel_in.envelope()}, "
f"gridsize={size_new} {size_new}, "
f"cornerofcorners=yes"
)
dfNum = ilwis.DataDefinition(ilwis.NumericDomain('code=value'), ilwis.NumericRange(0.0, 25000.0, 0))
rcNew = ilwis.RasterCoverage()
rcNew.setSize(ilwis.Size(size_new,size_new,S2_sel_in.size().zsize))
rcNew.setGeoReference(grf_new)
rcNew.setDataDef(dfNum)
512
In [32]:
data = np.array(output_array).flatten()
data.shape
Out[32]:
(1048576,)
In [33]:
#add the data to the raster
rcNew.array2raster(data)
print(rcNew.size())
Size(512, 512, 4)
In [34]:
#stretch the rasterbands using a loop
multiple_stretch = []
multiple_bands = ilwis.do('selection',rcNew,"rasterbands(0..3)")
ls = ilwis.do('linearstretch',multiple_bands, 1) #using an upper and lower data limit defined by the cumulative 1 and 99 % thresholds
mb_stretch = ilwis.do('setvaluerange', ls, 0, 255, 1) #set output to byte range
In [35]:
#store the results as an ILWIS maplist - display the map using the ilwis386 desktop software
mb_stretch.store('SR_out.mpl')
In [ ]:
In [ ]: