Oil spill detection using Sentinel 1 and Sentinel 2, example Al Qibiliya - Oman

Notebook prepared by Ben Maathuis. ITC-University of Twente, Enschede. The Netherlands

Microwave sensing:

• Satellites like Copernicus Sentinel-1 use Synthetic Aperture Radar to measure sea surface roughness. Floating oil dampens small wind-driven waves, causing the oiled zone to show up as a distinct dark smudge or streak compared to rougher, reflective clean water.

Multispectral sensing:

• in the Visible Light Bands (RGB): Blue, green, and red wavelengths detect high-contrast oil patches. Thick oil looks brown or black, while thin sheens create a silvery sunglint anomaly against the water.
• In the Near-Infrared (NIR): Water strongly absorbs NIR light, appearing almost black. Oil reflects NIR highly, making the spill stand out as a bright, high-contrast shape against the dark sea.
• In the Shortwave Infrared (SWIR): SWIR bands pinpoint specific chemical absorption features unique to hydrocarbons. This allows analysts to distinguish true crude oil from organic matter like algae or sargassum weed.
• In the Thermal Infrared (TIR): Oil absorbs solar radiation faster than water during the day, making thick oil patches appear much warmer than the surrounding sea in thermal images. At night, the trend reverses. Note that the spectral information from this part of the electromagntic spectrum is not processed within this notebook.

Review:

• https://blog.response.restoration.noaa.gov/how-thick-oil-slick
• https://www.sciencedirect.com/science/article/pii/S0034425719304407?via%3Dihub
• https://www.sciencedirect.com/science/article/pii/S2215016121001205

For additional background information on this oil spill incident see:

• https://www.nhregister.com/news/world/article/oil-spill-from-grounded-tanker-off-oman-expands-22373858.php
• https://www.khaleejtimes.com/world/gulf/oman-oil-spill-al-hallaniyat-islands-ship-grounds

Note: the stranded ship location is directly south of the south-west extension of the Al Qibiliya island in the Arabian Sea in front of the coast of South Oman. The prevailing wind direction at the surface during these days is from the south west to the north east and according the newspaper information additional oil spill exension was detected on images from 05 to 07 August 2026.

Within this Notebook your are going to use ILWISPy in conjunction with a number of common used Python libraries (like Numpy and Matplotlib) for data processing as well as OpenEO for initial data retrieval. Use is made of Sentinel1 and Sentinel2, first the images form both satellites are visualized seperately, then the data is merged / fused and a unbsupervised classification is subsequently performed to distinguish the identified features.

Download and pre-processing data

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 images check the Copernicus Browser, available at: https://browser.dataspace.copernicus.eu/

#import required libraries - eventually install if not available 
import os
import ilwis
import openeo
from openeo.rest.auth.config import RefreshTokenStore
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import matplotlib.ticker as mticker
import cartopy
from cartopy import crs as ccrs, feature as cfeature
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from sklearn.cluster import KMeans
from matplotlib.colors import BoundaryNorm

#suppress warnings
import warnings
warnings.filterwarnings('ignore')
warnings.simplefilter('ignore')

#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\oilspill_detection
current working directory is: d:\jupyter\notebook_scripts\Special\oilspill_detection/result
Folder exists

#set the working directory for ILWISPy
ilwis.setWorkingCatalog(work_dir)
print(work_dir)

d:\jupyter\notebook_scripts\Special\oilspill_detection/result

ilwis.version()

'1.0 build 20260708'

#remove previous login credentials, uncomment if running the notebook for the first time
#RefreshTokenStore().remove()

#ensure openeo is installed and you have registered
connection = openeo.connect("openeo.dataspace.copernicus.eu").authenticate_oidc()

Authenticated using refresh token.

Indentify the shipwreck location

#check the meta data information
connection.describe_collection("SENTINEL2_L2A")

Collect a Sentinel 2 image before the ship wreck event

#note the date 
t = ["2026-06-18", "2026-06-18"]
s2_cube = connection.load_collection("SENTINEL2_L2A",
    spatial_extent={'west': 56.31, 'east': 56.355, 'south': 17.485, 'north': 17.515,"crs": "EPSG:4326"},
    temporal_extent= t,
    bands=['B02', 'B03', 'B04', 'B08'],
    max_cloud_cover=100,
)

s2_cube.download(work_dir+"/s2_background.tif")

#read the S2 image in ilwispy
s2_background = ilwis.RasterCoverage ('s2_background.tif')
print(s2_background.size())

Size(480, 335, 4)

The initial Sentinel-2 data is in int16 format (data range from -32768 to 32767). To convert to reflectance the scaling factor (0.0001) should be used. Here the initial data is stretched to ubyte (data range from 0 to 255), as we are using it for visualization!

#stretch the rasterbands contained in s2_background using a loop
multiple_stretch = []
multiple_bands = ilwis.do('selection',s2_background,"rasterbands(0..2)") 
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

#store the results as an ILWIS maplist - display the map using the ilwis386 desktop software
mb_stretch.store('S2_background.mpl')

#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)")

#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))

#create a numpy 3D data stack
ncol = np.dstack((Reds_2np, Greens_2np, Blues_2np))

img_extent = (56.31, 56.355, 17.485, 17.515)

fig = plt.figure(figsize=(12, 10))
ax = plt.axes(projection=ccrs.PlateCarree())

plt.title('Sentinel-2 composite of 2026-06-18 of Al Qibiliya - Arabian Sea, Oman')

# Data raster
ax.imshow(
    ncol,
    origin='upper',
    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.01)
gl.ylocator = mticker.MultipleLocator(0.01)

# Labels
gl.top_labels = False
gl.right_labels = False

gl.xlabel_style = {
    "size": 10,
}
gl.ylabel_style = {
    "size": 10,
}

plt.show()

Process the same area of interest after ship stranded on the southwestern coast of the island

#note the date and spectral channels
t = ["2026-07-28", "2026-07-28"]
s2_cube = connection.load_collection("SENTINEL2_L2A",
    spatial_extent={'west': 56.31, 'east': 56.355, 'south': 17.485, 'north': 17.515,"crs": "EPSG:4326"},
    temporal_extent= t,
    bands=['B02', 'B03', 'B04', 'B08', 'B11', 'B12'],
    max_cloud_cover=100,
)

s2_cube.download(work_dir+"/s2_ship_stranded.tif")

#read the S2 image in ilwispy
s2_ship = ilwis.RasterCoverage ('s2_ship_stranded.tif')
print(s2_ship.size())

Size(480, 335, 6)

#stretch the rasterbands using a loop
multiple_stretch = []
multiple_bands = ilwis.do('selection',s2_ship,"rasterbands(0..5)") 
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

#store the results as an ILWIS maplist - display the map using the ilwis386 desktop software
mb_stretch.store('S2_ship_stranded.mpl')

#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)")
IRs = ilwis.do('selection',mb_stretch,"rasterbands(3)")
SWIR11s = ilwis.do('selection',mb_stretch,"rasterbands(4)")
SWIR12s = ilwis.do('selection',mb_stretch,"rasterbands(5)")

#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))

#create a numpy 3D data stack
ncol = np.dstack((Reds_2np, Greens_2np, Blues_2np))

img_extent = (56.31, 56.355, 17.485, 17.515)

fig = plt.figure(figsize=(12, 10))
ax = plt.axes(projection=ccrs.PlateCarree())

plt.title('Sentinel-2 natural color composite of 2026-07-28 showing location of shipwreck\n'
          'and initial oil spill features, southwest of Al Qibiliya - Arabian Sea, Oman')

# Data raster
ax.imshow(
    ncol,
    origin='upper',
    extent=img_extent,
    transform=ccrs.PlateCarree()
)

# Shipwreck location
lon = 56.3213
lat = 17.4938

circle = patches.Circle(
    (lon, lat),
    radius=0.0015,          #in degrees
    fill=False,
    edgecolor='red',
    linewidth=2,
    transform=ccrs.PlateCarree(),
    zorder=5
)

ax.add_patch(circle)

# 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.01)
gl.ylocator = mticker.MultipleLocator(0.01)

# Labels
gl.top_labels = False
gl.right_labels = False

gl.xlabel_style = {
    "size": 10,
}
gl.ylabel_style = {
    "size": 10,
}

plt.show()

Enhancement of the oil spill features, band transforms after Rajendran et al. and creating color composite. The following band ratios are used, note the Sentinel-2 spectral channel details from the meta data given above:

• Red: B3/B2
• Green: (B11+B12)/B8
• Blue: (B3+B4)/B2

#The band ratios are derived using the original data, so once more load the original individual Sentinel-2 spectral channels
B2_org = ilwis.do('selection',multiple_bands,"rasterbands(0)")
B3_org = ilwis.do('selection',multiple_bands,"rasterbands(1)")
B4_org = ilwis.do('selection',multiple_bands,"rasterbands(2)")
B8_org = ilwis.do('selection',multiple_bands,"rasterbands(3)")
SWIR11_org = ilwis.do('selection',multiple_bands,"rasterbands(4)")
SWIR12_org = ilwis.do('selection',multiple_bands,"rasterbands(5)")

#calculate the band ratios
Red = ilwis.do('mapcalc','(@2/@1)', B2_org,  B3_org)
Green = ilwis.do('mapcalc','(@2+@3)/@1', B2_org, B3_org, B4_org)
Blue = ilwis.do('mapcalc','(@2+@3)/@1', B8_org, SWIR11_org, SWIR12_org)

