Morphology: MSPA and SPA#

This notebook demonstrates the two morphological pattern tools of pyGuidos on a binary forest map:

``mspa()`` — Morphological Spatial Pattern Analysis: the full segmentation of the foreground into mutually exclusive morphological classes (Core, Islet, Edge, Perforation, Bridge, Loop, Branch and their variants). It is computed by the original miallib C engine of Soille and Vogt, embedded in pyGuidos, and reproduces the GuidosToolbox (GTB) output exactly.
``spa()`` — Simplified Pattern Analysis: a fast, pure-Python approximation of MSPA that classifies the foreground into up to 6 of the most used classes.

Input data: binary Forest / Non-Forest map derived from Corine Land Cover 2018 at 100m resolution, Corsica, France:

Pixel Value	Class
2	Foreground (Forest)
1	Background (Non-Forest)
0	NoData (water, outside extent)

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 path: binary Forest/Non-Forest map ---
fnf_tiff = pg.DATA_DIR / "CLC2018_corsica_FNF.tif"

[ ]:

from pathlib import Path
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from matplotlib.colors import ListedColormap, BoundaryNorm
import rasterio

[ ]:

# --- 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. Inspect Input Map#

Both mspa() and spa() write their result GeoTIFF with an embedded GTB colour palette. We render the maps with the built-in utils.get_tif_colormap() helper, which reads that palette back so the output displays with the exact GTB colours.

[ ]:

# Metadata
info = utils.get_raster_info(fnf_tiff)
print(f"Size      : {info['rows']} rows x {info['cols']} cols")
print(f"Dtype     : {info['dtype']}")
print(f"Resolution: {info['resX']} x {info['resY']} m")
print(f"EPSG      : {info['epsg']}")

# Pixel frequencies
with rasterio.open(fnf_tiff) as src:
    fnf_data = src.read(1)

fnf_freq = utils.get_pxl_freq(fnf_data)
tot = info['rows'] * info['cols']

print("\nPixel value distribution:")
labels = {0: 'NoData', 1: 'Background', 2: 'Foreground (Forest)'}
for val, name in labels.items():
    n = fnf_freq.get(val, 0)
    print(f"  Value {val} — {name:<20}: {n:>10} px  ({n/tot*100:5.2f}%)")

[ ]:

# Visualise input binary map
fnf_cmap = ListedColormap(['white', 'lightgrey', 'darkgreen'])
fnf_norm = BoundaryNorm([0, 1, 2, 3], fnf_cmap.N)

fig, ax = plt.subplots(figsize=(7, 9))
ax.imshow(fnf_data, cmap=fnf_cmap, norm=fnf_norm, interpolation='none')
ax.set_title('Forest / Non-Forest Input Map — Corsica\nCLC 2018, 100m resolution',
             fontsize=13, pad=15)
ax.axis('off')

legend_patches = [
    mpatches.Patch(facecolor='white', edgecolor='black', label='NoData (0)'),
    mpatches.Patch(color='lightgrey', label='Background (1)'),
    mpatches.Patch(color='darkgreen', label='Foreground / Forest (2)'),
]
ax.legend(handles=legend_patches, loc='upper left', fontsize=11, framealpha=0.9)
plt.tight_layout()
plt.show()

3. Morphological Spatial Pattern Analysis (MSPA)#

We run mspa() with the standard parameters:

connectivity=8 — 8-connected foreground
edge_width=1 — 1-pixel wide edge/transition zone
transition=True — keep Loop/Bridge pixels crossing an Edge/Perforation as their own class
intext=True — separate internal (hole-facing) features from external ones

The output GeoTIFF encodes each foreground pixel into one of the MSPA morphological classes and carries the GTB colour palette.

[ ]:

print("Running MSPA (connectivity=8, edge_width=1, transition=True, intext=True)...")
mspa_result = pg.mspa(
    in_tiff=fnf_tiff,
    connectivity=8,
    edge_width=1,
    transition=True,
    intext=True,
    outdir=OUT_DIR,
    statists=True,
    stat_files=True,
    verb=False
)
print("\nMSPA completed.")

[ ]:

# --- Output paths and input stats ---
print("Output paths:")
for key, path in mspa_result['output paths'].items():
    print(f"  {key:<10}: {path}")

print("\nInput pixel counts:")
for k, v in mspa_result['input stats'].items():
    print(f"  {k:<16}: {v:>10}")

[ ]:

# --- MSPA aggregated class statistics ---
out_stats = mspa_result['output stats']
agg = out_stats['class freq aggr']
fg = mspa_result['input stats']['foreground pxl']

print("Aggregated foreground classes:")
print(f"  {'Class':<14}{'Pixels':>12}{'% FG':>10}")
print("  " + "-" * 36)
for cls, n in agg.items():
    pct = n / fg * 100 if fg > 0 else 0
    print(f"  {cls:<14}{n:>12}{pct:>9.2f}%")

print(f"\n  Integral Foreground : {out_stats['integral foregr']:>12}")
print(f"  Porosity [%]        : {out_stats['porosity']:>12.4f}")

[ ]:

# --- Visualise MSPA output with its embedded GTB palette ---
mspa_tiff = mspa_result['output paths']['path tif']
cmap_mspa, norm_mspa = utils.get_tif_colormap(mspa_tiff)
with rasterio.open(mspa_tiff) as src:
    data_mspa = src.read(1)

# Legend patches for the main MSPA classes (GTB colours, transition palette)
mspa_legend = [
    mpatches.Patch(color='#00c800', label='Core (17)'),
    mpatches.Patch(color='#a03c00', label='Islet (9)'),
    mpatches.Patch(color='#000000', label='Edge (3)'),
    mpatches.Patch(color='#0000ff', label='Perforation (5)'),
    mpatches.Patch(color='#ff0000', label='Bridge (33)'),
    mpatches.Patch(color='#ffff00', label='Loop (65)'),
    mpatches.Patch(color='#ff8c00', label='Branch (1)'),
    mpatches.Patch(color='#dcdcdc', label='Background (0)'),
    mpatches.Patch(color='#C2C2C2', label='Border-opening (220)'),
    mpatches.Patch(color='#888888', label='Core-opening (100)'),
]

fig, ax = plt.subplots(figsize=(8, 10))
ax.imshow(data_mspa, cmap=cmap_mspa, norm=norm_mspa, interpolation='none')
ax.set_title('MSPA — Corsica\nCLC 2018, 100m  (8, 1, transition=True, intext=True)',
             fontsize=12, pad=15)
ax.axis('off')
ax.legend(handles=mspa_legend, loc='upper left', fontsize=9, framealpha=0.9)
plt.tight_layout()
plt.show()

4. Simplified Pattern Analysis (SPA)#

spa() is a fast, pure-Python approximation of MSPA. With classes=6 it classifies the foreground into Core, Edge, Perforation, Islet and Margin (Linear), plus Core-Opening. It is useful when the full MSPA segmentation is not required.

[ ]:

print("Running SPA (edge_width=1, classes=6)...")
spa_result = pg.spa(
    in_tiff=fnf_tiff,
    edge_width=1,
    classes=6,
    outdir=OUT_DIR,
    statists=True,
    stat_files=True,
    verb=False
)
print("\nSPA completed.")

[ ]:

# --- SPA class statistics ---
# 'class freq' is keyed by descriptive labels (e.g. '1 Core (17)'), and its
# exact set of keys depends on the number of classes requested.
spa_freq = spa_result['output stats']['class freq']
spa_fg = spa_result['input stats']['foreground pxl']

print("SPA class pixel counts:")
print(f"  {'Class':<22}{'Pixels':>12}{'% FG':>10}")
print("  " + "-" * 44)
for label, n in spa_freq.items():
    pct = n / spa_fg * 100 if spa_fg > 0 else 0
    print(f"  {label:<22}{n:>12}{pct:>9.2f}%")

[ ]:

# --- Visualise SPA output with its embedded GTB palette ---
spa_tiff = spa_result['output paths']['path tif']
cmap_spa, norm_spa = utils.get_tif_colormap(spa_tiff)
with rasterio.open(spa_tiff) as src:
    data_spa = src.read(1)

spa_legend = [
    mpatches.Patch(color='#00c800', label='Core (17)'),
    mpatches.Patch(color='#000000', label='Edge (3)'),
    mpatches.Patch(color='#0000ff', label='Perforation (5)'),
    mpatches.Patch(color='#a03c00', label='Islet (9)'),
    mpatches.Patch(color='#ff8c00', label='Margin (1)'),
    mpatches.Patch(color='#dcdcdc', label='Background (0)'),
    mpatches.Patch(color='#888888', label='Core-Opening (100)'),
]

fig, ax = plt.subplots(figsize=(8, 10))
ax.imshow(data_spa, cmap=cmap_spa, norm=norm_spa, interpolation='none')
ax.set_title('SPA (6 classes) — Corsica\nCLC 2018, 100m  (edge_width=1, classes=6)',
             fontsize=12, pad=15)
ax.axis('off')
ax.legend(handles=spa_legend, loc='upper left', fontsize=9, framealpha=0.9)
plt.tight_layout()
plt.show()

5. Zoomed MSPA vs SPA — Side by Side#

MSPA provides the complete morphological segmentation (including Bridge, Loop and Branch connectivity features), while SPA offers a lighter, pure-Python classification of the most used classes. Both share the same GTB colour convention for the classes they have in common.

[ ]:

# Define your region of interest (in pixel coordinates / array indices)
x_min, x_max = 620, 720  # Horizontal range (columns)
y_min, y_max = 620, 720  # Vertical range (rows)

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

# Subplot 1
axes[0].imshow(data_mspa, cmap=cmap_mspa, norm=norm_mspa, interpolation='none')
axes[0].set_xlim(x_min, x_max)
axes[0].set_ylim(y_max, y_min)  # Inverted Y-axis keeps origin at the top-left
axes[0].set_title('MSPA — full segmentation (Zoomed)', fontsize=12, pad=15)
axes[0].axis('off')

# Subplot 2
axes[1].imshow(data_spa, cmap=cmap_spa, norm=norm_spa, interpolation='none')
axes[1].set_xlim(x_min, x_max)
axes[1].set_ylim(y_max, y_min)
axes[1].set_title('SPA — simplified (Zoomed)', fontsize=12, pad=15)
axes[1].axis('off')

fig.suptitle('Morphology — Corsica, CLC 2018 (100m)', fontsize=14, y=1.01)
plt.tight_layout()
plt.show()

6. Summary#

In this notebook we applied both morphology tools to the Corsica Forest/Non-Forest map:

``mspa()`` produced the complete morphological segmentation, distinguishing the forest Core from its Edges, Perforations, Islets and the connecting Bridge/Loop/Branch features. The aggregated class statistics and the derived Porosity summarise the pattern quantitatively.
``spa()`` produced a simplified 6-class map covering the most used classes, at a lower computational cost and with no compiled dependency.

The MSPA output is bit-identical to the GuidosToolbox result for the same parameters, while SPA is a fast native-Python approximation. Standalone statistics can be recomputed from an existing output GeoTIFF at any time with pg.mspa_stats() or pg.spa_stats().

[ ]: