L39 地形傾斜図 2件統合分析 — DEM由来の傾斜角分布が描く地形リスク構造 と 土砂災害警戒区域の物理的論理

データ取得手順

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

実行コマンド:

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

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

学習目標と問い

本記事の独自性 — 砂防工学標準閾値 × 警戒区域 × L38 整合性

地形傾斜図はオープンデータとして広く公開されていますが、本記事は以下の 3 軸で独自分析を行います:

• 砂防工学標準 5 階級分類: 平地<5° / 緩 5-15° / 中 15-30° / 急 30-45° / 急峻>45° の 標準閾値で分類し、3 エリアの「地形リスク構造」を比較。30°は急傾斜地崩壊危険箇所の 法的指定基準でもある。
• 警戒区域との空間結合: 土砂災害警戒区域 (yellow) と特別警戒区域 (red) のポリゴンを 傾斜ラスタの座標系に rasterize し、区域内 vs 区域外で平均傾斜・階級構成を比較。 「警戒区域は急傾斜に指定されているか」を量的に検証する。
• L38 CS立体図との Pearson 相関: 本記事の傾斜 (L39) と L38 の RGB 派生 Saturation を 同自治体・同 overview shape で読み込み、ピクセル単位の Pearson 相関を計算。 両者は同じ DEM から派生した別表現なので、強い正相関が出ることを期待 (H7)。

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

このシリーズは 1m 版 / 50cm 版 の 2 解像度系統で、各自治体ごとに GeoTIFF (Float32 度数) を切出して公開する構造。 1m 版は 22 自治体、50cm 版は 12 自治体に対応。

• 1m 版 (dsid 1630): 22 自治体 GeoTIFF。 本記事は世羅町 (resource_id 177725)と熊野町 (176866)を採用 (= L36/L37/L38 と同自治体ペアで LiDAR ファミリの相補比較が可能)。
• 50cm 版 (dsid 1631): 12 自治体 GeoTIFF。 熊野町は 50cm 版に存在しないため、代わりに坂町 (176905)を採用。 坂町は平成 30 年 7 月豪雨の被災地で、土石流災害が多発した地域。 L38 でも同じく坂町 50cm 版を採用しており、本記事と L38 は完全な自治体ペアを構成する。
• 3 エリアは 1m × 2 + 50cm × 1 の構成で、解像度差と地形類型差を同時に観察可能。

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

研究の問い (RQ)

広島県の地形傾斜図 2 解像度 (1m / 50cm) は、それぞれ 世羅町 (中山間/林業地), 熊野町 (都市近郊/山岳混在), 坂町 (平成30豪雨被災地) という異なる地形リスク特性のエリアを公開している。 3 エリアの傾斜分布を 5 階級で量化し、地形リスク構造を比較する。 さらに土砂災害警戒区域との空間結合で「警戒区域 = 急傾斜地」という法的論理が データに刻まれているかを量的に検証する。最後にL38 CS立体図との物理的整合性を ピクセル相関で確認する。

1. 2 dataset の解像度・面積・有効データ率・ピクセル数の構造比較は?
2. 3 エリアの傾斜角ヒストグラム (2.5° bin) はどう違うか?
3. 砂防工学標準 5 階級の構成比は?
4. 30°超急傾斜セルは空間集塊するか? 8 近傍同種率は?
5. 土砂災害警戒区域の内 vs 外で平均傾斜・階級構成はどう違うか?
6. 解像度差 (1m vs 50cm) で急傾斜検出率はどう変わるか?
7. L38 CS立体図 Saturation と本記事の傾斜は正相関するか?

仮説 H1〜H7

• H1 (世羅町は緩傾斜丘陵): 平均傾斜 ≤ 12°, 30°超セル比率 ≤ 12%。 中山間の植林地は緩斜面が広い。
• H2 (熊野町は急峻地形): 平均傾斜 ≥ 18°, 30°超セル比率 ≥ 25%。 花崗岩マサ土地形で急傾斜面が多い。
• H3 (坂町は最急峻): 平均傾斜が 3 エリアで最大、急峻 (>45°) 比率も最大。 平成30豪雨で多数の土石流が発生した急斜面。
• H4 (警戒区域は急傾斜): 警戒区域内の平均傾斜は区域外の1.5 倍以上。 急傾斜地崩壊危険箇所の法的論理がデータに刻まれているはず。
• H5 (急傾斜は集塊): 30°超セルは空間的にランダムではなく、 8 近傍同種率 ≥ 0.85 で強集塊する。
• H6 (50cm 版は急傾斜検出率高): 坂町 (50cm) は世羅・熊野 (1m) より 30°超セル比率が高い (高解像度で局所急斜面を解像)。
• H7 (L38 CS Sat と正相関): 同自治体ピクセル比較で 傾斜 (L39) と CS立体図 Saturation (L38) の Pearson r ≥ 0.4。 両者は同じ DEM から派生しているので強い物理的整合性。

到達点

3 エリアの地形傾斜図を「3 自治体の代表エリア」として量的に解析し、 合計 633 M セルの Float32 傾斜データを overview ピラミッドで効率的に概観する。 さらに土砂災害警戒区域 Shapefile (~ 2,000 ポリゴン) との空間結合で 「法的指定論理の量的検証」を行い、L38 CS立体図とのピクセル相関で 「LiDAR ファミリ内の物理的整合性」を実証する。 これは傾斜単独・警戒区域単独・CS立体図単独では不可能な、3 軸統合分析である。

使用データ

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

データ仕様 (両 dataset 共通)

