Fragmentation Change#

This notebook demonstrates the frag_change() module, which computes multi-temporal fragmentation transitions between two time periods.

We will:

Compute fragmentation (FAD, window=27) for both the CLC 2000 and CLC 2018 Forest/Non-Forest maps
Run the frag_change() analysis to detect pixel-level connectivity changes
Visualise the spatial change map
Display the connectivity change histogram
Analyse the land and class transition matrices

Study area: Corsica, France — 100 m resolution, EPSG:3035

1. Import Libraries#

[ ]:

# Install pyguidos if running in Google Colab
try:
    import pyguidos as pg
except ImportError:
    !pip install -q pyguidos
    import pyguidos as pg
from pyguidos import utils
print(f"pyGuidos version: {pg.__version__}")

[ ]:

from pathlib import Path
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap, BoundaryNorm
import matplotlib.cm as cm
from matplotlib.cm import ScalarMappable
import rasterio
import pandas as pd

2. Define Paths#

[ ]:

# Input Forest/Non-Forest maps for two time periods
fnf_2000 = pg.DATA_DIR / "CLC2000_corsica_FNF.tif"
fnf_2018 = pg.DATA_DIR / "CLC2018_corsica_FNF.tif"

# Verify files exist
for f in [fnf_2000, fnf_2018]:
    status = "\u2713" if f.exists() else "\u2717 NOT FOUND"
    print(f"{status}  {f.name}")

[ ]:

# Output directory
CONFIG_PATH = pg.PROJECT_ROOT / ".notebook_config"

if CONFIG_PATH.exists():
    OUT_DIR = Path(CONFIG_PATH.read_text(encoding="utf-8").strip())
    print(f"Loaded existing workspace: {OUT_DIR}")
else:
    OUT_DIR = pg.PROJECT_ROOT / "output"
    print(f"Using default workspace: {OUT_DIR}")

OUT_DIR.mkdir(parents=True, exist_ok=True)

3. Compute Fragmentation for Both Time Periods#

We compute the Foreground Area Density (FAD) with a 27-pixel window for both the year 2000 and 2018 maps.

[ ]:

print("Computing fragmentation for CLC 2000...")
frag_2000 = pg.frag(
    in_tiff=fnf_2000,
    method='FAD',
    window_size=27,
    outdir=OUT_DIR,
    statists=True,
    stat_files=True,
    verb=True
)

print("Computing fragmentation for CLC 2018...")
frag_2018 = pg.frag(
    in_tiff=fnf_2018,
    method='FAD',
    window_size=27,
    outdir=OUT_DIR,
    statists=True,
    stat_files=True,
    verb=True
)

[ ]:

# Paths to the fragmentation outputs
frag_tiff_2000 = Path(frag_2000['output paths']['path tif'])
frag_tiff_2018 = Path(frag_2018['output paths']['path tif'])

print(f"T1 (2000): {frag_tiff_2000.name}")
print(f"T2 (2018): {frag_tiff_2018.name}")

4. Run Fragmentation Change Analysis#

frag_change() compares the two fragmentation rasters pixel by pixel, computing:

A categorical change map (spatial output)
A transition matrix tracking all class movements
Connectivity change statistics

[ ]:

print("Computing fragmentation change (2000 -> 2018)...")
fch_result = pg.frag_change(
    in_tiff_t1=frag_tiff_2000,
    in_tiff_t2=frag_tiff_2018,
    outdir=OUT_DIR,
    statists=True,
    stat_files=True,
    verb=True
)

5. Visualise the Change Map#

The output change map encodes connectivity change as pixel values 0-200, where:

Values < 99 = decrease in fragmentation in connectivity (T1 had higher FAD than T2)
Value [99-101] = no or insignificant change in fragmentation
Values > 101 = increase in fragmentation in connectivity (T2 has higher FAD than T1)
Values 250-254 = special codes (background dynamics, missing data)

[ ]:

