これは何？

Watershed（分水界）アルゴリズムは、地形に水を流す比喩で動く画像セグメンテーション手法だ。

森林LiDARでの使いどころ：

樹高マップ（CHM）から1本1本の樹冠ポリゴンを切り出す
樹冠面積・等価半径の計算（樹高だけでは出せない情報）
炭素量推定の精度向上（面積から断面積を推定）

局所最大値法が「木の頭の座標」を出すのに対し、Watershedは「木が占める面積」まで求められる。

局所最大値法との違い

比較項目	局所最大値法	Watershed
出力	点（樹頂点の座標）	ポリゴン（樹冠の輪郭）
樹高	◯	◯
樹冠面積	✕	◯
等価半径	✕（推定のみ）	◯
計算コスト	軽い	やや重い
過分割リスク	低い	sigmaが小さいと多め

同じマーカーを使えば検出本数は同じ。違いは「各木の領域」が分かるかどうかだ。

アルゴリズム：山に水を流す

Watershedの直感的なイメージ：

CHM（樹高マップ）を「山の地形」として見る

  高い ── ★         ★         ★     ← 樹頂点（局所最大値）
           |\       /|\       /|
  低い ──  | \─────/ | \─────/ |     ← 谷（樹冠の境界）

手順

CHMを反転する（高い木 → 深い谷になる）
各樹頂点にマーカー番号を振る（この谷は何番の木か）
Watershedが各マーカーから水を満たしながら領域を広げる
隣のマーカーの水域と交わった境界 = 樹冠の境界線

樹冠は「隣の木と競合するまで広がる」ので、自然な輪郭が得られる。

Pythonでの実装

1. マーカーの作成

局所最大値法とまったく同じパラメータで種点を作る。

import numpy as np
import rasterio
from scipy.ndimage import gaussian_filter, maximum_filter, label as nd_label

SIGMA       = 1.5   # ガウスフィルターのsigma
WINDOW_SIZE = 5     # 局所最大値の探索窓（px）
MIN_HEIGHT  = 5.0   # 最低樹高（m）

# CHM読み込み
with rasterio.open('../data/processed/height_map.tif') as src:
    chm       = src.read(1).astype('float32')
    transform = src.transform
    crs       = src.crs
    nodata    = src.nodata
    pixel_w   = transform.a
    pixel_h   = abs(transform.e)

if nodata is not None:
    chm[chm == nodata] = np.nan

pixel_area = pixel_w * pixel_h  # m²/pixel

# スムージング
chm_filled = np.where(np.isnan(chm), 0, chm)
chm_smooth = gaussian_filter(chm_filled, sigma=SIGMA)
chm_smooth[np.isnan(chm)] = np.nan
chm_detect = np.where(np.isnan(chm_smooth), 0, chm_smooth)

# 局所最大値マスク
local_max_mask = (
    (chm_detect == maximum_filter(chm_detect, size=WINDOW_SIZE)) &
    (chm_detect >= MIN_HEIGHT) &
    ~np.isnan(chm)
)

# 各局所最大値に一意の整数ラベルを付ける
# 局所最大値は孤立した1点なので、1マーカー = 1本の木
markers, n_markers = nd_label(local_max_mask)
print(f'マーカー数: {n_markers:,}本')
# → マーカー数: 944本

nd_label は連結した True 領域ごとに番号を振る。局所最大値は孤立点なので、1マーカー = 1本の木になる。

2. Watershed実行

from skimage.segmentation import watershed

# CHMを反転：高いところ（木のてっぺん）が低い値になる
# → watershedは低い値から広がるので、木のてっぺんからスタートする
chm_inv = -chm_detect

# マスク：CHM > 0の領域のみ（空やnodataを除外）
ws_mask = (chm_detect > 0) & ~np.isnan(chm)

# Watershed実行
ws_labels = watershed(chm_inv, markers=markers, mask=ws_mask)

n_segments = len(np.unique(ws_labels[ws_labels > 0]))
print(f'セグメント数: {n_segments:,}')
# → セグメント数: 944

マーカー数とセグメント数が一致していれば正常動作。

3. メトリクス計算

import pandas as pd

segment_ids = np.unique(ws_labels[ws_labels > 0])
records = []

for sid in segment_ids:
    seg_mask = ws_labels == sid
    heights  = chm[seg_mask]
    heights  = heights[~np.isnan(heights)]
    if len(heights) == 0:
        continue
    px_count     = int(seg_mask.sum())
    crown_area   = px_count * pixel_area          # m²
    crown_radius = float(np.sqrt(crown_area / np.pi))  # 等価半径
    records.append({
        'segment_id'    : int(sid),
        'height_m'      : float(heights.max()),
        'mean_height_m' : float(heights.mean()),
        'crown_area_m2' : float(crown_area),
        'crown_radius_m': crown_radius,
    })

df = pd.DataFrame(records)
print(df[['height_m', 'crown_area_m2', 'crown_radius_m']].describe().round(2))

       height_m  crown_area_m2  crown_radius_m
count    944.00         944.00          944.00
mean      30.44         101.54            5.32
std        8.43          80.90            2.00
min        8.48           6.00            1.38
25%       24.93          48.00            3.91
50%       30.36          80.00            5.05
75%       35.99         129.00            6.41
max       49.99         867.00           16.61

御岳山では平均樹冠面積 101.5 m²、等価半径 5.3 m。高木が多いエリアらしい結果だ。

sigmaが分割数に与える影響

sigmaを変えると検出本数がどう変わるか：

from scipy.ndimage import maximum_filter

sigmas = [0.5, 1.0, 1.5, 2.0, 3.0, 5.0]

for s in sigmas:
    s_smooth = gaussian_filter(chm_filled, sigma=s)
    s_lm = (
        (s_smooth == maximum_filter(s_smooth, size=WINDOW_SIZE)) &
        (s_smooth >= MIN_HEIGHT)
    )
    s_markers, s_n = nd_label(s_lm)
    s_ws = watershed(-s_smooth, markers=s_markers,
                     mask=(s_smooth > 0))
    s_segs = len(np.unique(s_ws[s_ws > 0]))
    mark = '← 採用' if s == 1.5 else ''
    print(f'sigma={s:.1f}: {s_segs:,}本  {mark}')

sigma=0.5:  1,843本  ← 過分割（1本が複数セグメントに）
sigma=1.0:  1,388本
sigma=1.5:    944本  ← 採用
sigma=2.0:    647本
sigma=3.0:    317本
sigma=5.0:    139本  ← 過統合（複数本が1セグメントに）

sigmaを大きくするほど分割数は単調に減る。「正解」は現地調査データとの照合で決める。

可視化：境界線をCHMに重ねる

import matplotlib.pyplot as plt
from skimage.segmentation import find_boundaries

boundaries = find_boundaries(ws_labels, mode='outer')

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

# 左: CHM + 境界線
axes[0].imshow(chm, cmap='YlGn', vmin=0, vmax=50)
bound_overlay = np.zeros((*chm.shape, 4), dtype=np.float32)
bound_overlay[boundaries] = [0.9, 0.1, 0.1, 0.85]  # 赤い境界線
axes[0].imshow(bound_overlay)
axes[0].set_title(f'CHM + Watershed境界線 ({n_segments:,}本)')

# 右: セグメントをランダムカラーで色分け
rng = np.random.default_rng(42)
n_max = int(ws_labels.max())
palette = np.zeros((n_max + 1, 4), dtype=np.float32)
palette[1:] = np.column_stack([rng.random((n_max, 3)), np.full(n_max, 0.75)])
axes[1].imshow(palette[ws_labels])
axes[1].set_title('セグメント色分け（各色 = 1本）')

plt.tight_layout()
plt.show()

GeoPackageとして保存（QGIS連携）

import geopandas as gpd
from shapely.geometry import shape
from rasterio.features import shapes as rio_shapes

metrics_dict = df.set_index('segment_id').to_dict('index')
geom_records = []

for geom_dict, seg_id_f in rio_shapes(
    ws_labels.astype(np.int32),
    mask=(ws_labels > 0).astype(np.uint8),
    transform=transform
):
    sid = int(seg_id_f)
    if sid == 0 or sid not in metrics_dict:
        continue
    m = metrics_dict[sid]
    geom_records.append({
        'geometry'      : shape(geom_dict),
        'segment_id'    : sid,
        'height_m'      : m['height_m'],
        'crown_area_m2' : m['crown_area_m2'],
        'crown_radius_m': m['crown_radius_m'],
    })

gdf = gpd.GeoDataFrame(geom_records, crs=crs)
gdf.to_file('../data/processed/mitake_watershed.gpkg', driver='GPKG')
print(f'保存完了: {len(gdf):,}ポリゴン')
# → 保存完了: 944ポリゴン

QGISで開くと、1本ごとの樹冠が面データとして表示される。シンボロジで crown_area_m2 を色分けすると、大木の分布が一目でわかる。

まとめ

項目	内容
何をするか	CHMを反転し、各樹頂点（マーカー）から領域を広げて樹冠を切り分ける
核心	`watershed(-chm, markers=markers, mask=...)` の3引数
マーカー	局所最大値法と同じ → 検出本数は変わらない
付加価値	樹冠面積・等価半径が計算できる
sigmaの影響	小さい → 過分割、大きい → 過統合。現地データで調整する

# 最小構成
chm_inv  = -chm_detect
markers, _ = nd_label(local_max_mask)
ws_labels  = watershed(chm_inv, markers=markers, mask=ws_mask)

このシリーズは、東京都多摩地域の森林LiDARデータをPythonで分析する中で出会ったアルゴリズムを1つずつ解説しています。
実装コードは tama-forest notebooks で公開中。

作者のその他の実験 → https://niikun.net