• 形式: GeoTIFF, 1 バンド Float32 = 傾斜角 (度, 0-90)
• CRS: EPSG:6671 (JGD2011 / Japan Plane Rectangular CS III) — m 単位の平面直角座標
• nodata: 1m 版は -9999.0, 50cm 版 (坂町) は -9999.0
• 値の意味: 各ピクセルの局所傾斜角 (度)。0=水平、90=垂直崖
• 派生過程: 標高 (DTM) → 3×3 セル窓で勾配計算 → arctan で度数化 → ラスタ出力
• 切出単位: 各自治体の行政界 + バッファで矩形切出 (L36/L37/L38 と同様)
• 用途: 土砂災害感受性評価、急傾斜地崩壊危険箇所の同定、農地適性評価、林道計画
• ライセンス: CC-BY 4.0 (測量法 44 条に基づく公共測量 2 次著作物)

L36/L37/L38 との関係 (本記事は単独, 参考併置のみ)

3 自治体の特性比較

データ品質メモ

• Float32 度数: 各セル 0-90 の浮動小数点。0 = 水平面、90 = 垂直崖。 実データ上は 90° 近傍は稀 (= 真の崖は地表で稀)。
• 1m 版は overview あり (2/4/8/16/32/64/128), 50cm 版も overview あり。 本記事は overview /1/32〜/1/64 で読み込み、メモリを抑える。
• 傾斜計算法は明示されていないが、結果分布から Horn 法 (3×3 窓) または類似と推定。
• 解像度差: 1m メッシュより 50cm メッシュの方が局所急傾斜を細かく解像する (= 同じ崖でも 50cm では真の傾斜に近い 50-70°が出るが、1m では平均化されて 30-50°になる)。
• nodata の扱い: 行政界外 = nodata。本記事の overview 読込では average resampling のため 端のセルは nodata と有効値の混合になり、わずかにバイアスが入る (ただし全体傾向には影響しない)。

ダウンロード

生データ (DoBoX 直リンク)

• 地形傾斜図（1m）_世羅町: 直 DL ZIP / DoBoX resource 177725 / dsid 1630
• 地形傾斜図（1m）_熊野町: 直 DL ZIP / DoBoX resource 176866 / dsid 1630
• 地形傾斜図（50cm）_坂町: 直 DL ZIP / DoBoX resource 176905 / dsid 1631

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

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

L38 CS立体図 TIFF (H7 の相関検証に必要, optional)

本記事は L38 と同じ自治体ペアを採用しています。 data/extras/L38_cs_terrain/ 配下の TIFF (CS立体図 RGB) があると自動的に ピクセル相関を計算します。L38 を先に実行することで全機能が有効になります。 ただし L38 TIFF が無くても L39 単独で全分析が走ります (H7 の判定のみ NA になる)。

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

• L39_series.csv — 3 dataset の一覧 (代表自治体・解像度・サイズ)
• L39_tiff_meta.csv — TIFF メタ (CRS, bounds, セル数, overview, dtype, nodata)
• L39_basic_stats.csv — 傾斜の平均・中央値・std・パーセンタイル・階級率
• L39_slope_classes.csv — 5 階級の構成比 (平地/緩/中/急/急峻)
• L39_slope_histogram.csv — 2.5° bin の傾斜ヒストグラム
• L39_sediment_slope_compare.csv — 警戒区域 内/外 の傾斜統計比較
• L39_zone_class_breakdown.csv — zone × 階級の細分構成
• L39_steep_clustering.csv — 30° 超急傾斜セルの集塊性 (8 近傍同種率)
• L39_l38_correlation.csv — L38 CS立体図 Saturation との Pearson 相関
• L39_hypothesis_results.csv / L39_hypothesis_results.json — H1〜H7 検証結果

図 (本記事生成 PNG)

• L39_fig1_dataset_overview.png — 3 dataset カード + 傾斜ヒスト
• L39_fig2_slope_maps.png — 傾斜ラスタ主題図 (3 エリア, グラディエント)
• L39_fig3_class_map.png — 5 階級分類マップ
• L39_fig4_class_composition.png — 5 階級構成比 (積み上げ)
• L39_fig5_sediment_overlay.png — 傾斜 × 土砂警戒区域 重ね合わせ
• L39_fig6_sediment_compare.png — 警戒区域 内/外 平均傾斜・p90・30°超比率
• L39_fig7_small_multiples.png — 階級別 small multiples (3×5)
• L39_fig8_center_samples.png — 中央 1024×1024 セル本来解像度サンプル
• L39_fig9_l38_correlation.png — L38 CS Sat vs 傾斜 散布 (Pearson)
• L39_fig10_concept.png — 傾斜計算と 5 階級の概念図
• L39_fig11_zone_class.png — zone別 5 階級構成 stack-bar

再現用 Python スクリプト

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

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

狙い

2 dataset (1m / 50cm) と 3 自治体 (世羅 / 熊野 / 坂) の規模・形状・傾斜分布の概形を 1 枚の絵で示す。傾斜ヒストグラムを重ねて地形リスク構造を最初に視覚化する。

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

L36-L38 と全く同じ rasterio スタイル。傾斜 TIFF は1 バンド Float32 で、各セルが 傾斜角 (度) を持つ。読込手順:

• overview ピラミッド: 大きい TIFF (数 GB) を直接読むとメモリが破裂するので、 ds.overviews(1) で利用可能な /1/2, /1/4, ... のピラミッドを取得し、 /1/32 〜 /1/64 を選んで読み込む。実効セルサイズ ≈ 32m。
• average resampling: 粗化時の補間方法。傾斜の平均をピクセル平均値とすることで 分布形状を保つ。最大値や最小値ではなく平均が適切 (極端値が紛れ込まない)。
• nodata マスク: 行政界外は ds.nodata 値もしくは NaN。 (arr != nodata) & (arr >= 0) & (arr <= 90.001) で有効セル抽出。

実装

↑ L39_terrain_slope.py 行 1768–1816

 1
 2
 3
 4
 5
 6
 7
 8
 9
1777
1778
1779
1780
1781
1782
1783
1784
1785
1786
1787
1788

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

        # 1 バンド読込 (傾斜は単一バンド) → shape (H, W)
        arr = ds.read(1, out_shape=(ov_h, ov_w),
                      resampling=Resampling.average).astype(np.float32)
        nodata = ds.nodata if ds.nodata is not None else -99999.0
        valid_mask = (~np.isnan(arr)) & (arr != nodata) & (arr >= 0) & (arr <= 90.001)
        area["arr"] = arr
        area["valid_mask"] = valid_mask

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

なぜこの図か: 3 dataset の規模と特性をカードで言語的に、 かつ 2.5° bin のヒストグラムで「傾斜分布の形」として量的に、同時に 1 枚で伝えるため。 階級境界 (5°/15°/30°/45°) を縦線で示し、5 階級分類の文脈を即座に与える。

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

読み取り:

• 3 エリアとも 0° 近傍にピーク (= 平地・緩傾斜が最頻ビン) を持つ右裾分布。 平均値は19.8° / 22.3° / 23.8° と 世羅 < 熊野 < 坂の順 (もしくは検証結果に応じて差がない場合もあり)。
• 世羅町は 30° 縦線より左側に分布が集中 = 緩傾斜が支配的。30° 超の右裾は薄い。
• 熊野町は世羅より右にシフトし、30° 縦線を跨ぐ分布が見られる。
• 坂町は 30° 〜 45° 帯に厚みがあり、急峻 (45° 超) も観測される ⇒ 平成30豪雨被災地の急斜面構造がデータに刻まれている可能性。
• p90 (90 percentile) は 3 エリアでばらつきがあり、地形リスク構造の差を端的に示す。

表 (基本統計)

読み取り: 平均・中央値・p90 が 3 エリアで明確に階層化される。 30° 超セル比率は世羅 14.6% / 熊野 32.2% / 坂 45.8%で、 急傾斜地崩壊危険箇所の法的指定基準 (30°) を超える面積比が地域差を端的に示す。 急峻 (>45°) 比率は世羅 0.09% / 熊野 0.19% / 坂 0.20%で、 極端な崖地形の存在を量化する。

分析 2: 傾斜ラスタ主題図化

狙い

傾斜の絶対値分布を 3 エリア並列で表示し、地形リスクの空間パターンを視覚的に理解する。 緑 → 黄 → 赤 のグラディエントで「平地は緑、急峻は赤」と直感的に読める。

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

連続値ラスタの主題図化は matplotlib.pyplot.imshow に vmin=0, vmax=60, cmap=... を渡すだけで実現できる。 本記事では5 階級境界 (5/15/30/45°) に対応する 6 色のグラディエントを LinearSegmentedColormap で構築し、緑 (平地) → 黄 (中傾斜) → 赤 (急峻) の 段彩を作る。これにより「色が傾斜カテゴリと対応」する視認性が得られる。

実装

↑ L39_terrain_slope.py 行 972–995

972
973
974
975
976
977
978
979
980

from matplotlib.colors import LinearSegmentedColormap
slope_cmap = LinearSegmentedColormap.from_list(
    "slope", ["#1a9850", "#a6d96a", "#fee08b", "#fdae61", "#d73027", "#7f0000"], N=256
)

# nodata は NaN 化して透明表示
arr_show = np.where(valid_mask, arr, np.nan)
im = ax.imshow(arr_show, cmap=slope_cmap, vmin=0, vmax=60,
               extent=area_bounds, origin="upper", interpolation="nearest")

図と読み取り (図 2: 傾斜主題図)

なぜこの図か: 5 階級分類 (図 3) の前に、まず連続値の生分布を見て 「どこが急斜面か」を直感的に把握する。グラディエント表示はカテゴリ化前の本来の信号。

図 2: 3 エリア 傾斜ラスタ主題図 (緑=平地 → 暗赤=急峻)

読み取り:

• 世羅町 (左): 全体に緑〜黄緑が支配的で、河川沿いの黒線 (谷底) と尾根筋の薄黄色が 見える程度。急斜面 (赤) は限定的な集塊として観察される。
• 熊野町 (中): 中央〜南東に黄〜橙の中傾斜域が広がる。 山地北側に赤い急傾斜帯が連続的に走り、花崗岩マサ土の急峻地形が可視化される。
• 坂町 (右): 全体に橙〜赤の急傾斜が支配的。南東部の住宅地と山地境界に 暗赤色 (急峻) が斑状に出現 ⇒ 平成30豪雨で土石流が発生した地形と整合。
• 3 エリアとも北 (上) が山地、南 (下) が平地・海岸寄りの傾向。 坂町は最も小さく、急峻地形の解像度が高い。

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

読み取り: overview /1/32 で実効セルサイズ 32-32m に揃え、 3 エリア間で比較可能な解像度に標準化している。坂町は 50cm × overview /64 で実効 32m となるため、解像度差は overview レベルでは消える (= 図 8 の生解像度サンプルでのみ見える)。

分析 3: 5 階級分類 (砂防工学標準閾値)

狙い

連続値の傾斜ラスタを砂防工学標準の 5 階級にカテゴリ化し、 3 エリアの「地形リスク構造」を比較可能な形にする。 階級境界は法的・経験的に意味のある値 (5° 平地基準, 15° 緩傾斜限界, 30° 急傾斜地法定基準, 45° 崖基準)。

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

連続値の傾斜を 5 種のラベル (1-5) に変換する。閾値は絶対値で固定するため、 3 エリア間で横断比較が可能になる (L38 のパーセンタイル適応閾値とは設計思想が異なる)。

限界と代替:

• 限界: 絶対値閾値はエリア横断比較に有利だが、地域特性 (花崗岩マサ土 vs 粘板岩等)に応じた 崩壊リスク評価には不十分。「30° 以上=リスク」という単純化はあくまで第一近似。
• 代替案: パーセンタイル相対閾値 (L38 で採用) で「そのエリア内の上位 25%」を急傾斜とする方法。 ただし横断比較は不可になる。
• 真の正解: 傾斜 + 方位 + 曲率 + 表層地質 + 降水量を統合した 土砂災害感受性指数 (Susceptibility Index)。発展課題で扱う。

実装

↑ L39_terrain_slope.py 行 1918–1951

 1
 2
 3
 4
 5
 6
 7
 8
 9
1927
1928
1929
1930
1931

SLOPE_CLASSES = [
    ("平地",   -0.001,  5.0),
    ("緩傾斜",   5.0, 15.0),
    ("中傾斜", 15.0, 30.0),
    ("急傾斜", 30.0, 45.0),
    ("急峻",   45.0, 90.001),
]

def classify_slope(arr, valid_mask):
    label = np.zeros(arr.shape, dtype=np.int8)
    for i, (_, lo, hi) in enumerate(SLOPE_CLASSES, start=1):
        m = (arr > lo) & (arr <= hi) & valid_mask
        label[m] = i
    return label

図と読み取り (図 3: 5 階級分類マップ)

なぜこの図か: 連続グラディエント (図 2) では「リスクのカテゴリ」が読みにくい。 階級分けすることで、「30°以上の急傾斜面はどこか」が即座に分かる。

図 3: 5 階級分類マップ

読み取り:

• 世羅町: 中傾斜 (黄) が 54% で支配的。 緑 (平地) は河川沿いの帯状、橙 (急傾斜) は山地に散在。 平地比率 4.1% で 3 エリア中で最低 = 全体的に丘陵地形。
• 熊野町: 中傾斜 (黄) と急傾斜 (橙) がほぼ同等の比率 (32% が 30° 超) で、 急傾斜面が広く分布する花崗岩マサ土山地の特徴を量化。
• 坂町: 急傾斜 (橙) が支配的で 46% に達する。 平地 (緑) は南西の住宅地に集塊し、その境界が急傾斜面に隣接する構造 = 平成30豪雨で土石流が住宅地を直撃した立地条件と整合。

図と読み取り (図 4: 階級構成比 stack-bar)

図 4: 5 階級構成比の積み上げ棒グラフ

読み取り: 30° 超 (急 + 急峻) の比率は世羅 14.6% / 熊野 32.2% / 坂 45.8%と 2-3 倍の階層差がある。逆に平地 (<5°)は世羅 4.1% / 熊野 6.9% / 坂 16.4% で、 坂町が最大 (= 海岸に近い扇状地・住宅地の平坦地が広い)。 これは「坂町=急峻」のステレオタイプ的予想とは逆方向で、平地と急傾斜が共存する地形構造が読み取れる。

表 (傾斜階級構成比 wide)

図と読み取り (図 7: 階級別 small multiples)

図 7: 階級別の空間分布 small multiples (3 エリア × 5 階級)

読み取り:

• 平地列 (緑): 世羅は河川沿いに細長く、熊野は南部市街地に集塊、坂は南西の住宅地と扇状地に集塊。 3 エリアで「平地の地理的位置」が産業構造と対応。
• 急傾斜列 (橙): 世羅は山地周辺に散在、熊野は山地北側に連続帯、坂は町域全体に支配的。 坂町で急傾斜面が住宅地のすぐ隣に存在することが視覚的に確認できる。
• 急峻列 (赤): 3 エリアで稀少 (0.1-0.2%) だが、それぞれ特定の谷頭・崖地形に集塊。 これらは古崩壊地形・現役崩壊地形の候補で、発展課題 Z3 の対象。

分析 4: 警戒区域 × 傾斜 — ポリゴン×ラスタの空間結合 (本記事の核心)

狙い (本記事の核心の 1 つ)

土砂災害警戒区域 (yellow) と特別警戒区域 (red) のポリゴンを傾斜ラスタの座標系に rasterize して、「警戒区域内 vs 区域外で平均傾斜・階級構成はどう違うか」を量的に比較する。 法的指定論理 (急傾斜地崩壊危険箇所基準) がデータに刻まれているかを検証する。

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

ポリゴン → ラスタ化 (rasterize): 不規則な形のポリゴン (警戒区域)を 「画素ごとに 0/1 のマスク」に変換する操作。L38 と同じ rasterio.features.rasterize を使用。

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

実装

 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 を傾斜ラスタの (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)]:
    v = arr[m]
    mean_slope = v.mean()
    p90 = np.percentile(v, 90)
    pct_ge_30 = (v >= 30).sum() / len(v) * 100

図と読み取り (図 5: 警戒区域 × 傾斜 重ね合わせ)

なぜこの図か: 警戒区域がどの傾斜域に重なるかを視覚的に確認。 傾斜段彩を下地、警戒区域を境界線で重ねる。

図 5: 傾斜 + 土砂災害警戒区域の重ね合わせ

読み取り:

• 3 エリアとも警戒区域 (黄/黒線) が町域に広く分散。
• 警戒区域の外縁が急斜面 (橙〜赤) に接しているように見えるが、 内部は中傾斜〜緩傾斜のセルも多く含む。
• これは警戒区域が「急斜面そのもの」ではなく、急斜面から土砂が到達する範囲 (扇状地・谷出口)を 含む形で指定されているためと考えられる ⇒ 後の量的比較で確認 (図 6)。

図と読み取り (図 6: 警戒区域 内/外 量的比較 grouped bar)

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

図 6: 警戒区域 内/外 の平均傾斜・p90・30°超比率

読み取り (重要発見):

• 仮説 H4 (警戒区域内の平均傾斜は外の 1.5 倍以上) は反証された! しかし反対方向の関係が 3 エリア共通で量的に確認された: 警戒区域内の平均傾斜は区域外より低い。
• 世羅 内 15.18° vs 外 20.05° (×0.76), 熊野 内 12.08° vs 外 24.47° (×0.49), 坂 内 17.64° vs 外 24.58° (×0.72)。すべてのエリアで内 < 外。
• これは何を意味するか? 警戒区域は急傾斜面そのものではなく、 急斜面から土砂が到達する範囲 (谷出口・扇状地・崖下)を含むように指定されているため。 = 「土石流が止まる場所」を保護するために指定されているので、地形上は緩傾斜に分類される。
• 特別警戒 (red) は警戒 (yellow) より平均傾斜が高い: 世羅 21.4 vs 15.2°, 坂 25.5 vs 17.6°。 これは特別警戒が「実際に建物倒壊リスクのある急斜面」に集中する結果と整合。
• 30° 超セル比率: 区域外 (世羅 15%, 熊野 39%, 坂 51%) が警戒区域内 (3-12%) より桁違いに高い。 ⇒ 「警戒区域 = 30° 超急斜面」という単純解釈は誤りであり、 警戒区域は影響範囲全体を含む点が量的に明確化された。
• これは「土砂災害警戒区域の物理的論理」がオープンデータの 2 dataset 結合で 再現できることを示す重要な発見。 仮説 H4 は反証されたが、反対方向の関係が新知見として得られた (反証も研究的価値あり)。

表 (警戒区域 内/外 の傾斜統計)

図と読み取り (図 11: zone別 5 階級構成 stack-bar)

なぜこの図か: 平均値だけでは分布の形が見えない。階級構成比を見れば 「警戒区域内は中傾斜が多いか平地が多いか」が分かる。

図 11: zone別 5 階級構成 stack-bar (3 エリア)

読み取り:

• 区域外 (左バー) は急傾斜 (橙) と中傾斜 (黄) の比率が高く、 とくに熊野・坂町は橙が支配的。
• 警戒 (中バー) は緩傾斜 (緑系) と中傾斜の比率が増え、急傾斜が減る。 = yellow 警戒区域は谷出口・扇状地・緩急遷移帯に多い。
• 特別警戒 (右バー) は中傾斜が支配的で、急傾斜の比率が警戒よりは増える。 特別警戒は実際の建物倒壊リスク場所 (= 急斜面下) に集中するため。
• 3 エリアで「区域外 → 警戒 → 特別警戒」と進むにつれて急傾斜比率が一度減って増える U字パターンが 見える ⇒ 警戒区域の指定論理 (土石流発生源 + 流下範囲 + 堆積範囲) を量化。

分析 5: 急傾斜地の空間集塊性 (8 近傍同種率)

狙い

30° 超急傾斜セルが空間的に集塊するか散在するかを量化する。 集塊性は地形リスク評価の重要指標で、集塊する=地続きの急斜面塊で 土砂災害リスクが連鎖的に高い、散在する=点状リスクで影響範囲が限定的、と読める。

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

本記事独自の「8 近傍同種率」を使う。これは Moran's I の計算量軽減版で、 全ペア相関を取らずに「8 近傍にも同種がいる確率」だけを集計する。

解釈:

• 1.0 に近い: 急傾斜セルは完全な塊で分布 (= 地続きの広域急斜面)
• 0.5 付近: 50% は単独、50% は塊 = 中程度の集塊
• 0.0 に近い: ランダム散在 (= 単独セルが散らばる)

実装

↑ L39_terrain_slope.py 行 2120–2280

 1
 2
 3
 4
 5
 6
 7
 8
 9
2129
2130
2131
2132
2133

def neighbor_same_rate(mask, valid):
    H, W = mask.shape
    pairs_total, pairs_same = 0, 0
    for dy in (-1, 0, 1):
        for dx in (-1, 0, 1):
            if dy == 0 and dx == 0: continue
            ma = mask[max(0,dy):H+min(0,dy), max(0,dx):W+min(0,dx)]
            mb = mask[max(0,-dy):H+min(0,-dy), max(0,-dx):W+min(0,-dx)]
            va = valid[max(0,dy):H+min(0,dy), max(0,dx):W+min(0,dx)]
            vb = valid[max(0,-dy):H+min(0,-dy), max(0,-dx):W+min(0,-dx)]
            both_valid = va & vb
            pairs_total += int(((ma | mb) & both_valid).sum())
            pairs_same  += int((ma & mb & both_valid).sum())
    return pairs_same / max(1, pairs_total)

図と読み取り

図 7 (small multiples の急傾斜列) を再度参照。各エリアの急傾斜塊の連続性を量化した結果が下の表。

表 (急傾斜地集塊性)

読み取り (重要発見):

• 仮説 H5 (8 近傍同種率 ≥ 0.85) は反証。実測値は世羅 0.601 / 熊野 0.696 / 坂 0.787。
• ただし解像度が高い (= 50cm 坂町) ほど同種率が高い傾向が明確 (世羅 1m 0.60 < 坂 50cm 0.79)。 = 高解像度ほど急斜面が連続的な塊として捉えられる。
• 0.6 程度の同種率は「中程度の集塊」に該当: 急傾斜セルは線状の崖や尾根筋に並ぶことが多く、 斜面方向には連続するが直交方向には不連続なため、隣接 8 セルすべてが急傾斜になる確率は限定的。
• 地理学的には、急傾斜は線形構造 (尾根線・崖線)として分布するため、 面的集塊性 (8 近傍) の指標では中程度になりやすい ⇒ 仮説 H5 の閾値設定が高すぎたという反省。
• これは「線状地形構造の集塊性は 8 近傍では過小評価される」という発見で、 発展課題 Z2 (方位を加味した集塊性) に繋がる。

分析 6: 中央サンプルと L38 CS立体図との Pearson 相関 (H7)

狙い

2 つの分析を並列で扱う:

1. 中央 1024×1024 セルを本来解像度で読込み、解像度差 (1m vs 50cm) を体感する
2. L38 CS立体図 Saturation との同自治体ピクセル相関を取り、両者の物理的整合性を検証 (H7)

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

本来解像度サンプル: rasterio.windows.Window で各 TIFF の中央 1024×1024 セルだけを 読込み、本来解像度で表示。1m なら 1024m × 1024m、50cm なら 512m × 512m の領域。

L38 との相関: 同自治体の L38 TIFF (CS立体図 RGB) を本記事の overview と同 shapeで 読み込み (rasterio の out_shape 機能で空間分解能を揃える)、 ピクセル単位で傾斜 (L39) vs CS立体図 Saturation (L38) の Pearson 相関を計算。 両者は同じ DEM から派生した別表現なので、強い正相関を期待 (H7)。

実装

↑ L39_terrain_slope.py 行 2178–2234

 1
 2
 3
 4
 5
 6
 7
 8
 9
2187
2188
2189
2190
2191
2192
2193
2194
2195
2196
2197
2198
2199
2200
2201
2202

import rasterio
from rasterio.windows import Window
from rasterio.enums import Resampling

# (1) 中央サンプル
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(1, window=win).astype(np.float32)

# (2) L38 との Pearson 相関
for area in AREAS:
    with rasterio.open(area["l38_tif"]) as ds_cs:
        rgb = ds_cs.read(out_shape=(3, area["arr_h"], area["arr_w"]),
                         resampling=Resampling.average).astype(np.float32)
        # CS Saturation 計算
        maxv = np.maximum.reduce(rgb)
        minv = np.minimum.reduce(rgb)
        sat = np.where(maxv > 0, (maxv - minv) / np.maximum(maxv, 1), 0)
    # 共通有効マスクで Pearson 相関
    both = area["valid_mask"] & (rgb.sum(axis=0) > 0)
    r = np.corrcoef(area["arr"][both], sat[both])[0, 1]

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

なぜこの図か: overview /1/32-/1/64 で粗化された 分析と異なり、本来解像度で見ると急斜面の細部構造が分かる。 1m と 50cm の差はここで初めて視認可能。

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

読み取り:

• 世羅町中央 (左): 1024m × 1024m の本来解像度。緩傾斜 (緑) が支配的で、 平均 15.1°。河川沿いに薄い橙 (中傾斜) が線状に並ぶ。
• 熊野町中央 (中): 同 1024m × 1024m だが平均 9.1°と 世羅より低い (= サンプル位置による偶然性)。中央サンプルは町中央 (山地中) なので 分析全域 (overview) の傾向と乖離することがある。
• 坂町中央 (右): 50cm × 1024 = 512m × 512m。狭い範囲だが急傾斜 (橙) が支配的で 平均 27.0°。50cm 解像度で細かい崖線・凹凸が拾われている。
• 3 エリアの中央サンプル平均は必ずしも全域平均と一致しない = サンプリングの代表性の問題。 全域分析 (overview) との対比で、教材的に重要な GIS リテラシ。

図と読み取り (図 9: L38 CS立体図 Saturation との Pearson 相関)

なぜこの図か: H7 (傾斜 vs CS Sat の正相関) を視覚的・量的に検証する。 散布図に線形回帰線を重ね、Pearson r を表示。

図 9: L39 傾斜 × L38 CS Saturation 散布図

読み取り (重要発見):

• 仮説 H7 (Pearson r ≥ 0.4 が 2 エリア以上) は強支持された: 世羅 r = +0.796, 熊野 r = +0.669, 坂 r = +0.824。 3 エリアすべてで r >= 0.6 の強い正相関。
• これは「CS立体図の Saturation は傾斜の代理指標として機能する」ことの量的実証。 両者は同じ DEM から派生しているため物理的に整合するのは当然だが、 実データでこれだけ強い相関が出るのはDoBoX の傾斜計算と CS立体図生成が同一 DEM で行われている証拠でもある。
• 線形回帰の傾きは傾斜 1° 増加で Sat が 0.001-0.003 増加程度 (= 傾斜 30° で Sat が 0.05 程度上がる)。これが「急斜面ほど CS図が彩度高」の量化。
• L38 と本記事は同自治体ペアを意図的に選択しているため、 この相関は本記事と L38 の記事間整合性を裏付ける研究的意義もある。
• 参考: slope vs darkness (CS図の暗さ): 世羅 r=+0.944, 熊野 r=+0.913, 坂 r=+0.321。 Saturation よりさらに強い相関 (世羅・熊野で r>0.9!) = CS立体図の暗さ (G 低)は 傾斜の最良の代理指標。これは L38 の概念図 (G=傾斜・明度) と整合。

表 (L38 相関)

図と読み取り (図 10: 傾斜計算と 5 階級の概念図)

図 10: 傾斜計算 (DEM → 勾配 → arctan) と 5 階級の概念図

読み取り:

• 左図: 地形断面に接線を貼り、各点の傾斜角を arctan で計算するイメージ。 緩 (10°)、中 (20°)、急 (40°) の 3 例を色で示す。
• 右凡例: 5 階級の閾値と地形的意味。30° = 急傾斜地崩壊危険箇所基準は法的閾値で、 45° = 経験的人家不可上限。
• 傾斜は単一指標 (例: 標高) より「土砂災害感受性」に直結する点が独特。 学習者には「傾斜は土砂災害評価の最も基本的な物理量」という概念の良い例。

仮説検証と考察

本記事で立てた 7 つの仮説 (H1〜H7) と、3 エリアの傾斜ラスタ解析 + 警戒区域空間結合 + L38 ピクセル相関の照合結果:

判定サマリ: 支持 4 / 部分支持 1 / 反証 2 (全 7 仮説中)

主な発見 (5 点)

1. H4 (警戒区域は急傾斜) は反証されたが、反対方向の関係が 3 エリア共通で量的に確認された: 警戒区域内の平均傾斜は区域外より低い (世羅 ×0.76, 熊野 ×0.49, 坂 ×0.72)。 これは警戒区域が「急傾斜面そのもの」ではなく「急斜面から土砂が到達する範囲 (谷出口・扇状地)を 含むように指定されているためで、土砂災害警戒区域の物理的指定論理を量的に明らかにした。 特別警戒 (red) は警戒 (yellow) より平均傾斜が高く、3 ゾーン間に明確な階層。 これは本記事の最も重要な発見。
2. H7 (L38 CS Sat との正相関) が強支持された: Pearson r = +0.796 / +0.669 / +0.824 (世羅/熊野/坂)。 3 エリアすべてで強相関 (r ≥ 0.6)。さらに slope vs CS darkness では r >= 0.9 (世羅・熊野) と 極めて強い相関が出た ⇒ L38 CS立体図の暗さ (G低) は傾斜の最良の代理指標。 これは記事間整合性 (L38 と L39 が同自治体・同 DEM 由来であること) の量的裏付け。
3. H3 (坂町は最急峻) が支持された: 平均傾斜 23.78° (3エリア最大), 30°超セル比率 45.8% (3 エリア最大)。50cm 解像度の効果も含めて、 平成30豪雨被災地の急斜面構造が量的に確認された。 ただし急峻 (>45°) 比率は 0.2% で熊野と同等 = 極端な崖は 3 エリアで稀。
4. H6 (50cm 版は急傾斜検出率高) が支持された: 30°超比率は世羅 14.6% / 熊野 32.2% / 坂 45.8% で 坂町 (50cm) が圧倒。ただし「解像度差」と「地形類型差」が交絡しており、 50cm の効果単独を分離するには同自治体 1m vs 50cmの比較が必要 (発展課題)。 集塊性 (8 近傍同種率) も坂町 0.79 で最大 = 高解像度ほど急斜面塊が連続的に解像される。
5. H5 (急傾斜の集塊性 ≥ 0.85) は反証: 実測 0.60-0.79。 ただし解像度に応じて単調増加する傾向は確認 ⇒ 仮説の閾値設定が高すぎた。 急傾斜は線形構造 (尾根線・崖線)として分布するため、面的指標 (8 近傍) では 過小評価される。これは「線状地形の集塊性は方向性を加味すべき」という方法論的発見。

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

考察

本記事の主たる発見は、「土砂災害警戒区域は急傾斜面そのものではなく、その影響範囲を含む」 ことが 3 エリア共通で量的に確認されたことである。仮説 H4 (警戒区域内の方が急傾斜) は反証され、 反対方向の関係 (警戒区域内 < 区域外) が一貫して観察された。これは 土砂災害警戒区域の指定目的 (= 土石流到達範囲・堆積範囲の保護) と整合する物理的論理であり、 オープンデータの 2 dataset 結合だけで再現可能なことを示している。 特別警戒 (red) は警戒 (yellow) より平均傾斜が高く、両ゾーンの機能差 (建物倒壊リスク vs 注意喚起)が 傾斜分布に明確に刻まれている。

L38 CS立体図 Saturation との Pearson 相関 (r=0.67-0.82) は記事間整合性の量的実証であり、 DoBoX の同シリーズ (LiDAR ファミリ) が同一 DEM から派生していることを裏付ける。 特に slope vs CS darkness の相関は r=0.9+ でCS立体図の暗さ成分は傾斜の最良代理と判明。 これは L38 が視覚解析向きの圧縮表現として有効であることの定量裏付けでもある。

解像度差 (1m vs 50cm) は overview レベルでは認識困難だが、本来解像度のサンプル (図 8) で 50cm の優位性が見え、また 30°超比率と集塊性の解像度依存性も確認された。 ただし同自治体での 1m vs 50cm 比較は本データセットの制約 (各自治体は 1 解像度のみ公開) で 不可能で、解像度差と地形類型差が交絡する点には留意が必要 (発展課題)。

急傾斜地の集塊性 (8 近傍同種率) は仮説 0.85 を下回り 0.6-0.79 だったが、 解像度に応じて単調増加する傾向は確認された。これは「線状地形構造の集塊性は面的指標で過小評価される」 という方法論的発見で、方位 (aspect) を加味した方向性集塊性を計算する発展課題に繋がる。

発展課題

結果 X1 → 仮説 Y1 → 課題 Z1: 警戒区域内外の傾斜分布の物理的解釈深掘り (本記事の発見の深掘り)

結果 X1: 3 エリア共通で警戒区域内 < 区域外の平均傾斜。特別警戒は警戒より急傾斜。 H4 を反対方向で支持する量的傾向。

新仮説 Y1: 警戒区域は「土石流の発生源 (急傾斜面 30°超) + 流下範囲 (15-30°) + 堆積範囲 (扇状地 5°以下)」 の 3 帯を組み合わせて指定されており、平均値だけでは見えない3 段の幾何構造を持つ。

課題 Z1: 各警戒区域ポリゴンを個別に傾斜分布解析。最大傾斜・最小傾斜・分布形状で 3 タイプ (急斜面卓越型 / 流下範囲型 / 堆積範囲型) に分類。 さらに土石流・急傾斜・地すべりの 3 種別ごとに比較すると、各災害タイプの幾何構造が見える。 これにより警戒区域の地形リスク類型を定量化できる。

結果 X2 → 仮説 Y2 → 課題 Z2: 方位を加味した方向性集塊性

結果 X2: 急傾斜の 8 近傍同種率は 0.6-0.79 で仮説 0.85 を下回ったが、 線状地形 (尾根線・崖線) の存在から方向性が示唆される。

新仮説 Y2: 急傾斜セルは等高線に沿った方向 (= 斜面の走向)には連続するが、 直交方向 (斜面方向) には不連続。方位を加味すれば線状集塊度は 0.9+ になる。

課題 Z2: 各セルで方位 (aspect) を計算 (DEM 勾配ベクトルの方向)、 等高線方向 (= 方位 ⊥) のセル間の連続性のみを評価する方向性集塊度を計算。 標準の Moran's I を方位重み付き隣接行列で実装。 これは「線状地形の集塊性」を初めて量化する方法論的貢献になる。

結果 X3 → 仮説 Y3 → 課題 Z3: 急峻地形と古崩壊地の同定

結果 X3: 45° 超急峻セルは 3 エリアで 0.1-0.2% と稀少。これらは特定の崖・崩壊跡に 集塊する可能性。

新仮説 Y3: 急峻セル + その周辺 (50m バッファ) の中傾斜セル混在域は、 古崩壊地形 (paleolandslide) の候補。過去の崩壊跡が地形に残る。

課題 Z3: 各エリアの 45° 超セルのクラスタリング (DBSCAN) で連続塊を抽出。 各クラスタの周辺バッファに含まれる傾斜分布を解析し、古崩壊地形類型 (新鮮崩壊・古崩壊・侵食地形) に分類。 H30 豪雨の崩壊実績 (国交省データ) と空間照合して精度評価。 高精度なら「傾斜ラスタからの古崩壊地形自動検出」を提案できる ⇒ 砂防研究の応用可能。

結果 X4 → 仮説 Y4 → 課題 Z4: 標高ラスタからの傾斜再計算と本データとの照合

結果 X4: L38 CS立体図 Saturation との r=0.67-0.82 は強相関だが完全ではない。 他の派生過程 (例えば曲率と傾斜の重ね合わせ) も影響している。

新仮説 Y4: 標高ラスタ (DTM, dsid 1635/1636) から直接 Horn 法で傾斜を再計算すると、 DoBoX 公式傾斜と Pearson r ≥ 0.99 の極めて強い一致が出る ⇒ 公式が Horn 法を使っている裏付け。

課題 Z4: 標高 GeoTIFF を読み込んで scipy.ndimage.sobel や手計算で 3×3 Horn 法を実装し、本データと相関係数で照合。さらに計算法の違い (Horn vs Zevenbergen) で結果がどう変わるかも比較。教材としても優れた発展課題。

結果 X5 → 仮説 Y5 → 課題 Z5: 多自治体一括の傾斜指標化と地形類型クラスタリング

結果 X5: 本記事は 3 自治体のみ。30°超比率で世羅 < 熊野 < 坂 の階層が見えたが、 サンプル数が少なすぎて地形類型の自然分類は仮説止まり。

新仮説 Y5: 全 22 自治体 (1m 版) で平均傾斜・p90・30°超比率を計算すると、 中山間 / 都市近郊 / 沿岸被災地 / 島嶼 などの地形類型クラスタが自然分類される。

課題 Z5: 22 自治体すべての 1m 版を順次 DL → overview 読込 → 5 階級構成計算で 22 個の特徴ベクトル (5次元) を作る。k-means や階層クラスタリングで類型化。 さらに各自治体の警戒区域率と平均傾斜を散布図にプロット ⇒ 「平均傾斜 vs 警戒区域率の正相関」が量的に確認されれば 「傾斜分布は土砂災害リスクの代理指標」という新提案ができる。

結果 X6 → 仮説 Y6 → 課題 Z6: 1m vs 50cm の同自治体比較 (解像度効果の純粋分離)

結果 X6: 30°超比率は坂町 (50cm) で最大だが、地形類型差と解像度差が交絡。

新仮説 Y6: 同自治体で 1m と 50cm の両方が公開されている自治体 (例: 三原市・呉市・尾道市等) を使えば、 解像度差の純粋効果が分離できる。同自治体での 50cm vs 1m は、急傾斜比率を 1.2-1.5 倍に増やす。

課題 Z6: 同自治体で両解像度公開の市町村を選び、ピクセル単位で空間整合した上で 分布差を計算。さらに 1m を 50cm に再投影 (resampling) して計算上の同 shape にしてから比較すると、 「解像度差の数値効果」が分離可能になる。これは GIS リテラシ教材として優れた題材。

結果 X7 → 仮説 Y7 → 課題 Z7: 過去被災地のセル単位検出 (ML 応用)

結果 X7: 坂町の特別警戒区域 (0.40 km²) と急傾斜セル (45.8%) は 相互に重なるが、過去の被災実績との照合は未実施。

新仮説 Y7: 平成30年7月豪雨で実際に崩壊が発生したセルは傾斜・方位・曲率の 特定パターン (傾斜 30-50° + 凹型曲率 + 北向き等) を持つ。これらの特徴量で機械学習による被災予測が可能。

課題 Z7: 国交省の H30 豪雨被災調査データ (発生位置 GIS) を入手し、 傾斜 + 方位 + 曲率 + L38 CS Sat の特徴でロジスティック回帰や Random Forest による 被災予測モデルを訓練。AUC 0.8 以上が出れば傾斜ベースの土砂災害発生予測モデルとして実装可能。 これは社会実装に直結する研究的発展。