L38 CS立体図 2件統合分析 — 曲率×傾斜の RGB 段彩が描く地形微細構造 と 土砂災害警戒区域の物理的論理

データ取得手順

✅ このスクリプトは初回実行時にデータを自動取得します（DoBoX からの直接ダウンロード）。

実行コマンド:

cd "2026 DoBoX 教材"
python -X utf8 lessons/L38_cs_terrain.py

DoBoX のオープンデータは申請不要・商用/非商用とも利用可。 data/extras/ は .gitignore 対象（約 57 GB のキャッシュ）。 スクリプト実行で自動再生成されます。

学習目標と問い

本記事の独自性 — RGB 3バンド × 警戒区域の量的結合

CS立体図は専門家向けの視覚的地形可視化として広く使われていますが、 本記事はそのRGB 数値そのものを解析対象として扱う点で他と異なります:

• 視覚→数値: 通常 CS立体図は「目で見て地形を読む」ものですが、本記事では R/G/B の各バンド値を統計的に集計し、3 エリアで比較します。
• 地形タイプ自動分類: パーセンタイル適応閾値で「尾根 / 谷 / 急傾斜 / 平地 / 複雑」の 5 クラスにラベル付けし、どのクラスがどこに集塊するかを地図化します。
• 警戒区域との空間結合: 土砂災害警戒区域 (yellow) と特別警戒区域 (red) の ポリゴンを overview ラスタの座標系に rasterize し、各 zone 内の RGB 平均値を比較。 「警戒区域は CS図でどんな色をしているか」を量的に語る。

シリーズ構造判定: (a) 同一手法・解像度違い + (b) 自治体ごとの分割

このシリーズは 1m 版 (区分A) / 50cm 版 (区分B) の 2 解像度系統で、 各自治体ごとに GeoTIFF (RGB 3バンド uint8) を切出して公開する構造。 1m 版は 22 自治体、50cm 版は 12 自治体に対応。

• 1m 版 (dsid 1628): 県内 22 自治体 GeoTIFF。 本記事は世羅町 (resource_id 177706)と熊野町 (176860)を採用 (= L36/L37 と同自治体ペアで LiDAR ファミリの相補比較が可能)。
• 50cm 版 (dsid 1629): 県内 12 自治体 GeoTIFF (1m 版より少ない)。 熊野町は 50cm 版に存在しないため、代わりに坂町 (176895)を採用。 坂町は平成 30 年 7 月豪雨の被災地で、土石流災害が多発した地域 ⇒ 「CS立体図が崩壊地形をどう描くか」を見る研究的価値が高い。
• 3 エリアは 1m × 2 + 50cm × 1 の構成で、解像度差と地形類型差を同時に観察可能。

独自用語の定義 (本記事内で固定)

研究の問い (RQ)

広島県の CS立体図 2 解像度 (1m / 50cm) は、それぞれ 世羅町 (中山間/林業地), 熊野町 (都市近郊/山岳混在), 坂町 (平成30豪雨被災地) という異なる地形リスク特性のエリアを公開している。 3 エリアの CS立体図を RGB 3 軸で比較し、地形微細構造の構成比を量化する。 さらに土砂災害警戒区域との空間結合で、CS図のどの色域が警戒域に多いかを定量化する。

1. 3 dataset の解像度・面積・有効データ率・ピクセル数の構造比較は?
2. R, G, B 各バンドの値分布はエリアごとにどう違うか?
3. 適応閾値による地形タイプ自動分類の構成比は?
4. HSV 散布から見える 3 エリアの色相・彩度パターン差は?
5. 土砂災害警戒区域内 vs 区域外で R, G, B, Saturation はどう違うか?
6. 解像度差 (1m vs 50cm) で微地形コントラストはどう変わるか?

仮説 H1〜H6

• H1 (世羅町は緩傾斜丘陵): 世羅の CS図 Saturation 平均 ≤ 0.10。 中山間の植林地は緩い丘陵地形が中心で、極端な凹凸が少ない。
• H2 (熊野町は急峻地形): 熊野の Sat 平均 > 世羅。 熊野は花崗岩マサ土地形で急傾斜面が多い。
• H3 (坂町は崩壊跡を含む急傾斜): 坂町の Sat 平均が 3 エリアで最大。 平成30豪雨の崩壊地形を CS図が刻み込んでいる。
• H4 (警戒区域内の CS図は微地形強): 特別警戒区域内の Saturation は 区域外より大きい。急斜面で凹凸が強い場所は警戒区域に指定されているはず。
• H5 (色相は地形タイプを区別): HSV の Hue 標準偏差 > 0.05 で、 赤系 (尾根) と青系 (谷) の双方が混在することを確認。
• H6 (50cm 版は微細構造を捉える): 坂町 (50cm) の Sat 標準偏差は 他 1m 版より大きい。高解像度ほど微地形のばらつきを拾う。

到達点

3 エリアの CS立体図を「3 自治体の代表エリア」として読み解き、 合計 633 M セル × 3 バンド ≈ 1,900 M 数値の RGB 段彩データを overview ピラミッドで効率的に概観する。 さらに土砂災害警戒区域 Shapefile (~ 2,000 ポリゴン) との空間結合で 「CS図と警戒区域の関係」を量的に検証する。 これは CS立体図単独・警戒区域単独では不可能な、両者の同時空間整合分析である。

使用データ

本記事が使用する 2 dataset の一覧。両者は同一手法 (CS立体図 = 曲率×傾斜の RGB 段彩) の 解像度違いで、各自治体ごとに GeoTIFF を切り出して公開する構造。本記事は 2 dataset から 合計 3 自治体を代表として用います。

データ仕様 (両 dataset 共通)

• 形式: GeoTIFF, 3 バンド uint8 (0-255) = RGB 段彩図
• CRS: EPSG:6671 (JGD2011 / Japan Plane Rectangular CS III) — m 単位の平面直角座標
• nodata: 1m 版は 0 (3バンド全て 0 = 無効), 50cm 版 (坂町) は明示なし (推定: 全バンド 0 が境界)
• 値の意味: R/G/B 各バンド 0-255。本記事の解釈では: R = 尾根強度 (凸地形), G = 傾斜 (明度), B = 谷強度 (凹地形)
• 派生過程: 標高 (DTM) → 曲率計算 + 傾斜計算 → RGB 合成 → 段彩図
• 切出単位: 各自治体の行政界 + 100m バッファで矩形切出 (L36/L37 と同様)
• 用途: 地形微細構造の視覚判読 (尾根/谷/段丘/崩壊地形/古砂防 等)
• ライセンス: CC-BY 4.0 (測量法 44 条に基づく公共測量 2 次著作物)

L36 樹高 / L37 点群密度 との関係 (本記事は単独, 参考併置のみ)

3 自治体の特性比較

データ品質メモ

• 3 バンド uint8: 各バンド 0-255 の整数値。RGB 3 チャンネルが それぞれ独立した地形指標 (尾根/傾斜/谷) を表現。
• 1m 版は overview あり (2/4/8/16/32/64/128), 熊野町は overview なし (= 内部に ピラミッドが構築されていない)。坂町 (50cm) は overview あり。
• nodata の扱い: 1m 版は (0,0,0) = 無効。50cm 版 (坂町) は nodata 未指定だが、 実データでは (0,0,0) を境界マスクとして同様に扱う。
• 値の典型範囲: 全エリアで R/G/B 平均が 160-200 の中間色帯に集中 (= パステル系の地図画像)。極端な R 高 (>240) や B 高 (>240) は少数。
• 解像度差は overview で見えにくい: 全エリアを overview /1/32〜64 で読込むため、 解像度差 (1m vs 50cm) は生解像度サンプルでのみ確認できる (図 8)。

ダウンロード

生データ (DoBoX 直リンク)

• CS立体図（1m）_世羅町: 直 DL ZIP / DoBoX resource 177706 / dsid 1628
• CS立体図（1m）_熊野町: 直 DL ZIP / DoBoX resource 176860 / dsid 1628
• CS立体図（50cm）_坂町: 直 DL ZIP / DoBoX resource 176895 / dsid 1629

土砂災害警戒区域データ (本記事の空間結合に必要)

本記事はL10 土砂災害警戒区域 × 用途地域と 同じ Shapefile (DoBoX dsid 48) を使用します。 data/extras/sediment_shp/ に展開済みであれば自動利用。

中間データ (本記事生成 CSV)

• L38_series.csv — 3 dataset の一覧 (代表自治体・解像度・サイズ)
• L38_tiff_meta.csv — TIFF メタ (CRS, bounds, セル数, overview, dtype, nodata, 3 バンド)
• L38_basic_stats.csv — RGB 各バンド統計 (平均・中央値・std) + Saturation
• L38_terrain_classes.csv — 5 種地形タイプの構成比 (尾根/谷/急/平地/複雑)
• L38_rgb_histogram.csv — R/G/B 各バンドの 32 ビン ヒストグラム
• L38_sediment_rgb_compare.csv — 警戒区域 内/外 の RGB 平均比較
• L38_hypothesis_results.csv / L38_hypothesis_results.json — H1〜H6 検証結果

図 (本記事生成 PNG)

• L38_fig1_dataset_overview.png — 3 dataset カード + Saturation ヒスト
• L38_fig2_raster_maps.png — CS立体図 RGB 主題図 (3 エリア)
• L38_fig3_terrain_classification.png — 地形タイプ自動分類マップ
• L38_fig4_rgb_histograms.png — R/G/B 各バンドのヒストグラム比較
• L38_fig5_terrain_composition.png — 5 種地形タイプ構成比 (積み上げ)
• L38_fig6_sediment_overlay.png — CS立体図 × 土砂警戒区域 重ね合わせ
• L38_fig7_sediment_rgb_compare.png — 警戒区域 内/外 RGB 平均比較
• L38_fig8_center_samples.png — 中央 1024×1024 セル本来解像度サンプル
• L38_fig9_concept.png — CS立体図のしくみ概念図
• L38_fig10_hsv_scatter.png — HSV 散布図 (色相 × 彩度)
• L38_fig11_small_multiples.png — 地形タイプ別 small multiples (3×4)

再現用 Python スクリプト

L38_cs_terrain.py を取得して プロジェクトルートで py -X utf8 lessons/L38_cs_terrain.py を実行。 TIFF が無ければ ZIP 自動 DL → 自動展開 → overview 読込で分析が走ります (初回 DL を除けば約 13 秒で完走)。

分析 1: 3 dataset とエリアの構造を可視化

狙い

2 dataset (1m / 50cm) と 3 自治体 (世羅 / 熊野 / 坂) の規模・形状・データ量を 1 枚の絵で示す。Saturation ヒストグラムを重ねて地形微細構造の強さを最初に視覚化する。

手法 (リテラシレベル解説)

L36 / L37 と全く同じ rasterio スタイルだが、CS立体図には重要な違いがある:

• 3 バンド読込: ds.read() でバンド指定なしに呼ぶと (3, H, W) の 3D 配列が返る。 各バンドは uint8 (0-255) の整数値。
• RGB ではなく数値: バンド 1 = R, バンド 2 = G, バンド 3 = B。各セルは 3 つの 0-255 整数で「色」が指定される。本記事はこの 3 値を独立した地形指標として扱う。
• Saturation 算出: HSV 変換で「彩度」を計算。S = (max(R,G,B) - min(R,G,B)) / max。 微地形コントラストが強い場所ほど Sat が高い。

実装

↑ L38_cs_terrain.py 行 1659–1721

 1
 2
 3
 4
 5
 6
 7
 8
 9
1668
1669
1670
1671
1672
1673
1674
1675
1676
1677
1678
1679
1680
1681
1682
1683
1684
1685
1686

import rasterio
from rasterio.enums import Resampling
import numpy as np

TARGET_CELL_M = 32   # overview 目標セルサイズ
for area in AREAS:
    with rasterio.open(area["tif"]) as ds:
        ovr_choices = ds.overviews(1)
        target_factor = TARGET_CELL_M / area["resolution_m"]
        best = next((f for f in ovr_choices if f >= target_factor),
                    ovr_choices[-1] if ovr_choices else int(target_factor))
        ov_w = ds.width // best
        ov_h = ds.height // best

        # 3 バンドを 1 度に読込 → shape (3, H, W)
        rgb = ds.read(out_shape=(3, ov_h, ov_w),
                      resampling=Resampling.average).astype(np.float32)
        area["rgb"] = rgb

        # nodata マスク (全バンド 0 = 無効)
        valid_mask = rgb.sum(axis=0) > 0
        area["valid_mask"] = valid_mask

        # HSV の Saturation 計算
        maxv = np.maximum.reduce(rgb)
        minv = np.minimum.reduce(rgb)
        sat = np.where(maxv > 0, (maxv - minv) / np.maximum(maxv, 1), 0)
        area["sat"] = sat

図と読み取り (図 1: dataset カード + Saturation ヒスト)

なぜこの図か: 3 dataset の規模と特性をカードで言語的に、 かつ Saturation ヒストグラムで「微地形コントラスト分布」として量的に、 同時に 1 枚で伝えるため。

図 1: 3 dataset カード (左) + Saturation 分布 (右)

読み取り:

• 3 エリアとも RGB 平均が R > G > B の順で減衰する典型的なパステル色域に集中。 Sera 202.4 / 192.5 / 187.9, Kumano 196.3 / 185.2 / 180.8, Saka 175.5 / 162.2 / 159.4。
• Saturation ヒスト (右) は3 エリアともほぼ 0 にピークを持つ右裾分布。 多くのセルが「ほぼ無彩色」 (中間グレー)、少数のセルが「彩度高め」 (強い凹凸)。
• Sat 平均は世羅 0.074 < 熊野 0.083 < 坂 0.096。順位は仮説 (世羅 < 熊野 < 坂) と一致。 特に坂町は明確に右にずれており、急傾斜・崩壊地形の影響が見える。
• 坂町の R 平均が顕著に低い (175.5) のは「全体的に暗い色 = 急斜面影が多い」と読める。

表 (基本統計)

読み取り: R/G/B 平均は 3 エリアで似た値域 (160-202) だが、標準偏差は明確に違う ⇒ 世羅 R-std=15.7, 熊野 19.7, 坂 42.0。坂町は標準偏差が約 3 倍も大きい = 「明るい部分と暗い部分が広く分布」 ⇒ 急傾斜面と平地の混在による陰影コントラストの強さを示す。 Sat 平均と Sat 標準偏差はどちらも坂町が最大で、仮説 H3 (坂町は崩壊跡を含む急傾斜) を支持する強い量的証拠。

分析 2: CS立体図 RGB 主題図化

狙い

CS立体図のRGB 元画像を 3 エリア並列で表示し、 「地形微細構造はどう見えるか」を視覚的に理解する。これは CS立体図本来の使い方 (目で見て地形を読む) を踏襲したセクション。

手法 (リテラシレベル解説)

matplotlib.pyplot.imshow に (H, W, 3) 形状の uint8 配列を渡すと RGB 画像として表示される。rasterio の read() 結果は (3, H, W) なので np.stack([R, G, B], axis=-1) で軸を入れ替える必要がある。

• 赤い領域 = R が突出 = 尾根線 (凸地形)。標高ラスタの正の曲率部分。
• 青い領域 = B が突出 = 谷線・崩壊地形 (凹地形)。標高ラスタの負の曲率部分。
• 暗い領域 = 全バンド低 = 急斜面の影。標高ラスタの傾斜が強い場所。
• パステルピンク〜灰 = 緩傾斜の平地・台地・尾根丘陵。CS立体図の「背景色」。

実装

↑ L38_cs_terrain.py 行 1738–1753

1738
1739
1740
1741
1742

# 3 バンドを RGB 画像形式に整形
rgb_disp = np.stack([rgb[0], rgb[1], rgb[2]], axis=-1).astype(np.uint8)
# nodata セルを白に置換
rgb_disp_show = np.where(valid_mask[..., None], rgb_disp, 255)
ax.imshow(rgb_disp_show, extent=area_bounds, origin="upper", interpolation="nearest")

図と読み取り (図 2: CS立体図 RGB 主題図)

なぜこの図か: 3 エリアの地形微細構造を視覚的に並べ、 専門家のように「赤=尾根、青=谷、暗=急斜面」を読み取る練習をするため。

図 2: 3 エリア CS立体図 RGB 元画像

読み取り:

• 世羅町 (左): 全体に淡いピンク〜灰色のパステル基調。河川沿いに細い青線 (= 谷地形) が走り、それ以外は緩やかな起伏。中山間の典型的な丘陵地形。
• 熊野町 (中): 中央〜南東に暗い色域 (= 急斜面が広い山地)。 谷部の青色が筋状に強く出るのが特徴的で、花崗岩マサ土地形特有の急峻な谷。
• 坂町 (右): 全体に色のコントラストが強く、暗い影 (急斜面) と明るい尾根が 斑模様に交互する。南東の住宅地付近に明確な凹凸地形 (= 平成30豪雨で土石流が 発生した山地と扇状地の境界域) が見える。
• 3 エリアとも北 (上) が山地、南 (下) が平地・海岸寄りの傾向。 坂町は最も小さい (4.4 × 6.9 km) ため細部構造が見える。

表 (TIFF メタ + 行政界面積)

読み取り: 世羅町は 27 × 16 km と広大で 446 M セル × 3 = 約 1.3 GB の生 TIFF。 熊野町は内部 overview を持たない (= ピラミッド未構築) ため out_shape で擬似粗化。 坂町は 50cm 解像度のため 30 km² で 122 M セル相当。解像度差 = ピクセル数 4 倍がはっきり分かる。

分析 3: 地形タイプ自動分類 (パーセンタイル適応閾値, 5 クラス)

狙い

専門家の目で「赤い場所は尾根、青い場所は谷」と判読する代わりに、 RGB 値からアルゴリズム的に地形タイプを自動分類する。 これにより 3 エリアの地形タイプ構成比を量的に比較可能になる。

手法 (リテラシレベル解説)

本記事独自のパーセンタイル適応閾値分類を使う。CS立体図の RGB は パステル色 (中間値) が大多数なので、絶対値閾値ではほぼ全セルが同じクラスに振られる。 そこで各エリアの分布から相対的に閾値を決める。

限界と代替:

• 限界: 適応閾値はエリア横断で「絶対値が違う」ので 「世羅の尾根」と「坂町の尾根」は厳密には同じ閾値ではない。 ただし各エリア内での相対ランキングとして意味がある。
• 代替案: 絶対値閾値 (例 R-B >= 30) で揃える。ただしこれだと CS立体図の中間色域では ほぼ全セルが「平地」になってしまう (最初の試行では 100% 平地)。
• 真の正解: 標高ラスタ (DTM) から曲率・傾斜を直接計算する (発展課題 Z2)。

実装

↑ L38_cs_terrain.py 行 1816–1873

 1
 2
 3
 4
 5
 6
 7
 8
 9
1825
1826
1827
1828
1829
1830
1831
1832
1833
1834
1835
1836
1837
1838
1839
1840
1841

def classify_terrain(rgb, valid_mask):
    """RGB を 5 種地形タイプにラベル付け (パーセンタイル適応閾値)"""
    R, G, B = rgb[0], rgb[1], rgb[2]
    diff_RB = R - B
    bright = (R + G + B) / 3.0
    band_std = np.stack([R, G, B], axis=0).std(axis=0)

    # 各エリアの分布から閾値を決める
    p75 = np.percentile(diff_RB[valid_mask], 75)
    p25 = np.percentile(diff_RB[valid_mask], 25)
    bp15 = np.percentile(bright[valid_mask], 15)
    sp25 = np.percentile(band_std[valid_mask], 25)

    label = np.zeros(R.shape, dtype=np.int8)
    is_ridge  = (diff_RB >= p75) & valid_mask
    is_valley = (diff_RB <= p25) & valid_mask
    is_steep  = (bright <= bp15) & ~is_ridge & ~is_valley & valid_mask
    is_flat   = (band_std <= sp25) & ~is_ridge & ~is_valley & ~is_steep & valid_mask
    is_complex = valid_mask & ~is_ridge & ~is_valley & ~is_steep & ~is_flat

    label[is_ridge]   = 1
    label[is_valley]  = 2
    label[is_steep]   = 3
    label[is_flat]    = 4
    label[is_complex] = 5
    return label

図と読み取り (図 3: 地形タイプ分類マップ)

なぜこの図か: RGB 元画像 (図 2) の隣に分類結果を並べて、 「色がどの分類クラスに対応するか」を視覚的に確認する。

図 3: 地形タイプ自動分類マップ (5 クラス)

読み取り:

• 3 エリアとも赤 (尾根) と青 (谷) が拮抗して 25-30% ずつ分布。 これはパーセンタイル閾値が「上位 25% 以上」「下位 25% 以下」で対称定義のため。
• 急傾斜 (濃灰) は 3 エリアで 4-6% と少なめ。CS立体図の「明度」軸が圧縮されているため。
• 複雑地形 (紫) は 38-42% と最大の比率を占める。 = 中間色だが std や bright で例外的なセル群 = 細かい起伏が混在する典型的な山岳地形。
• 世羅町は河川沿いに青い「谷」帯が連続的に走り、その周囲に赤い「尾根」が伴う。 = 水系構造が CS立体図に明確に刻まれている。
• 坂町は各タイプが斑状に細かく混在。50cm 解像度で微地形が拾えていることが原因。

図と読み取り (図 5: 地形タイプ構成比 stack-bar)

図 5: 5 種地形タイプの構成比 (3 エリア)

読み取り: 3 エリアの構成比は非常に似ている (尾根 25-28%, 谷 25-30%, 急 4-6%, 平地 0-1%, 複雑 38-42%)。 これは「パーセンタイル適応閾値」の必然的帰結 — 各エリア内で上位/下位 25% を取るので、 クラス間バランスはエリア横断でほぼ均一になる。 意味のある違いは「どの場所がそのクラスに分類されているか」 (= 空間分布) にあり、 それは図 3 と次節 (small multiples) で見る。

図と読み取り (図 11: 地形タイプ別 small multiples)

図 11: 地形タイプ別の空間分布 small multiples (3 エリア × 4 タイプ)

読み取り:

• 尾根 (赤列): 3 エリアとも町域全体に網目状に広がる。世羅は線状に連続、 熊野は山地北側に集塊、坂は短い尾根が斑状に分散。
• 谷 (青列): 水系のネットワーク構造がはっきり可視化される。 特に世羅町は河川沿いに谷が連続的に走り、樹枝状の流域パターンが読める。
• 急傾斜 (黒列): 全エリアで少数だが、熊野町は山地に集塊、 坂町は南東部 (= 豪雨被災地) に集塊する明瞭なパターン。
• 平地 (灰列): 3 エリアでほぼ 0% に近い (band_std P25 がそもそも低いため)。 この分類軸は本データセットでは情報量が少ない。

表 (地形タイプ構成比 wide)

分析 4: RGB バンド分布と HSV 散布

狙い

3 エリアのR, G, B 各バンドの値分布を直接比較し、 HSV 散布図で「色相 × 彩度」の 2 次元パターンを観察する。 バンドごとの違いは「どの地形指標がエリアごとに強いか」を語る。

手法 (リテラシレベル解説)

2 つの可視化を併用:

• バンドごとのヒスト: R, G, B それぞれを 0-255 で 4 刻みのヒストにし、3 エリア重ね。 各バンドが「どこにピークを持ち、どんな裾を持つか」を確認。
• HSV 散布: RGB を HSV (Hue, Saturation, Value) に変換し、Hue × Sat 平面に散布。 各点を実際の RGB 色で描画する (= 各セルが実物のどの色かを視認可能)。

HSV とは: RGB の代替表現で、 H = 色相 (0=赤, 0.33=緑, 0.66=青), S = 彩度, V = 明度。 S = (max - min) / max の式で「色がどれだけ偏っているか」を 0-1 で表現する。 グレー (R=G=B) は S=0、純赤 (R=255, G=B=0) は S=1。

実装

1
2
3
4
5
6
7
8
9

import matplotlib

# RGB を HSV に変換 (各値 0-1 にスケール)
rgb_norm = (rgb_arr / 255.0).reshape(-1, 3)
hsv = matplotlib.colors.rgb_to_hsv(rgb_norm.reshape(1, -1, 3)).reshape(-1, 3)
hue, sat, val = hsv[:, 0], hsv[:, 1], hsv[:, 2]

# 散布図: 各点の色は実際の RGB そのまま
ax.scatter(hue, sat, c=rgb_norm, s=2, alpha=0.5)

図と読み取り (図 4: RGB バンドヒスト)

なぜこの図か: 3 バンドを別々に並べて、各エリアの値域シフトを見るため。 散布図ではなくヒストにすることで、各バンドの中央値・裾の重さが分かりやすい。

図 4: R/G/B 各バンドの 3 エリア比較ヒストグラム

読み取り:

• R バンド (左): 3 エリアとも 180-220 にピーク。世羅 (μ=202) が最右、坂 (μ=176) が最左。 R が低い = 全体的に赤の成分が弱い = 尾根の凸成分が抑えられている。
• G バンド (中): ほぼ同形だがピーク位置が R より約 10 だけ低い。 傾斜による明度低下を表現する G の役割が見える。
• B バンド (右): 値域は R, G より幅広い (std が大きい)。 = 谷の凹成分は地形によって大きく変動する。
• 坂町だけ裾が左に伸びる (R, G, B とも 0-100 域にも分布) ⇒ 暗い色 (急斜面の影) を多く含む = 急傾斜地形が広がる証拠。

図と読み取り (図 10: HSV 散布)

なぜこの図か: RGB の 3 軸では分布パターンが見えづらい。 HSV に変換して 2 次元の Hue × Sat 平面で見ると、 「赤系の凸地形」「青系の凹地形」「灰系の平地」のクラスタリング構造がよく見える。

図 10: HSV 散布図 — Hue × Sat 平面 (各点を実 RGB 色で描画)

読み取り:

• 3 エリアとも Hue が 0 (赤系) と 0.6-0.7 (青系) の 2 帯に分かれて分布。 仮説 H5 (色相は 2 峰性) を視覚的に支持。
• Sat 軸は 3 エリアとも 0 - 0.2 の狭い範囲に圧縮 (= パステル色)。 坂町だけ 0.25 付近に裾を持つ (= 微地形コントラスト強)。
• Hue 0.0 (赤) のクラスタは「尾根」、0.66 (青) は「谷」、0.33 (緑) は「傾斜」、 中間 (灰) は「平地」と読める。CS立体図の配色設計が散布図上で 3 主帯に分かれるのが視覚化される。
• Hue 標準偏差は世羅 0.057 < 熊野 0.136 < 坂 0.252と大幅に増加。 坂町は色相が広く分散 = 地形タイプが多様。

分析 5: 警戒区域 × CS立体図 — ポリゴン×ラスタの空間結合 (本記事の核心)

狙い (本記事の核心)

土砂災害警戒区域 (yellow) と特別警戒区域 (red) のポリゴンを CS立体図のラスタ座標系に rasterize して、 「警戒区域内 vs 区域外で CS立体図の RGB 平均値はどう違うか」を量的に比較する。 これは CS立体図単独・警戒区域単独では出せない本記事の独自貢献。

手法 (リテラシレベル解説)

「ポリゴン → ラスタ化 (rasterize)」とは: 不規則な形のポリゴン (警戒区域) を「画素ごとに 0/1 のマスク」に変換する操作。具体的には:

1. 警戒区域ポリゴンを rasterio の features.rasterize に渡し、 対象ラスタと同じ shape・transform で描画。
2. 得られたバイナリマスク (1=区域内, 0=外) と CS立体図ラスタをセル単位でペア化。
3. 各 zone (特別警戒 / 警戒 / 区域外) に分けて R, G, B, Sat の平均を計算。
4. 3 ゾーン × 3 エリア × 4 指標 (R, G, B, Sat) のクロス表を作る。

パッケージ: rasterio.features.rasterize。geopandas.overlay や shapely.contains よりもラスタ向けに最適化されており、 60K セル × 数千ポリゴンでも数秒で完了。

実装

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

from rasterio.features import rasterize
from rasterio.transform import from_bounds

def rasterize_polys(polys_gdf, area):
    """警戒区域 polygons をエリア overview の (H, W) shape にラスタ化"""
    bb = area["bounds"]
    transform = from_bounds(bb.left, bb.bottom, bb.right, bb.top,
                            area["arr_w"], area["arr_h"])
    shapes = ((geom, 1) for geom in polys_gdf.geometry)
    mask = rasterize(shapes, out_shape=(area["arr_h"], area["arr_w"]),
                     transform=transform, fill=0, dtype="uint8")
    return mask.astype(bool)

# 警戒区域 / 特別警戒 / 区域外 で各バンド平均
yellow_mask = rasterize_polys(area["sed_yellow"], area)
red_mask = rasterize_polys(area["sed_red"], area)
in_yellow = yellow_mask & area["valid_mask"]
in_red    = red_mask    & area["valid_mask"]
out_warn  = ~yellow_mask & ~red_mask & area["valid_mask"]

for zone, m in [("特別警戒", in_red), ("警戒", in_yellow), ("区域外", out_warn)]:
    R_mean = rgb[0][m].mean()
    G_mean = rgb[1][m].mean()
    B_mean = rgb[2][m].mean()
    Sat_mean = sat[m].mean()

図と読み取り (図 6: 警戒区域 × CS overlay)

なぜこの図か: 警戒区域がCS立体図のどの色域に集塊しているかを視覚的に確認。 ポリゴンの透明度を調整して下地の RGB 段彩を見せつつ、警戒域 (yellow/red) を半透明で重ねる。

図 6: CS立体図 + 土砂災害警戒区域 (yellow/red) の重ね合わせ

読み取り:

• 世羅町: 警戒区域 (黄) が町域全体に分散。地すべり警戒区域が広い。 特別警戒 (赤) は山地の急斜面に集塊。
• 熊野町: 警戒区域は山地に集中。平地住宅地はほぼ警戒区域外。 花崗岩マサ土地形の急傾斜面が警戒域に明確に対応。
• 坂町: 警戒区域 (黄) と特別警戒 (赤) が町域の山地側を広く覆う。 南西の住宅地と山地の境界に特別警戒が集中 = 平成30豪雨で土石流が 発生した地形と完全に対応。
• 3 エリアともに警戒区域は急斜面 (暗い色域) や谷 (青色域) に重なる傾向が見える。

図と読み取り (図 7: RGB 平均比較 grouped bar)

なぜこの図か: 視覚的な印象 (図 6) を量的に裏付ける。 3 エリア × 3 ゾーン × 3 指標 (R, R-B, Sat) のクロス比較。

図 7: 警戒区域 内/外 の RGB / R-B / Sat 平均比較

読み取り (重要発見):

• 仮説 H4 (警戒区域内の Sat は外より大きい) は反証された! 3 エリアとも警戒区域内の Sat 平均は区域外より低い: 世羅 内 0.068 vs 外 0.075, 熊野 0.050 vs 0.091, 坂 0.090 vs 0.100。
• 同様に R-B 平均も警戒区域内が低い: 世羅 内 +13.1 vs 外 +14.7, 熊野 +9.3 vs +17.0, 坂 +15.5 vs +16.7。
• これは何を意味するか? 警戒区域は谷底や急斜面の影に多く重なるため、 R (尾根) より B (谷)が強く出る = R-B 差が小さくなる。 尾根の凸地形は警戒区域として指定されにくい (= 雨水が集まりにくいため)。
• これは土砂災害指定の物理的論理 (= 谷頭・急斜面で土石流リスク大) が CS立体図の RGB に量的に反映されていることを示す重要な発見。
• 仮説 H4 は反証されたが、反対方向の関係が量的に確認できたため、 研究としては新しい知見が得られた (反証も研究的価値あり)。

表 (警戒区域 内/外 の RGB 平均)

読み取り: 警戒区域内のセル数は世羅町で 14,179, 熊野町で 6,445, 坂町で 2,197。 特別警戒区域は各エリアで 390-2,503 セルと小規模だが、量的比較は十分可能。 diff_RB_mean (R-B 平均)を見ると: 区域内は +9〜+15、区域外は +14〜+17 と一貫して 区域外の方が高い (= 区域外は尾根成分が強い、区域内は谷成分が強い)。

分析 6: 中央サンプルと CS立体図のしくみ概念

狙い

各エリアの中央 1024×1024 セルを本来解像度で読込み、CS立体図の細かい構造を体感する。 さらに概念図で「CS立体図がどうやって作られているか」を学習者に理解させる。

手法 (リテラシレベル解説)

rasterio.windows.Window で各 TIFF の中央 1024×1024 セルだけを読込み、 本来解像度で表示。1m なら 1024m × 1024m の正方領域、50cm なら 512m × 512m の領域。 3 バンド同時読込で RGB 画像として可視化。

実装

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

from rasterio.windows import Window

RAW_SAMPLE = 1024
for area in AREAS:
    with rasterio.open(area["tif"]) as ds:
        cx, cy = ds.width//2, ds.height//2
        win = Window(cx-RAW_SAMPLE//2, cy-RAW_SAMPLE//2,
                     RAW_SAMPLE, RAW_SAMPLE)
        raw = ds.read(window=win).astype(np.uint8)  # (3, H, W)
        # imshow に渡せる (H, W, 3) に転置
        rgb_disp = np.stack([raw[0], raw[1], raw[2]], axis=-1)
        ax.imshow(rgb_disp)

図と読み取り (図 8: 中央サンプル)

なぜこの図か: overview 表示では微地形の細部が見えない。 本来解像度で見ると CS立体図の真価 (尾根線・谷線・崩壊跡地が筋状に明瞭) が分かる。

図 8: 中央 1024×1024 セルの本来解像度サンプル

読み取り:

• 世羅町中央 (左): 1024m × 1024m = 1km 四方の本来解像度。 樹枝状の谷ネットワークが淡い青で連続的に走り、その間に薄ピンクの尾根が網目状。 = 中山間の典型的な地形パターン。
• 熊野町中央 (中): 同じ 1024m × 1024m。急斜面の暗い影と谷頭の青い凹みが 明瞭に分かれる。建物 (=四角い影) も視認可能。
• 坂町中央 (右): 50cm 解像度なので 512m × 512m の狭い範囲だが 空間解像度は約 4 倍。微細な凹凸構造 (= 古砂防地形 / 平成30崩壊跡 の可能性) が 1m 版より明確に拾われている。
• これがCS立体図が「専門家の地形読み」に有効な理由 ⇒ 標高ラスタや傾斜図と異なり、 凹凸の符号が色で直感的に区別できる。

図と読み取り (図 9: CS立体図のしくみ概念図)

図 9: CS立体図の概念図 — 標高 → 曲率 + 傾斜 → RGB 合成

読み取り:

• 左図: 地形断面 (尾根 + 谷 + 平地) に CS立体図の配色を重ねた模式。 尾根 = 赤丸, 谷 = 青丸, 平地 = 中間色, 急斜面 = 暗い影 (G 低)。
• 右凡例: 5 種の RGB パターンが対応する地形タイプ。本記事の「地形タイプ自動分類」の根拠。
• CS立体図は単一指標 (例: 傾斜のみ, 標高のみ) よりも3 指標を 1 枚に圧縮している点が 独特。学習者には「色が複数の地形情報を担う」概念の良い例。

表 (中央サンプル統計)

読み取り: 中央サンプルの RGB 平均は 3 エリアで類似 (160-200) だが、 サイズ差 (= 実距離差)が解像度差を反映する。坂町は 50cm 解像度のため 同じ 1024×1024 でも 512m 四方しかカバーしない。 これは「ピクセル数 = 情報密度ではない」という重要な GIS リテラシ。 解像度で正規化して比較すべき。

仮説検証と考察

本記事で立てた 6 つの仮説 (H1〜H6) と、3 エリアの CS立体図 RGB 解析 + 警戒区域空間結合の照合結果:

判定サマリ: 支持 5 / 部分支持 0 / 反証 1 (全 6 仮説中)

主な発見 (3-5 点)

1. H4 (警戒区域内の Sat は外より大きい) は反証されたが、その反対方向が 3 エリア共通で量的に確認された: 警戒区域は CS立体図で「より低い R-B 差・より低い Saturation」を持つ。 = 警戒区域は谷底や急斜面の影に多く重なり、尾根の凸地形は警戒区域として 指定されにくいという物理的論理が CS立体図の RGB に明瞭に刻まれている。 これは本記事の最も重要な発見であり、 「土砂災害指定の地形物理学」がオープンデータの 2 dataset 結合で再現できることを示す。
2. H3 (坂町は崩壊跡を含む急傾斜) が強支持された: 坂町の Saturation 平均 0.096 は 3 エリアで最大。 さらに R/G/B の標準偏差も最大 (R-std=42.0 vs 世羅 15.7)。 平成30年7月豪雨の被災地形が CS立体図に量的に保存されていることを示す。
3. H2 (熊野町は急峻地形) は支持されたが差は微小 (Sat 0.083 vs 世羅 0.074)。 熊野は花崗岩マサ土地形で急傾斜面が多いと予想したが、実データ上は世羅との 差は限定的。これは「都市近郊の山地」の急峻度が中山間の丘陵地形と大差ない という意外な発見 (発展課題 Z3 で再検証)。
4. H5 (色相は 2 峰性) が支持された: 全エリアで Hue が 赤帯 (0.0) と青帯 (0.66) の双方に分布。Hue 標準偏差は世羅 0.057 < 熊野 0.136 < 坂 0.252と 段階的に増加 = エリアの地形多様度を量的に区別する指標として機能。
5. H6 (50cm 版は微細構造強) が支持された: 坂町 (50cm) の Sat 標準偏差 0.060 は 1m 版 (世羅 0.043, 熊野 0.060) より大きい。 高解像度ほど微地形のばらつきを拾うことを実証。 ただし熊野 (1m) も同等の Sat std を持つため、解像度差より地形類型差の方が Sat std に効く可能性も残る。

3 エリアの地形特性比較 (本記事の量的サマリ)

考察

本記事の主たる発見は、「警戒区域の指定論理が CS立体図の RGB に量的に刻まれている」 ことが 3 エリア共通で確認されたことである。仮説 H4 は反証されたが、 反対方向の関係 (警戒区域は谷底寄り = R-B 平均が低い) が一貫して観察された。 これは「土砂災害警戒区域 = 谷頭 / 急斜面の影」という地形物理的論理の量的証拠であり、 オープンデータの 2 dataset 結合だけで再現可能なことを示している。

同時に、解像度差 (1m vs 50cm)は overview レベルでは認識困難だが、 本来解像度のサンプル (図 8) では 50cm の優位性が見える。微細な凹凸 (古砂防地形・ 古崩壊跡など) は 50cm 版でしか拾えない可能性が高く、地形リテラシ研究には 50cm 版が必須。

3 エリアの絶対値の RGB スケールには大差はないが、標準偏差・彩度には明確な 階層 (世羅 < 熊野 ≦ 坂) が見える。これは「地形微細構造の強さ」を CS立体図 RGB から 量的に評価する手法として有効。仮説を立てる前にデータの値域を観察する重要性を示す事例。

パーセンタイル適応閾値での地形タイプ分類 (尾根/谷/急/平地/複雑) は各エリア内の相対分類として 意味があるが、絶対値での横断比較には不向き。絶対閾値 vs 相対閾値の使い分けは リテラシ上の重要論点 (発展課題 Z2)。

発展課題

結果 X1 → 仮説 Y1 → 課題 Z1: 警戒区域 ⇔ 谷地形の量的検証 (本記事の発見の深掘り)

結果 X1: 3 エリア共通で警戒区域の R-B 平均は区域外より低い (= 谷成分強)。 H4 を反対方向で支持する量的傾向。

新仮説 Y1: 警戒区域内の「谷タイプ (B>R)」セルの比率は区域外より明確に高い。 50%以上の警戒域が谷タイプに分類される。

課題 Z1: 警戒区域内の地形タイプ構成比 (尾根/谷/急/平地) を本記事のラスタ化マスクで集計し、 区域外と並列比較。クロス表で「警戒域は谷タイプが何 % か」を量化。 さらに土石流・急傾斜・地すべりの 3 種別ごとに比較すると、各災害タイプの 地形特性が CS立体図でどう違うか分かる。

結果 X2 → 仮説 Y2 → 課題 Z2: 標高ラスタからの直接的な曲率計算

結果 X2: 本記事の「地形タイプ自動分類」は CS立体図の RGB から逆推定したものだが、 適応閾値が必要なため絶対比較が困難。

新仮説 Y2: 標高ラスタ (DTM, dsid 1635/1636) から直接 曲率・傾斜を計算した方が 量的に精確で、絶対閾値での横断比較が可能。

課題 Z2: 標高 GeoTIFF を読み込んで scipy.ndimage.laplace で曲率を、 numpy.gradient で傾斜を直接計算。本記事の RGB 推定と相関係数で照合。 高い相関が出れば「CS立体図の配色は曲率・傾斜の忠実な可視化」と確認できる。

結果 X3 → 仮説 Y3 → 課題 Z3: 多自治体一括の地形指標化

結果 X3: 本記事は 3 自治体のみ。Sat 平均 (= 微地形コントラスト) で 世羅 < 熊野 < 坂 の階層が見えたが、サンプル数が少なすぎて地形類型と Sat の関係が 仮説止まり。

新仮説 Y3: 全 22 自治体 (1m 版) で Sat 平均を計算すると、 中山間 / 都市近郊 / 沿岸被災地 / 島嶼 などの地形類型クラスタが自然分類できる。

課題 Z3: 22 自治体すべての CS立体図を順次 DL → overview 読込 → Sat 平均計算で 22 値の分布を作る。k-means や階層クラスタリングで類型化。 さらに各自治体の警戒区域率と Sat 平均を散布図にプロット ⇒ 「Sat 高 = 警戒区域率高」の正相関が出れば 「CS立体図の Sat は土砂災害リスクの代理指標」という新提案ができる。

結果 X4 → 仮説 Y4 → 課題 Z4: 古崩壊地形の自動検出

結果 X4: 坂町 (H30豪雨被災地) は 3 エリアで最高の Sat 標準偏差を示した。 微細な凹凸 (= 過去の崩壊地形 / 古砂防地形) が CS立体図に保存されていると示唆。

新仮説 Y4: 「Sat の局所変動が大きい場所」は古崩壊地形の可能性が高い。 過去の航空写真や砂防履歴と照合可能。

課題 Z4: 各セルの周囲 5×5 セルでの Sat 局所標準偏差を計算し、 高い局所 std セルを抽出。これが既知の H30 崩壊跡 (地理院地図の被災地マップ) と 空間的に一致するかを精度評価。一致が高ければ「Sat 局所 std による古崩壊地形自動検出」を提案できる。 これは砂防研究・防災マップ更新に直接応用可能。

結果 X5 → 仮説 Y5 → 課題 Z5: L36/L37 との 3 軸結合 (LiDAR 派生物の三脚)

結果 X5: 本記事は CS立体図のみ単独扱い (合体禁止)。 L36 樹高 / L37 点群密度との関係は未検証。

新仮説 Y5: 同じ航空レーザ点群から派生した 3 ラスタ (樹高 / 点群密度 / CS立体図) は 3 軸の地形・植生指標として機能。 3 軸の組合せで「裸地崩壊」「森林急斜面」「都市平地」等の 9-12 クラス自動分類が可能。

課題 Z5: X 系の横断記事として、3 ラスタを同一座標系で空間結合し、 GMM (Gaussian Mixture Model) で 9-12 クラスタリング。各クラスタが 何の地形/植生に対応するかを衛星 NDVI や建物 footprint で検証。 これは LiDAR ファミリ全体を統合する集大成記事になる。

結果 X6 → 仮説 Y6 → 課題 Z6: 過去被災地のセル単位検出

結果 X6: 坂町の特別警戒区域 (0.40 km²) と CS立体図の急斜面 (G 低) は 重なるが、過去の被災実績との照合は未実施。

新仮説 Y6: 平成30年7月豪雨で実際に崩壊が発生したセルは CS立体図で 特定の RGB 特徴 (Sat 中, 谷タイプ, 急傾斜) を持つ。RGB 特徴量で機械学習による被災予測が可能。

課題 Z6: 国交省の H30 豪雨被災調査データ (発生位置 GIS) を入手し、 CS立体図の RGB 特徴でロジスティック回帰や Random Forest による被災予測モデルを訓練。 特徴量重要度から「どのバンドが被災予測に効くか」を判定。 高い AUC が出ればCS立体図ベースの土砂災害発生予測モデルとして実装可能。