# Helper function for rendering the map
def frag_change_legend(ax):
    """Adds a fragmentation change legend to the given axes."""
    import matplotlib.patches as mpatches

    legend_items = [
        mpatches.Patch(color='#CC5500', label='High decrease [-100, -21]'),
        mpatches.Patch(color='#FF8C00', label='Medium decrease [-20, -11]'),
        mpatches.Patch(color='#FFB732', label='Low decrease [-10, -2]'),
        mpatches.Patch(color='#FFFFCC', label='Stable/insignificant [-1, 1]'),
        mpatches.Patch(color='#66CCAA', label='Low increase [2, 10]'),
        mpatches.Patch(color='#339966', label='Medium increase [11, 20]'),
        mpatches.Patch(color='#006633', label='High increase [21, 100]'),
        mpatches.Patch(color='#00FF00', label='Forest gain [250]'),
        mpatches.Patch(color='#000000', label='Forest loss [251]'),
        mpatches.Patch(color='#D3D3D3', label='Nonforest [252]'),
        mpatches.Patch(color='#0066FF', label='Water [253]'),
        mpatches.Patch(facecolor='#FFFFFF', edgecolor='black', label='Missing [254]'),
    ]

    ax.legend(handles=legend_items, loc='upper left', bbox_to_anchor=(1.02, 1),
          fontsize=9, framealpha=0.9, title='ΔFAD', borderaxespad=0)

[ ]:

# Read the change map
change_tiff = Path(fch_result['output paths']['path tif'])
with rasterio.open(change_tiff) as src:
    data_change = src.read(1)

# Visualise island-level fragmentation
cmap_change, norm_change = pg.utils.get_tif_colormap(change_tiff)

fig, ax = plt.subplots(figsize=(7, 9))
ax.imshow(data_change, cmap=cmap_change, norm=norm_change, interpolation='none')
ax.set_title('Fragmentation Change Map — Corsica\nCLC 2000 \u2192 CLC 2018 (FAD, window=27)',
             fontsize=13, pad=15)
ax.axis('off')
frag_change_legend(ax)

plt.tight_layout()
plt.show()

6. Display the Connectivity Change Histogram#

The frag_change() module automatically generates a histogram PNG. Let’s also display the change distribution from the statistics.

[ ]:

# Display the auto-generated histogram
from IPython.display import Image, display

png_path = Path(fch_result['output paths']['path png'])
if png_path.exists():
    display(Image(filename=str(png_path), width=600))
else:
    print("Histogram PNG not found.")

[ ]:

# Print the 7-class change summary
change_freq = fch_result['output stats']['Frag change freq']

print("="*60)
print("Fragmentation Change Classes (2000 \u2192 2018)")
print("="*60)

total_fg = sum(change_freq.values())
for cls_name, count in change_freq.items():
    pct = count / total_fg * 100 if total_fg > 0 else 0
    print(f"  {cls_name:<30} : {count:>10} px  ({pct:6.2f}%)")
print(f"{'':30}   {'':>10}     ------")
print(f"  {'Total common foreground':<30} : {total_fg:>10} px  (100.00%)")

7. Land Cover Transition Matrix#

The land transition matrix shows pixel movements between three land categories:

Foreground (forest pixels with FAD values 0-100)
Background (non-forest pixels: values 101, 105, 106)
Missing (NoData: value 102)

[ ]:

# Land cover transition matrix
land_matrix = fch_result['output stats']['Land change matrix']

labels = ['Foreground', 'Background', 'Missing']
df_land = pd.DataFrame(land_matrix, index=[f'T1: {l}' for l in labels],
                        columns=[f'T2: {l}' for l in labels])

print("Land Cover Transition Matrix (pixels)")
print("Rows = Time 1 (2000), Columns = Time 2 (2018)")
print("="*60)
print(df_land.to_string())
print()

[ ]:

# Visualise land transition matrix as heatmap
fig, ax = plt.subplots(figsize=(7, 5))
im = ax.imshow(land_matrix, cmap='YlOrRd', interpolation='none')

# Labels
ax.set_xticks(range(len(labels)))
ax.set_yticks(range(len(labels)))
ax.set_xticklabels([f'T2: {l}' for l in labels], fontsize=11)
ax.set_yticklabels([f'T1: {l}' for l in labels], fontsize=11)

# Annotate cells
for i in range(len(labels)):
    for j in range(len(labels)):
        val = land_matrix[i, j]
        color = 'white' if val > land_matrix.max() * 0.5 else 'black'
        ax.text(j, i, f'{val:,.0f}', ha='center', va='center', fontsize=11, color=color)

ax.set_title('Land Cover Transition Matrix\n(2000 \u2192 2018)', fontsize=13, pad=15)
plt.colorbar(im, ax=ax, fraction=0.046, pad=0.04, label='Pixel count')
plt.tight_layout()
plt.show()

8. Fragmentation Class Transition Matrix#

This matrix shows transitions between the 5 fragmentation classes plus background:

Background (non-forest)
Rare (FAD 0-9%)
Patchy (FAD 10-39%)
Transitional (FAD 40-59%)
Dominant (FAD 60-89%)
Interior (FAD 90-100%)

[ ]:

# Class transition matrix
class_matrix = fch_result['output stats']['Class change matrix']

class_labels = ['Background', 'Rare', 'Patchy', 'Transitional', 'Dominant', 'Interior']
df_class = pd.DataFrame(class_matrix,
                         index=[f'T1: {l}' for l in class_labels],
                         columns=[f'T2: {l}' for l in class_labels])

print("Fragmentation Class Transition Matrix (pixels)")
print("Rows = Time 1 (2000), Columns = Time 2 (2018)")
print("="*75)
print(df_class.to_string())

[ ]:

# Visualise class transition matrix as heatmap
fig, ax = plt.subplots(figsize=(9, 7))

# Use log scale for better visibility of small transitions
from matplotlib.colors import LogNorm
matrix_plot = class_matrix.astype(float)
matrix_plot[matrix_plot == 0] = np.nan  # Hide zero cells

im = ax.imshow(matrix_plot, cmap='YlOrRd', norm=LogNorm(vmin=1, vmax=matrix_plot[~np.isnan(matrix_plot)].max()),
               interpolation='none')

# Labels
ax.set_xticks(range(len(class_labels)))
ax.set_yticks(range(len(class_labels)))
ax.set_xticklabels([f'T2:\n{l}' for l in class_labels], fontsize=10)
ax.set_yticklabels([f'T1: {l}' for l in class_labels], fontsize=10)

# Annotate cells with values
for i in range(len(class_labels)):
    for j in range(len(class_labels)):
        val = class_matrix[i, j]
        if val > 0:
            color = 'white' if val > class_matrix.max() * 0.1 else 'black'
            ax.text(j, i, f'{val:,.0f}', ha='center', va='center', fontsize=9, color=color)

ax.set_title('Fragmentation Class Transition Matrix\n(2000 \u2192 2018)', fontsize=13, pad=15)
plt.colorbar(im, ax=ax, fraction=0.046, pad=0.04, label='Pixel count (log scale)')
plt.tight_layout()
plt.show()

9. Connectivity Indices Comparison#

Compare the overall connectivity indices between the two time periods.

[ ]:

# Connectivity indices
avcon_A = fch_result['output stats']['A avcon']
avcon_B = fch_result['output stats']['B avcon']
fad_av_A = fch_result['output stats']['A fad_av']
fad_av_B = fch_result['output stats']['B fad_av']

print("="*60)
print("Connectivity Indices Comparison")
print("="*60)
print(f"{'Index':<20} {'T1 (2000)':>12} {'T2 (2018)':>12} {'Change':>12} {'Rel. [%]':>10}")
print("-"*60)
print(f"{'AVcon':<20} {avcon_A:>12.4f} {avcon_B:>12.4f} {avcon_B-avcon_A:>+12.4f} {(avcon_B-avcon_A)/avcon_A*100:>+10.2f}")
print(f"{'FAD_av':<20} {fad_av_A:>12.4f} {fad_av_B:>12.4f} {fad_av_B-fad_av_A:>+12.4f} {(fad_av_B-fad_av_A)/fad_av_A*100:>+10.2f}")
print()
print("AVcon  = Average connectivity (over reporting unit area)")
print("FAD_av = Average FAD (over foreground pixels only)")

10. Summary of Output Files#

[ ]:

print("Generated output files:")
print("="*60)
for key, path in fch_result['output paths'].items():
    p = Path(path)
    size_kb = p.stat().st_size / 1024 if p.exists() else 0
    print(f"  {key:<12} : {p.name:<30} ({size_kb:.1f} KB)")

11. Summary#

In this notebook we have:

Computed FAD fragmentation for both the CLC 2000 and CLC 2018 forest maps
Performed a multi-temporal fragmentation change analysis with pg.frag_change()
Visualised the spatial change map showing areas of connectivity increase and decrease
Examined the auto-generated connectivity change histogram
Analysed the land cover and fragmentation class transition matrices
Compared connectivity indices (AVcon and FAD_av) between the two time periods

The frag_change() module directly supports the EU Nature Restoration Regulation temporal reporting requirement by quantifying changes in forest connectivity status across monitoring periods.

[ ]: