{ "cells": [ { "cell_type": "markdown", "id": "144c0f32-91c6-447e-8171-e7ee6d15be1b", "metadata": {}, "source": [ "# Oil spill detection using Sentinel 1 and Sentinel 2, example Al Qibiliya - Oman" ] }, { "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": [ "Microwave sensing:\n", "+ 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. \n", "\n", "Multispectral sensing:\n", "+ 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.\n", "+ 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.\n", "+ 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.\n", "+ 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.\n", "\n", "Review:\n", "+ https://blog.response.restoration.noaa.gov/how-thick-oil-slick\n", "+ https://www.sciencedirect.com/science/article/pii/S0034425719304407?via%3Dihub\n", "+ https://www.sciencedirect.com/science/article/pii/S2215016121001205" ] }, { "cell_type": "markdown", "id": "395cd765-b118-4225-abf5-e2bea35a07ac", "metadata": {}, "source": [ "For additional background information on this oil spill incident see:\n", "+ https://www.nhregister.com/news/world/article/oil-spill-from-grounded-tanker-off-oman-expands-22373858.php\n", "+ https://www.khaleejtimes.com/world/gulf/oman-oil-spill-al-hallaniyat-islands-ship-grounds\n", "\n", "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." ] }, { "cell_type": "markdown", "id": "bb6851d7-5a50-4802-a773-aeee0f07c807", "metadata": {}, "source": [ "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." ] }, { "cell_type": "markdown", "id": "e01b103f-d0b1-40b9-bfd3-b7934fc54f21", "metadata": {}, "source": [ "### Download and pre-processing data\n", "\n", "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/" ] }, { "cell_type": "code", "execution_count": null, "id": "7c4ac941-af9c-423c-ba33-93e37a8db12a", "metadata": {}, "outputs": [], "source": [ "#import required libraries - eventually install if not available \n", "import os\n", "import ilwis\n", "import openeo\n", "from openeo.rest.auth.config import RefreshTokenStore\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import matplotlib.patches as patches\n", "import matplotlib.ticker as mticker\n", "import cartopy\n", "from cartopy import crs as ccrs, feature as cfeature\n", "from sklearn.decomposition import PCA\n", "from sklearn.preprocessing import StandardScaler\n", "from sklearn.cluster import KMeans\n", "from matplotlib.colors import BoundaryNorm\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": [ "#remove previous login credentials, uncomment if running the notebook for the first time\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": [ "### Indentify the shipwreck location" ] }, { "cell_type": "code", "execution_count": null, "id": "db0f248c-7518-490c-aa4f-6ff127c62fb7", "metadata": {}, "outputs": [], "source": [ "#check the meta data information\n", "connection.describe_collection(\"SENTINEL2_L2A\")" ] }, { "cell_type": "markdown", "id": "906ade48-6afd-4eb2-9d7e-064666a998c0", "metadata": {}, "source": [ "Collect a Sentinel 2 image before the ship wreck event" ] }, { "cell_type": "code", "execution_count": null, "id": "8b474eaf-3424-40a7-b37a-6addc4756027", "metadata": {}, "outputs": [], "source": [ "#note the date \n", "t = [\"2026-06-18\", \"2026-06-18\"]\n", "s2_cube = connection.load_collection(\"SENTINEL2_L2A\",\n", " spatial_extent={'west': 56.31, 'east': 56.355, 'south': 17.485, 'north': 17.515,\"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": "62ae0e7d-1f2d-4cd7-b522-61f09f3440a9", "metadata": {}, "outputs": [], "source": [ "s2_cube.download(work_dir+\"/s2_background.tif\")" ] }, { "cell_type": "code", "execution_count": null, "id": "6e8b41c8-d6d5-456d-bafc-6170771bcb20", "metadata": {}, "outputs": [], "source": [ "#read the S2 image in ilwispy\n", "s2_background = ilwis.RasterCoverage ('s2_background.tif')\n", "print(s2_background.size())" ] }, { "cell_type": "markdown", "id": "2d48cd2e-fc81-43b3-bc8b-cf0e93ff2563", "metadata": {}, "source": [ "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!" ] }, { "cell_type": "code", "execution_count": null, "id": "4f9b71e3-5767-4242-9dcb-9d548720ab22", "metadata": {}, "outputs": [], "source": [ "#stretch the rasterbands contained in s2_background using a loop\n", "multiple_stretch = []\n", "multiple_bands = ilwis.do('selection',s2_background,\"rasterbands(0..2)\") \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": "8ae73649-b5ff-4022-b358-fc81d0b39201", "metadata": {}, "outputs": [], "source": [ "#store the results as an ILWIS maplist - display the map using the ilwis386 desktop software\n", "mb_stretch.store('S2_background.mpl')" ] }, { "cell_type": "code", "execution_count": null, "id": "fd0d3e34-5582-4fee-a55d-50ff9d01cbe3", "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": "94e65ae2-07b9-4283-97e4-14bf050bb477", "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": "e647e332-fdca-42ed-8125-5cd59bbfe187", "metadata": {}, "outputs": [], "source": [ "#create a numpy 3D data stack\n", "ncol = np.dstack((Reds_2np, Greens_2np, Blues_2np))" ] }, { "cell_type": "code", "execution_count": null, "id": "15059802-54ae-492f-b8c9-7b23b7bf373d", "metadata": {}, "outputs": [], "source": [ "img_extent = (56.31, 56.355, 17.485, 17.515)\n", "\n", "fig = plt.figure(figsize=(12, 10))\n", "ax = plt.axes(projection=ccrs.PlateCarree())\n", "\n", "plt.title('Sentinel-2 composite of 2026-06-18 of Al Qibiliya - Arabian Sea, Oman')\n", "\n", "# Data raster\n", "ax.imshow(\n", " ncol,\n", " origin='upper',\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.01)\n", "gl.ylocator = mticker.MultipleLocator(0.01)\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": "e6163bf1-88cd-4d00-b5c8-a363b84bf98d", "metadata": {}, "source": [ "Process the same area of interest after ship stranded on the southwestern coast of the island" ] }, { "cell_type": "code", "execution_count": null, "id": "e3967531-8833-4123-96bd-620de8fb4fc1", "metadata": {}, "outputs": [], "source": [ "#note the date and spectral channels\n", "t = [\"2026-07-28\", \"2026-07-28\"]\n", "s2_cube = connection.load_collection(\"SENTINEL2_L2A\",\n", " spatial_extent={'west': 56.31, 'east': 56.355, 'south': 17.485, 'north': 17.515,\"crs\": \"EPSG:4326\"},\n", " temporal_extent= t,\n", " bands=['B02', 'B03', 'B04', 'B08', 'B11', 'B12'],\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_ship_stranded.tif\")" ] }, { "cell_type": "code", "execution_count": null, "id": "a2145e03-ec90-48c0-9de4-7d86dfe3a5b6", "metadata": {}, "outputs": [], "source": [ "#read the S2 image in ilwispy\n", "s2_ship = ilwis.RasterCoverage ('s2_ship_stranded.tif')\n", "print(s2_ship.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_ship,\"rasterbands(0..5)\") \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_ship_stranded.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)\")\n", "IRs = ilwis.do('selection',mb_stretch,\"rasterbands(3)\")\n", "SWIR11s = ilwis.do('selection',mb_stretch,\"rasterbands(4)\")\n", "SWIR12s = ilwis.do('selection',mb_stretch,\"rasterbands(5)\")" ] }, { "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 3D data stack\n", "ncol = np.dstack((Reds_2np, Greens_2np, Blues_2np))" ] }, { "cell_type": "code", "execution_count": null, "id": "31e93451-dde6-4a35-ac06-2a5aeaec175c", "metadata": {}, "outputs": [], "source": [ "img_extent = (56.31, 56.355, 17.485, 17.515)\n", "\n", "fig = plt.figure(figsize=(12, 10))\n", "ax = plt.axes(projection=ccrs.PlateCarree())\n", "\n", "plt.title('Sentinel-2 natural color composite of 2026-07-28 showing location of shipwreck\\n'\n", " 'and initial oil spill features, southwest of Al Qibiliya - Arabian Sea, Oman')\n", "\n", "# Data raster\n", "ax.imshow(\n", " ncol,\n", " origin='upper',\n", " extent=img_extent,\n", " transform=ccrs.PlateCarree()\n", ")\n", "\n", "# Shipwreck location\n", "lon = 56.3213\n", "lat = 17.4938\n", "\n", "circle = patches.Circle(\n", " (lon, lat),\n", " radius=0.0015, #in degrees\n", " fill=False,\n", " edgecolor='red',\n", " linewidth=2,\n", " transform=ccrs.PlateCarree(),\n", " zorder=5\n", ")\n", "\n", "ax.add_patch(circle)\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.01)\n", "gl.ylocator = mticker.MultipleLocator(0.01)\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": "99166c18-b192-44bf-b615-7263e6260d3c", "metadata": {}, "source": [ "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:\n", "\n", "+ Red: B3/B2\n", "+ Green: (B11+B12)/B8\n", "+ Blue: (B3+B4)/B2" ] }, { "cell_type": "code", "execution_count": null, "id": "1440496f-7b91-43c7-876e-84ba275dcdab", "metadata": {}, "outputs": [], "source": [ "#The band ratios are derived using the original data, so once more load the original individual Sentinel-2 spectral channels\n", "B2_org = ilwis.do('selection',multiple_bands,\"rasterbands(0)\")\n", "B3_org = ilwis.do('selection',multiple_bands,\"rasterbands(1)\")\n", "B4_org = ilwis.do('selection',multiple_bands,\"rasterbands(2)\")\n", "B8_org = ilwis.do('selection',multiple_bands,\"rasterbands(3)\")\n", "SWIR11_org = ilwis.do('selection',multiple_bands,\"rasterbands(4)\")\n", "SWIR12_org = ilwis.do('selection',multiple_bands,\"rasterbands(5)\")" ] }, { "cell_type": "code", "execution_count": null, "id": "b9025122-f2c3-4a6d-932f-27585ee2e654", "metadata": {}, "outputs": [], "source": [ "#calculate the band ratios\n", "Red = ilwis.do('mapcalc','(@2/@1)', B2_org, B3_org)\n", "Green = ilwis.do('mapcalc','(@2+@3)/@1', B2_org, B3_org, B4_org)\n", "Blue = ilwis.do('mapcalc','(@2+@3)/@1', B8_org, SWIR11_org, SWIR12_org)" ] }, { "cell_type": "code", "execution_count": null, "id": "f9b0efc3-4d73-4b33-851c-2a7867cd2748", "metadata": {}, "outputs": [], "source": [ "#perform the image stretch and transform to numpy array\n", "stat_Red = Red.statistics(ilwis.PropertySets.pHISTOGRAM, 65535)\n", "minPerc1, maxPerc1 = stat_Red.calcStretchRange(1)\n", "Reds = ilwis.do('linearstretch',Red, minPerc1, maxPerc1)\n", "Reds = ilwis.do('setvaluerange', Reds, 0, 255, 1)\n", "\n", "R_2np = np.fromiter(iter(Reds), np.ubyte, Red.size().linearSize()) \n", "R_2np = R_2np.reshape((Red.size().ysize, Red.size().xsize))\n", "\n", "\n", "stat_Green = Green.statistics(ilwis.PropertySets.pHISTOGRAM, 65535)\n", "minPerc1, maxPerc1 = stat_Green.calcStretchRange(1)\n", "Greens = ilwis.do('linearstretch',Green, minPerc1, maxPerc1)\n", "Greens = ilwis.do('setvaluerange', Greens, 0, 255, 1)\n", "\n", "G_2np = np.fromiter(iter(Greens), np.ubyte, Red.size().linearSize()) \n", "G_2np = G_2np.reshape((Red.size().ysize, Red.size().xsize))\n", "\n", "stat_Blue = Blue.statistics(ilwis.PropertySets.pHISTOGRAM, 65535)\n", "minPerc1, maxPerc1 = stat_Blue.calcStretchRange(1)\n", "Blues = ilwis.do('linearstretch',Blue, minPerc1, maxPerc1)\n", "Blues = ilwis.do('setvaluerange', Blues, 0, 255, 1)\n", "\n", "B_2np = np.fromiter(iter(Blues), np.ubyte, Red.size().linearSize()) \n", "B_2np = B_2np.reshape((Red.size().ysize, Red.size().xsize))" ] }, { "cell_type": "code", "execution_count": null, "id": "24bad390-ace1-486f-ac13-12aa1c53c929", "metadata": {}, "outputs": [], "source": [ "# create color composite (in RGB)\n", "RGB_S2 = np.dstack((R_2np, G_2np, B_2np))" ] }, { "cell_type": "code", "execution_count": null, "id": "38f3a07d-4e68-4713-8dd7-640c6822e548", "metadata": {}, "outputs": [], "source": [ "img_extent = (56.31, 56.355, 17.485, 17.515)\n", "\n", "fig = plt.figure(figsize=(12, 10))\n", "\n", "ax = plt.axes(projection=ccrs.PlateCarree())\n", "\n", "# Data raster\n", "ax.imshow(\n", " RGB_S2,\n", " origin='upper',\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.01)\n", "gl.ylocator = mticker.MultipleLocator(0.01)\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.title('RGB composite of 20260728 using band transforms for oil spill enhancement')\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "12b8a22c-2b33-4a82-befd-f3f78dbcac0a", "metadata": {}, "source": [ "### Retrieve Sentinel 1 GRD image from the oil spill event on 20260726" ] }, { "cell_type": "code", "execution_count": null, "id": "212c5049-065d-4e15-a2fc-3ce1f1928a6f", "metadata": {}, "outputs": [], "source": [ "connection.describe_collection(\"SENTINEL1_GRD\")" ] }, { "cell_type": "code", "execution_count": null, "id": "19d70cac-a477-4b2e-a882-09216a88ceea", "metadata": {}, "outputs": [], "source": [ "#note the date first time step\n", "t = [\"2026-07-26\", \"2026-07-26\"]\n", "s1_cube = connection.load_collection(\n", " \"SENTINEL1_GRD\",\n", " spatial_extent={'west': 56.30, 'east': 56.90, 'south': 17.46, 'north': 17.90,\"crs\": \"EPSG:4326\"},\n", " temporal_extent=t,\n", " bands=[\"VV\", \"VH\"]\n", ")" ] }, { "cell_type": "code", "execution_count": null, "id": "0c35a19e-b4c8-42b4-98fe-3374681f870e", "metadata": {}, "outputs": [], "source": [ "s1_cube.download(work_dir+'/S1_spill_20260726.tif')" ] }, { "cell_type": "markdown", "id": "68d19b86-6254-46cd-93a6-dfa399321370", "metadata": {}, "source": [ "#### Create radar RGB composite \n", "+ VV - red\n", "+ VH - green\n", "+ VC - blue (cross polarized)" ] }, { "cell_type": "code", "execution_count": null, "id": "597488f8-d09a-4aaa-9ec2-64d4543be524", "metadata": {}, "outputs": [], "source": [ "s1 = ilwis.RasterCoverage('S1_spill_20260726.tif')\n", "print(s1.size())\n", "print(s1.envelope())\n", "coordSys = s1.coordinateSystem()\n", "coordSys.toWKT()" ] }, { "cell_type": "code", "execution_count": null, "id": "afbb1aab-ef6d-4f81-b687-2a1d9f68f157", "metadata": {}, "outputs": [], "source": [ "VV = ilwis.do('selection',s1,\"rasterbands(0)\")" ] }, { "cell_type": "code", "execution_count": null, "id": "58c5fb51-c600-418c-a49e-db1928b87cff", "metadata": { "scrolled": true }, "outputs": [], "source": [ "#create histogram using large number of bins, given the fact that the input image is a float64, with both positive and negative numbers\n", "hist = VV.statistics(ilwis.PropertySets.pHISTOGRAM, 65535)\n", "minPerc1, maxPerc1 = hist.calcStretchRange(1)\n", "VVs = ilwis.do('linearstretch',VV, minPerc1, maxPerc1)\n", "VVs = ilwis.do('mapcalc','iff(@1==?,255,@1)', VVs) #modify pixels over water\n", "VVs = ilwis.do('setvaluerange', VVs, 0, 255, 1)\n", "#uncomment lines below if you want to store the stretched image\n", "#VVs.store('VVs.mpr')\n", "\n", "VV_2np = np.fromiter(iter(VVs), np.ubyte, VVs.size().linearSize()) \n", "VV_2np = VV_2np.reshape((VVs.size().ysize, VVs.size().xsize))" ] }, { "cell_type": "code", "execution_count": null, "id": "1c19886c-815d-475a-a63c-ff83a3e70bb3", "metadata": {}, "outputs": [], "source": [ "fig = plt.figure(figsize=(15, 10))\n", "\n", "plt.imshow(VV_2np, vmin=0, vmax=255, cmap= \"turbo\")\n", "plt.axis(\"off\")\n", "plt.title('VV - S1 image from 20260726 - retrieved from OpenEO service provider');" ] }, { "cell_type": "code", "execution_count": null, "id": "7432ad9c-218b-4548-af39-a9cb8f516196", "metadata": {}, "outputs": [], "source": [ "VH = ilwis.do('selection',s1,\"rasterbands(1)\")" ] }, { "cell_type": "code", "execution_count": null, "id": "f74ff83b-6679-4426-97f8-1529ad0f4d6c", "metadata": {}, "outputs": [], "source": [ "#stretch the VH band\n", "hist = VH.statistics(ilwis.PropertySets.pHISTOGRAM, 65535)\n", "minPerc1, maxPerc1 = hist.calcStretchRange(1)\n", "VHs = ilwis.do('linearstretch',VH, minPerc1, maxPerc1)\n", "VHs = ilwis.do('mapcalc','iff(@1==?,0,@1)', VHs) #remove bad pixels over water\n", "VHs = ilwis.do('setvaluerange', VHs, 0, 255, 1)\n", "#uncomment lines below if you want to store the stretched image\n", "#VHs.store('VHs.mpr')\n", "\n", "VH_2np = np.fromiter(iter(VHs), np.ubyte, VHs.size().linearSize()) \n", "VH_2np = VH_2np.reshape((VHs.size().ysize, VHs.size().xsize))" ] }, { "cell_type": "code", "execution_count": null, "id": "b18a1af0-08a1-41f2-bcda-7836099152ab", "metadata": {}, "outputs": [], "source": [ "#calculate the cross polarized image (VV/VH)\n", "VC = ilwis.do('mapcalc','(@1/@2)', VV, VH)\n", "#uncomment lines below if you want to store the S1 cross polazization image\n", "#VC.store('VC.mpr')" ] }, { "cell_type": "code", "execution_count": null, "id": "0fa6dbe6-96a7-4e97-81b7-b34a25e52c52", "metadata": {}, "outputs": [], "source": [ "#create histogram using large number of bins, given the fact that the input image is a float64, with both positive and negative numbers\n", "hist = VC.statistics(ilwis.PropertySets.pHISTOGRAM, 65535)\n", "minPerc1, maxPerc1 = hist.calcStretchRange(1)\n", "VCs = ilwis.do('linearstretch',VC, minPerc1, maxPerc1)\n", "VCs = ilwis.do('mapcalc','iff(@1==?,255,@1)', VCs) #modify pixels over water\n", "VCs = ilwis.do('setvaluerange', VCs, 0, 255, 1)\n", "#uncomment lines below if you want to store the stretched image\n", "#VCs.store('VCs.mpr')\n", "\n", "VC_2np = np.fromiter(iter(VCs), np.ubyte, VCs.size().linearSize()) \n", "VC_2np = VC_2np.reshape((VCs.size().ysize, VCs.size().xsize))" ] }, { "cell_type": "code", "execution_count": null, "id": "f140ce80-8883-4d5a-8a5a-6bab1055ff3a", "metadata": {}, "outputs": [], "source": [ "# create color composite (in RGB)\n", "cc = np.dstack((VV_2np, VH_2np, VC_2np))\n", "\n", "#band interleaved\n", "S1ALL = np.array([VC_2np, VH_2np, VV_2np])" ] }, { "cell_type": "code", "execution_count": null, "id": "8f2b38cb-c398-45d2-a1ff-3d82608ab65c", "metadata": {}, "outputs": [], "source": [ "fig, axes = plt.subplots(1, 2, figsize=(14, 7))\n", "\n", "# Full image\n", "axes[0].imshow(cc)\n", "axes[0].axis(\"off\")\n", "axes[0].set_title(\"S1 polarization composite\")\n", "\n", "# South-west 1000 x 1000 pixels\n", "sw = cc[-1000:, :1000]\n", "\n", "axes[1].imshow(sw)\n", "axes[1].axis(\"off\")\n", "axes[1].set_title(\"South-west image corner enlargement\")\n", "\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "id": "44281f2c-b9d5-4954-8490-18013b521f19", "metadata": {}, "outputs": [], "source": [ "#create empty ilwis raster \n", "S1_ilw = ilwis.RasterCoverage()\n", "defNumr = ilwis.DataDefinition(ilwis.NumericDomain('code=value'), ilwis.NumericRange(0, 255, 1))\n", "S1_ilw.setDataDef(defNumr)\n", "S1_ilw.setSize(ilwis.Size(745, 446, 3)) \n", "S1_ilw.setGeoReference(s1.geoReference())" ] }, { "cell_type": "code", "execution_count": null, "id": "7e0e583f-a619-4aa3-a797-00f02b6c9f72", "metadata": {}, "outputs": [], "source": [ "#read the numpy array\n", "S1_ilw.array2raster(S1ALL.flatten())" ] }, { "cell_type": "code", "execution_count": null, "id": "bd2e0a89-5f3e-4345-9652-40008f2a6d30", "metadata": {}, "outputs": [], "source": [ "#comment if you don't want to store the image in an ILWIS maplist format\n", "S1_ilw.store('S1_spill_20260726.mpl')\n", "\n", "#uncomment if you want to store the image in geotif format\n", "#S1_ilw.store(\"S1_spill_20260726.tif\", \"GTiff\", \"gdal\")" ] }, { "cell_type": "markdown", "id": "708569f6-06c9-46e8-b190-ec7c0c49c45d", "metadata": {}, "source": [ "### Retrieve S1 image from second event (20260807)" ] }, { "cell_type": "code", "execution_count": null, "id": "4ea7af45-3ded-4b25-8162-0db1fff2fd99", "metadata": {}, "outputs": [], "source": [ "#note the date second time step\n", "t = [\"2026-08-07\", \"2026-08-07\"]\n", "s1_cube = connection.load_collection(\n", " \"SENTINEL1_GRD\",\n", " spatial_extent={'west': 56.30, 'east': 56.90, 'south': 17.46, 'north': 17.90,\"crs\": \"EPSG:4326\"},\n", " temporal_extent=t,\n", " bands=[\"VV\", \"VH\"]\n", ")" ] }, { "cell_type": "code", "execution_count": null, "id": "6c8d7072-4f3b-44b7-918a-8a27c7fc2dd0", "metadata": {}, "outputs": [], "source": [ "s1_cube.download(work_dir+'/S1_spill_20260807.tif')" ] }, { "cell_type": "code", "execution_count": null, "id": "b15bbba4-b589-4b3b-9e11-feb3ad3fb429", "metadata": {}, "outputs": [], "source": [ "s1_new = ilwis.RasterCoverage('S1_spill_20260807.tif')\n", "print(s1_new.size())\n", "print(s1_new.envelope())\n", "coordSys = s1_new.coordinateSystem()\n", "coordSys.toWKT()" ] }, { "cell_type": "code", "execution_count": null, "id": "9afc5aa6-5285-4feb-b298-d58eef5b2705", "metadata": {}, "outputs": [], "source": [ "VV_new = ilwis.do('selection',s1_new,\"rasterbands(0)\")" ] }, { "cell_type": "code", "execution_count": null, "id": "2ab0456b-e45a-4f41-abde-f3d497533442", "metadata": {}, "outputs": [], "source": [ "#create histogram using large number of bins, given the fact that the input image is a float64, with both positive and negative numbers\n", "hist = VV_new.statistics(ilwis.PropertySets.pHISTOGRAM, 65535)\n", "minPerc1, maxPerc1 = hist.calcStretchRange(1)\n", "VVs_new = ilwis.do('linearstretch',VV_new, minPerc1, maxPerc1)\n", "VVs_new = ilwis.do('mapcalc','iff(@1==?,255,@1)', VVs_new) #modify pixels over water\n", "VVs_new = ilwis.do('setvaluerange', VVs_new, 0, 255, 1)\n", "#uncomment lines below if you want to store the stretched image\n", "#VVs_new.store('VVs_new.mpr')\n", "\n", "VVnew_2np = np.fromiter(iter(VVs_new), np.ubyte, VVs_new.size().linearSize()) \n", "VVnew_2np = VVnew_2np.reshape((VVs_new.size().ysize, VVs_new.size().xsize))" ] }, { "cell_type": "code", "execution_count": null, "id": "72271088-de88-4a06-ac20-96665132877c", "metadata": {}, "outputs": [], "source": [ "fig = plt.figure(figsize=(15, 10))\n", "\n", "plt.imshow(VVnew_2np, vmin=0, vmax=255, cmap= \"turbo\")\n", "plt.axis(\"off\")\n", "plt.title('VV - S1 image from 20260807 - retrieved from OpenEO service provider');" ] }, { "cell_type": "markdown", "id": "a63b69be-01d7-45d2-9db8-f22b32727db2", "metadata": {}, "source": [ "To see the temporal developments, a difference image is calculated" ] }, { "cell_type": "code", "execution_count": null, "id": "a2ce8c03-60c5-49a9-a38d-03713068184a", "metadata": {}, "outputs": [], "source": [ "vv_dif = ilwis.do('mapcalc','(@2 - @1)', VVs, VVs_new)\n", "#uncomment lines below if you want to store the difference image\n", "#vv_dif.store('vv_dif.mpr')" ] }, { "cell_type": "code", "execution_count": null, "id": "3345cb61-2ed0-49a4-a660-1b232404b5c0", "metadata": {}, "outputs": [], "source": [ "vv_dif_2np = np.fromiter(iter(vv_dif), np.float64, vv_dif.size().linearSize()) \n", "vv_dif_2np = vv_dif_2np.reshape((vv_dif.size().ysize, vv_dif.size().xsize))" ] }, { "cell_type": "code", "execution_count": null, "id": "6db57592-3a39-4b37-bde6-e5f7d14eb38d", "metadata": {}, "outputs": [], "source": [ "fig = plt.figure(figsize=(15, 10))\n", "\n", "plt.imshow(vv_dif_2np, vmin=-100, vmax=100, cmap= \"RdBu_r\")\n", "plt.colorbar(shrink=0.45, extend = 'both')\n", "plt.axis(\"off\")\n", "plt.title('VV difference - S1 images from 20260726 and 20260807');" ] }, { "cell_type": "markdown", "id": "04eab66a-c089-4d9e-862e-0c01781e47cc", "metadata": {}, "source": [ "Interpretation of the difference image:\n", "\n", "+ vv_diff < 0 = backscatter decreased\n", "+ vv_diff > 0 = backscatter increased\n", "+ vv_diff ≈ 0 = little change\n", "\n", "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." ] }, { "cell_type": "markdown", "id": "66661470-e9ad-46aa-ab23-b7d838806123", "metadata": {}, "source": [ "### Retrieve Sentinel 2 L2A image from the second oil spill event" ] }, { "cell_type": "code", "execution_count": null, "id": "dbed769b-e33f-4d96-aecb-72719222b3aa", "metadata": {}, "outputs": [], "source": [ "#note the date - second event\n", "t = [\"2026-08-07\", \"2026-08-07\"]\n", "s2_cube = connection.load_collection(\"SENTINEL2_L2A\",\n", " spatial_extent={'west': 56.30, 'east': 56.90, 'south': 17.46, 'north': 17.90,\"crs\": \"EPSG:4326\"},\n", " temporal_extent= t,\n", " bands=['B02', 'B03', 'B04', 'B08', 'B11', 'B12'],\n", " max_cloud_cover=100,\n", ")" ] }, { "cell_type": "code", "execution_count": null, "id": "f7ac5f68-3d16-4931-b275-b858c33308af", "metadata": {}, "outputs": [], "source": [ "s2_cube.download(work_dir+'/S2_spill_20260807.tif')" ] }, { "cell_type": "code", "execution_count": null, "id": "5bf75db0-4394-42db-a687-726610d2d8c1", "metadata": {}, "outputs": [], "source": [ "#read the S2 image in ilwispy\n", "s2_date2 = ilwis.RasterCoverage ('S2_spill_20260807.tif')\n", "print(s2_date2.size())" ] }, { "cell_type": "markdown", "id": "33b748ef-c3b9-42fa-9867-e0b7e377dd13", "metadata": {}, "source": [ "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" ] }, { "cell_type": "code", "execution_count": null, "id": "85d4fb30-ef3d-4727-aca9-36d1dc9f5e0e", "metadata": {}, "outputs": [], "source": [ "#stretch the rasterbands using a loop\n", "multiple_stretch = []\n", "multiple_bands = ilwis.do('selection',s2_date2,\"rasterbands(0..5)\") \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": "050c66cc-c5a4-4723-91cc-6cc1b7daae49", "metadata": {}, "outputs": [], "source": [ "#note index starts from 0\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)\")\n", "IRs = ilwis.do('selection',mb_stretch,\"rasterbands(3)\")\n", "SWIR11s = ilwis.do('selection',mb_stretch,\"rasterbands(4)\")\n", "SWIR12s = ilwis.do('selection',mb_stretch,\"rasterbands(5)\")" ] }, { "cell_type": "code", "execution_count": null, "id": "7e4414b4-c428-4315-a3e4-fb5cd15c1a61", "metadata": {}, "outputs": [], "source": [ "#Transform to a numpy array for visualization using imshow\n", "row = IRs.size().ysize \n", "col = IRs.size().xsize\n", "dim = IRs.size().linearSize()\n", "print(row, col, dim)\n", "\n", "Blues_2np = np.fromiter(iter(Blues), np.ubyte, dim) \n", "Blues_2np = Blues_2np.reshape((row, col))\n", "\n", "Greens_2np = np.fromiter(iter(Greens), np.ubyte, dim) \n", "Greens_2np = Greens_2np.reshape((row, col))\n", "\n", "Reds_2np = np.fromiter(iter(Reds), np.ubyte, dim) \n", "Reds_2np = Reds_2np.reshape((row, col))\n", "\n", "IRs_2np = np.fromiter(iter(IRs), np.ubyte, dim) \n", "IRs_2np = IRs_2np.reshape((row, col))\n", "\n", "SWIR11s_2np = np.fromiter(iter(SWIR11s), np.ubyte, dim) \n", "SWIR11s_2np = SWIR11s_2np.reshape((row, col))\n", "\n", "SWIR12s_2np = np.fromiter(iter(SWIR12s), np.ubyte, dim) \n", "SWIR12s_2np = SWIR12s_2np.reshape((row, col))" ] }, { "cell_type": "markdown", "id": "18125cc1-3a47-4d18-bc9d-c113d90f79d5", "metadata": {}, "source": [ "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" ] }, { "cell_type": "code", "execution_count": null, "id": "6cae9bea-8edd-455d-ad20-3bb04bc73b65", "metadata": {}, "outputs": [], "source": [ "# create color composite (in RGB)\n", "#pixel interleaved\n", "VIS_NIR_S2 = np.dstack((IRs_2np, Reds_2np, Greens_2np))\n", "NIR_SWIR_S2 = np.dstack((Reds_2np, SWIR11s_2np, SWIR12s_2np))\n", "\n", "#band interleaved\n", "S2ALL = np.array([Blues_2np, Greens_2np, Reds_2np, IRs_2np, SWIR11s_2np, SWIR12s_2np])" ] }, { "cell_type": "markdown", "id": "8a93dc0f-88b1-4221-8b0f-0c3bb903f4d2", "metadata": {}, "source": [ "Create a RGB plot showing a VIS-NIR color composite as well as a NIR-SWIR colour composite" ] }, { "cell_type": "code", "execution_count": null, "id": "e674757c-2917-41fc-9354-28887c04067d", "metadata": {}, "outputs": [], "source": [ "fig1 = plt.figure(figsize=(15, 10))\n", "\n", "plt.subplot(1, 2, 1)\n", "plt.imshow(VIS_NIR_S2)\n", "plt.axis(\"off\")\n", "plt.title('VIS-NIR image retrieved from 2026-08-07')\n", "\n", "plt.subplot(1, 2, 2)\n", "plt.imshow(NIR_SWIR_S2)\n", "plt.axis(\"off\")\n", "plt.title('NIR-SWIR image retrieved from 2026-08-07');" ] }, { "cell_type": "code", "execution_count": null, "id": "424b57c3-16eb-493c-93d2-e3db0fa262f1", "metadata": {}, "outputs": [], "source": [ "#create empty ilwis raster \n", "S2_ilw = ilwis.RasterCoverage()\n", "defNumr = ilwis.DataDefinition(ilwis.NumericDomain('code=value'), ilwis.NumericRange(0, 255, 1))\n", "S2_ilw.setDataDef(defNumr)\n", "S2_ilw.setSize(ilwis.Size(6374, 4883, 6)) \n", "S2_ilw.setGeoReference(s2_date2.geoReference())" ] }, { "cell_type": "code", "execution_count": null, "id": "aa0827fb-6579-43a5-9fae-d64bad6b042b", "metadata": {}, "outputs": [], "source": [ "#read the numpy array\n", "S2_ilw.array2raster(S2ALL.flatten())" ] }, { "cell_type": "code", "execution_count": null, "id": "524fea5a-ecf0-4045-96bc-bb7993bdddbf", "metadata": {}, "outputs": [], "source": [ "#comment if you don't want to store the image in an ILWIS maplist format\n", "S2_ilw.store('S2_spill_20260807.mpl')\n", "\n", "#uncomment if you want to store the image in geotif format\n", "#S2_ilw.store(\"S2_spill_20260807.tif\", \"GTiff\", \"gdal\")" ] }, { "cell_type": "markdown", "id": "27071974-aac6-42bd-813d-2364a820ef2b", "metadata": {}, "source": [ "From the image select once more the original data values and assign the spectral bands as seperate variables" ] }, { "cell_type": "code", "execution_count": null, "id": "7dfdc0fd-1070-43e7-99de-5bdb7101035c", "metadata": {}, "outputs": [], "source": [ "#note index starts from 0\n", "Blue = ilwis.do('selection',s2_date2,\"rasterbands(0)\")\n", "Green = ilwis.do('selection',s2_date2,\"rasterbands(1)\")\n", "Red = ilwis.do('selection',s2_date2,\"rasterbands(2)\")\n", "IR = ilwis.do('selection',s2_date2,\"rasterbands(3)\")\n", "SWIR11 = ilwis.do('selection',s2_date2,\"rasterbands(4)\")\n", "SWIR12 = ilwis.do('selection',s2_date2,\"rasterbands(5)\")" ] }, { "cell_type": "markdown", "id": "7887724e-f393-4bae-84d2-a7cb735fecc0", "metadata": {}, "source": [ "Calculation of the band ratios in a similar manner as before." ] }, { "cell_type": "code", "execution_count": null, "id": "d1ee4c46-c9b0-4ed5-9caa-0431843ca856", "metadata": {}, "outputs": [], "source": [ "Red = ilwis.do('mapcalc','(@2/@1)', Blue, Green)\n", "Green = ilwis.do('mapcalc','(@2+@3)/@1', Blue, Green, Red)\n", "Blue = ilwis.do('mapcalc','(@2+@3)/@1', IR, SWIR11, SWIR12)" ] }, { "cell_type": "code", "execution_count": null, "id": "3a07c563-a17e-45d4-bd48-cefef95c0d7b", "metadata": {}, "outputs": [], "source": [ "#perform the image stretch and transform to numpy array\n", "stat_Red = Red.statistics(ilwis.PropertySets.pHISTOGRAM, 65535)\n", "minPerc1, maxPerc1 = stat_Red.calcStretchRange(1)\n", "Reds = ilwis.do('linearstretch',Red, minPerc1, maxPerc1)\n", "Reds = ilwis.do('setvaluerange', Reds, 0, 255, 1)\n", "\n", "R_2np = np.fromiter(iter(Reds), np.ubyte, Red.size().linearSize()) \n", "R_2np = R_2np.reshape((Red.size().ysize, Red.size().xsize))\n", "\n", "stat_Green = Green.statistics(ilwis.PropertySets.pHISTOGRAM, 65535)\n", "minPerc1, maxPerc1 = stat_Green.calcStretchRange(1)\n", "Greens = ilwis.do('linearstretch',Green, minPerc1, maxPerc1)\n", "Greens = ilwis.do('setvaluerange', Greens, 0, 255, 1)\n", "\n", "G_2np = np.fromiter(iter(Greens), np.ubyte, Red.size().linearSize()) \n", "G_2np = G_2np.reshape((Red.size().ysize, Red.size().xsize))\n", "\n", "stat_Blue = Blue.statistics(ilwis.PropertySets.pHISTOGRAM, 65535)\n", "minPerc1, maxPerc1 = stat_Blue.calcStretchRange(1)\n", "Blues = ilwis.do('linearstretch',Blue, minPerc1, maxPerc1)\n", "Blues = ilwis.do('setvaluerange', Blues, 0, 255, 1)\n", "\n", "B_2np = np.fromiter(iter(Blues), np.ubyte, Red.size().linearSize()) \n", "B_2np = B_2np.reshape((Red.size().ysize, Red.size().xsize))" ] }, { "cell_type": "code", "execution_count": null, "id": "f9944f43-e5ad-4111-901a-d96c6cb59d1e", "metadata": {}, "outputs": [], "source": [ "# create color composite (in RGB)\n", "RGB_S2 = np.dstack((R_2np, G_2np, B_2np))" ] }, { "cell_type": "code", "execution_count": null, "id": "01cb7bad-c39c-4efd-a91b-4c5ffcb3d051", "metadata": {}, "outputs": [], "source": [ "fig = plt.figure(figsize=(15, 10))\n", "\n", "plt.imshow(RGB_S2)\n", "plt.axis(\"off\")\n", "plt.title('S2 band ratio composite from 20260807');" ] }, { "cell_type": "markdown", "id": "984849b6-9f58-45f6-a7e7-366a0b47fe81", "metadata": {}, "source": [ "Integration of Sentinel 1 and 2" ] }, { "cell_type": "code", "execution_count": null, "id": "b3c4615f-6e40-4518-8a4d-b3df8c09572b", "metadata": {}, "outputs": [], "source": [ "fused_rgb = np.stack([\n", " VV_2np,\n", " Greens_2np,\n", " SWIR12s_2np \n", "], axis=-1)\n", "\n", "plt.figure(figsize=(10, 10))\n", "plt.imshow(fused_rgb)\n", "plt.axis(\"off\")\n", "plt.title(\"Sentinel-1 VV / Sentinel-2 Green / SWIR\")\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "b97d8da1-2fc7-445e-b74c-15ae8ae769bf", "metadata": {}, "source": [ "### Data reduction and unsupervised image classification\n", "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" ] }, { "cell_type": "code", "execution_count": null, "id": "93cb3c93-d59a-4505-86af-3ab5af2bac8d", "metadata": {}, "outputs": [], "source": [ "data_stack = np.dstack((VV_2np, VVnew_2np, vv_dif_2np, Blues_2np, Greens_2np, Reds_2np, IRs_2np, SWIR11s_2np, SWIR12s_2np))" ] }, { "cell_type": "markdown", "id": "a4ef1cd0-e131-48c1-8c44-1f68360fc085", "metadata": {}, "source": [ "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" ] }, { "cell_type": "code", "execution_count": null, "id": "6a6dac6c-d1dc-4a6d-ad35-30b9b7a47db1", "metadata": {}, "outputs": [], "source": [ "print(data_stack.shape)" ] }, { "cell_type": "code", "execution_count": null, "id": "e505b1b2-96b7-4051-86d7-4d9d5dcbcc1a", "metadata": {}, "outputs": [], "source": [ "nrows, ncols, nbands = data_stack.shape\n", "\n", "X = data_stack.reshape(-1, nbands).astype(np.float32)\n", "\n", "print(X.shape)" ] }, { "cell_type": "markdown", "id": "d4b5893b-62ac-411d-8d45-7881705f248d", "metadata": {}, "source": [ "Handle NaN/Inf values, this is important for satellite imagery" ] }, { "cell_type": "code", "execution_count": null, "id": "70d63ee6-7dca-47c3-9079-809948d04482", "metadata": {}, "outputs": [], "source": [ "valid = np.all(np.isfinite(X), axis=1)\n", "\n", "X_valid = X[valid]\n", "\n", "print(\"Original pixels:\", X.shape[0])\n", "print(\"Valid pixels:\", X_valid.shape[0])" ] }, { "cell_type": "markdown", "id": "cf70bbc5-227a-4fe7-aa60-a0291839bd0f", "metadata": {}, "source": [ "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." ] }, { "cell_type": "code", "execution_count": null, "id": "5b3e54a1-2322-4d22-9bcd-c6081fc4b06b", "metadata": {}, "outputs": [], "source": [ "scaler = StandardScaler()\n", "\n", "X_scaled = scaler.fit_transform(X_valid)" ] }, { "cell_type": "markdown", "id": "6c578df3-a859-4bf4-9575-13bf4314d55f", "metadata": {}, "source": [ "Perform the PCA. Initially calculate all 9 components." ] }, { "cell_type": "code", "execution_count": null, "id": "725a37bf-2b21-4eef-99a3-0f0910c23769", "metadata": {}, "outputs": [], "source": [ "pca = PCA(n_components=9)\n", "\n", "X_pca = pca.fit_transform(X_scaled)" ] }, { "cell_type": "markdown", "id": "508d1918-50e3-4cdf-8e08-ce21f28c6ff3", "metadata": {}, "source": [ "Check the explained variance" ] }, { "cell_type": "code", "execution_count": null, "id": "cbc8a6bc-f9e4-491a-b018-672a7039c750", "metadata": {}, "outputs": [], "source": [ "for i, variance in enumerate(pca.explained_variance_ratio_, start=1):\n", " print(f\"PC{i}: {variance*100:.2f}%\")" ] }, { "cell_type": "markdown", "id": "5ab62ed8-fd3a-4db1-a2d1-21053a793697", "metadata": {}, "source": [ "Show the cumulative variance" ] }, { "cell_type": "code", "execution_count": null, "id": "2bd94d07-14a3-40a7-85eb-0ae98dc8ebcc", "metadata": {}, "outputs": [], "source": [ "print(\n", " \"Cumulative variance:\",\n", " np.cumsum(pca.explained_variance_ratio_) * 100\n", ")" ] }, { "cell_type": "markdown", "id": "a0698eed-531c-4b50-af32-2f9e6db8b80f", "metadata": {}, "source": [ "Plot the PCA components, only the first 5 components seem te be relevant, given the variance explained" ] }, { "cell_type": "code", "execution_count": null, "id": "41601a36-8d45-4ccb-a744-9d2549699b70", "metadata": {}, "outputs": [], "source": [ "#reshape the 1-D arrays\n", "n_pc = 5\n", "\n", "pca_image = np.full(\n", " (X.shape[0], n_pc),\n", " np.nan\n", ")\n", "\n", "pca_image[valid] = X_pca[:, :n_pc]\n", "\n", "pca_image = pca_image.reshape(\n", " nrows,\n", " ncols,\n", " n_pc\n", ")" ] }, { "cell_type": "code", "execution_count": null, "id": "97b55b44-c26b-4dba-8637-208d0f71f1b1", "metadata": {}, "outputs": [], "source": [ "#plot individual PCAs\n", "fig, axes = plt.subplots(\n", " 1, 5,\n", " figsize=(20, 5)\n", ")\n", "\n", "for i, ax in enumerate(axes):\n", "\n", " im = ax.imshow(\n", " pca_image[:, :, i],\n", " cmap=\"turbo\"\n", " )\n", "\n", " ax.set_title(f\"PC{i+1}\")\n", " ax.axis(\"off\")\n", "\n", " plt.colorbar(\n", " im,\n", " ax=ax,\n", " shrink=0.5\n", " )\n", "\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "762b2f00-afbc-42d2-b2a9-f7287eec775d", "metadata": {}, "source": [ "Create a RGB composite. Because PCA components can have negative values and different ranges, normalize each component before displaying " ] }, { "cell_type": "code", "execution_count": null, "id": "fd0aa5ce-b123-4d80-9e8f-fc2414b5abd3", "metadata": {}, "outputs": [], "source": [ "def normalize_pc(arr, pmin=2, pmax=98):\n", " \"\"\"Percentile stretch to 0-1 for RGB visualization.\"\"\"\n", " vmin, vmax = np.nanpercentile(arr, [pmin, pmax])\n", " return np.clip((arr - vmin) / (vmax - vmin), 0, 1)\n", "\n", "# -------------------------------------------------------\n", "# Composite 1: PC1 = Red, PC2 = Green, PC5 = Blue\n", "# -------------------------------------------------------\n", "pc1 = normalize_pc(pca_image[:, :, 0])\n", "pc2 = normalize_pc(pca_image[:, :, 1])\n", "pc5 = normalize_pc(pca_image[:, :, 4])\n", "\n", "rgb_125 = np.dstack((pc1, pc2, pc5))\n", "\n", "\n", "# -------------------------------------------------------\n", "# Composite 2: PC2 = Red, PC3 = Green, PC4 = Blue\n", "# -------------------------------------------------------\n", "pc2 = normalize_pc(pca_image[:, :, 1])\n", "pc3 = normalize_pc(pca_image[:, :, 2])\n", "pc4 = normalize_pc(pca_image[:, :, 3])\n", "\n", "rgb_234 = np.dstack((pc2, pc3, pc4))" ] }, { "cell_type": "code", "execution_count": null, "id": "8fc4768d-7009-4f08-a866-8496ef35f61f", "metadata": {}, "outputs": [], "source": [ "# -------------------------------------------------------\n", "# Plot\n", "# -------------------------------------------------------\n", "fig, axes = plt.subplots(1, 2, figsize=(18, 8))\n", "\n", "axes[0].imshow(rgb_125)\n", "axes[0].set_title(\"PCA composite: PC1 - PC2 - PC5\")\n", "axes[0].axis(\"off\")\n", "\n", "axes[1].imshow(rgb_234)\n", "axes[1].set_title(\"PCA composite: PC2 - PC3 - PC4\")\n", "axes[1].axis(\"off\")\n", "\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "b4476063-52c5-486d-825a-ebd7aa429912", "metadata": {}, "source": [ "### K-means clustering\n", "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?" ] }, { "cell_type": "code", "execution_count": null, "id": "5e729dd4-f4e4-44bd-a1cf-23cf9bb1e4fc", "metadata": {}, "outputs": [], "source": [ "#using the higher order PCA components only - see variance explained, cumulative variance\n", "n_pc = 5\n", "\n", "X_pca_reduced = X_pca[:, :n_pc]" ] }, { "cell_type": "markdown", "id": "9c3e17bf-3c9a-4824-9d11-83ed329e0f08", "metadata": {}, "source": [ "Perform K-means clustering, for example, start with 10 classes" ] }, { "cell_type": "code", "execution_count": null, "id": "bc8bfc75-b527-4050-8e88-a8096f989ad7", "metadata": {}, "outputs": [], "source": [ "n_classes = 20\n", "\n", "kmeans = KMeans(\n", " n_clusters=n_classes,\n", " random_state=42,\n", " n_init=10\n", ")\n", "\n", "clusters = kmeans.fit_predict(X_pca_reduced)" ] }, { "cell_type": "markdown", "id": "1b93464e-dcf3-4c55-9a35-7d71ae0240cc", "metadata": {}, "source": [ "Convert the classification back to an image" ] }, { "cell_type": "code", "execution_count": null, "id": "45f62940-7548-4eec-a6dc-5757945f6683", "metadata": {}, "outputs": [], "source": [ "classification = np.full(\n", " X.shape[0],\n", " -1,\n", " dtype=np.int16\n", ")\n", "\n", "classification[valid] = clusters\n", "\n", "classification = classification.reshape(\n", " nrows,\n", " ncols\n", ")" ] }, { "cell_type": "markdown", "id": "c50d7462-ab12-4137-a316-fd39ab396692", "metadata": {}, "source": [ "Plot the result" ] }, { "cell_type": "code", "execution_count": null, "id": "3519eb5c-c52d-43ca-b18e-c5c4c26d9033", "metadata": {}, "outputs": [], "source": [ "# -------------------------------------------------------\n", "# Create a discrete colormap\n", "# -------------------------------------------------------\n", "cmap = plt.get_cmap(\"tab20\", n_classes)\n", "\n", "# Boundaries centered around integer class numbers\n", "bounds = np.arange(-0.5, n_classes + 0.5, 1)\n", "\n", "norm = BoundaryNorm(\n", " bounds,\n", " cmap.N\n", ")" ] }, { "cell_type": "code", "execution_count": null, "id": "b84118df-804a-4913-a8c2-2d8af2cd508f", "metadata": {}, "outputs": [], "source": [ "# -------------------------------------------------------\n", "# Plot\n", "# -------------------------------------------------------\n", "fig, ax = plt.subplots(figsize=(12, 10))\n", "\n", "im = ax.imshow(\n", " classification,\n", " cmap=cmap,\n", " norm=norm\n", ")\n", "\n", "cbar = plt.colorbar(\n", " im,\n", " ax=ax,\n", " boundaries=bounds,\n", " ticks=np.arange(n_classes),\n", " spacing=\"proportional\", shrink=0.5,\n", ")\n", "\n", "cbar.set_label(\"K-means cluster / class\")\n", "\n", "ax.set_title(f\"Oil spill K-means classification - {n_classes} classes\")\n", "ax.axis(\"off\")\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "ffe0d53d-8a61-444d-83fb-9f4a1a8699d8", "metadata": {}, "source": [ "Export the results to an ilwis or geotif raster" ] }, { "cell_type": "code", "execution_count": null, "id": "eeb80627-1042-4a26-afc0-25af9be6fbc5", "metadata": {}, "outputs": [], "source": [ "#create empty ilwis raster \n", "Km_ilw = ilwis.RasterCoverage()\n", "defNumr = ilwis.DataDefinition(ilwis.NumericDomain('code=value'), ilwis.NumericRange(0, 255, 1))\n", "Km_ilw.setDataDef(defNumr)\n", "Km_ilw.setSize(ilwis.Size(6374, 4883, 1)) \n", "Km_ilw.setGeoReference(s2_date2.geoReference())" ] }, { "cell_type": "code", "execution_count": null, "id": "ab376946-ca16-435b-8dad-930e71e73465", "metadata": {}, "outputs": [], "source": [ "#read the numpy array\n", "Km_ilw.array2raster(classification.flatten())" ] }, { "cell_type": "code", "execution_count": null, "id": "3b3397f2-53df-4277-baf2-08e1c8302c58", "metadata": {}, "outputs": [], "source": [ "#comment if you don't want to store the image in an ILWIS maplist format\n", "Km_ilw.store('oil_spill_classification.mpr')\n", "\n", "#uncomment if you want to store the image in geotif format\n", "#Km_ilw.store(\"oil_spill_classification.tif\", \"GTiff\", \"gdal\")" ] }, { "cell_type": "markdown", "id": "b7a15412-528b-41dc-bc27-ea9ff9d0f370", "metadata": {}, "source": [ "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." ] }, { "cell_type": "code", "execution_count": null, "id": "d85b7009-aa30-4583-a2a9-78a4f91bda6b", "metadata": {}, "outputs": [], "source": [ "X_class = np.column_stack([\n", " X_pca[:, :5],\n", " X_scaled[:, 2] # VV difference (0-based indexing)\n", "])\n", "\n", "kmeans = KMeans(\n", " n_clusters=15,\n", " random_state=42,\n", " n_init=10\n", ")\n", "\n", "clusters = kmeans.fit_predict(X_class)" ] }, { "cell_type": "code", "execution_count": null, "id": "89acdb09-31d0-4d20-a450-cfb3feb4cff0", "metadata": {}, "outputs": [], "source": [ "classification = np.full(\n", " X.shape[0],\n", " -1,\n", " dtype=np.int16\n", ")\n", "\n", "classification[valid] = clusters\n", "\n", "classification1 = classification.reshape(\n", " nrows,\n", " ncols\n", ")" ] }, { "cell_type": "code", "execution_count": null, "id": "d6820e0e-8ab5-47f0-8b08-12762001f268", "metadata": {}, "outputs": [], "source": [ "# -------------------------------------------------------\n", "# Plot\n", "# -------------------------------------------------------\n", "fig, ax = plt.subplots(figsize=(12, 10))\n", "\n", "im = ax.imshow(\n", " classification1,\n", " cmap=cmap,\n", " norm=norm\n", ")\n", "\n", "cbar = plt.colorbar(\n", " im,\n", " ax=ax,\n", " boundaries=bounds,\n", " ticks=np.arange(n_classes),\n", " spacing=\"proportional\", shrink=0.5,\n", ")\n", "\n", "cbar.set_label(\"K-means cluster / class\")\n", "\n", "ax.set_title(f\"Oil spill K-means classification - {n_classes} classes - VV-difference enforced\")\n", "ax.axis(\"off\")\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "412a8972-e97b-4d13-98eb-9c7e836d458d", "metadata": {}, "source": [ "Export the results to an ilwis or geotif raster" ] }, { "cell_type": "code", "execution_count": null, "id": "eef5e65d-dc33-4d50-8812-ef2abff59693", "metadata": {}, "outputs": [], "source": [ "#read the numpy array\n", "Km_ilw.array2raster(classification1.flatten())" ] }, { "cell_type": "code", "execution_count": null, "id": "2f70272f-45cc-442a-8dfd-ce0838fff92b", "metadata": {}, "outputs": [], "source": [ "#comment if you don't want to store the image in an ILWIS maplist format\n", "Km_ilw.store('oil_spill_classification1.mpr')\n", "\n", "#uncomment if you want to store the image in geotif format\n", "#Km_ilw.store(\"oil_spill_classification1.tif\", \"GTiff\", \"gdal\")" ] }, { "cell_type": "markdown", "id": "568bd330-189d-445b-a181-f92509945025", "metadata": {}, "source": [ "Main conclusions:\n", "+ Sentinel-1 VV is highly useful for detecting the oil-spill signature.\n", "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.\n", "+ The temporal Sentinel-1 information is more informative than a single VV image.\n", "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.\n", "+ Sentinel-2 provides valuable complementary information.\n", "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.\n", "+ PCA reduces the dimensionality of the combined dataset.\n", "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.\n", "+ The PCA composites can enhance the visual separation of surface features.\n", "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.\n", "+ K-means can separate several surface classes, but oil and some other ocean / cloud features may remain confused.\n", "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.\n", "+ The VV temporal difference is particularly important for resolving the oil/other class ambiguity.\n", "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." ] }, { "cell_type": "code", "execution_count": null, "id": "854e869c-14f3-4f1e-b40f-e3910c0791b4", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "id": "c9cc6849-5fcb-43af-b5a6-3e7878526284", "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 }