{ "cells": [ { "cell_type": "markdown", "id": "144c0f32-91c6-447e-8171-e7ee6d15be1b", "metadata": {}, "source": [ "## Super-resolve Sentinel-2 image resolution while ensuring spatial and radiometric consistency" ] }, { "cell_type": "markdown", "id": "13045c13-4d10-4230-9dd2-3020ea2a7f13", "metadata": {}, "source": [ "Notebook prepared by Ben Maathuis. ITC-University of Twente, Enschede. The Netherlands" ] }, { "cell_type": "markdown", "id": "7d6e1607-c0ca-4fc9-ae39-988eef0f8e32", "metadata": {}, "source": [ "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/" ] }, { "cell_type": "markdown", "id": "395cd765-b118-4225-abf5-e2bea35a07ac", "metadata": {}, "source": [ "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." ] }, { "cell_type": "markdown", "id": "81159c61-d960-4981-b557-852360895f00", "metadata": {}, "source": [ "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/" ] }, { "cell_type": "markdown", "id": "39493f02-4d83-456f-8b37-ba93af045161", "metadata": {}, "source": [ "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." ] }, { "cell_type": "markdown", "id": "70ee055f-5e99-43dd-bcff-018530d5117e", "metadata": {}, "source": [ "Ensure that the required packages are installed, especially torch, OmegaConf and opensr_model" ] }, { "cell_type": "markdown", "id": "81d25d96-1ccf-4de2-bd5e-ac6699365db9", "metadata": {}, "source": [ "### Installing the packages, setting the folders and connecting to the Copernicus data space ecosystem" ] }, { "cell_type": "code", "execution_count": null, "id": "ee8f95ef-78bd-4ac8-8ea6-ba430d001e8a", "metadata": {}, "outputs": [], "source": [ "#import required libraries - eventually install if not available \n", "import torch\n", "import rasterio\n", "from omegaconf import OmegaConf\n", "import opensr_model\n", "import numpy as np\n", "import os\n", "from io import StringIO\n", "import requests\n", "import matplotlib.pyplot as plt\n", "import matplotlib.ticker as mticker\n", "from matplotlib.ticker import FormatStrFormatter\n", "import cartopy\n", "from cartopy import crs as ccrs, feature as cfeature\n", "import ilwis\n", "import openeo\n", "from openeo.rest.auth.config import RefreshTokenStore\n", "import math\n", "\n", "#suppress warnings\n", "import warnings\n", "warnings.filterwarnings('ignore')\n", "warnings.simplefilter('ignore')" ] }, { "cell_type": "code", "execution_count": null, "id": "f91300e3-977c-4541-ba85-adb8fea0334b", "metadata": {}, "outputs": [], "source": [ "#within the notebook folder a directory is created to store your processing results\n", "work_dir = os.getcwd()+'/result'\n", "\n", "print(\"current dir is: %s\" % (os.getcwd()))\n", "print(\"current working directory is:\",work_dir) \n", "\n", "if os.path.isdir(work_dir):\n", " print(\"Folder exists\")\n", "else:\n", " print(\"Folder doesn't exists\")\n", " os.mkdir(work_dir)" ] }, { "cell_type": "code", "execution_count": null, "id": "f3512453-d4b1-40a6-a329-0d7e84773f75", "metadata": {}, "outputs": [], "source": [ "#set the working directory for ILWISPy\n", "ilwis.setWorkingCatalog(work_dir)\n", "print(work_dir)" ] }, { "cell_type": "code", "execution_count": null, "id": "7ec744dd-7a7d-4837-8b19-6b1d330eb074", "metadata": {}, "outputs": [], "source": [ "ilwis.version()" ] }, { "cell_type": "code", "execution_count": null, "id": "ecc691d2-31ad-45c1-834b-32836e60e3ed", "metadata": {}, "outputs": [], "source": [ "#uncomment line below to remove previous login credentials\n", "#RefreshTokenStore().remove()" ] }, { "cell_type": "code", "execution_count": null, "id": "7f08fa67-eebe-4833-841b-fad346ea79a7", "metadata": {}, "outputs": [], "source": [ "#ensure openeo is installed and you have registered\n", "connection = openeo.connect(\"openeo.dataspace.copernicus.eu\").authenticate_oidc()" ] }, { "cell_type": "markdown", "id": "bb7e5c40-fcae-4308-a2af-a36d1c1f75a8", "metadata": {}, "source": [ "### Download and visualize selected Sentinel L2A image subset " ] }, { "cell_type": "code", "execution_count": null, "id": "db0f248c-7518-490c-aa4f-6ff127c62fb7", "metadata": {}, "outputs": [], "source": [ "connection.describe_collection(\"SENTINEL2_L2A\")" ] }, { "cell_type": "markdown", "id": "5593ea9d-b08d-478c-b638-992d6ba79ddc", "metadata": {}, "source": [ "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" ] }, { "cell_type": "code", "execution_count": null, "id": "5cb1fcab-7212-4eb6-b07e-10a1386d860c", "metadata": {}, "outputs": [], "source": [ "# 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\n", "lon_ll = 6.75479\n", "lat_ll = 52.29420" ] }, { "cell_type": "code", "execution_count": null, "id": "7f595cb1-a417-434f-a993-cb6ba608b409", "metadata": {}, "outputs": [], "source": [ "#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\n", "date_start =\"2026-08-04\"\n", "date_end = date_start" ] }, { "cell_type": "code", "execution_count": null, "id": "8fdc77bf-22f8-413c-9c4c-1bc03ba9746f", "metadata": {}, "outputs": [], "source": [ "# Sentinel-2 MSI Resolution\n", "resolution = 10 # metres/pixel\n", "\n", "# Image dimensions of selected window\n", "ncols = 128\n", "nrows = 128\n", "\n", "# Approximate conversion from metres to degrees\n", "meters_per_degree_lat = 111320\n", "meters_per_degree_lon = (\n", " 111320 * np.cos(np.deg2rad(lat_ll))\n", ")\n", "\n", "dlon = resolution / meters_per_degree_lon\n", "dlat = resolution / meters_per_degree_lat\n", "\n", "# Upper-right pixel centre\n", "lon_ur = lon_ll + (ncols - 1) * dlon\n", "lat_ur = lat_ll + (nrows - 1) * dlat\n", "\n", "print(f\"Pixel size: {dlon:.8f}° lon × {dlat:.8f}° lat\")\n", "print(f\"Lower-left: ({lon_ll:.8f}, {lat_ll:.8f})\")\n", "print(f\"Upper-right: ({lon_ur:.8f}, {lat_ur:.8f})\")" ] }, { "cell_type": "code", "execution_count": null, "id": "e3967531-8833-4123-96bd-620de8fb4fc1", "metadata": {}, "outputs": [], "source": [ "#note the date \n", "t = [date_start, date_end]\n", "s2_cube = connection.load_collection(\"SENTINEL2_L2A\",\n", " spatial_extent={'west': lon_ll, 'east': lon_ur, 'south': lat_ll, 'north': lat_ur, \"crs\": \"EPSG:4326\" }, \n", " temporal_extent= t,\n", " bands=['B02', 'B03', 'B04', 'B08'],\n", " max_cloud_cover=100,\n", ")" ] }, { "cell_type": "code", "execution_count": null, "id": "8896b864-a8c1-4c35-ad1c-298d9745bd1e", "metadata": {}, "outputs": [], "source": [ "s2_cube.download(work_dir+\"/S2_selected.tif\")" ] }, { "cell_type": "code", "execution_count": null, "id": "a2145e03-ec90-48c0-9de4-7d86dfe3a5b6", "metadata": {}, "outputs": [], "source": [ "#read the S2 image in ilwispy and note the number of spectral bands and image size\n", "S2_in = ilwis.RasterCoverage ('S2_selected.tif')\n", "print(S2_in.size())" ] }, { "cell_type": "code", "execution_count": null, "id": "a344126f-aca4-4ec3-83e6-43cb558fe4f1", "metadata": {}, "outputs": [], "source": [ "#stretch the rasterbands using a loop\n", "multiple_stretch = []\n", "multiple_bands = ilwis.do('selection',S2_in,\"rasterbands(0..3)\") \n", "ls = ilwis.do('linearstretch',multiple_bands, 1) #using an upper and lower data limit defined by the cumulative 1 and 99 % thresholds\n", "mb_stretch = ilwis.do('setvaluerange', ls, 0, 255, 1) #set ouput to byte range" ] }, { "cell_type": "code", "execution_count": null, "id": "e8b337d6-97af-49cb-b8dc-d1cbffc16011", "metadata": {}, "outputs": [], "source": [ "#store the results as an ILWIS maplist - display the map using the ilwis386 desktop software\n", "mb_stretch.store('S2_selected_stretch.mpl')" ] }, { "cell_type": "code", "execution_count": null, "id": "816cad61-84b0-4ace-820e-d0793562659b", "metadata": {}, "outputs": [], "source": [ "#load the individual spectral channels\n", "Blues = ilwis.do('selection',mb_stretch,\"rasterbands(0)\")\n", "Greens = ilwis.do('selection',mb_stretch,\"rasterbands(1)\")\n", "Reds = ilwis.do('selection',mb_stretch,\"rasterbands(2)\")" ] }, { "cell_type": "code", "execution_count": null, "id": "c13cef55-7300-4b44-b732-d5eb6b574204", "metadata": {}, "outputs": [], "source": [ "#transform the spectral channels from ilwis format to a numpy array using the iterator\n", "Blues_2np = np.fromiter(iter(Blues), np.ubyte, Blues.size().linearSize()) \n", "Blues_2np = Blues_2np.reshape((Blues.size().ysize, Blues.size().xsize))\n", "\n", "Greens_2np = np.fromiter(iter(Greens), np.ubyte, Blues.size().linearSize()) \n", "Greens_2np = Greens_2np.reshape((Blues.size().ysize, Blues.size().xsize))\n", "\n", "Reds_2np = np.fromiter(iter(Reds), np.ubyte, Blues.size().linearSize()) \n", "Reds_2np = Reds_2np.reshape((Blues.size().ysize, Blues.size().xsize))" ] }, { "cell_type": "code", "execution_count": null, "id": "f533c5ad-3c41-4eb3-a449-33f140d81564", "metadata": {}, "outputs": [], "source": [ "#create a numpy natural color 3D data stack\n", "ncol = np.dstack((Reds_2np, Greens_2np, Blues_2np))" ] }, { "cell_type": "markdown", "id": "3ea54ac9-bedd-40ab-b338-6f21158a9c5a", "metadata": {}, "source": [ "Check the image obtained" ] }, { "cell_type": "code", "execution_count": null, "id": "31e93451-dde6-4a35-ac06-2a5aeaec175c", "metadata": {}, "outputs": [], "source": [ "img_extent = (lon_ll, lon_ur, lat_ll, lat_ur)\n", "\n", "fig = plt.figure(figsize=(10, 14))\n", "ax = plt.axes(projection=ccrs.PlateCarree())\n", "\n", "plt.title('Selected Sentinel-2 composite of Area of Interest of ' + str(date_start))\n", "\n", "# Data raster\n", "ax.imshow(\n", " ncol,\n", " extent=img_extent,\n", " transform=ccrs.PlateCarree()\n", ")\n", "\n", "# Coordinate grid\n", "gl = ax.gridlines(\n", " crs=ccrs.PlateCarree(),\n", " draw_labels=True,\n", " linewidth=0.6,\n", " color=\"gray\",\n", " alpha=0.7,\n", " linestyle=\"--\",\n", " zorder=5\n", ")\n", "\n", "# Grid spacing\n", "gl.xlocator = mticker.MultipleLocator(0.005)\n", "gl.ylocator = mticker.MultipleLocator(0.005)\n", "\n", "# Labels\n", "gl.top_labels = False\n", "gl.right_labels = False\n", "\n", "gl.xlabel_style = {\n", " \"size\": 10,\n", "}\n", "gl.ylabel_style = {\n", " \"size\": 10,\n", "}\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "3219821b-8c7c-42a2-9711-07faa78ba382", "metadata": {}, "source": [ "Ensure the submap has 128 lines and 128 columns for further processing" ] }, { "cell_type": "code", "execution_count": null, "id": "3613cf5b-70c1-4c5b-9e65-ab3b6d7b2cbd", "metadata": {}, "outputs": [], "source": [ "# sub map creation to ensure that the final AoI has 128 lines by 128 columns\n", "rcSelect = ilwis.do('selection',S2_in,'boundingbox(1 1, 128 128)')\n", "print(rcSelect.size().xsize)\n", "print(rcSelect.size().ysize)\n", "print(rcSelect.size().zsize)" ] }, { "cell_type": "markdown", "id": "7272b74b-7451-4a6a-aa41-a3759be5a5db", "metadata": {}, "source": [ "Save the results, as a tiff file" ] }, { "cell_type": "code", "execution_count": null, "id": "456c3bcd-9777-4ddc-a431-0b101e566829", "metadata": {}, "outputs": [], "source": [ "rcSelect.store(\"AoI.tif\", \"GTiff\", \"gdal\")" ] }, { "cell_type": "markdown", "id": "ace4ad1a-965a-46b7-a5c7-6be9b514b3c8", "metadata": {}, "source": [ "### Start the Spatial Super Resolution Processing procedure" ] }, { "cell_type": "code", "execution_count": null, "id": "0dd666fb-8190-4f64-82c8-4dfba4e02ade", "metadata": {}, "outputs": [], "source": [ "# -------------------------------------------------------------\n", "# 1. Environment Setup & Configuration Setup\n", "# -------------------------------------------------------------\n", "# Determine device acceleration (GPU is heavily recommended for Latent Diffusion)\n", "device = \"cuda\" if torch.cuda.is_available() else \"cpu\"\n", "print(f\"Using processing engine: {device.upper()}\")\n", "\n", "# Fetch the official neural network configuration mapping from ESA OpenSR\n", "config_url = (\n", " \"https://raw.githubusercontent.com/\"\n", " \"ESAOpenSR/opensr-model/refs/heads/main/\"\n", " \"opensr_model/configs/config_10m.yaml\"\n", ")\n", "response = requests.get(config_url)\n", "config = OmegaConf.load(StringIO(response.text))\n", "\n", "# Instantiate the model and load its pretrained 4x upscaling weights\n", "print(\"Downloading and preparing the Latent Diffusion model...\")\n", "model = opensr_model.SRLatentDiffusion(config, device=device)\n", "model.load_pretrained(config.ckpt_version)\n", "model.eval() # Put neural network into evaluation mode" ] }, { "cell_type": "code", "execution_count": null, "id": "71e0f4d8-b98e-4e74-9e2b-fbcdc7d8a035", "metadata": {}, "outputs": [], "source": [ "# -------------------------------------------------------------\n", "# Loading and Standardizing the 10m Sentinel-2 AoI\n", "# -------------------------------------------------------------\n", "input_path = work_dir+\"/AoI.tif\" # 4-band file: Blue, Green, Red, NIR\n", "output_path = work_dir+\"/AoI_2_5m_upscaled.tif\"\n", "\n", "# with rasterio.open(input_path) as src:\n", "# meta = src.meta.copy()\n", "# # Read bands 1, 2, 3, 4 corresponding to Blue, Green, Red, NIR\n", "# img_array = src.read([1, 2, 3, 4]) \n", "\n", "\n", "S2_sel_in = ilwis.RasterCoverage (input_path)\n", "print(S2_sel_in.size())\n", "#print(rcSelect.size().zsize)\n", "\n", "#load the individual spectral channels 1, 2, 3, 4 corresponding to Blue, Green, Red, NIR\n", "B = ilwis.do('selection',S2_sel_in,\"rasterbands(0)\")\n", "G = ilwis.do('selection',S2_sel_in,\"rasterbands(1)\")\n", "R = ilwis.do('selection',S2_sel_in,\"rasterbands(2)\")\n", "N = ilwis.do('selection',S2_sel_in,\"rasterbands(3)\")\n", "\n", "#transform the spectral channels from ilwis format to a numpy array using the iterator\n", "B_2np = np.fromiter(iter(B), np.uint32, B.size().linearSize()) \n", "B_2np = B_2np.reshape((B.size().ysize, B.size().xsize))\n", "\n", "G_2np = np.fromiter(iter(G), np.uint32, B.size().linearSize()) \n", "G_2np = G_2np.reshape((B.size().ysize, B.size().xsize))\n", "\n", "R_2np = np.fromiter(iter(R), np.uint32, B.size().linearSize()) \n", "R_2np = R_2np.reshape((B.size().ysize, B.size().xsize))\n", "\n", "N_2np = np.fromiter(iter(N), np.uint32, B.size().linearSize()) \n", "N_2np = N_2np.reshape((B.size().ysize, B.size().xsize))\n", "\n", "img_array = np.dstack((B_2np, G_2np, R_2np, N_2np))\n", "# rasters_np[name] = img_array\n", "# Convert list of arrays to one NumPy array\n", "#img_array = np.stack(img_array, axis=0)\n", "\n", "# Move bands from last axis to first\n", "img_array = np.moveaxis(img_array, -1, 0)\n", "\n", "print(\"input img_array shape:\", img_array.shape)\n", "\n", "# Transform spatial array to a PyTorch tensor scaled between 0.0 and 1.0\n", "# The model expects dimensions shaped as: (1, Bands, Height, Width)\n", "img_tensor = torch.from_numpy(img_array).float() / 10000.0 # Normalized from S2 L2A SR scale\n", "img_tensor = img_tensor.unsqueeze(0).to(device)\n", "\n", "print(\"Model input:\", img_tensor.shape)" ] }, { "cell_type": "code", "execution_count": null, "id": "a508fb08-cd37-4929-b838-303cbcaef950", "metadata": {}, "outputs": [], "source": [ "#img_array" ] }, { "cell_type": "code", "execution_count": null, "id": "bb0573e3-1393-42ed-b966-139b50b2d01d", "metadata": {}, "outputs": [], "source": [ "# -------------------------------------------------------------\n", "# Executing Super-Resolution (Processing)\n", "# -------------------------------------------------------------\n", "print(\"Processing imagery from 10m down to 2.5m resolution...\")\n", "with torch.no_grad():\n", " # Pass tensor through the model pipeline\n", " output_tensor = model.forward(img_tensor)\n", " \n", " # Strip batch dimension, move back to CPU, and scale back to reflectance values\n", " output_array = output_tensor.squeeze(0).cpu().numpy()\n", " output_array = (output_array * 10000.0).astype(\"uint16\")" ] }, { "cell_type": "code", "execution_count": null, "id": "65910b96-ab5f-421d-8fd0-91e09aaa788c", "metadata": {}, "outputs": [], "source": [ "#note the new size obtained\n", "output_array.shape" ] }, { "cell_type": "code", "execution_count": null, "id": "71c816c0-25da-410b-8897-8750218c268b", "metadata": {}, "outputs": [], "source": [ "# -------------------------------------------------------------\n", "# 1. load the selected 128 lines by 128 columns original 10 m image submap\n", "# Input order: Blue, Green, Red, NIR\n", "# -------------------------------------------------------------\n", "print(\"Original image size = \",rcSelect.size())\n", "\n", "#load the individual spectral channels and convert to numpy array\n", "Borg = ilwis.do('selection',rcSelect,\"rasterbands(0)\")\n", "b2_2np = np.fromiter(iter(Borg), np.uint16, Borg.size().linearSize()) \n", "B2org_2np = b2_2np.reshape((Borg.size().ysize, Borg.size().xsize))\n", "\n", "Gorg = ilwis.do('selection',rcSelect,\"rasterbands(1)\")\n", "b3_2np = np.fromiter(iter(Gorg), np.uint16, Borg.size().linearSize()) \n", "B3org_2np = b3_2np.reshape((Borg.size().ysize, Borg.size().xsize))\n", "\n", "Rorg = ilwis.do('selection',rcSelect,\"rasterbands(2)\")\n", "b4_2np = np.fromiter(iter(Rorg), np.uint16, Borg.size().linearSize()) \n", "B4org_2np = b4_2np.reshape((Borg.size().ysize, Borg.size().xsize))\n", "\n", "Norg = ilwis.do('selection',rcSelect,\"rasterbands(3)\")\n", "b8_2np = np.fromiter(iter(Norg), np.uint16, Borg.size().linearSize()) \n", "B8org_2np = b8_2np.reshape((Borg.size().ysize, Borg.size().xsize))\n", "\n", "#create a numpy 3D data stack\n", "original_rgb = np.dstack((B4org_2np, B3org_2np, B2org_2np))#note B8org_2np not used, change accordingly for false color composite\n", "\n", "scaling = 1500 #change the value is visualization is too light or dark\n", "\n", "# Visualization applying the scaling factor\n", "original_rgb_display = np.clip(\n", " original_rgb / scaling,\n", " 0,\n", " 1\n", ")\n", "\n", "# -------------------------------------------------------------\n", "# 2. Create RGB from OpenSR output\n", "# output_array shape: (4, 512, 512) - see above\n", "# Input order: Blue, Green, Red, NIR\n", "# -------------------------------------------------------------\n", "rgb = np.moveaxis(\n", " output_array[[2, 1, 0]],\n", " 0,\n", " -1\n", ")\n", "\n", "rgb_display = np.clip(\n", " rgb / scaling,\n", " 0,\n", " 1\n", ")\n", "\n", "output_array_display = np.moveaxis(output_array, 0, -1)\n", "print(\"SR enhanced image size = Size\",output_array_display.shape)\n", "# -------------------------------------------------------------\n", "# 3. Plot side by side\n", "# -------------------------------------------------------------\n", "fig, axes = plt.subplots(\n", " 1, 2,\n", " figsize=(14, 7)\n", ")\n", "\n", "# Original\n", "axes[0].imshow(original_rgb_display)\n", "axes[0].set_title(\"Original Sentinel-2 (10 m)\")\n", "axes[0].axis(\"off\")\n", "\n", "# OpenSR\n", "axes[1].imshow(rgb_display)\n", "axes[1].set_title(\"OpenSR Super-Resolved (2.5 m)\")\n", "axes[1].axis(\"off\")\n", "\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "id": "b0f4287e-f637-4480-bf6f-51bb116273e0", "metadata": {}, "outputs": [], "source": [ "# Use cartopy and show with geographical coordinates\n", "# Approximate extent based on pixel centres\n", "img_extent = (lon_ll, lon_ur, lat_ll, lat_ur)\n", "fig, axes = plt.subplots(1, 2, figsize=(14, 7))\n", "\n", "axes[0].imshow(\n", " original_rgb_display,\n", " extent=img_extent,\n", ")\n", "\n", "axes[1].imshow(\n", " rgb_display,\n", " extent=img_extent,\n", ")\n", "\n", "# for ax in axes:\n", "# ax.ticklabel_format(\n", "# axis=\"y\",\n", "# style=\"plain\",\n", "# useOffset=False\n", "# )\n", "# ax.tick_params(\n", "# axis=\"both\",\n", "# labelsize=8\n", "# )\n", "\n", "\n", "\n", "for ax in axes:\n", " ax.xaxis.set_major_formatter(FormatStrFormatter(\"%.4f\"))\n", " ax.yaxis.set_major_formatter(FormatStrFormatter(\"%.4f\"))\n", " ax.tick_params(axis=\"both\", labelsize=8)\n", "\n", "axes[0].set_title(\"Original Sentinel-2 (10 m)\")\n", "axes[1].set_title(\"OpenSR Super-Resolved (2.5 m)\")\n", "\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "eb9289e4-4323-446d-bfaa-e7df465a8861", "metadata": {}, "source": [ "#### Save the resulting super-resolution image as geotif" ] }, { "cell_type": "code", "execution_count": null, "id": "1789dece-0a3e-4b3a-8659-15c8a76b0469", "metadata": {}, "outputs": [], "source": [ "#Calculate the UTM zone from the lower left coordinate provided\n", "\n", "def get_utm_epsg(longitude, latitude):\n", " zone = math.floor((longitude + 180) / 6) + 1\n", "\n", " if latitude >= 0:\n", " epsg = 32600 + zone\n", " else:\n", " epsg = 32700 + zone\n", "\n", " return zone, epsg\n", "\n", "\n", "zone, epsg = get_utm_epsg(lon_ll, lat_ll)\n", "\n", "print(f\"UTM zone: {zone}\")\n", "print(f\"EPSG: {epsg}\")" ] }, { "cell_type": "code", "execution_count": null, "id": "212cf493-4224-40b2-a4da-bfb9922d941d", "metadata": {}, "outputs": [], "source": [ "#get info on the original image sizes, image extent and the coordinate system used\n", "print(S2_sel_in.size())\n", "print(S2_sel_in.envelope())\n", "coordSys = S2_sel_in.coordinateSystem()\n", "coordSys.toWKT()" ] }, { "cell_type": "markdown", "id": "559520b8-721a-4d72-ba96-95000195cefc", "metadata": {}, "source": [ "Create empty raster using the epsg zone and the extent from the image evelope information, note the resolution improvement is a factor of 4 " ] }, { "cell_type": "code", "execution_count": null, "id": "1f70e7ac-05fa-4c83-8f71-5888931ebc2c", "metadata": {}, "outputs": [], "source": [ "size_new = S2_sel_in.size().xsize * 4\n", "print(size_new)\n", "\n", "\n", "grf_new = ilwis.GeoReference(\n", " f\"code=georef:type=corners, \"\n", " f\"csy=epsg:{epsg}, \"\n", " f\"envelope={S2_sel_in.envelope()}, \"\n", " f\"gridsize={size_new} {size_new}, \"\n", " f\"cornerofcorners=yes\"\n", ")\n", "dfNum = ilwis.DataDefinition(ilwis.NumericDomain('code=value'), ilwis.NumericRange(0.0, 25000.0, 0))\n", "rcNew = ilwis.RasterCoverage()\n", "rcNew.setSize(ilwis.Size(size_new,size_new,S2_sel_in.size().zsize))\n", "rcNew.setGeoReference(grf_new)\n", "rcNew.setDataDef(dfNum)" ] }, { "cell_type": "code", "execution_count": null, "id": "35e5054b-55c6-49fb-b4b1-a35a697fe815", "metadata": {}, "outputs": [], "source": [ "data = np.array(output_array).flatten()\n", "data.shape" ] }, { "cell_type": "code", "execution_count": null, "id": "4d6084bc-d166-4115-9a2c-c1417e764da1", "metadata": {}, "outputs": [], "source": [ "#add the data to the raster\n", "rcNew.array2raster(data)\n", "print(rcNew.size())" ] }, { "cell_type": "code", "execution_count": null, "id": "0482101c-f669-4e39-b225-66649e79fe29", "metadata": {}, "outputs": [], "source": [ "#stretch the rasterbands using a loop\n", "multiple_stretch = []\n", "multiple_bands = ilwis.do('selection',rcNew,\"rasterbands(0..3)\") \n", "ls = ilwis.do('linearstretch',multiple_bands, 1) #using an upper and lower data limit defined by the cumulative 1 and 99 % thresholds\n", "mb_stretch = ilwis.do('setvaluerange', ls, 0, 255, 1) #set output to byte range" ] }, { "cell_type": "code", "execution_count": null, "id": "83a19173-5074-45f6-a15f-7e04758d1a15", "metadata": {}, "outputs": [], "source": [ "#store the results as an ILWIS maplist - display the map using the ilwis386 desktop software\n", "mb_stretch.store('SR_out.mpl')" ] }, { "cell_type": "code", "execution_count": null, "id": "05563d23-51eb-4d6b-ac6b-d9e85b3b08cc", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "id": "6ba519fa-f48b-477e-a4f5-c54db4940501", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.13.3" } }, "nbformat": 4, "nbformat_minor": 5 }