これは何？
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