{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "fd4bf47a-83d4-4160-8c31-dec059c1fe54",
   "metadata": {},
   "source": [
    "## Monitoring of deforestation\n",
    "\n",
    "Notebook prepared by Ben Maathuis. ITC-University of Twente, Enschede. The Netherlands. \n",
    "\n",
    "In order to monitor the deforestation over the last 30 years of a selected area at the edge of the Amazon forest in Brazil use is made of time series of Landsat mosaic images (available till end 2025), combined with more recently acquired Sentinel-2 MSI. Satellite data is obtained using OpenEO and for each time step a NDVI is calculated which is subsequently combined to see the temporal developments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca52e35e-6088-41c3-8a45-81c48f3ccdc9",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Uncomment line below if these packages are not installed\n",
    "#!pip install contextily shapely pyproj rasterio"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca0b1722-d791-4f2f-949e-1d223756bc48",
   "metadata": {},
   "outputs": [],
   "source": [
    "#import required libraries - eventually install if not available \n",
    "import os\n",
    "import glob\n",
    "import ilwis\n",
    "import openeo\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.colors as mcolors\n",
    "from matplotlib.colors import ListedColormap\n",
    "import matplotlib.ticker as mticker\n",
    "from matplotlib.ticker import MultipleLocator\n",
    "import contextily as cx\n",
    "from shapely.geometry import box\n",
    "from pyproj import Transformer\n",
    "from rasterio.warp import reproject, calculate_default_transform, Resampling\n",
    "from rasterio.transform import from_bounds\n",
    "from rasterio.crs import CRS"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3bbf9b45-cdd2-4629-8064-e6dfe00d7df4",
   "metadata": {},
   "source": [
    "Area of interest selected (Amazon region Brazil - around Boca do Acre, situated in the southern part of the Brazilian state of Amazonas)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4322a4fe-415d-4d90-9b12-5e33e9438e41",
   "metadata": {},
   "outputs": [],
   "source": [
    "# AOI bounding box (EPSG:4326)\n",
    "xmin = -66.85\n",
    "xmax = -66.50\n",
    "ymin = -8.70\n",
    "ymax = -8.48"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "259730d9-10d8-401f-9e02-9f5903fb1c77",
   "metadata": {},
   "source": [
    "Display openstreetmap to see the location details of the area of interest specified"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c6bf5070-d42f-4cc9-89be-1290d0136f7a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Transform longitude/latitude to Web Mercator (EPSG:3857)\n",
    "transformer = Transformer.from_crs(\n",
    "    \"EPSG:4326\", \"EPSG:3857\", always_xy=True\n",
    ")\n",
    "\n",
    "xmin_m, ymin_m = transformer.transform(xmin, ymin)\n",
    "xmax_m, ymax_m = transformer.transform(xmax, ymax)\n",
    "\n",
    "# Create AOI rectangle\n",
    "aoi = box(xmin_m, ymin_m, xmax_m, ymax_m)\n",
    "\n",
    "# Plot map\n",
    "fig, ax = plt.subplots(figsize=(10, 8))\n",
    "\n",
    "# Set map extent with a small margin\n",
    "margin_x = (xmax_m - xmin_m) * 3.5\n",
    "margin_y = (ymax_m - ymin_m) * 3.5\n",
    "\n",
    "ax.set_xlim(xmin_m - margin_x, xmax_m + margin_x)\n",
    "ax.set_ylim(ymin_m - margin_y, ymax_m + margin_y)\n",
    "\n",
    "cx.add_basemap(\n",
    "    ax,\n",
    "    source=\"https://a.tile.openstreetmap.fr/hot/{z}/{x}/{y}.png\",\n",
    "    crs=\"EPSG:3857\"\n",
    ")\n",
    "\n",
    "# Draw AOI bounding box\n",
    "x, y = aoi.exterior.xy\n",
    "ax.plot(x, y, color=\"red\", linewidth=2, label=\"AOI boundary\")\n",
    "\n",
    "ax.set_title(\"Selected Amazon study area (AOI)\")\n",
    "ax.set_xlabel(\"Easting (m)\")\n",
    "ax.set_ylabel(\"Northing (m)\")\n",
    "ax.legend()\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3811758c-cb07-4618-9508-78e7e41c3437",
   "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": "fce686d3-1f94-4698-8536-a3d2c11cab3e",
   "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": "78000fc3-9d00-4059-862b-9bf5471212d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "ilwis.version()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f028fa24-0b06-4552-9062-2582dabb9dca",
   "metadata": {},
   "outputs": [],
   "source": [
    "#uncomment last line to remove any previous login credentials\n",
    "from openeo.rest.auth.config import RefreshTokenStore\n",
    "#RefreshTokenStore().remove()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c79fa580-f94b-4c66-866e-9fa45071c1c0",
   "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": "d5df3169-3dc4-4937-b378-908aa4994512",
   "metadata": {},
   "source": [
    "Get a listing of all collections available, display the id, title and description of the first 20 records, store full collection description details in csv file. Note the record numbers 5 and 17 and check their full content in the downloaded csv file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "35042062-ef28-49aa-a2d3-598d97d3e447",
   "metadata": {},
   "outputs": [],
   "source": [
    "collections = connection.list_collections()\n",
    "\n",
    "df_collections = pd.DataFrame(collections)\n",
    "display(df_collections[[\"id\", \"title\", \"description\"]].head(20))\n",
    "\n",
    "# Save to CSV\n",
    "outfile = work_dir + \"/openeo_collections.csv\"\n",
    "df_collections.to_csv(outfile, index=False, encoding=\"utf-8\")\n",
    "\n",
    "print(f\"Collections saved to: {outfile}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b60c023-3e52-49ad-ae58-58582aa63dc6",
   "metadata": {},
   "source": [
    "Get metadata description details for selected collection"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ff3480a-8cbe-4538-8249-5837f455b372",
   "metadata": {},
   "source": [
    "### The current situation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d1b336d7-4730-4fce-99b9-be0f27ebcd58",
   "metadata": {},
   "outputs": [],
   "source": [
    "connection.describe_collection(\"SENTINEL2_L2A\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bc5c3dc4-1f22-4c14-91f1-198b395c4898",
   "metadata": {},
   "source": [
    "To get an idea of the current situation an image from Sentinel-2 is going to be retrieved for the area of interest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd876109-5b04-47ae-8eb2-b1de0c36df88",
   "metadata": {},
   "outputs": [],
   "source": [
    "#note the date \n",
    "t = [\"2026-05-01\", \"2026-10-01\"]\n",
    "s2_cube = connection.load_collection(\"SENTINEL2_L2A\",\n",
    "    spatial_extent={'west': xmin, 'east': xmax, 'south': ymin, 'north': ymax,\"crs\": \"EPSG:4326\"},\n",
    "    temporal_extent= t,\n",
    "    bands=['B02', 'B03', 'B04', 'B08'],\n",
    "    max_cloud_cover=10,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0bac51a0-b5bf-4ffa-91da-aa11199a9f9e",
   "metadata": {},
   "source": [
    "Temporal aggregation - note that t = [\"2026-05-01\", \"2026-10-01\"] which is representing the dry season. \n",
    "From the selected images, perform a temporal aggregation, using the 'median_time' operation which calculates the median value for each pixel across all time steps in a data cube. It reduces the temporal dimension to produce a single image for each band. For example, if a pixel has digital numbers of 25, 60, 40, 90, 35, the median is 40. Outliers caused by clouds or shadow in your time series are likely excluded using such an approach and produces a representative composite from multiple observations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "04970cf5-016b-4bbc-80ae-5c82d2dad9d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "s2_cube_red = s2_cube.median_time()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d81313a0-9d8e-4b9b-b0ed-e7407a3a8b4a",
   "metadata": {},
   "outputs": [],
   "source": [
    "s2_cube_red.download(work_dir+\"/s2_2026.tif\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bfeb42b2-5b5b-40fe-b220-5238a787eda0",
   "metadata": {},
   "outputs": [],
   "source": [
    "#read the S2 image in ilwispy\n",
    "s2_2026 = ilwis.RasterCoverage ('s2_2026.tif')\n",
    "print(s2_2026.size())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d91df20a-669f-4036-aa9a-5412796c3dd6",
   "metadata": {},
   "outputs": [],
   "source": [
    "#stretch the rasterbands contained in s2_2026 using a loop\n",
    "multiple_stretch = []\n",
    "multiple_bands = ilwis.do('selection',s2_2026,\"rasterbands(0..3)\") \n",
    "ls = ilwis.do('linearstretch',multiple_bands, 1) #using an upper and lower data limit defined by the cumulative 1 and 99 % thresholds\n",
    "mb_stretch = ilwis.do('setvaluerange', ls, 0, 255, 1) #set ouput to byte range"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "837a2470-7983-463f-adf8-e64b9f9f1ecd",
   "metadata": {},
   "outputs": [],
   "source": [
    "#store the results as an ILWIS maplist - display the image as a color composite using the ilwis386 desktop software\n",
    "mb_stretch.store('S2_2026.mpl')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "580b508c-78cc-4c4a-bcd6-ac11ccd7f732",
   "metadata": {},
   "source": [
    "Calculate the NDVI from the Sentinel-2 spectral channels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ff5ff42e-d293-438b-a631-69d75d65d682",
   "metadata": {},
   "outputs": [],
   "source": [
    "# select required spectral channels from the initial image stack, here the visible red and NIR channels\n",
    "Nir_in = ilwis.do('selection',s2_2026,\"rasterbands(3)\")\n",
    "Red_in = ilwis.do('selection',s2_2026,\"rasterbands(2)\")\n",
    "\n",
    "S2_2026_ndvi = ilwis.do('mapcalc','(@1 - @2)/(@1 + @2)', Nir_in, Red_in)\n",
    "S2_2026_ndvi = ilwis.do('setvaluerange', S2_2026_ndvi, -1, 1, 0.001)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "78c5979b-daee-4a13-8bd5-9c392961004f",
   "metadata": {},
   "outputs": [],
   "source": [
    "#store the results as an ILWIS raster map - display the NDVI using the ilwis386 desktop software \n",
    "#check the map values, especially of the forest and the defeorested areas\n",
    "S2_2026_ndvi.store('S2_2026_ndvi.mpr')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0664d656-728e-4547-91c1-c856158b4d9c",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Transform to a numpy array for visualization using imshow\n",
    "row = S2_2026_ndvi.size().ysize \n",
    "col = S2_2026_ndvi.size().xsize\n",
    "dim = S2_2026_ndvi.size().linearSize()\n",
    "print(row, col, dim)\n",
    "\n",
    "raster_2np = np.fromiter(iter(S2_2026_ndvi), np.float64, dim) \n",
    "raster_2np = raster_2np.reshape((row, col))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0507317a-e204-4ce0-9448-74222c6b6bfe",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the NDVI derived from Sentinel-2 MSI\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(12, 8))\n",
    "\n",
    "im = ax.imshow(\n",
    "    raster_2np,\n",
    "    extent=[xmin, xmax, ymin, ymax],\n",
    "    cmap=\"RdYlGn\",\n",
    "    vmin=0,\n",
    "    vmax=1\n",
    ")\n",
    "\n",
    "# Colorbar\n",
    "cbar = plt.colorbar(im, ax=ax, shrink=0.5)\n",
    "cbar.set_label(\"NDVI\")\n",
    "\n",
    "# Axes\n",
    "ax.set_xlabel(\"Longitude (°E)\")\n",
    "ax.set_ylabel(\"Latitude (°N)\")\n",
    "ax.set_title(\"NDVI from Sentinel-2 Mosaic of \" + (t[0])+\" to \" + (t[1]))\n",
    "\n",
    "# Limit to the AOI\n",
    "ax.set_xlim(xmin, xmax)\n",
    "ax.set_ylim(ymin, ymax)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5f479366-9a09-4cec-94f7-bf9dff78e308",
   "metadata": {},
   "source": [
    "Visually compare the NDVI plot with the Sentinel-2 color composite using the ilwis386 desktop software"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b0f71f01-5d99-425e-8ca5-1353832403e7",
   "metadata": {},
   "source": [
    "### The historical deforestation developments"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "846f5e94-b00d-4d54-bda4-a8b4ceef1034",
   "metadata": {},
   "source": [
    "In order to get an idea of the deforestation developments over the last 3 decades use can be made of the Landsat image archive, especially the Landsat bimonthly mosaic. First we review the meta-data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3272059a-5ece-4ace-8c87-5927ac42c0da",
   "metadata": {},
   "outputs": [],
   "source": [
    "connection.describe_collection(\"LANDSAT_BIMONTHLY_MOSAIC\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7258a94-42c7-42f2-a475-2d4093d56338",
   "metadata": {},
   "source": [
    "As the Landsat bimonthly mosaic is available from 1997 to 2025, a timeseries with a temporal interval of 3 years, would provide a good impression of the temporal deforestation developments. Instead of downloading the individual spectral channels and calculate the NDVI on the local repository, within OpenEO the function cube.ndvi is used, which calculates the Normalized Difference Vegetation Index (NDVI) from the near-infrared (NIR) and red spectral bands of the selected satellite data cube. The median_time function is used to retrieve the final ndvi, which is subsequently downloaded for each of the temporal intervals selected."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "67912ab7-8734-448d-9a92-4b444850252e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Target years: 1997, 2000, ..., 2024 - note temporal extent from Landsatbi-monthly mosaic collection\n",
    "years = range(1997, 2025, 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2ce8788-ac5a-475a-b36c-139d0b3b34f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "#download derived NDVI (only) from Landsat mosaic collection from 1997 to 2025, interval of 3 years\n",
    "\n",
    "collection_id = \"LANDSAT_BIMONTHLY_MOSAIC\"\n",
    "\n",
    "for year in years:\n",
    "    print(f\"Processing NDVI for {year}\")\n",
    "\n",
    "    # Load May–October observations\n",
    "    cube = connection.load_collection(\n",
    "        collection_id,\n",
    "        spatial_extent={'west': (xmin), 'east': (xmax), 'south': (ymin), 'north': (ymax),\"crs\": \"EPSG:4326\"},\n",
    "        temporal_extent=[\n",
    "            f\"{year}-05-01\",\n",
    "            f\"{year}-10-01\"\n",
    "        ],\n",
    "        bands=[\"B03\", \"B04\"]\n",
    "    )\n",
    "\n",
    "    # NDVI = (NIR - Red) / (NIR + Red)\n",
    "    ndvi = cube.ndvi(\n",
    "        nir=\"B04\",\n",
    "        red=\"B03\"\n",
    "    )\n",
    "\n",
    "    # Calculate median across the selected period\n",
    "    ndvi_median = ndvi.median_time()\n",
    "\n",
    "    # Download GeoTIFF\n",
    "    outfile = os.path.join(\n",
    "        work_dir,\n",
    "        f\"NDVI_{year}_May_Oct.tif\"\n",
    "    )\n",
    "\n",
    "    ndvi_median.download(outfile)\n",
    "\n",
    "    print(f\"Saved: {outfile}\")\n",
    "\n",
    "print(\"All years processed.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "657bbd80-1c0f-41f2-a09d-1d90a2098a14",
   "metadata": {},
   "source": [
    "Load the data so it can be used in ilwispy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a4d70f29-0afd-48ab-8fae-2ba7383cbc1c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find all NDVI images\n",
    "ndvi_files = sorted(glob.glob(os.path.join(work_dir, \"NDVI_*_May_Oct.tif\")))\n",
    "\n",
    "# Load images into a dictionary\n",
    "ndvi = {}\n",
    "\n",
    "for filename in ndvi_files:\n",
    "    year = os.path.basename(filename).split(\"_\")[1]\n",
    "\n",
    "    ndvi[year] = ilwis.RasterCoverage(filename)\n",
    "\n",
    "    print(f\"Loaded NDVI image for {year}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "505af184-73a1-4ba4-bb59-37aadb864bd3",
   "metadata": {},
   "outputs": [],
   "source": [
    "#assign appropriate name \n",
    "ndvi_1997 = ndvi[\"1997\"]\n",
    "ndvi_2000 = ndvi[\"2000\"]\n",
    "ndvi_2003 = ndvi[\"2003\"]\n",
    "ndvi_2006 = ndvi[\"2006\"]\n",
    "ndvi_2009 = ndvi[\"2009\"]\n",
    "ndvi_2012 = ndvi[\"2012\"]\n",
    "ndvi_2015 = ndvi[\"2015\"]\n",
    "ndvi_2018 = ndvi[\"2018\"]\n",
    "ndvi_2021 = ndvi[\"2021\"]\n",
    "ndvi_2024 = ndvi[\"2024\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19d41104-ebd5-4f10-ace2-bcc0fcee34ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(ndvi_1997.size())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a8b26efd-1d0f-4b38-9064-9668bf3d48f7",
   "metadata": {},
   "outputs": [],
   "source": [
    "#store the NDVI as an ILWIS map - display the map using the ilwis386 desktop software\n",
    "ndvi_1997.store('landsat_ndvi1997.mpr')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c5b33f26-830f-4597-9190-1676f9958c24",
   "metadata": {},
   "source": [
    "Inspect the NDVI map and determine a NDVI value which would separate the forest from the cleared areas. Take special care of mixed pisels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "266333d7-f75e-4647-9330-3dd3af8490bd",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Transform to a numpy array for visualization using imshow for selected interval\n",
    "\n",
    "interval = ndvi_1997\n",
    "\n",
    "row = interval.size().ysize \n",
    "col = interval.size().xsize\n",
    "dim = interval.size().linearSize()\n",
    "print(row, col, dim)\n",
    "\n",
    "raster_2npl = np.fromiter(iter(interval), np.float64, dim) \n",
    "raster_2npl = raster_2npl.reshape((row, col))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "973e0c78-6477-4875-b256-c2dca3ed83a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a plot of the NDVI from the Landsat mosaic of 1997\n",
    "fig, ax = plt.subplots(figsize=(12, 8))\n",
    "\n",
    "im = ax.imshow(\n",
    "    raster_2npl,\n",
    "    extent=[xmin, xmax, ymin, ymax],\n",
    "    cmap=\"RdYlGn\",\n",
    "    vmin=0.5,\n",
    "    vmax=1\n",
    ")\n",
    "\n",
    "# Colorbar\n",
    "cbar = plt.colorbar(im, ax=ax, shrink=0.5)\n",
    "cbar.set_label(\"NDVI\")\n",
    "\n",
    "fname = os.path.basename(str(interval))\n",
    "year_title = fname.split(\"_\")[1]\n",
    "\n",
    "\n",
    "# Axes\n",
    "ax.set_xlabel(\"Longitude (°E)\")\n",
    "ax.set_ylabel(\"Latitude (°N)\")\n",
    "ax.set_title(\"NDVI from Landsat Mosaic of \" + year_title) \n",
    "\n",
    "# Limit to the AOI\n",
    "ax.set_xlim(xmin, xmax)\n",
    "ax.set_ylim(ymin, ymax)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dcdc314d-9280-4fe6-ac99-b70633093283",
   "metadata": {},
   "source": [
    "In the equations used below the deforestation for each time step is retrieved and combined into a single map. Here a cutoff value as derived from visual inspection of the Landsat ndvi map to separate forest from cleared - non forest pixels used is 0.815, if less or equal then this value, the pixel is assigned deforested. If you have defined another cutoff value it can be replaced in the field below.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e58539e-1513-4b5d-85c0-191b11f8f6cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "cutoff = 0.815 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f713a42c-ab23-407a-b2a3-7fed53b53a6c",
   "metadata": {},
   "outputs": [],
   "source": [
    "#decade-1\n",
    "#Note the use of an f-string to insert the Python variable into the ilwispy map calculation formula below\n",
    "change1 = ilwis.do('mapcalc', f'iff(@1<={cutoff}, 1, iff(@2<={cutoff},2, iff(@3<={cutoff}, 3, iff(@4<={cutoff}, 4, 5))))',ndvi_1997, ndvi_2000, ndvi_2003, ndvi_2006) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2c381bec-a799-46d5-a9f0-b657f7ac5e21",
   "metadata": {},
   "outputs": [],
   "source": [
    "#to check your results display the map using the ilwis386 desktop software\n",
    "change1.store('change_1997_2006_cl.mpr')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c50da117-0490-4fde-8e7c-28d39c8ecb87",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Transform to a numpy array for visualization using imshow for selected interval\n",
    "\n",
    "row = change1.size().ysize \n",
    "col = change1.size().xsize\n",
    "dim = change1.size().linearSize()\n",
    "print(row, col, dim)\n",
    "\n",
    "raster_2npl = np.fromiter(iter(change1), np.float64, dim) \n",
    "raster_2npl = raster_2npl.reshape((row, col))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c9117491-ca83-484d-a0e7-57b9fbba2aec",
   "metadata": {},
   "outputs": [],
   "source": [
    "#get the class statistics\n",
    "\n",
    "# Count pixels for each class\n",
    "values, counts = np.unique(\n",
    "    raster_2npl[np.isfinite(raster_2npl)],\n",
    "    return_counts=True\n",
    ")\n",
    "\n",
    "# Create pandas DataFrame\n",
    "df_classes = pd.DataFrame({\n",
    "    \"Class\": values,\n",
    "    \"Pixel_count\": counts\n",
    "})\n",
    "\n",
    "# Display without index\n",
    "display(df_classes.style.hide(axis=\"index\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8ba3ecab-7dcd-45f1-bc5f-2a1399e108ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a plot\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(12, 8))\n",
    "\n",
    "# Show every other coordinate tick\n",
    "ax.xaxis.set_major_locator(MultipleLocator(0.05)) #unit is degree\n",
    "ax.yaxis.set_major_locator(MultipleLocator(0.05))\n",
    "\n",
    "\n",
    "# Five classes: 1, 2, 3, 4, 5\n",
    "cmap = plt.get_cmap(\"jet\", 5)\n",
    "\n",
    "# Each integer class is centered in its color interval\n",
    "bounds = [0.5, 1.5, 2.5, 3.5, 4.5, 5.5]\n",
    "norm = mcolors.BoundaryNorm(bounds, cmap.N)\n",
    "\n",
    "im = ax.imshow(\n",
    "    raster_2npl,\n",
    "    extent=[xmin, xmax, ymin, ymax],\n",
    "    cmap=cmap,\n",
    "    norm=norm,\n",
    "    interpolation=\"nearest\"\n",
    ")\n",
    "\n",
    "\n",
    "# Color bar\n",
    "cbar = plt.colorbar(\n",
    "    im,\n",
    "    ax=ax,\n",
    "    boundaries=bounds,\n",
    "    ticks=[1, 2, 3, 4, 5],\n",
    "    spacing=\"uniform\",\n",
    "    shrink=0.4\n",
    ")\n",
    "\n",
    "cbar.set_ticklabels([\n",
    "    \"till-1997\",\n",
    "    \"till-2001\",\n",
    "    \"till-2003\",\n",
    "    \"till-2006\",\n",
    "    \"forest\"\n",
    "])\n",
    "cbar.set_label(\"Deforestation progress\")\n",
    "\n",
    "# Axes\n",
    "ax.set_xlabel(\"Longitude (°E)\")\n",
    "ax.set_ylabel(\"Latitude (°N)\")\n",
    "ax.set_title(\"Deforestation from Landsat Mosaic between 1997–2006\")\n",
    "\n",
    "ax.set_xlim(xmin, xmax)\n",
    "ax.set_ylim(ymin, ymax)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5df7d07a-4c63-438f-818c-9877947a083f",
   "metadata": {},
   "outputs": [],
   "source": [
    "#decade-2\n",
    "#calculate the deforestation change for the second time period\n",
    "change2 = ilwis.do('mapcalc', f'iff(@1<={cutoff}, 1, iff(@2<={cutoff},2, iff(@3<={cutoff}, 3, iff(@4<={cutoff}, 4, 5))))',ndvi_2009, ndvi_2012, ndvi_2015, ndvi_2018) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0d6f15b-e139-4d4b-a424-bd2c1e2b290b",
   "metadata": {},
   "outputs": [],
   "source": [
    "change2.store('change_2009_2018_cl.mpr')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0fa1fa1e-a552-422f-b3a3-5309e6e835c6",
   "metadata": {},
   "source": [
    "Review the change map of 2009-2018 and note in the equations below that the two change maps are combined. Check the change map of to whole period as well"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a44c158a-ad99-48a2-899f-948ed159994b",
   "metadata": {},
   "outputs": [],
   "source": [
    "combine = ilwis.do('mapcalc', 'iff(@1<=4, @1, @2+4)',change1, change2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ea3bdd63-f3cb-47b7-989a-e94e1f332586",
   "metadata": {},
   "outputs": [],
   "source": [
    "combine.store('change_1997_2018_cl.mpr')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e6863929-000a-42d2-825d-dc7ddff73f4a",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Transform to a numpy array for visualization using imshow for selected interval\n",
    "\n",
    "row = combine.size().ysize \n",
    "col = combine.size().xsize\n",
    "dim = combine.size().linearSize()\n",
    "print(row, col, dim)\n",
    "\n",
    "raster_2npl = np.fromiter(iter(combine), np.float64, dim) \n",
    "raster_2npl = raster_2npl.reshape((row, col))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "647a2593-8470-4e71-a96e-b2da0d19d98d",
   "metadata": {},
   "outputs": [],
   "source": [
    "#get the class statistics\n",
    "\n",
    "# Count pixels for each class\n",
    "values, counts = np.unique(\n",
    "    raster_2npl[np.isfinite(raster_2npl)],\n",
    "    return_counts=True\n",
    ")\n",
    "\n",
    "# Create pandas DataFrame\n",
    "df_classes = pd.DataFrame({\n",
    "    \"Class\": values,\n",
    "    \"Pixel_count\": counts\n",
    "})\n",
    "\n",
    "# Display without index\n",
    "display(df_classes.style.hide(axis=\"index\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "edce779f-81be-4370-869d-d20443bfc084",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a plot\n",
    "fig, ax = plt.subplots(figsize=(12, 8))\n",
    "\n",
    "# Show every other coordinate tick\n",
    "ax.xaxis.set_major_locator(MultipleLocator(0.05)) #unit is degree\n",
    "ax.yaxis.set_major_locator(MultipleLocator(0.05))\n",
    "\n",
    "\n",
    "# Five classes: 1, 2, 3, 4, 5, 6, 7, 8, 9 \n",
    "cmap = plt.get_cmap(\"jet\", 9)\n",
    "\n",
    "# Each integer class is centered in its color interval\n",
    "bounds = [0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5]\n",
    "norm = mcolors.BoundaryNorm(bounds, cmap.N)\n",
    "\n",
    "im = ax.imshow(\n",
    "    raster_2npl,\n",
    "    extent=[xmin, xmax, ymin, ymax],\n",
    "    cmap=cmap,\n",
    "    norm=norm,\n",
    "    interpolation=\"nearest\"\n",
    ")\n",
    "\n",
    "\n",
    "# Color bar\n",
    "cbar = plt.colorbar(\n",
    "    im,\n",
    "    ax=ax,\n",
    "    boundaries=bounds,\n",
    "    ticks=[1, 2, 3, 4, 5, 6, 7, 8, 9],\n",
    "    spacing=\"uniform\",\n",
    "    shrink=0.4\n",
    ")\n",
    "\n",
    "cbar.set_ticklabels([\n",
    "    \"till-1997\",\n",
    "    \"till-2001\",\n",
    "    \"till-2003\",\n",
    "    \"till-2006\",\n",
    "    \"till-2009\",\n",
    "    \"till-2012\",\n",
    "    \"till-2015\",\n",
    "    \"till-2018\",\n",
    "    \"forest\"\n",
    "])\n",
    "cbar.set_label(\"Deforestation progress\")\n",
    "\n",
    "# Axes\n",
    "ax.set_xlabel(\"Longitude (°E)\")\n",
    "ax.set_ylabel(\"Latitude (°N)\")\n",
    "ax.set_title(\"Deforestation from Landsat Mosaic between 1997–2018\")\n",
    "\n",
    "ax.set_xlim(xmin, xmax)\n",
    "ax.set_ylim(ymin, ymax)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd02584f-23ed-417a-b682-eae3ecd9aef8",
   "metadata": {},
   "source": [
    "In order to calculate the change till now, the recent image from Sentinel-2 has to be resampled to Landsat - review also the meta-data given above. Note that Sentinel-2 image retrieved has a spatial resolution of 10 meters, different number of rows and columns, other projection, etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2892a586-6549-4dd1-aec1-fa1b7be7a4f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Retrieve the Sentinel-2 derived NDVI georeference information \n",
    "print(S2_2026_ndvi.size())\n",
    "print(S2_2026_ndvi.envelope())\n",
    "\n",
    "print()\n",
    "\n",
    "csy = S2_2026_ndvi.coordinateSystem()\n",
    "print(\"Coordinate system:\", csy.toWKT())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c9c499c8-3991-4216-92a4-1f99dfaf9cfc",
   "metadata": {},
   "source": [
    "Resample Sentinel-2 image to Landsat - also have a look at the georeference and projection information"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b3ae7f9-a2a2-4d9f-a052-1772ec050f70",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Retrieve Landsat georeference information \n",
    "print(ndvi_1997.size())\n",
    "print(ndvi_1997.envelope())\n",
    "grfLandsat = ndvi_1997.geoReference()\n",
    "print(grfLandsat)\n",
    "csy = grfLandsat.coordinateSystem()\n",
    "\n",
    "print(\"Coordinate system:\", csy.toWKT())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d885702a-9835-473f-897c-ca26d34a94e7",
   "metadata": {},
   "outputs": [],
   "source": [
    "#conduct the resampling operation for the Sentinel-2 NDVI\n",
    "s2ndvi2026_res = ilwis.do('resample', S2_2026_ndvi, grfLandsat,'nearestneighbour')\n",
    "\n",
    "print(s2ndvi2026_res.size())\n",
    "print(s2ndvi2026_res.coordinateSystem())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2631d03c-e0d1-4b88-a21a-0f6d691e2beb",
   "metadata": {},
   "outputs": [],
   "source": [
    "#check resulting resampled ndvi map derived from Sentinel-2\n",
    "s2ndvi2026_res.store(\"s2_ndvi2026_res.mpr\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1e500338-0a5e-4b9e-8f87-64248edb9599",
   "metadata": {},
   "outputs": [],
   "source": [
    "#calculate the changes during the final decade, also using the resamped Sentinel-2 ndvi from 2026\n",
    "change3 = ilwis.do('mapcalc', f'iff(@1<={cutoff}, 1, iff(@2<={cutoff},2, iff(@3<={cutoff}, 3, 4)))',ndvi_2021, ndvi_2024, s2ndvi2026_res) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3f61548e-de9a-491f-8bce-18e22090e712",
   "metadata": {},
   "outputs": [],
   "source": [
    "change3.store('change_2021_2026_cl.mpr')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "63ab97c0-9d6f-4596-8158-674ef7d859db",
   "metadata": {},
   "source": [
    "Combine all timesteps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "823b2b2d-1adb-400e-a746-5f1ea8fc1999",
   "metadata": {},
   "outputs": [],
   "source": [
    "combine_all = ilwis.do('mapcalc', 'iff(@1<=8, @1, @2+8)',combine, change3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9011e678-9f8c-4e8a-81f2-a0092bbe4fb7",
   "metadata": {},
   "outputs": [],
   "source": [
    "#also check the results obtained using ilwis desktop\n",
    "combine_all.store('change_1997_2026_cl.mpr')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03954ea6-14fc-47d4-b62e-576968d21ee9",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Transform to a numpy array for visualization using imshow for selected interval\n",
    "\n",
    "row = combine_all.size().ysize \n",
    "col = combine_all.size().xsize\n",
    "dim = combine_all.size().linearSize()\n",
    "print(row, col, dim)\n",
    "\n",
    "raster_2npl = np.fromiter(iter(combine_all), np.float64, dim) \n",
    "raster_2npl = raster_2npl.reshape((row, col))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5bcbe77c-8867-4bc5-a956-d09e6201ff28",
   "metadata": {},
   "outputs": [],
   "source": [
    "#get the class statistics for the whole time series\n",
    "\n",
    "# Count pixels for each class\n",
    "values, counts = np.unique(\n",
    "    raster_2npl[np.isfinite(raster_2npl)],\n",
    "    return_counts=True\n",
    ")\n",
    "\n",
    "# Create pandas DataFrame\n",
    "df_classes = pd.DataFrame({\n",
    "    \"Class\": values,\n",
    "    \"Pixel_count\": counts\n",
    "})\n",
    "\n",
    "# Assuming the pixel counts are in the second column\n",
    "df_classes[\"Year\"] = 1997 + (df_classes.index) * 3\n",
    "\n",
    "# Set class 11 to 2026\n",
    "df_classes.loc[df_classes[\"Class\"] == 11, \"Year\"] = 2026\n",
    "\n",
    "# Remove class 12 from the dataframe\n",
    "df_classes = df_classes[df_classes[\"Class\"] != 12].copy()\n",
    "\n",
    "# Set Year to integer\n",
    "df_classes[\"Year\"] = df_classes[\"Year\"].astype(int)\n",
    "\n",
    "\n",
    "# Convert pixel size of 30 by 30 meters to hectares / km2\n",
    "# Area of one pixel: 30 m * 30 m = 900 m}2\n",
    "# Total area in square meters:  1354 pixels = 1218600 m^2\n",
    "# Conversion to Hectares: Divide by 10000 (1218600 / 10000) = 121.86 ha\n",
    "# Conversion to km²: Divide by 1000000 (1218600 / 1000000) = 1.2186 km²\n",
    "# or\n",
    "# 1 pixel = 900 m^2 and 1 hectare = 10,000 m^2\n",
    "# 900/10,000 = 0.09 hectares per pixel\n",
    "\n",
    "df_classes[\"Area_ha\"] = df_classes[\"Pixel_count\"] * 0.09\n",
    "df_classes[\"Area_km2\"] = df_classes[\"Pixel_count\"] * 0.0009\n",
    "\n",
    "# Display without index\n",
    "display(df_classes.style.hide(axis=\"index\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "87caeb90-83ea-4818-99a6-2fff1b40986a",
   "metadata": {},
   "outputs": [],
   "source": [
    "#get the total area deforested over time\n",
    "deforested = df_classes[\"Area_km2\"].sum()\n",
    "\n",
    "print(f\"Area deforested : {deforested:,.2f} km²\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "57b9f179-74f3-4256-940f-a07ad99379e7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sort by year\n",
    "df_plot = df_classes.sort_values(\"Year\")\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10, 5))\n",
    "\n",
    "ax.bar(df_plot[\"Year\"], df_plot[\"Area_km2\"], width=1.0)\n",
    "\n",
    "ax.set_xlabel(\"Year\")\n",
    "ax.set_ylabel(\"Area (km²)\")\n",
    "ax.set_title(\"Deforestation from 1997 to 2026\")\n",
    "ax.set_xticks(df_plot[\"Year\"])\n",
    "ax.tick_params(axis=\"x\", rotation=45)\n",
    "\n",
    "ax.grid(axis=\"y\", linestyle=\"--\", alpha=0.5)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "adb7f8d8-a22e-456a-a148-9f908ff8a981",
   "metadata": {},
   "source": [
    "Calculate percentage of deforestation over the time series period for the area of interest "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "38fb7b08-28ee-445b-ad08-536e8ba853ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "# AoI according to Landsat subset\n",
    "pixels = 1440 * 880 #number of columns / rows of AoI\n",
    "pixel_area = 900  # m² per pixel for Landsat\n",
    "total_area_m2 = pixels * pixel_area\n",
    "total_area_km2 = total_area_m2 / 1_000_000\n",
    "\n",
    "print(f\"Total Area: {total_area_km2:,.2f} km²\")\n",
    "print(f\"Area Deforested : {deforested:,.2f} km²\")\n",
    "\n",
    "percentage_deforest = (deforested / total_area_km2)*100\n",
    "print(f\"Percentage Deforestation in AoI : {percentage_deforest:,.2f} %\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "57c9a19b-d7c4-4494-b8c1-b96ff337f3cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot final results\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(12, 8))\n",
    "\n",
    "# Show every other coordinate tick\n",
    "ax.xaxis.set_major_locator(MultipleLocator(0.05)) #unit is degree\n",
    "ax.yaxis.set_major_locator(MultipleLocator(0.05))\n",
    "\n",
    "# twelve classes: 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12 \n",
    "cmap = plt.get_cmap(\"jet\", 12)\n",
    "\n",
    "# Each integer class is centered in its color interval\n",
    "bounds = [0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5]\n",
    "norm = mcolors.BoundaryNorm(bounds, cmap.N)\n",
    "\n",
    "im = ax.imshow(\n",
    "    raster_2npl,\n",
    "    extent=[xmin, xmax, ymin, ymax],\n",
    "    cmap=cmap,\n",
    "    norm=norm,\n",
    "    interpolation=\"nearest\"\n",
    ")\n",
    "\n",
    "# Color bar\n",
    "cbar = plt.colorbar(\n",
    "    im,\n",
    "    ax=ax,\n",
    "    boundaries=bounds,\n",
    "    ticks=[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12],\n",
    "    spacing=\"uniform\",\n",
    "    shrink=0.4\n",
    ")\n",
    "\n",
    "cbar.set_ticklabels([\n",
    "    \"till-1997\",\n",
    "    \"till-2001\",\n",
    "    \"till-2003\",\n",
    "    \"till-2006\",\n",
    "    \"till-2009\",\n",
    "    \"till-2012\",\n",
    "    \"till-2015\",\n",
    "    \"till-2018\",\n",
    "    \"till-2021\",\n",
    "    \"till-2024\",\n",
    "    \"till-2026\",\n",
    "    \"forest\"\n",
    "])\n",
    "cbar.set_label(\"Deforestation progress\")\n",
    "\n",
    "# Axes\n",
    "ax.set_xlabel(\"Longitude (°E)\")\n",
    "ax.set_ylabel(\"Latitude (°N)\")\n",
    "ax.set_title(\"Deforestation as observed from Landsat Mosaic and Sentinel-2 between 1997–2026\")\n",
    "\n",
    "ax.set_xlim(xmin, xmax)\n",
    "ax.set_ylim(ymin, ymax)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "469897e1-5986-42c1-a4fe-9cea1df8c417",
   "metadata": {},
   "source": [
    "### Checking your results using openstreetmap\n",
    "\n",
    "+ assuming that all deforested areas are depicted on the final Sentinel-2 image\n",
    "+ then in the final binary raster assign 0 to forest and 1 to deforested pixels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c96fb5f9-0607-49fb-968b-c9cf814667be",
   "metadata": {},
   "outputs": [],
   "source": [
    "change_mask = ilwis.do('mapcalc', f'iff(@1<={cutoff}, 1, 0)',s2ndvi2026_res) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf4e1e0b-addd-4a02-9aea-263ce5361a65",
   "metadata": {},
   "source": [
    "As can be observed from the deforestation time series map displayed above in the forest and along the clearing edges there is an issue of mixture. In order to deal with this mixture we apply a 3 by 3 mean filter and set a threshold if we still regard this pixel as forested, e.g. if within the filter kernel 6 out of the 9 pixels (the majority) are still marked as 1, then the average is 6/9=0.666, so if we set a threshold on 0.6 then based upon the majority of surrounding pixels composted of forest, the pixel will be assignd 1, else 0. This will substantially reduce the spickled appearance on the mask, especially caused by the 2026 (Sentinel-2) timestep as the cufoff value is likely not representative as it was based on the Landsat bi-monthly mosiac."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31525895-5cd4-4f60-b163-4037a4bb9410",
   "metadata": {},
   "outputs": [],
   "source": [
    "filter_avg = ilwis.do('linearrasterfilter', change_mask, 'average')\n",
    "filter_mask = ilwis.do('mapcalc', 'iff(@1>=0.75, 1, 0)',filter_avg) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f0f658e5-a08d-444d-a035-7fd9c29e6061",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Transform to a numpy array for visualization using imshow for selected interval\n",
    "\n",
    "row = filter_mask.size().ysize \n",
    "col = filter_mask.size().xsize\n",
    "dim = filter_mask.size().linearSize()\n",
    "print(row, col, dim)\n",
    "\n",
    "filter_mask_2npl = np.fromiter(iter(filter_mask), np.float64, dim) \n",
    "filter_mask_2npl = filter_mask_2npl.reshape((row, col))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "583a040c-e33f-4719-a6b4-e5e08d6f6143",
   "metadata": {},
   "outputs": [],
   "source": [
    "#display the numpy array created in the previous step using matplotlib \n",
    "fig = plt.figure(figsize =(12, 8))\n",
    "plt.imshow(filter_mask_2npl)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7c0396b2-59e4-4af2-866a-fe7179d7e14c",
   "metadata": {},
   "outputs": [],
   "source": [
    "#get the statistics for the filtered total deforested area\n",
    "\n",
    "# Count pixels for each class\n",
    "values, counts = np.unique(\n",
    "    filter_mask_2npl[np.isfinite(filter_mask_2npl)],\n",
    "    return_counts=True\n",
    ")\n",
    "\n",
    "# Create pandas DataFrame\n",
    "df_2classes_only = pd.DataFrame({\n",
    "    \"Class\": values,\n",
    "    \"Pixel_count\": counts\n",
    "})\n",
    "\n",
    "# Display without index\n",
    "display(df_2classes_only.style.hide(axis=\"index\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "74c10384-5675-467f-9d06-6eed9b4e75dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "#get the filtered / corrected total area deforested over time\n",
    "total_pixels = df_2classes_only[\"Pixel_count\"].sum()\n",
    "\n",
    "class1_percentage = (\n",
    "    df_2classes_only.loc[df_2classes_only[\"Class\"] == 1, \"Pixel_count\"].sum()\n",
    "    / total_pixels\n",
    ") * 100\n",
    "\n",
    "print(f\"Percentage of class 1: {class1_percentage:.2f}%\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6f2381aa-7743-4da3-8d86-36c631b35d7d",
   "metadata": {},
   "outputs": [],
   "source": [
    "#display results over openstreetsmap\n",
    "\n",
    "# AOI bounds (EPSG:4326)\n",
    "xmin = -66.85\n",
    "xmax = -66.50\n",
    "ymin = -8.70\n",
    "ymax = -8.48\n",
    "\n",
    "# Final binary raster: 0 = forest, 1 = deforested\n",
    "# final_map should have shape (rows, columns)\n",
    "height, width = raster_2npl.shape\n",
    "\n",
    "# Source raster georeferencing\n",
    "src_transform = from_bounds(\n",
    "    xmin, ymin, xmax, ymax, width, height\n",
    ")\n",
    "src_crs = CRS.from_epsg(4326)\n",
    "dst_crs = CRS.from_epsg(3857)\n",
    "\n",
    "# Calculate destination grid\n",
    "dst_transform, dst_width, dst_height = calculate_default_transform(\n",
    "    src_crs, dst_crs, width, height,\n",
    "    xmin, ymin, xmax, ymax\n",
    ")\n",
    "\n",
    "# Reproject categorical raster using nearest-neighbour resampling\n",
    "map_3857 = np.full((dst_height, dst_width), -1, dtype=np.int16)\n",
    "\n",
    "reproject(\n",
    "    source=filter_mask_2npl.astype(np.int16),\n",
    "    destination=map_3857,\n",
    "    src_transform=src_transform,\n",
    "    src_crs=src_crs,\n",
    "    dst_transform=dst_transform,\n",
    "    dst_crs=dst_crs,\n",
    "    src_nodata=-1,\n",
    "    dst_nodata=-1,\n",
    "    resampling=Resampling.nearest\n",
    ")\n",
    "\n",
    "# Transform AOI bounding box to Web Mercator\n",
    "transformer = Transformer.from_crs(\n",
    "    \"EPSG:4326\", \"EPSG:3857\", always_xy=True\n",
    ")\n",
    "\n",
    "xmin_m, ymin_m = transformer.transform(xmin, ymin)\n",
    "xmax_m, ymax_m = transformer.transform(xmax, ymax)\n",
    "\n",
    "aoi = box(xmin_m, ymin_m, xmax_m, ymax_m)\n",
    "\n",
    "# Plot\n",
    "fig, ax = plt.subplots(figsize=(12, 9))\n",
    "\n",
    "margin_x = (xmax_m - xmin_m) * 0.15\n",
    "margin_y = (ymax_m - ymin_m) * 0.15\n",
    "\n",
    "ax.set_xlim(xmin_m - margin_x, xmax_m + margin_x)\n",
    "ax.set_ylim(ymin_m - margin_y, ymax_m + margin_y)\n",
    "\n",
    "# OpenStreetMap basemap\n",
    "cx.add_basemap(\n",
    "    ax,\n",
    "    source=\"https://a.tile.openstreetmap.fr/hot/{z}/{x}/{y}.png\",\n",
    "    crs=\"EPSG:3857\"\n",
    ")\n",
    "\n",
    "# Transparent forest, red deforested areas, transparent outside raster\n",
    "cmap = ListedColormap([\"none\", \"red\"])\n",
    "cmap.set_bad(\"none\")\n",
    "\n",
    "plot_map = np.ma.masked_where(\n",
    "    (map_3857 != 0) & (map_3857 != 1),\n",
    "    map_3857\n",
    ")\n",
    "\n",
    "# Make forest pixels transparent so the OSM map remains visible\n",
    "plot_map = np.ma.masked_where(plot_map == 0, plot_map)\n",
    "\n",
    "left = dst_transform.c\n",
    "top = dst_transform.f\n",
    "right = left + dst_width * dst_transform.a\n",
    "bottom = top + dst_height * dst_transform.e\n",
    "\n",
    "ax.imshow(\n",
    "    plot_map,\n",
    "    extent=[left, right, bottom, top],\n",
    "    origin=\"upper\",\n",
    "    cmap=ListedColormap([\"red\"]),\n",
    "    vmin=0,\n",
    "    vmax=1,\n",
    "    interpolation=\"nearest\",\n",
    "    alpha=0.7,\n",
    "    zorder=2\n",
    ")\n",
    "\n",
    "# AOI boundary\n",
    "x, y = aoi.exterior.xy\n",
    "ax.plot(x, y, color=\"yellow\", linewidth=2, label=\"AOI boundary\")\n",
    "\n",
    "ax.set_title(\"Amazon deforestation map\")\n",
    "ax.set_xlabel(\"Easting (m)\")\n",
    "ax.set_ylabel(\"Northing (m)\")\n",
    "ax.legend()\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d1165102-a9f3-4a83-a21d-f55c73f06e68",
   "metadata": {},
   "source": [
    "Some stricking differences can be observed, therefor also compare your results with the (latest) Sentinel-2 satellite image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f877218f-0929-4d3e-a86b-cdb60dc5ef01",
   "metadata": {},
   "outputs": [],
   "source": [
    "#conduct the resampling operation for the Sentinel-2 MSI images\n",
    "s2_2026_res = ilwis.do('resample', mb_stretch, grfLandsat,'nearestneighbour')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ffaf552a-e93b-43b3-8381-e82d3102e8fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "#load the individual spectral channels of the resampled Sentinel-2 image\n",
    "Blues = ilwis.do('selection',s2_2026_res,\"rasterbands(0)\")\n",
    "Greens = ilwis.do('selection',s2_2026_res,\"rasterbands(1)\")\n",
    "Reds = ilwis.do('selection',s2_2026_res,\"rasterbands(2)\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "af904d80-b73f-4b5e-80a8-baf98223f709",
   "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": "bb147cf9-bb12-472c-8430-48fa0448d2ea",
   "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": "123de1a9-44f5-4caf-a025-c281e0fa678a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot Sentinel-2 composite and deforestation mask\n",
    "fig, ax = plt.subplots(figsize=(12, 8))\n",
    "\n",
    "# Background: Sentinel-2 colour composite\n",
    "ax.imshow(\n",
    "    ncol,\n",
    "    extent=[xmin, xmax, ymin, ymax],\n",
    "    origin=\"upper\"\n",
    ")\n",
    "\n",
    "# Mask forest pixels (0), retaining only deforestation pixels (1)\n",
    "deforestation = np.ma.masked_equal(filter_mask_2npl, 0)\n",
    "\n",
    "# Transparent background for masked pixels\n",
    "cmap = ListedColormap([\"yellow\"])\n",
    "cmap.set_bad((0, 0, 0, 0))\n",
    "\n",
    "# Overlay deforestation in semi-transparent yellow\n",
    "ax.imshow(\n",
    "    deforestation,\n",
    "    extent=[xmin, xmax, ymin, ymax],\n",
    "    origin=\"upper\",\n",
    "    cmap=cmap,\n",
    "    vmin=1,\n",
    "    vmax=1,\n",
    "    alpha=0.5,\n",
    "    interpolation=\"nearest\"\n",
    ")\n",
    "\n",
    "# Axes and title\n",
    "ax.set_xlabel(\"Longitude (°E)\")\n",
    "ax.set_ylabel(\"Latitude (°N)\")\n",
    "ax.set_title(\"Resampled Sentinel-2 composite – 2026 with deforestation mask\")\n",
    "\n",
    "ax.set_xlim(xmin, xmax)\n",
    "ax.set_ylim(ymin, ymax)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9e85a8f3-0146-4346-8adc-9c481883cbf7",
   "metadata": {},
   "source": [
    "For the overlay of the deforestation map in semi-transparent yellow you can modify the alpha value to make it more or less transparent to check the results obtained. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f47c63bf-af4d-4643-a7ed-461e6727b582",
   "metadata": {},
   "source": [
    "The Landsat bi-monthly mosaic time series in conjunction with the more recent Sentinel-2 capability show a lot of potential to conduct change detection studies over the last 3 decades. But is should be noted that spectral response can / is different between the sensors so a single threshold to define the ndvi threshold which is used to determine if a pixel is deforested might not be the correct procedure!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4177571e-fdf1-4ae6-971e-db6f0a877b43",
   "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
}