Regional Analysis#

This notebook demonstrates how to combine pyGuidos analysis tools with extract_by_polygon() to perform regional landscape analysis across multiple administrative units.

The workflow follows three steps:

Island-level fragmentation: Run FAD fragmentation on the full Corsica Forest/Non-Forest map. Computing fragmentation at island level is essential to preserve the spatial context of the moving window — if fragmentation were computed independently per region, pixels near region boundaries would have incomplete neighbourhoods, introducing artificial edge effects.
Extract by region: Clip the island-level fragmentation GeoTIFF to each administrative subdivision using extract_by_polygon(). The original FAD pixel values are fully preserved.
Regional statistics: Run frag_stats() on each extracted region to compute the AVcon index and fragmentation class distribution per region.

AVcon (Average Connectivity) is computed as:

\[AVcon = \frac{\sum_{v=0}^{100} v \cdot n_v}{N_{ru}}\]

where \(n_v\) is the number of pixels with FAD value \(v\) and \(N_{ru}\) is the total number of pixels in the reporting unit (foreground + background). Unlike FAD_av which considers only forested pixels, AVcon accounts for the entire landscape, making it a better indicator of overall forest connectivity at the regional scale.

Input data:

Forest/Non-Forest map: CLC 2018, 100m resolution, Corsica
Administrative subdivisions: GISCO Communes database, 5 regions (EPSG:3035)

Analysis parameters: FAD method, window size 27x27 pixels (2.7km x 2.7km, ~729 ha)

1. Import Libraries and Define Paths#

[ ]:

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

# --- Data paths ---
fnf_tiff    = pg.DATA_DIR / "CLC2018_corsica_FNF.tif"
vector_file = pg.DATA_DIR / "GISCO_adm_corsica.gpkg"

[ ]:

from pathlib import Path
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from matplotlib.colors import ListedColormap
import rasterio
import pyogrio
import geopandas as gp
import pandas as pd

[ ]:

# --- Load the Output directory ---
CONFIG_PATH = pg.PROJECT_ROOT / ".notebook_config"
if CONFIG_PATH.exists():
    try:
        OUT_DIR = Path(CONFIG_PATH.read_text(encoding="utf-8").strip())
        print(f"Workspace synced: {OUT_DIR}")
    except Exception as e:
        print(f"Error reading config, using default. {e}")
        OUT_DIR = pg.PROJECT_ROOT / "output"
else:
    # Fallback if the user skipped Notebook 1
    OUT_DIR = pg.PROJECT_ROOT / "output"
    print(f"Config not found. Using default: {OUT_DIR}")

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

2. Step 1 — Island-level Fragmentation Analysis#

We first run FAD fragmentation on the full Corsica island map. This ensures that the moving window always has access to the complete neighbourhood for every pixel, regardless of administrative boundaries.

Note: If you have already run fragmentation in Notebook 3, you can skip this cell and point frag_island_tiff directly to the existing output file CLC2018_corsica_FNF_frag_fad_27.tif.

[ ]:

print("Running island-level fragmentation (FAD, window=27)...")
frag_island = pg.frag(
    in_tiff=fnf_tiff,
    method='FAD',
    window_size=27,
    outdir=OUT_DIR,
    statists=True,
    stat_files=True,
    verb=True
)

[ ]:

# Visualise island-level fragmentation
frag_island_tiff = Path(frag_island['output paths']['path tif'])
with rasterio.open(frag_island_tiff) as src:
    data_island = src.read(1)
cmap_island, norm_island = pg.utils.get_tif_colormap(frag_island_tiff)

fig, ax = plt.subplots(figsize=(7, 9))
ax.imshow(data_island, cmap=cmap_island, norm=norm_island, interpolation='none')
ax.set_title('Fragmentation (FAD) — Corsica island\nCLC 2018, 100m, Window 27x27',
             fontsize=13, pad=15)
ax.axis('off')

# Legend
legend_patches_frag = [
    mpatches.Patch(color='#d73228', label='Rare (0-10%)'),
    mpatches.Patch(color='#fa8c5a', label='Patchy (10-40%)'),
    mpatches.Patch(color='#ffc800', label='Transitional (40-60%)'),
    mpatches.Patch(color='#8cc864', label='Dominant (60-90%)'),
    mpatches.Patch(color='#00af00', label='Interior (90-100%)'),
    mpatches.Patch(color='#afafaf', label='Non-forest'),
]
ax.legend(handles=legend_patches_frag, loc='upper left', fontsize=10, framealpha=0.9)

# Show island-level indices
avcon_island = frag_island['output stats']['avcon']
fad_island   = frag_island['output stats']['fad_av']
ax.text(0.02, 0.02, f'AVcon  = {avcon_island:.2f}\nFAD_av = {fad_island:.2f}',
        transform=ax.transAxes, fontsize=11, verticalalignment='bottom',
        bbox=dict(boxstyle='round', facecolor='white', alpha=0.8))

plt.tight_layout()
plt.show()

3. Step 2 — Extract Fragmentation Output by Region#

We now clip the island-level fragmentation GeoTIFF to each administrative subdivision. The original FAD pixel values [0-100] and fragmentation class encoding are fully preserved in each extracted file.

[ ]:

print("Extracting fragmentation map by region...")
pg.extract_by_polygon(
    vector_path=str(vector_file),
    geotiff_path=str(frag_island_tiff),
    output_dir=str(OUT_DIR),
    id_field='ADM_ID',
    name_prefix='FRAG_'
)
print("\nExtraction completed. Output files:")
for f in sorted(OUT_DIR.glob('FRAG_*.tif')):
    info = utils.get_raster_info(f)
    print(f"  {f.name:<45} {info['rows']} x {info['cols']} px")

[ ]:

# Visualise extracted fragmentation maps side by side
region_frag_tiffs = sorted(OUT_DIR.glob('FRAG_*.tif'))
n_regions = len(region_frag_tiffs)

fig, axes = plt.subplots(1, n_regions, figsize=(4 * n_regions, 8))
for ax, tif in zip(axes, region_frag_tiffs):
    with rasterio.open(tif) as src:
        data = src.read(1)
    cmap, norm = pg.utils.get_tif_colormap(tif)
    region_name = tif.stem.replace('FRAG_', '')
    ax.imshow(data, cmap=cmap, norm=norm, interpolation='none')
    ax.set_title(region_name, fontsize=11)
    ax.axis('off')

fig.suptitle('Fragmentation (FAD) — Extracted by Region\n(island-level values preserved)',
             fontsize=13, y=1.01)
plt.tight_layout()
plt.show()

4. Step 3 — Regional Statistics with frag_stats()#

We use frag_stats() to compute statistics for each extracted region. Since the extracted files retain the GTB metadata tag from the island-level analysis, frag_stats() correctly identifies the tool and parameters without requiring any additional input.

[ ]:

# Compute statistics for each region
region_stats = {}

for tif in sorted(OUT_DIR.glob('FRAG_*.tif')):
    region_name = tif.stem.replace('FRAG_', '')
    stats = pg.frag_stats(
        frag_tiff=tif,
        stat_files=True,
        outdir=OUT_DIR
    )
    region_stats[region_name] = stats
    avcon  = stats['output stats']['avcon']
    fad_av = stats['output stats']['fad_av']
    print(f"  AVcon={avcon:.2f}  FAD_av={fad_av:.2f}")

print("\nStatistics completed for all regions.")

5. Regional Comparison — Summary Table and Bar Charts#

5.1 Summary Table#

[ ]:

# Build summary table
print(f"{'Region':<12} {'Forest px':>12} {'Non-For px':>12} {'Forest %':>10} "
      f"{'FAD_av':>8} {'AVcon':>8}")
print("-" * 68)

# Island-level reference row
fg_isl  = frag_island['input stats']['foreground pxl']
bg_isl  = frag_island['input stats']['background pxl']
tot_isl = fg_isl + bg_isl
print(f"{'CORSICA':<12} {fg_isl:>12} {bg_isl:>12} "
      f"{fg_isl/tot_isl*100:>10.1f} "
      f"{frag_island['output stats']['fad_av']:>8.2f} "
      f"{frag_island['output stats']['avcon']:>8.2f}  <- Map reference")
print("-" * 68)

summary = []
for region_name, stats in region_stats.items():
    inp   = stats['input stats']
    out   = stats['output stats']
    fg    = inp['foreground pxl']
    bg    = inp['background pxl']
    tot   = fg + bg
    pct   = fg / tot * 100 if tot > 0 else 0
    fad   = out['fad_av']
    avcon = out['avcon']
    print(f"{region_name:<12} {fg:>12} {bg:>12} {pct:>10.1f} {fad:>8.2f} {avcon:>8.2f}")
    summary.append({
        'region'       : region_name,
        'forest_pxl'   : fg,
        'nonforest_pxl': bg,
        'forest_pct'   : pct,
        'fad_av'       : fad,
        'avcon'        : avcon
    })

5.2 Bar Charts#

[ ]:

regions_list = [s['region']     for s in summary]
avcon_list   = [s['avcon']      for s in summary]
fad_list     = [s['fad_av']     for s in summary]
forest_pct   = [s['forest_pct'] for s in summary]

x = np.arange(len(regions_list))
width = 0.35

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

# Left — AVcon per region with island and mean reference lines
bars1 = axes[0].bar(x, avcon_list, width=0.5, color='#4e79a7',
                    edgecolor='black', linewidth=0.7)
axes[0].set_ylabel('AVcon', fontsize=12)
axes[0].set_title('Average Connectivity (AVcon)\nper Administrative Region',
                  fontsize=12, pad=12)
axes[0].set_xticks(x)
axes[0].set_xticklabels(regions_list, rotation=0, ha='center', fontsize=10)
axes[0].set_ylim(0, 100)
axes[0].axhline(y=avcon_island, color='red', linestyle='--', linewidth=1.5,
                label=f'Island AVcon: {avcon_island:.2f}')
axes[0].axhline(y=np.mean(avcon_list), color='orange', linestyle=':', linewidth=1.5,
                label=f'Regional AVcon mean: {np.mean(avcon_list):.2f}')
axes[0].legend(fontsize=10)
for bar, val in zip(bars1, avcon_list):
    axes[0].text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.8,
                 f'{val:.1f}', ha='center', va='bottom', fontsize=10)

# Right — Forest % vs FAD_av
axes[1].bar(x - width/2, forest_pct, width=width, color='darkgreen',
            edgecolor='black', linewidth=0.7, label='Forest %')
axes[1].bar(x + width/2, fad_list,   width=width, color='#59a14f',
            edgecolor='black', linewidth=0.7, label='FAD_av')
axes[1].set_ylabel('%', fontsize=12)
axes[1].set_title('Forest Cover (%) vs FAD_av Index\nper Administrative Region',
                  fontsize=12, pad=12)
axes[1].set_xticks(x)
axes[1].set_xticklabels(regions_list, rotation=0, ha='center', fontsize=10)
axes[1].set_ylim(0, 100)
axes[1].legend(fontsize=10)

fig.suptitle('Fragmentation Analysis — Regional Comparison\nCorsica CLC 2018 (100m), Window 27x27',
             fontsize=13, y=1.02)
plt.tight_layout()
plt.show()

6. Interpretation#

The regional comparison reveals differences in forest connectivity across the 5 administrative subdivisions of Corsica:

The island-level AVcon (red dashed line) serves as the reference for overall forest connectivity of Corsica. Regions above this line have higher than average connectivity; regions below have lower connectivity.
AVcon accounts for the entire reporting unit area, rewarding regions where forest patches are large and well-connected relative to total landscape
FAD_av only considers forested pixels — a region with the same forest cover but more fragmented patches will show a lower FAD_av
The gap between Forest % and FAD_av indicates the degree of fragmentation: a small gap suggests compact, well-connected forest; a large gap suggests fragmented patches with strong edge effects

Key methodological note: because fragmentation was computed at island level before extraction, the AVcon values per region reflect the true connectivity context within the broader Corsican landscape — not an isolated regional analysis. This is scientifically more meaningful for comparing regions within a continuous landscape.

7. Summary#

In this notebook we demonstrated a complete regional analysis workflow following best practice for landscape fragmentation analysis:

``frag()`` computed island-level FAD fragmentation, preserving the full spatial context of the moving window
``extract_by_polygon()`` clipped the fragmentation output to each of the 5 administrative regions, preserving original pixel values and GTB tags
``frag_stats()`` computed regional statistics directly from the extracted fragmentation GeoTIFFs, enabling meaningful cross-regional comparison

This three-step workflow — analyse at landscape level, extract by unit, compute statistics per unit — is the recommended approach for all pyGuidos tools when performing regional or comparative analysis.

The same workflow can be applied to any spatial unit — countries, protected areas, watershed boundaries — simply by replacing the input vector file.

[ ]: