#Uncomment line below if these packages are not installed
#!pip install contextily shapely pyproj rasterio

#import required libraries - eventually install if not available 
import os
import glob
import ilwis
import openeo
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
from matplotlib.colors import ListedColormap
import matplotlib.ticker as mticker
from matplotlib.ticker import MultipleLocator
import contextily as cx
from shapely.geometry import box
from pyproj import Transformer
from rasterio.warp import reproject, calculate_default_transform, Resampling
from rasterio.transform import from_bounds
from rasterio.crs import CRS

# AOI bounding box (EPSG:4326)
xmin = -66.85
xmax = -66.50
ymin = -8.70
ymax = -8.48

# Transform longitude/latitude to Web Mercator (EPSG:3857)
transformer = Transformer.from_crs(
    "EPSG:4326", "EPSG:3857", always_xy=True
)

xmin_m, ymin_m = transformer.transform(xmin, ymin)
xmax_m, ymax_m = transformer.transform(xmax, ymax)

# Create AOI rectangle
aoi = box(xmin_m, ymin_m, xmax_m, ymax_m)

# Plot map
fig, ax = plt.subplots(figsize=(10, 8))

# Set map extent with a small margin
margin_x = (xmax_m - xmin_m) * 3.5
margin_y = (ymax_m - ymin_m) * 3.5

ax.set_xlim(xmin_m - margin_x, xmax_m + margin_x)
ax.set_ylim(ymin_m - margin_y, ymax_m + margin_y)

cx.add_basemap(
    ax,
    source="https://a.tile.openstreetmap.fr/hot/{z}/{x}/{y}.png",
    crs="EPSG:3857"
)

# Draw AOI bounding box
x, y = aoi.exterior.xy
ax.plot(x, y, color="red", linewidth=2, label="AOI boundary")

ax.set_title("Selected Amazon study area (AOI)")
ax.set_xlabel("Easting (m)")
ax.set_ylabel("Northing (m)")
ax.legend()
plt.tight_layout()
plt.show()

#within the notebook folder a directory is created to store your processing results
work_dir = os.getcwd()+'/result'

print("current dir is: %s" % (os.getcwd()))
print("current working directory is:",work_dir) 

if os.path.isdir(work_dir):
    print("Folder exists")
else:
    print("Folder doesn't exists")
    os.mkdir(work_dir)

current dir is: d:\jupyter\notebook_scripts\Special\Deforestation
current working directory is: d:\jupyter\notebook_scripts\Special\Deforestation/result
Folder exists

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

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

ilwis.version()

'1.0 build 20260708'

#uncomment last line to remove any previous login credentials
from openeo.rest.auth.config import RefreshTokenStore
#RefreshTokenStore().remove()

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

Authenticated using refresh token.

collections = connection.list_collections()

df_collections = pd.DataFrame(collections)
display(df_collections[["id", "title", "description"]].head(20))

# Save to CSV
outfile = work_dir + "/openeo_collections.csv"
df_collections.to_csv(outfile, index=False, encoding="utf-8")

print(f"Collections saved to: {outfile}")

Collections saved to: d:\jupyter\notebook_scripts\Special\Deforestation/result/openeo_collections.csv

connection.describe_collection("SENTINEL2_L2A")

#note the date 
t = ["2026-05-01", "2026-10-01"]
s2_cube = connection.load_collection("SENTINEL2_L2A",
    spatial_extent={'west': xmin, 'east': xmax, 'south': ymin, 'north': ymax,"crs": "EPSG:4326"},
    temporal_extent= t,
    bands=['B02', 'B03', 'B04', 'B08'],
    max_cloud_cover=10,
)

s2_cube_red = s2_cube.median_time()

s2_cube_red.download(work_dir+"/s2_2026.tif")

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

Size(3871, 2459, 4)

#stretch the rasterbands contained in s2_2026 using a loop
multiple_stretch = []
multiple_bands = ilwis.do('selection',s2_2026,"rasterbands(0..3)") 
ls = ilwis.do('linearstretch',multiple_bands, 1) #using an upper and lower data limit defined by the cumulative 1 and 99 % thresholds
mb_stretch = ilwis.do('setvaluerange', ls, 0, 255, 1) #set ouput to byte range

#store the results as an ILWIS maplist - display the image as a color composite using the ilwis386 desktop software
mb_stretch.store('S2_2026.mpl')

# select required spectral channels from the initial image stack, here the visible red and NIR channels
Nir_in = ilwis.do('selection',s2_2026,"rasterbands(3)")
Red_in = ilwis.do('selection',s2_2026,"rasterbands(2)")

S2_2026_ndvi = ilwis.do('mapcalc','(@1 - @2)/(@1 + @2)', Nir_in, Red_in)
S2_2026_ndvi = ilwis.do('setvaluerange', S2_2026_ndvi, -1, 1, 0.001)

#store the results as an ILWIS raster map - display the NDVI using the ilwis386 desktop software 
#check the map values, especially of the forest and the defeorested areas
S2_2026_ndvi.store('S2_2026_ndvi.mpr')

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

raster_2np = np.fromiter(iter(S2_2026_ndvi), np.float64, dim) 
raster_2np = raster_2np.reshape((row, col))

2459 3871 9518789

# Plot the NDVI derived from Sentinel-2 MSI

fig, ax = plt.subplots(figsize=(12, 8))

im = ax.imshow(
    raster_2np,
    extent=[xmin, xmax, ymin, ymax],
    cmap="RdYlGn",
    vmin=0,
    vmax=1
)

# Colorbar
cbar = plt.colorbar(im, ax=ax, shrink=0.5)
cbar.set_label("NDVI")

# Axes
ax.set_xlabel("Longitude (°E)")
ax.set_ylabel("Latitude (°N)")
ax.set_title("NDVI from Sentinel-2 Mosaic of " + (t[0])+" to " + (t[1]))

# Limit to the AOI
ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)

plt.tight_layout()
plt.show()

connection.describe_collection("LANDSAT_BIMONTHLY_MOSAIC")

# Target years: 1997, 2000, ..., 2024 - note temporal extent from Landsatbi-monthly mosaic collection
years = range(1997, 2025, 3)

#download derived NDVI (only) from Landsat mosaic collection from 1997 to 2025, interval of 3 years

collection_id = "LANDSAT_BIMONTHLY_MOSAIC"

for year in years:
    print(f"Processing NDVI for {year}")

    # Load May–October observations
    cube = connection.load_collection(
        collection_id,
        spatial_extent={'west': (xmin), 'east': (xmax), 'south': (ymin), 'north': (ymax),"crs": "EPSG:4326"},
        temporal_extent=[
            f"{year}-05-01",
            f"{year}-10-01"
        ],
        bands=["B03", "B04"]
    )

    # NDVI = (NIR - Red) / (NIR + Red)
    ndvi = cube.ndvi(
        nir="B04",
        red="B03"
    )

    # Calculate median across the selected period
    ndvi_median = ndvi.median_time()

    # Download GeoTIFF
    outfile = os.path.join(
        work_dir,
        f"NDVI_{year}_May_Oct.tif"
    )

    ndvi_median.download(outfile)

    print(f"Saved: {outfile}")

print("All years processed.")

Processing NDVI for 1997
Saved: d:\jupyter\notebook_scripts\Special\Deforestation/result\NDVI_1997_May_Oct.tif
Processing NDVI for 2000
Saved: d:\jupyter\notebook_scripts\Special\Deforestation/result\NDVI_2000_May_Oct.tif
Processing NDVI for 2003
Saved: d:\jupyter\notebook_scripts\Special\Deforestation/result\NDVI_2003_May_Oct.tif
Processing NDVI for 2006
Saved: d:\jupyter\notebook_scripts\Special\Deforestation/result\NDVI_2006_May_Oct.tif
Processing NDVI for 2009
Saved: d:\jupyter\notebook_scripts\Special\Deforestation/result\NDVI_2009_May_Oct.tif
Processing NDVI for 2012
Saved: d:\jupyter\notebook_scripts\Special\Deforestation/result\NDVI_2012_May_Oct.tif
Processing NDVI for 2015
Saved: d:\jupyter\notebook_scripts\Special\Deforestation/result\NDVI_2015_May_Oct.tif
Processing NDVI for 2018
Saved: d:\jupyter\notebook_scripts\Special\Deforestation/result\NDVI_2018_May_Oct.tif
Processing NDVI for 2021
Saved: d:\jupyter\notebook_scripts\Special\Deforestation/result\NDVI_2021_May_Oct.tif
Processing NDVI for 2024
Saved: d:\jupyter\notebook_scripts\Special\Deforestation/result\NDVI_2024_May_Oct.tif
All years processed.

# Find all NDVI images
ndvi_files = sorted(glob.glob(os.path.join(work_dir, "NDVI_*_May_Oct.tif")))

# Load images into a dictionary
ndvi = {}

for filename in ndvi_files:
    year = os.path.basename(filename).split("_")[1]

    ndvi[year] = ilwis.RasterCoverage(filename)

    print(f"Loaded NDVI image for {year}")

Loaded NDVI image for 1997
Loaded NDVI image for 2000
Loaded NDVI image for 2003
Loaded NDVI image for 2006
Loaded NDVI image for 2009
Loaded NDVI image for 2012
Loaded NDVI image for 2015
Loaded NDVI image for 2018
Loaded NDVI image for 2021
Loaded NDVI image for 2024

#assign appropriate name 
ndvi_1997 = ndvi["1997"]
ndvi_2000 = ndvi["2000"]
ndvi_2003 = ndvi["2003"]
ndvi_2006 = ndvi["2006"]
ndvi_2009 = ndvi["2009"]
ndvi_2012 = ndvi["2012"]
ndvi_2015 = ndvi["2015"]
ndvi_2018 = ndvi["2018"]
ndvi_2021 = ndvi["2021"]
ndvi_2024 = ndvi["2024"]

print(ndvi_1997.size())

Size(1400, 880, 1)

#store the NDVI as an ILWIS map - display the map using the ilwis386 desktop software
ndvi_1997.store('landsat_ndvi1997.mpr')

#Transform to a numpy array for visualization using imshow for selected interval

interval = ndvi_1997

row = interval.size().ysize 
col = interval.size().xsize
dim = interval.size().linearSize()
print(row, col, dim)

raster_2npl = np.fromiter(iter(interval), np.float64, dim) 
raster_2npl = raster_2npl.reshape((row, col))

880 1400 1232000

# Create a plot of the NDVI from the Landsat mosaic of 1997
fig, ax = plt.subplots(figsize=(12, 8))

im = ax.imshow(
    raster_2npl,
    extent=[xmin, xmax, ymin, ymax],
    cmap="RdYlGn",
    vmin=0.5,
    vmax=1
)

# Colorbar
cbar = plt.colorbar(im, ax=ax, shrink=0.5)
cbar.set_label("NDVI")

fname = os.path.basename(str(interval))
year_title = fname.split("_")[1]


# Axes
ax.set_xlabel("Longitude (°E)")
ax.set_ylabel("Latitude (°N)")
ax.set_title("NDVI from Landsat Mosaic of " + year_title) 

# Limit to the AOI
ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)

plt.tight_layout()
plt.show()

cutoff = 0.815

#decade-1
#Note the use of an f-string to insert the Python variable into the ilwispy map calculation formula below
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)

#to check your results display the map using the ilwis386 desktop software
change1.store('change_1997_2006_cl.mpr')

#Transform to a numpy array for visualization using imshow for selected interval

row = change1.size().ysize 
col = change1.size().xsize
dim = change1.size().linearSize()
print(row, col, dim)

raster_2npl = np.fromiter(iter(change1), np.float64, dim) 
raster_2npl = raster_2npl.reshape((row, col))

880 1400 1232000

#get the class statistics

# Count pixels for each class
values, counts = np.unique(
    raster_2npl[np.isfinite(raster_2npl)],
    return_counts=True
)

# Create pandas DataFrame
df_classes = pd.DataFrame({
    "Class": values,
    "Pixel_count": counts
})

# Display without index
display(df_classes.style.hide(axis="index"))

# Create a plot

fig, ax = plt.subplots(figsize=(12, 8))

# Show every other coordinate tick
ax.xaxis.set_major_locator(MultipleLocator(0.05)) #unit is degree
ax.yaxis.set_major_locator(MultipleLocator(0.05))


# Five classes: 1, 2, 3, 4, 5
cmap = plt.get_cmap("jet", 5)

# Each integer class is centered in its color interval
bounds = [0.5, 1.5, 2.5, 3.5, 4.5, 5.5]
norm = mcolors.BoundaryNorm(bounds, cmap.N)

im = ax.imshow(
    raster_2npl,
    extent=[xmin, xmax, ymin, ymax],
    cmap=cmap,
    norm=norm,
    interpolation="nearest"
)


# Color bar
cbar = plt.colorbar(
    im,
    ax=ax,
    boundaries=bounds,
    ticks=[1, 2, 3, 4, 5],
    spacing="uniform",
    shrink=0.4
)

cbar.set_ticklabels([
    "till-1997",
    "till-2001",
    "till-2003",
    "till-2006",
    "forest"
])
cbar.set_label("Deforestation progress")

# Axes
ax.set_xlabel("Longitude (°E)")
ax.set_ylabel("Latitude (°N)")
ax.set_title("Deforestation from Landsat Mosaic between 1997–2006")

ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)

plt.tight_layout()
plt.show()

#decade-2
#calculate the deforestation change for the second time period
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)

change2.store('change_2009_2018_cl.mpr')

combine = ilwis.do('mapcalc', 'iff(@1<=4, @1, @2+4)',change1, change2)

combine.store('change_1997_2018_cl.mpr')

#Transform to a numpy array for visualization using imshow for selected interval

row = combine.size().ysize 
col = combine.size().xsize
dim = combine.size().linearSize()
print(row, col, dim)

raster_2npl = np.fromiter(iter(combine), np.float64, dim) 
raster_2npl = raster_2npl.reshape((row, col))

880 1400 1232000

#get the class statistics

# Count pixels for each class
values, counts = np.unique(
    raster_2npl[np.isfinite(raster_2npl)],
    return_counts=True
)

# Create pandas DataFrame
df_classes = pd.DataFrame({
    "Class": values,
    "Pixel_count": counts
})

# Display without index
display(df_classes.style.hide(axis="index"))

# Create a plot
fig, ax = plt.subplots(figsize=(12, 8))

# Show every other coordinate tick
ax.xaxis.set_major_locator(MultipleLocator(0.05)) #unit is degree
ax.yaxis.set_major_locator(MultipleLocator(0.05))


# Five classes: 1, 2, 3, 4, 5, 6, 7, 8, 9 
cmap = plt.get_cmap("jet", 9)

# Each integer class is centered in its color interval
bounds = [0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5]
norm = mcolors.BoundaryNorm(bounds, cmap.N)

im = ax.imshow(
    raster_2npl,
    extent=[xmin, xmax, ymin, ymax],
    cmap=cmap,
    norm=norm,
    interpolation="nearest"
)


# Color bar
cbar = plt.colorbar(
    im,
    ax=ax,
    boundaries=bounds,
    ticks=[1, 2, 3, 4, 5, 6, 7, 8, 9],
    spacing="uniform",
    shrink=0.4
)

cbar.set_ticklabels([
    "till-1997",
    "till-2001",
    "till-2003",
    "till-2006",
    "till-2009",
    "till-2012",
    "till-2015",
    "till-2018",
    "forest"
])
cbar.set_label("Deforestation progress")

# Axes
ax.set_xlabel("Longitude (°E)")
ax.set_ylabel("Latitude (°N)")
ax.set_title("Deforestation from Landsat Mosaic between 1997–2018")

ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)

plt.tight_layout()
plt.show()

# Retrieve the Sentinel-2 derived NDVI georeference information 
print(S2_2026_ndvi.size())
print(S2_2026_ndvi.envelope())

print()

csy = S2_2026_ndvi.coordinateSystem()
print("Coordinate system:", csy.toWKT())

Size(3871, 2459, 1)
736550.000000 9037400.000000 775260.000000 9061990.000000

Coordinate system: PROJCS["s2_2026.tif",GEOCS["s2_2026.tif",DATUM[" WGS 84",[DWGS84],ELLIPSOID["WGS 84",6378137.000000000000,298.257223563000],PRIMEM["Greenwich",0, AUTHORITY["EPSG",8901"]]],PROJECTION["Transverse_Mercator"],PARAMETER[,No]PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],PARAMETER["scale",0.9996],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-69],UNIT[meter,1.0]]

# Retrieve Landsat georeference information 
print(ndvi_1997.size())
print(ndvi_1997.envelope())
grfLandsat = ndvi_1997.geoReference()
print(grfLandsat)
csy = grfLandsat.coordinateSystem()

print("Coordinate system:", csy.toWKT())

Size(1400, 880, 1)
-66.850000 -8.700000 -66.500000 -8.480000
NDVI_1997_May_Oct.tif
Coordinate system: PROJCS["ndvi_1997_may_oct.tif",GEOCS["ndvi_1997_may_oct.tif",DATUM[" WGS 84",[DWGS84],ELLIPSOID["WGS 84",6378137.000000000000,298.257223563000],PRIMEM["Greenwich",0, AUTHORITY["EPSG",8901"]]],PROJECTION[""],,UNIT[degrees,1.0]]

#conduct the resampling operation for the Sentinel-2 NDVI
s2ndvi2026_res = ilwis.do('resample', S2_2026_ndvi, grfLandsat,'nearestneighbour')

print(s2ndvi2026_res.size())
print(s2ndvi2026_res.coordinateSystem())

Size(1400, 880, 1)
ndvi_1997_may_oct.tif

#check resulting resampled ndvi map derived from Sentinel-2
s2ndvi2026_res.store("s2_ndvi2026_res.mpr")

#calculate the changes during the final decade, also using the resamped Sentinel-2 ndvi from 2026
change3 = ilwis.do('mapcalc', f'iff(@1<={cutoff}, 1, iff(@2<={cutoff},2, iff(@3<={cutoff}, 3, 4)))',ndvi_2021, ndvi_2024, s2ndvi2026_res)

change3.store('change_2021_2026_cl.mpr')

combine_all = ilwis.do('mapcalc', 'iff(@1<=8, @1, @2+8)',combine, change3)

#also check the results obtained using ilwis desktop
combine_all.store('change_1997_2026_cl.mpr')

#Transform to a numpy array for visualization using imshow for selected interval

row = combine_all.size().ysize 
col = combine_all.size().xsize
dim = combine_all.size().linearSize()
print(row, col, dim)

raster_2npl = np.fromiter(iter(combine_all), np.float64, dim) 
raster_2npl = raster_2npl.reshape((row, col))

880 1400 1232000

#get the class statistics for the whole time series

# Count pixels for each class
values, counts = np.unique(
    raster_2npl[np.isfinite(raster_2npl)],
    return_counts=True
)

# Create pandas DataFrame
df_classes = pd.DataFrame({
    "Class": values,
    "Pixel_count": counts
})

# Assuming the pixel counts are in the second column
df_classes["Year"] = 1997 + (df_classes.index) * 3

# Set class 11 to 2026
df_classes.loc[df_classes["Class"] == 11, "Year"] = 2026

# Remove class 12 from the dataframe
df_classes = df_classes[df_classes["Class"] != 12].copy()

# Set Year to integer
df_classes["Year"] = df_classes["Year"].astype(int)


# Convert pixel size of 30 by 30 meters to hectares / km2
# Area of one pixel: 30 m * 30 m = 900 m}2
# Total area in square meters:  1354 pixels = 1218600 m^2
# Conversion to Hectares: Divide by 10000 (1218600 / 10000) = 121.86 ha
# Conversion to km²: Divide by 1000000 (1218600 / 1000000) = 1.2186 km²
# or
# 1 pixel = 900 m^2 and 1 hectare = 10,000 m^2
# 900/10,000 = 0.09 hectares per pixel

df_classes["Area_ha"] = df_classes["Pixel_count"] * 0.09
df_classes["Area_km2"] = df_classes["Pixel_count"] * 0.0009

# Display without index
display(df_classes.style.hide(axis="index"))

#get the total area deforested over time
deforested = df_classes["Area_km2"].sum()

print(f"Area deforested : {deforested:,.2f} km²")

Area deforested : 496.71 km²

# Sort by year
df_plot = df_classes.sort_values("Year")

fig, ax = plt.subplots(figsize=(10, 5))

ax.bar(df_plot["Year"], df_plot["Area_km2"], width=1.0)

ax.set_xlabel("Year")
ax.set_ylabel("Area (km²)")
ax.set_title("Deforestation from 1997 to 2026")
ax.set_xticks(df_plot["Year"])
ax.tick_params(axis="x", rotation=45)

ax.grid(axis="y", linestyle="--", alpha=0.5)
plt.tight_layout()
plt.show()

# AoI according to Landsat subset
pixels = 1440 * 880 #number of columns / rows of AoI
pixel_area = 900  # m² per pixel for Landsat
total_area_m2 = pixels * pixel_area
total_area_km2 = total_area_m2 / 1_000_000

print(f"Total Area: {total_area_km2:,.2f} km²")
print(f"Area Deforested : {deforested:,.2f} km²")

percentage_deforest = (deforested / total_area_km2)*100
print(f"Percentage Deforestation in AoI : {percentage_deforest:,.2f} %")

Total Area: 1,140.48 km²
Area Deforested : 496.71 km²
Percentage Deforestation in AoI : 43.55 %

# Plot final results

fig, ax = plt.subplots(figsize=(12, 8))

# Show every other coordinate tick
ax.xaxis.set_major_locator(MultipleLocator(0.05)) #unit is degree
ax.yaxis.set_major_locator(MultipleLocator(0.05))

# twelve classes: 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12 
cmap = plt.get_cmap("jet", 12)

# Each integer class is centered in its color interval
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]
norm = mcolors.BoundaryNorm(bounds, cmap.N)

im = ax.imshow(
    raster_2npl,
    extent=[xmin, xmax, ymin, ymax],
    cmap=cmap,
    norm=norm,
    interpolation="nearest"
)

# Color bar
cbar = plt.colorbar(
    im,
    ax=ax,
    boundaries=bounds,
    ticks=[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12],
    spacing="uniform",
    shrink=0.4
)

cbar.set_ticklabels([
    "till-1997",
    "till-2001",
    "till-2003",
    "till-2006",
    "till-2009",
    "till-2012",
    "till-2015",
    "till-2018",
    "till-2021",
    "till-2024",
    "till-2026",
    "forest"
])
cbar.set_label("Deforestation progress")

# Axes
ax.set_xlabel("Longitude (°E)")
ax.set_ylabel("Latitude (°N)")
ax.set_title("Deforestation as observed from Landsat Mosaic and Sentinel-2 between 1997–2026")

ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)

plt.tight_layout()
plt.show()

change_mask = ilwis.do('mapcalc', f'iff(@1<={cutoff}, 1, 0)',s2ndvi2026_res)

filter_avg = ilwis.do('linearrasterfilter', change_mask, 'average')
filter_mask = ilwis.do('mapcalc', 'iff(@1>=0.75, 1, 0)',filter_avg)

#Transform to a numpy array for visualization using imshow for selected interval

row = filter_mask.size().ysize 
col = filter_mask.size().xsize
dim = filter_mask.size().linearSize()
print(row, col, dim)

filter_mask_2npl = np.fromiter(iter(filter_mask), np.float64, dim) 
filter_mask_2npl = filter_mask_2npl.reshape((row, col))

880 1400 1232000

#display the numpy array created in the previous step using matplotlib 
fig = plt.figure(figsize =(12, 8))
plt.imshow(filter_mask_2npl)

<matplotlib.image.AxesImage at 0x18988f2ccd0>

#get the statistics for the filtered total deforested area

# Count pixels for each class
values, counts = np.unique(
    filter_mask_2npl[np.isfinite(filter_mask_2npl)],
    return_counts=True
)

# Create pandas DataFrame
df_2classes_only = pd.DataFrame({
    "Class": values,
    "Pixel_count": counts
})

# Display without index
display(df_2classes_only.style.hide(axis="index"))

#get the filtered / corrected total area deforested over time
total_pixels = df_2classes_only["Pixel_count"].sum()

class1_percentage = (
    df_2classes_only.loc[df_2classes_only["Class"] == 1, "Pixel_count"].sum()
    / total_pixels
) * 100

print(f"Percentage of class 1: {class1_percentage:.2f}%")

Percentage of class 1: 35.07%

#display results over openstreetsmap

# AOI bounds (EPSG:4326)
xmin = -66.85
xmax = -66.50
ymin = -8.70
ymax = -8.48

# Final binary raster: 0 = forest, 1 = deforested
# final_map should have shape (rows, columns)
height, width = raster_2npl.shape

# Source raster georeferencing
src_transform = from_bounds(
    xmin, ymin, xmax, ymax, width, height
)
src_crs = CRS.from_epsg(4326)
dst_crs = CRS.from_epsg(3857)

# Calculate destination grid
dst_transform, dst_width, dst_height = calculate_default_transform(
    src_crs, dst_crs, width, height,
    xmin, ymin, xmax, ymax
)

# Reproject categorical raster using nearest-neighbour resampling
map_3857 = np.full((dst_height, dst_width), -1, dtype=np.int16)

reproject(
    source=filter_mask_2npl.astype(np.int16),
    destination=map_3857,
    src_transform=src_transform,
    src_crs=src_crs,
    dst_transform=dst_transform,
    dst_crs=dst_crs,
    src_nodata=-1,
    dst_nodata=-1,
    resampling=Resampling.nearest
)

# Transform AOI bounding box to Web Mercator
transformer = Transformer.from_crs(
    "EPSG:4326", "EPSG:3857", always_xy=True
)

xmin_m, ymin_m = transformer.transform(xmin, ymin)
xmax_m, ymax_m = transformer.transform(xmax, ymax)

aoi = box(xmin_m, ymin_m, xmax_m, ymax_m)

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

margin_x = (xmax_m - xmin_m) * 0.15
margin_y = (ymax_m - ymin_m) * 0.15

ax.set_xlim(xmin_m - margin_x, xmax_m + margin_x)
ax.set_ylim(ymin_m - margin_y, ymax_m + margin_y)

# OpenStreetMap basemap
cx.add_basemap(
    ax,
    source="https://a.tile.openstreetmap.fr/hot/{z}/{x}/{y}.png",
    crs="EPSG:3857"
)

# Transparent forest, red deforested areas, transparent outside raster
cmap = ListedColormap(["none", "red"])
cmap.set_bad("none")

plot_map = np.ma.masked_where(
    (map_3857 != 0) & (map_3857 != 1),
    map_3857
)

# Make forest pixels transparent so the OSM map remains visible
plot_map = np.ma.masked_where(plot_map == 0, plot_map)

left = dst_transform.c
top = dst_transform.f
right = left + dst_width * dst_transform.a
bottom = top + dst_height * dst_transform.e

ax.imshow(
    plot_map,
    extent=[left, right, bottom, top],
    origin="upper",
    cmap=ListedColormap(["red"]),
    vmin=0,
    vmax=1,
    interpolation="nearest",
    alpha=0.7,
    zorder=2
)

# AOI boundary
x, y = aoi.exterior.xy
ax.plot(x, y, color="yellow", linewidth=2, label="AOI boundary")

ax.set_title("Amazon deforestation map")
ax.set_xlabel("Easting (m)")
ax.set_ylabel("Northing (m)")
ax.legend()
plt.tight_layout()
plt.show()

#conduct the resampling operation for the Sentinel-2 MSI images
s2_2026_res = ilwis.do('resample', mb_stretch, grfLandsat,'nearestneighbour')

#load the individual spectral channels of the resampled Sentinel-2 image
Blues = ilwis.do('selection',s2_2026_res,"rasterbands(0)")
Greens = ilwis.do('selection',s2_2026_res,"rasterbands(1)")
Reds = ilwis.do('selection',s2_2026_res,"rasterbands(2)")

#transform the spectral channels from ilwis format to a numpy array using the iterator
Blues_2np = np.fromiter(iter(Blues), np.ubyte, Blues.size().linearSize()) 
Blues_2np = Blues_2np.reshape((Blues.size().ysize, Blues.size().xsize))

Greens_2np = np.fromiter(iter(Greens), np.ubyte, Blues.size().linearSize()) 
Greens_2np = Greens_2np.reshape((Blues.size().ysize, Blues.size().xsize))

Reds_2np = np.fromiter(iter(Reds), np.ubyte, Blues.size().linearSize()) 
Reds_2np = Reds_2np.reshape((Blues.size().ysize, Blues.size().xsize))

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

# Plot Sentinel-2 composite and deforestation mask
fig, ax = plt.subplots(figsize=(12, 8))

# Background: Sentinel-2 colour composite
ax.imshow(
    ncol,
    extent=[xmin, xmax, ymin, ymax],
    origin="upper"
)

# Mask forest pixels (0), retaining only deforestation pixels (1)
deforestation = np.ma.masked_equal(filter_mask_2npl, 0)

# Transparent background for masked pixels
cmap = ListedColormap(["yellow"])
cmap.set_bad((0, 0, 0, 0))

# Overlay deforestation in semi-transparent yellow
ax.imshow(
    deforestation,
    extent=[xmin, xmax, ymin, ymax],
    origin="upper",
    cmap=cmap,
    vmin=1,
    vmax=1,
    alpha=0.5,
    interpolation="nearest"
)

# Axes and title
ax.set_xlabel("Longitude (°E)")
ax.set_ylabel("Latitude (°N)")
ax.set_title("Resampled Sentinel-2 composite – 2026 with deforestation mask")

ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)

plt.tight_layout()
plt.show()

	id	title	description
0	SENTINEL1_GRD	Sentinel-1 SAR GRD: C-band Synthetic Aperture ...	The Sentinel-1 mission provides data from a du...
1	COPERNICUS_30	Copernicus Global 30 meter Digital Elevation M...	Copernicus Global 30 meter Digital Elevation M...
2	SENTINEL3_OLCI_L1B	Sentinel 3 OLCI	Sentinel 3 imagery captured by OLCI sensor
3	SENTINEL3_SLSTR	Sentinel 3 SLSTR	Sentinel 3 imagery captured by SLSTR sensor
4	SENTINEL_5P_L2	Sentinel 5 Precursor	Sentinel 5 Precursor imagery captured by TROPO...
5	LANDSAT_BIMONTHLY_MOSAIC	Landsat mosaic (bi-monthly, v1.0.1)	Global, cloud-free, and reconstructed historic...
6	SENTINEL3_SYN_L2_SYN	Daily global Sentinel-3 Synergy products at 30...	[Experimental] Sentinel-3 SYNERGY products ar...
7	SENTINEL3_SLSTR_L2_LST	Daily global Sentinel-3 Land surface temperatu...	Experimental Daily global Sentinel-3 land su...
8	SENTINEL5P_L2_CO	Sentinel-5P TROPOMI Carbon Monoxide L2	This Collection provides Sentinel-5P Level-2 C...
9	SENTINEL5P_L2_NO2	Sentinel-5P TROPOMI Nitrogendioxide L2	This Collection provides Sentinel-5P Level-2 N...
10	SENTINEL5P_L2_CH4	Sentinel-5P TROPOMI Methane L2	This Collection provides Sentinel-5P Level-2 C...
11	COPERNICUS_VEGETATION_PHENOLOGY_PRODUCTIVITY_1...	Vegetation Phenology and Productivity paramete...	[EXPERIMENTAL] Copernicus High Resolution Land...
12	COPERNICUS_VEGETATION_PHENOLOGY_PRODUCTIVITY_1...	Vegetation Phenology and Productivity paramete...	[EXPERIMENTAL] Copernicus High Resolution Land...
13	COPERNICUS_PLANT_PHENOLOGY_INDEX	Seasonal PPI Trajectories, 10-daily, UTM proje...	[EXPERIMENTAL] Plant phenology index (PPI) at ...
14	ESA_WORLDCOVER_10M_2020_V1	Global landcover at 10m resolution for 2020.	WorldCover provides global land cover products...
15	ESA_WORLDCOVER_10M_2021_V2	Global landcover at 10m resolution for 2021.	WorldCover provides global land cover products...
16	SENTINEL1_GLOBAL_MOSAICS	Sentinel 1 Monthly Mosaics	[EXPERIMENTAL] Sentinel-1 monthly mosaics are ...
17	SENTINEL2_L2A	Sentinel-2 L2A	SENTINEL-2 is a wide-swath, high-resolution, m...
18	SENTINEL2_L1C	Sentinel-2 L1C	SENTINEL-2 is a wide-swath, high-resolution, m...
19	SENTINEL3_OLCI_L2_LAND	Sentinel-3 OLCI Land Full Resolution (NRT)	This collection provides Sentinel-3 OLCI Level...

Class	Pixel_count
1.000000	7163
2.000000	9864
3.000000	11805
4.000000	65228
5.000000	1137940

Class	Pixel_count
1.000000	7163
2.000000	9864
3.000000	11805
4.000000	65228
5.000000	35024
6.000000	49499
7.000000	91889
8.000000	85051
9.000000	876477

Class	Pixel_count	Year	Area_ha	Area_km2
1.000000	7163	1997	644.670000	6.446700
2.000000	9864	2000	887.760000	8.877600
3.000000	11805	2003	1062.450000	10.624500
4.000000	65228	2006	5870.520000	58.705200
5.000000	35024	2009	3152.160000	31.521600
6.000000	49499	2012	4454.910000	44.549100
7.000000	91889	2015	8270.010000	82.700100
8.000000	85051	2018	7654.590000	76.545900
9.000000	66964	2021	6026.760000	60.267600
10.000000	71883	2024	6469.470000	64.694700
11.000000	57533	2026	5177.970000	51.779700

Monitoring of deforestation¶

The current situation¶

The historical deforestation developments¶

Checking your results using openstreetmap¶

Class	Pixel_count
0.000000	799955
1.000000	432045