#perform the image stretch and transform to numpy array
stat_Red = Red.statistics(ilwis.PropertySets.pHISTOGRAM, 65535)
minPerc1, maxPerc1 = stat_Red.calcStretchRange(1)
Reds = ilwis.do('linearstretch',Red, minPerc1, maxPerc1)
Reds = ilwis.do('setvaluerange', Reds, 0, 255, 1)

R_2np = np.fromiter(iter(Reds), np.ubyte, Red.size().linearSize()) 
R_2np = R_2np.reshape((Red.size().ysize, Red.size().xsize))


stat_Green = Green.statistics(ilwis.PropertySets.pHISTOGRAM, 65535)
minPerc1, maxPerc1 = stat_Green.calcStretchRange(1)
Greens = ilwis.do('linearstretch',Green, minPerc1, maxPerc1)
Greens = ilwis.do('setvaluerange', Greens, 0, 255, 1)

G_2np = np.fromiter(iter(Greens), np.ubyte, Red.size().linearSize()) 
G_2np = G_2np.reshape((Red.size().ysize, Red.size().xsize))

stat_Blue = Blue.statistics(ilwis.PropertySets.pHISTOGRAM, 65535)
minPerc1, maxPerc1 = stat_Blue.calcStretchRange(1)
Blues = ilwis.do('linearstretch',Blue, minPerc1, maxPerc1)
Blues = ilwis.do('setvaluerange', Blues, 0, 255, 1)

B_2np = np.fromiter(iter(Blues), np.ubyte, Red.size().linearSize()) 
B_2np = B_2np.reshape((Red.size().ysize, Red.size().xsize))

# create color composite (in RGB)
RGB_S2 = np.dstack((R_2np, G_2np, B_2np))

img_extent = (56.31, 56.355, 17.485, 17.515)

fig = plt.figure(figsize=(12, 10))

ax = plt.axes(projection=ccrs.PlateCarree())

# Data raster
ax.imshow(
    RGB_S2,
    origin='upper',
    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.01)
gl.ylocator = mticker.MultipleLocator(0.01)

# Labels
gl.top_labels = False
gl.right_labels = False

gl.xlabel_style = {
    "size": 10,
}
gl.ylabel_style = {
    "size": 10,
}

plt.title('RGB composite of 20260728 using band transforms for oil spill enhancement')

plt.show()

Retrieve Sentinel 1 GRD image from the oil spill event on 20260726

connection.describe_collection("SENTINEL1_GRD")

#note the date first time step
t = ["2026-07-26", "2026-07-26"]
s1_cube = connection.load_collection(
  "SENTINEL1_GRD",
  spatial_extent={'west': 56.30, 'east': 56.90, 'south': 17.46, 'north': 17.90,"crs": "EPSG:4326"},
  temporal_extent=t,
  bands=["VV", "VH"]
)

s1_cube.download(work_dir+'/S1_spill_20260726.tif')

Create radar RGB composite

• VV - red
• VH - green
• VC - blue (cross polarized)

s1 = ilwis.RasterCoverage('S1_spill_20260726.tif')
print(s1.size())
print(s1.envelope())
coordSys = s1.coordinateSystem()
coordSys.toWKT()

Size(6374, 4883, 2)
425670.000000 1930440.000000 489410.000000 1979270.000000

'PROJCS["s1_spill_20260726.tif",GEOCS["s1_spill_20260726.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",57],UNIT[meter,1.0]]'

VV = ilwis.do('selection',s1,"rasterbands(0)")

#create histogram using large number of bins, given the fact that the input image is a float64, with both positive and negative numbers
hist = VV.statistics(ilwis.PropertySets.pHISTOGRAM, 65535)
minPerc1, maxPerc1 = hist.calcStretchRange(1)
VVs = ilwis.do('linearstretch',VV, minPerc1, maxPerc1)
VVs = ilwis.do('mapcalc','iff(@1==?,255,@1)', VVs) #modify pixels over water
VVs = ilwis.do('setvaluerange', VVs, 0, 255, 1)
#uncomment lines below if you want to store the stretched image
#VVs.store('VVs.mpr')

VV_2np = np.fromiter(iter(VVs), np.ubyte, VVs.size().linearSize()) 
VV_2np = VV_2np.reshape((VVs.size().ysize, VVs.size().xsize))

fig = plt.figure(figsize=(15, 10))

plt.imshow(VV_2np, vmin=0, vmax=255, cmap= "turbo")
plt.axis("off")
plt.title('VV - S1 image from 20260726 - retrieved from OpenEO service provider');

VH = ilwis.do('selection',s1,"rasterbands(1)")

#stretch the VH band
hist = VH.statistics(ilwis.PropertySets.pHISTOGRAM, 65535)
minPerc1, maxPerc1 = hist.calcStretchRange(1)
VHs = ilwis.do('linearstretch',VH, minPerc1, maxPerc1)
VHs = ilwis.do('mapcalc','iff(@1==?,0,@1)', VHs) #remove bad pixels over water
VHs = ilwis.do('setvaluerange', VHs, 0, 255, 1)
#uncomment lines below if you want to store the stretched image
#VHs.store('VHs.mpr')

VH_2np = np.fromiter(iter(VHs), np.ubyte, VHs.size().linearSize()) 
VH_2np = VH_2np.reshape((VHs.size().ysize, VHs.size().xsize))

#calculate the cross polarized image (VV/VH)
VC = ilwis.do('mapcalc','(@1/@2)', VV, VH)
#uncomment lines below if you want to store the S1 cross polazization image
#VC.store('VC.mpr')

#create histogram using large number of bins, given the fact that the input image is a float64, with both positive and negative numbers
hist = VC.statistics(ilwis.PropertySets.pHISTOGRAM, 65535)
minPerc1, maxPerc1 = hist.calcStretchRange(1)
VCs = ilwis.do('linearstretch',VC, minPerc1, maxPerc1)
VCs = ilwis.do('mapcalc','iff(@1==?,255,@1)', VCs) #modify pixels over water
VCs = ilwis.do('setvaluerange', VCs, 0, 255, 1)
#uncomment lines below if you want to store the stretched image
#VCs.store('VCs.mpr')

VC_2np = np.fromiter(iter(VCs), np.ubyte, VCs.size().linearSize()) 
VC_2np = VC_2np.reshape((VCs.size().ysize, VCs.size().xsize))

# create  color composite (in RGB)
cc = np.dstack((VV_2np, VH_2np, VC_2np))

#band interleaved
S1ALL = np.array([VC_2np, VH_2np, VV_2np])

fig, axes = plt.subplots(1, 2, figsize=(14, 7))

# Full image
axes[0].imshow(cc)
axes[0].axis("off")
axes[0].set_title("S1 polarization composite")

# South-west 1000 x 1000 pixels
sw = cc[-1000:, :1000]

axes[1].imshow(sw)
axes[1].axis("off")
axes[1].set_title("South-west image corner enlargement")

plt.tight_layout()
plt.show()

#create empty ilwis raster 
S1_ilw = ilwis.RasterCoverage()
defNumr = ilwis.DataDefinition(ilwis.NumericDomain('code=value'), ilwis.NumericRange(0, 255, 1))
S1_ilw.setDataDef(defNumr)
S1_ilw.setSize(ilwis.Size(745, 446, 3)) 
S1_ilw.setGeoReference(s1.geoReference())

#read the numpy array
S1_ilw.array2raster(S1ALL.flatten())

#comment if you don't want to store the image in an ILWIS maplist format
S1_ilw.store('S1_spill_20260726.mpl')

#uncomment if you want to store the image in geotif format
#S1_ilw.store("S1_spill_20260726.tif", "GTiff", "gdal")

Retrieve S1 image from second event (20260807)

#note the date second time step
t = ["2026-08-07", "2026-08-07"]
s1_cube = connection.load_collection(
  "SENTINEL1_GRD",
  spatial_extent={'west': 56.30, 'east': 56.90, 'south': 17.46, 'north': 17.90,"crs": "EPSG:4326"},
  temporal_extent=t,
  bands=["VV", "VH"]
)

s1_cube.download(work_dir+'/S1_spill_20260807.tif')

s1_new = ilwis.RasterCoverage('S1_spill_20260807.tif')
print(s1_new.size())
print(s1_new.envelope())
coordSys = s1_new.coordinateSystem()
coordSys.toWKT()

Size(6374, 4883, 2)
425670.000000 1930440.000000 489410.000000 1979270.000000

'PROJCS["s1_spill_20260807.tif",GEOCS["s1_spill_20260807.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",57],UNIT[meter,1.0]]'

VV_new = ilwis.do('selection',s1_new,"rasterbands(0)")

#create histogram using large number of bins, given the fact that the input image is a float64, with both positive and negative numbers
hist = VV_new.statistics(ilwis.PropertySets.pHISTOGRAM, 65535)
minPerc1, maxPerc1 = hist.calcStretchRange(1)
VVs_new = ilwis.do('linearstretch',VV_new, minPerc1, maxPerc1)
VVs_new = ilwis.do('mapcalc','iff(@1==?,255,@1)', VVs_new) #modify pixels over water
VVs_new = ilwis.do('setvaluerange', VVs_new, 0, 255, 1)
#uncomment lines below if you want to store the stretched image
#VVs_new.store('VVs_new.mpr')

VVnew_2np = np.fromiter(iter(VVs_new), np.ubyte, VVs_new.size().linearSize()) 
VVnew_2np = VVnew_2np.reshape((VVs_new.size().ysize, VVs_new.size().xsize))

fig = plt.figure(figsize=(15, 10))

plt.imshow(VVnew_2np, vmin=0, vmax=255, cmap= "turbo")
plt.axis("off")
plt.title('VV - S1 image from 20260807 - retrieved from OpenEO service provider');

To see the temporal developments, a difference image is calculated

vv_dif = ilwis.do('mapcalc','(@2 - @1)', VVs, VVs_new)
#uncomment lines below if you want to store the difference image
#vv_dif.store('vv_dif.mpr')

vv_dif_2np = np.fromiter(iter(vv_dif), np.float64, vv_dif.size().linearSize()) 
vv_dif_2np = vv_dif_2np.reshape((vv_dif.size().ysize, vv_dif.size().xsize))

fig = plt.figure(figsize=(15, 10))

plt.imshow(vv_dif_2np, vmin=-100, vmax=100, cmap= "RdBu_r")
plt.colorbar(shrink=0.45, extend = 'both')
plt.axis("off")
plt.title('VV difference - S1 images from 20260726 and 20260807');

Interpretation of the difference image:

• vv_diff < 0 = backscatter decreased
• vv_diff > 0 = backscatter increased
• vv_diff ≈ 0 = little change

For an oil spill, a negative change in VV backscatter may be observed where the oil dampens short surface waves. This is clearly visible in the extended darker bluish areas. The image also shows the disappearance of the dark reddish linear feature, which likely represented an oil sheen surface pollution feature.

Retrieve Sentinel 2 L2A image from the second oil spill event

#note the date - second event
t = ["2026-08-07", "2026-08-07"]
s2_cube = connection.load_collection("SENTINEL2_L2A",
    spatial_extent={'west': 56.30, 'east': 56.90, 'south': 17.46, 'north': 17.90,"crs": "EPSG:4326"},
    temporal_extent= t,
    bands=['B02', 'B03', 'B04', 'B08', 'B11', 'B12'],
    max_cloud_cover=100,
)

s2_cube.download(work_dir+'/S2_spill_20260807.tif')

#read the S2 image in ilwispy
s2_date2 = ilwis.RasterCoverage ('S2_spill_20260807.tif')
print(s2_date2.size())

Size(6374, 4883, 6)

For visualization in Matplotlib (using Imshow) the data is stretched, transformed into a numpy array, using the ILWISPY operation 'iter', each pixel is assigned into a 1 D array, which by numpy is reshaped back into a 2-D array - for each of the spectral channels

#stretch the rasterbands using a loop
multiple_stretch = []
multiple_bands = ilwis.do('selection',s2_date2,"rasterbands(0..5)") 
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

#note index starts from 0
Blues = ilwis.do('selection',mb_stretch,"rasterbands(0)")
Greens = ilwis.do('selection',mb_stretch,"rasterbands(1)")
Reds = ilwis.do('selection',mb_stretch,"rasterbands(2)")
IRs = ilwis.do('selection',mb_stretch,"rasterbands(3)")
SWIR11s = ilwis.do('selection',mb_stretch,"rasterbands(4)")
SWIR12s = ilwis.do('selection',mb_stretch,"rasterbands(5)")

#Transform to a numpy array for visualization using imshow
row = IRs.size().ysize 
col = IRs.size().xsize
dim = IRs.size().linearSize()
print(row, col, dim)

Blues_2np = np.fromiter(iter(Blues), np.ubyte, dim) 
Blues_2np = Blues_2np.reshape((row, col))

Greens_2np = np.fromiter(iter(Greens), np.ubyte, dim) 
Greens_2np = Greens_2np.reshape((row, col))

Reds_2np = np.fromiter(iter(Reds), np.ubyte, dim) 
Reds_2np = Reds_2np.reshape((row, col))

IRs_2np = np.fromiter(iter(IRs), np.ubyte, dim) 
IRs_2np = IRs_2np.reshape((row, col))

SWIR11s_2np = np.fromiter(iter(SWIR11s), np.ubyte, dim) 
SWIR11s_2np = SWIR11s_2np.reshape((row, col))

SWIR12s_2np = np.fromiter(iter(SWIR12s), np.ubyte, dim) 
SWIR12s_2np = SWIR12s_2np.reshape((row, col))

4883 6374 31124242

A numpy multi dimensional data stack is created, for visualization in Matplotlib the data is organized pixel interleaved, to use the results later as an ILWIS Maplist, the data stack is transformed into a multi-dimensional array which is band interleaved

# create  color composite (in RGB)
#pixel interleaved
VIS_NIR_S2 = np.dstack((IRs_2np, Reds_2np, Greens_2np))
NIR_SWIR_S2 = np.dstack((Reds_2np, SWIR11s_2np, SWIR12s_2np))

#band interleaved
S2ALL = np.array([Blues_2np, Greens_2np, Reds_2np, IRs_2np, SWIR11s_2np, SWIR12s_2np])

Create a RGB plot showing a VIS-NIR color composite as well as a NIR-SWIR colour composite

fig1 = plt.figure(figsize=(15, 10))

plt.subplot(1, 2, 1)
plt.imshow(VIS_NIR_S2)
plt.axis("off")
plt.title('VIS-NIR image retrieved from 2026-08-07')

plt.subplot(1, 2, 2)
plt.imshow(NIR_SWIR_S2)
plt.axis("off")
plt.title('NIR-SWIR image retrieved from 2026-08-07');

#create empty ilwis raster 
S2_ilw = ilwis.RasterCoverage()
defNumr = ilwis.DataDefinition(ilwis.NumericDomain('code=value'), ilwis.NumericRange(0, 255, 1))
S2_ilw.setDataDef(defNumr)
S2_ilw.setSize(ilwis.Size(6374, 4883, 6)) 
S2_ilw.setGeoReference(s2_date2.geoReference())

#read the numpy array
S2_ilw.array2raster(S2ALL.flatten())

#comment if you don't want to store the image in an ILWIS maplist format
S2_ilw.store('S2_spill_20260807.mpl')

#uncomment if you want to store the image in geotif format
#S2_ilw.store("S2_spill_20260807.tif", "GTiff", "gdal")

From the image select once more the original data values and assign the spectral bands as seperate variables

#note index starts from 0
Blue = ilwis.do('selection',s2_date2,"rasterbands(0)")
Green = ilwis.do('selection',s2_date2,"rasterbands(1)")
Red = ilwis.do('selection',s2_date2,"rasterbands(2)")
IR = ilwis.do('selection',s2_date2,"rasterbands(3)")
SWIR11 = ilwis.do('selection',s2_date2,"rasterbands(4)")
SWIR12 = ilwis.do('selection',s2_date2,"rasterbands(5)")

Calculation of the band ratios in a similar manner as before.

Red = ilwis.do('mapcalc','(@2/@1)', Blue, Green)
Green = ilwis.do('mapcalc','(@2+@3)/@1', Blue, Green, Red)
Blue = ilwis.do('mapcalc','(@2+@3)/@1', IR, SWIR11, SWIR12)

#perform the image stretch and transform to numpy array
stat_Red = Red.statistics(ilwis.PropertySets.pHISTOGRAM, 65535)
minPerc1, maxPerc1 = stat_Red.calcStretchRange(1)
Reds = ilwis.do('linearstretch',Red, minPerc1, maxPerc1)
Reds = ilwis.do('setvaluerange', Reds, 0, 255, 1)

R_2np = np.fromiter(iter(Reds), np.ubyte, Red.size().linearSize()) 
R_2np = R_2np.reshape((Red.size().ysize, Red.size().xsize))

stat_Green = Green.statistics(ilwis.PropertySets.pHISTOGRAM, 65535)
minPerc1, maxPerc1 = stat_Green.calcStretchRange(1)
Greens = ilwis.do('linearstretch',Green, minPerc1, maxPerc1)
Greens = ilwis.do('setvaluerange', Greens, 0, 255, 1)

G_2np = np.fromiter(iter(Greens), np.ubyte, Red.size().linearSize()) 
G_2np = G_2np.reshape((Red.size().ysize, Red.size().xsize))

stat_Blue = Blue.statistics(ilwis.PropertySets.pHISTOGRAM, 65535)
minPerc1, maxPerc1 = stat_Blue.calcStretchRange(1)
Blues = ilwis.do('linearstretch',Blue, minPerc1, maxPerc1)
Blues = ilwis.do('setvaluerange', Blues, 0, 255, 1)

B_2np = np.fromiter(iter(Blues), np.ubyte, Red.size().linearSize()) 
B_2np = B_2np.reshape((Red.size().ysize, Red.size().xsize))

# create  color composite (in RGB)
RGB_S2 = np.dstack((R_2np, G_2np, B_2np))

fig = plt.figure(figsize=(15, 10))

plt.imshow(RGB_S2)
plt.axis("off")
plt.title('S2 band ratio composite from 20260807');

Integration of Sentinel 1 and 2

fused_rgb = np.stack([
    VV_2np,
    Greens_2np,
    SWIR12s_2np    
], axis=-1)

plt.figure(figsize=(10, 10))
plt.imshow(fused_rgb)
plt.axis("off")
plt.title("Sentinel-1 VV / Sentinel-2  Green / SWIR")
plt.show()

Data reduction and unsupervised image classification

With a 9-layer stack, one can perform PCA first, then use the principal components as input into an unsupervised classification such as K-means

data_stack = np.dstack((VV_2np, VVnew_2np, vv_dif_2np, Blues_2np, Greens_2np, Reds_2np, IRs_2np, SWIR11s_2np, SWIR12s_2np))

PCA expects observations × variables, check the shape and then reshape the data stack. While reshaping the stack the data is transformed from float64 to float32, to conserve memory

print(data_stack.shape)

(4883, 6374, 9)

nrows, ncols, nbands = data_stack.shape

X = data_stack.reshape(-1, nbands).astype(np.float32)

print(X.shape)

(31124242, 9)

Handle NaN/Inf values, this is important for satellite imagery

valid = np.all(np.isfinite(X), axis=1)

X_valid = X[valid]

print("Original pixels:", X.shape[0])
print("Valid pixels:", X_valid.shape[0])

Original pixels: 31124242
Valid pixels: 31124242

Standardize the 9 bands. This is particularly important in this case because layers can have very different units/ranges. For example, Sentinel-1 VV may be in dB, whereas Sentinel-2 bands may have DN/reflectance values. This gives each input layer approximately: mean = 0 and standard deviation = 1. Without this step, bands with larger numerical ranges can dominate the PCA.

scaler = StandardScaler()

X_scaled = scaler.fit_transform(X_valid)

Perform the PCA. Initially calculate all 9 components.

pca = PCA(n_components=9)

X_pca = pca.fit_transform(X_scaled)

Check the explained variance

for i, variance in enumerate(pca.explained_variance_ratio_, start=1):
    print(f"PC{i}: {variance*100:.2f}%")

PC1: 59.65%
PC2: 19.65%
PC3: 12.65%
PC4: 5.20%
PC5: 2.05%
PC6: 0.39%
PC7: 0.31%
PC8: 0.11%
PC9: 0.00%

Show the cumulative variance

print(
    "Cumulative variance:",
    np.cumsum(pca.explained_variance_ratio_) * 100
)

Cumulative variance: [59.651993 79.301094 91.94712  97.14799  99.196884 99.583885 99.891685
 99.99999  99.99999 ]

Plot the PCA components, only the first 5 components seem te be relevant, given the variance explained

#reshape the 1-D arrays
n_pc = 5

pca_image = np.full(
    (X.shape[0], n_pc),
    np.nan
)

pca_image[valid] = X_pca[:, :n_pc]

pca_image = pca_image.reshape(
    nrows,
    ncols,
    n_pc
)

#plot individual PCAs
fig, axes = plt.subplots(
    1, 5,
    figsize=(20, 5)
)

for i, ax in enumerate(axes):

    im = ax.imshow(
        pca_image[:, :, i],
        cmap="turbo"
    )

    ax.set_title(f"PC{i+1}")
    ax.axis("off")

    plt.colorbar(
        im,
        ax=ax,
        shrink=0.5
    )

plt.tight_layout()
plt.show()

Create a RGB composite. Because PCA components can have negative values and different ranges, normalize each component before displaying

def normalize_pc(arr, pmin=2, pmax=98):
    """Percentile stretch to 0-1 for RGB visualization."""
    vmin, vmax = np.nanpercentile(arr, [pmin, pmax])
    return np.clip((arr - vmin) / (vmax - vmin), 0, 1)

# -------------------------------------------------------
# Composite 1: PC1 = Red, PC2 = Green, PC5 = Blue
# -------------------------------------------------------
pc1 = normalize_pc(pca_image[:, :, 0])
pc2 = normalize_pc(pca_image[:, :, 1])
pc5 = normalize_pc(pca_image[:, :, 4])

rgb_125 = np.dstack((pc1, pc2, pc5))


# -------------------------------------------------------
# Composite 2: PC2 = Red, PC3 = Green, PC4 = Blue
# -------------------------------------------------------
pc2 = normalize_pc(pca_image[:, :, 1])
pc3 = normalize_pc(pca_image[:, :, 2])
pc4 = normalize_pc(pca_image[:, :, 3])

rgb_234 = np.dstack((pc2, pc3, pc4))

# -------------------------------------------------------
# Plot
# -------------------------------------------------------
fig, axes = plt.subplots(1, 2, figsize=(18, 8))

axes[0].imshow(rgb_125)
axes[0].set_title("PCA composite: PC1 - PC2 - PC5")
axes[0].axis("off")

axes[1].imshow(rgb_234)
axes[1].set_title("PCA composite: PC2 - PC3 - PC4")
axes[1].axis("off")

plt.tight_layout()
plt.show()

K-means clustering

From the above PCA visualizations a lot of oil spill details and temporal changes can be observed. Can these also be classified using an unsupervised image classification routine?

#using the higher order PCA components only - see variance explained, cumulative variance
n_pc = 5

X_pca_reduced = X_pca[:, :n_pc]

Perform K-means clustering, for example, start with 10 classes

n_classes = 20

kmeans = KMeans(
    n_clusters=n_classes,
    random_state=42,
    n_init=10
)

clusters = kmeans.fit_predict(X_pca_reduced)

Convert the classification back to an image

classification = np.full(
    X.shape[0],
    -1,
    dtype=np.int16
)

classification[valid] = clusters

classification = classification.reshape(
    nrows,
    ncols
)

Plot the result

# -------------------------------------------------------
# Create a discrete colormap
# -------------------------------------------------------
cmap = plt.get_cmap("tab20", n_classes)

# Boundaries centered around integer class numbers
bounds = np.arange(-0.5, n_classes + 0.5, 1)

norm = BoundaryNorm(
    bounds,
    cmap.N
)

# -------------------------------------------------------
# Plot
# -------------------------------------------------------
fig, ax = plt.subplots(figsize=(12, 10))

im = ax.imshow(
    classification,
    cmap=cmap,
    norm=norm
)

cbar = plt.colorbar(
    im,
    ax=ax,
    boundaries=bounds,
    ticks=np.arange(n_classes),
    spacing="proportional", shrink=0.5,
)

cbar.set_label("K-means cluster / class")

ax.set_title(f"Oil spill K-means classification - {n_classes} classes")
ax.axis("off")

plt.show()

Export the results to an ilwis or geotif raster

#create empty ilwis raster 
Km_ilw = ilwis.RasterCoverage()
defNumr = ilwis.DataDefinition(ilwis.NumericDomain('code=value'), ilwis.NumericRange(0, 255, 1))
Km_ilw.setDataDef(defNumr)
Km_ilw.setSize(ilwis.Size(6374, 4883, 1)) 
Km_ilw.setGeoReference(s2_date2.geoReference())

#read the numpy array
Km_ilw.array2raster(classification.flatten())

#comment if you don't want to store the image in an ILWIS maplist format
Km_ilw.store('oil_spill_classification.mpr')

#uncomment if you want to store the image in geotif format
#Km_ilw.store("oil_spill_classification.tif", "GTiff", "gdal")

As there is some confusion between the clouds and oil spill, include the Sentinel-1 VV difference image explicitly in the classification, this gives K-means the first five PCA dimensions plus the original standardized VV difference.

X_class = np.column_stack([
    X_pca[:, :5],
    X_scaled[:, 2]       # VV difference (0-based indexing)
])

kmeans = KMeans(
    n_clusters=15,
    random_state=42,
    n_init=10
)

clusters = kmeans.fit_predict(X_class)

classification = np.full(
    X.shape[0],
    -1,
    dtype=np.int16
)

classification[valid] = clusters

classification1 = classification.reshape(
    nrows,
    ncols
)

# -------------------------------------------------------
# Plot
# -------------------------------------------------------
fig, ax = plt.subplots(figsize=(12, 10))

im = ax.imshow(
    classification1,
    cmap=cmap,
    norm=norm
)

cbar = plt.colorbar(
    im,
    ax=ax,
    boundaries=bounds,
    ticks=np.arange(n_classes),
    spacing="proportional", shrink=0.5,
)

cbar.set_label("K-means cluster / class")

ax.set_title(f"Oil spill K-means classification - {n_classes} classes - VV-difference enforced")
ax.axis("off")

plt.show()

Export the results to an ilwis or geotif raster

#read the numpy array
Km_ilw.array2raster(classification1.flatten())

#comment if you don't want to store the image in an ILWIS maplist format
Km_ilw.store('oil_spill_classification1.mpr')

#uncomment if you want to store the image in geotif format
#Km_ilw.store("oil_spill_classification1.tif", "GTiff", "gdal")

Main conclusions:

• Sentinel-1 VV is highly useful for detecting the oil-spill signature. The negative VV difference between the two acquisition dates indicates a reduction in radar backscatter. This is consistent with an oil film damping short surface waves, causing the sea surface to appear darker in Sentinel-1 imagery.
• The temporal Sentinel-1 information is more informative than a single VV image. Comparing the two VV observations allows you to identify areas that changed between the acquisitions. This helps distinguish persistent dark-water areas from newly developed features potentially associated with the spill.
• Sentinel-2 provides valuable complementary information. The optical bands help distinguish water, land, clouds and other surface features. The combination of Sentinel-1 and Sentinel-2 therefore provides more information than either sensor alone.
• PCA reduces the dimensionality of the combined dataset. The 9-layer stack contains correlated information. PCA transforms these variables into a smaller number of components while retaining most of the variance. The PCA composites can reveal spatial patterns that are difficult to see in individual bands.
• The PCA composites can enhance the visual separation of surface features. The PC1-PC2-PC5 and PC2-PC3-PC4 composites provide different representations of the scene. If the potential oil-spill region has a distinct colour or tonal response in one of these composites, this suggests that the combined spectral and SAR information contains useful discriminatory information.
• K-means can separate several surface classes, but oil and some other ocean / cloud features may remain confused. This is an important result rather than necessarily a failure. Increasing the number of clusters allows the algorithm to subdivide the feature space and potentially separate different types of water, cloud, cloud shadow, oil sheen and other features. Remember that K-means is unsupervised. It identifies statistically similar groups of pixels, not physical oil classes.
• The VV temporal difference is particularly important for resolving the oil/other class ambiguity. e.g. clouds can be very similar to some oil-spill features in optical imagery, but they should not produce the same temporal SAR backscatter response. Including vv_dif_2np in the classification therefore provides an important additional discriminator.