地球を5mごとに切ると、隣町はどこへ行くのか

世界地図を、5mごとの標高で切ったらどうなるだろう。

海抜0m以上5m未満、5m以上10m未満、10m以上15m未満。同じ高さの土地だけを残して眺めたら、ふだん見ている地図とは違うまとまりが現れるかもしれない。

僕が気になったのは、色分けしたあとのつながり方だった。同じ標高帯の土地は、どれくらい連続しているのか。すぐ近くに見える二つの場所も、間に谷や丘があれば、別々のまとまりになるのだろうか。

そこで、標高帯を「層」、その中で連続する土地を「飛び地」として数えるコードをLayered Earthに実装した。最初の対象は関東周辺。NOAAのETOPO 2022を使い、実際に5mごとの地図と統計を出してみた。

高さを切ってから、つながりを数える

標高を持つ地図をDEM、数値標高モデルとして扱う。地図上の有効なセルの集合を $M$ 、セル $p$ の標高を $h(p)$ m、層幅を $\Delta=5$ mとする。

層番号は、標高を5で割って切り下げればよい。

L(p)=\left\lfloor\frac{h(p)}{\Delta}\right\rfloor

第 $k$ 層の領域は、

R_k=\{p\in M\mid k\Delta\le h(p)<(k+1)\Delta\}

となる。これは今回の実装で採用した定義だ。15mちょうどは第3層、20mちょうどは第4層に入る。上端を含めないので、境界上の土地を二重に数えずに済む。海抜−0.1mなら第−1層であり、0へ丸めない。今回の実データ実験では標高0m以上を対象にしたので、負の層は集計していない。数理モデル v0.1

この定義には、国境が出てこない。日本にも米国にも同じ第3層がある。その層の中で、隣り合うセルをたどって到達できるまとまりを抽出する。これが地図上の「飛び地」、数学では連結成分になる。

ただし、ここでいう飛び地は行政上の飛び地や実際の島ではない。「選んだ標高帯の中で、指定した規則に従ってつながるセルの集合」という、この計算の中での呼び名だ。

角で触れた土地を、つながったと呼ぶか

ラスター上の連結には選択がある。

4近傍なら上下左右のセルだけをつなぐ。8近傍なら斜めのセルもつなぐ。SciPyの連結成分ラベリングも、この接続規則を指定できる。

実地形ではない、2×2セルの小さな例を置いてみる。

標高・m	西側	東側
北側	1	6
南側	6	1

層幅5mなら、1mのセル同士と6mのセル同士が、それぞれ斜めに向き合う。4近傍では全体で四つの飛び地、8近傍では二つになる。地形の数字を一つも変えずに、飛び地数が半分になった。

この例は実装のテストにも入っている。同じDEM・同じ層幅・同じ対象セルなら、8近傍では接続を追加するだけなので、飛び地数は4近傍以下になる。

ただし、角で接するセルを通れるとみなすことは、幅のある道路が存在することを保証しない。逆に4近傍でも、セル内部の細い溝や段差は見えていない。

だから4近傍を基本にしつつ、8近傍との差も測る。この差が大きい地域は、生活圏が多いと即断するより先に、接続の定義や地図の解像度を疑いたくなる。

上へ向かう矢印は、まだ道路ではない

飛び地を点として、それぞれが地図上で接しているかを調べる。異なる層のセルが辺を共有していれば、低い層の飛び地から高い層の飛び地へ矢印を引く。

現在の実装では、成分抽出に8近傍を選んでも、この接触判定は辺の共有だけで行う。対角接触は接触辺に数えない。接触グラフの実装

すべての矢印で層番号が増えるので、矢印をたどって元の点へ戻ることはない。有向非巡回グラフ、DAGになる。もし一周できるなら層番号は一周して増加しながら元の値に戻らなければならず、矛盾するからだ。

ただし、この矢印は登れることを意味しない。崖でも地図上では接するし、接点まで飛び地内を移動する際に下りが必要かもしれない。粗いDEMでは隣接セルの標高差が大きく、層をいくつも飛び越す矢印もできる。コードではその差を記録し、間の土地を存在したことにはしない。

道路や橋による接続を加えるなら、別種の辺が必要になる。往復できる道を含む交通グラフは一般にDAGではない。低い層から高い層へ向けたのは今回の解析上の規則で、現実の移動を一方向に制限しているわけではない。

上向き接触グラフは、地形についての薄い骨格だ。そこへ移動条件を足して、初めて交通のモデルになる。

面積を測るときにも、何の面積かを決める

飛び地の広さは、セル数だけでは比較しづらい。緯度経度の格子では、同じ角度幅でも高緯度ほど東西方向が短くなる。

現在の実装では地球を半径 $R=6371.0088$ kmの球と近似し、一つのセルの面積を、

a(p)=R^2 \left|\sin\phi_N-\sin\phi_S\right| \left|\lambda_E-\lambda_W\right|

で計算する。 $\phi_N,\phi_S$ は北端と南端の緯度、 $\lambda_E,\lambda_W$ は東端と西端の経度で、いずれもラジアン。面積はkm²になる。球面の面積要素 $R^2\cos\phi\,d\phi\,d\lambda$ を、セルの範囲で積分した式だ。

飛び地 $C_{k,j}$ の面積は、そのセル面積の和になる。

A_{k,j}=\sum_{p\in C_{k,j}}a(p)

全飛び地の面積を足せば、有効セル全体の面積と一致するはずだ。

\sum_k\sum_j A_{k,j}=\sum_{p\in M}a(p)

これは分類の二重計上や取りこぼしを見つける検証に使える。ただし測っているのは球面上の面積で、起伏に沿った地表面積でも、耕作可能面積でもない。

同じ層の総面積が100km²でも、一つの100km²の飛び地なのか、1km²が百個なのかで景色は変わる。そこで次の集計案として、層内の最大飛び地が占める比率を置ける。

q_k=\frac{\max_j A_{k,j}}{\sum_j A_{k,j}}

占有面積が正の層だけで定義する。前者なら $q_k=1$ 、後者なら $q_k=0.01$ 。同じ総面積でもまとまり方を区別できる。これは本文で追加した指標案で、現在のコードの出力項目ではない。

この値もあくまで地形のまとまり方を測る指標で、道路の接続や土地利用は表していない。

関東周辺のデータで実行してみる

使ったのは、NOAA ETOPO 2022 Surfaceの60秒角データ。東経138〜141度、北緯34〜37度の矩形を、元のセルを補間せずに切り出した。行政上の関東地方そのものではなく、周辺も含む範囲になる。

項目	今回の条件
入力	ETOPO 2022 Surface、60秒角
切り出し範囲	東経138〜141度、北緯34〜37度
格子	180×180セル
対象セル	有限の標高値を持ち、標高0m以上の18,350セル
層幅・近傍	5m、4近傍と8近傍
リサンプリング	なし

解析スクリプトは入力ファイルのSHA-256を照合し、同じデータから結果を再現できるようにしてある。標高帯、連結成分、面積、接触辺を計算し、CSVと図へ出力した。実行条件と再現手順

関東周辺の0〜5m、5〜10m、10〜15m、20〜25mの標高帯を青で示した4枚の地図

青が各パネルで指定した標高帯、灰色がそれ以外の標高0m以上のセル、白が今回の集計から除外したセル。白をそのまま海と読むことはできない。経緯度で表示し、縦横比は北緯35.5度の距離比に合わせている。図は0〜5m、5〜10m、10〜15m、20〜25mの四つを抜粋したもので、15〜20mを含むほかの占有層も計算には入っている。

同じ範囲でも、取り出す標高帯が変わると、青い部分の広がり方が変わる。隣のパネルへ移るだけで、残る場所も消える場所もある。連続した地形を高さで切るという操作が、まず一枚の図として見えるようになった。

接続規則だけで、成分数は1,088変わった

全占有層を集計した結果はこうなった。

近傍	連結成分数	占有層数	対象面積 km²	接触辺数	1セル成分数
4近傍	14,563	511	50,970.97	29,663	13,345
8近傍	13,475	511	50,970.97	27,542	11,941

出典：実験のパラメータと集計結果。

対象セルも総面積も変えていない。斜めの接続を認めるだけで、成分数は1,088、4近傍を基準に約7.47％減った。小さな2×2の例で見た性質が、実データの集計でも現れている。

200m未満の各標高帯の4近傍・8近傍成分数と、4近傍で面積が大きい上位6成分の棒グラフ

層別CSVと面積ランキングCSVから再描画した図。左は下端が200m未満の層についての成分数、右は全占有層から選んだ4近傍の面積上位6成分。左のグラフだけで全511層を表しているわけではない。面積は球面セル面積の合計だ。

気になるのは、4近傍の成分のうち13,345、約91.64％が1セルだけだったことだ。8近傍でも11,941、約88.62％が1セルだった。この割合の分母は成分数であって、対象面積ではない。

つまり、今回の条件では細かく分かれた成分が数としては大半を占める。ただし、これを「関東には一万以上の独立した土地がある」とは読めない。約1.5〜1.9km間隔の格子を、鉛直方向には5mという細い帯で分類している。隣接セルが別の帯に入れば、その間に現実には同じ高さの細い通路があっても、この格子ではつながらない。

ここからは推測だが、多数の1セル成分には、実際の地形の分断だけでなく、水平解像度と層幅の組み合わせが強く影響している可能性がある。それを切り分けるには、同じ範囲を異なる水平解像度で解析して比べたい。

5mに切れることと、5mを見分けられること

地図を5m刻みに色分けする計算は簡単だ。でも、その色の境界がどれだけ信頼できるかは別の問題になる。

層幅の5mは鉛直方向の分類幅であり、水平のセル間隔でも、鉛直精度の保証でもない。今回使った60秒角のデータは、この緯度では概ね南北1.85km・東西1.51km間隔だ。それを5m帯に分けても、近所の階段や狭い道路が見えるようにはならない。

標高の推定値を $\hat h$ m、仮に誤差の絶対値が $\varepsilon$ m以下と分かっているとする。区間 $[\hat h-\varepsilon,\hat h+\varepsilon]$ が一つの層に完全に収まるセルなら、その誤差範囲で分類は変わらない。層境界をまたぐセルは、所属層が変わり得る。

さらに、分類の揺れるセルが二つの大きな領域をつなぐ細い通路にあると、その一セルの変化で一つの飛び地が二つになる。面積の誤差は小さくても、接続の誤差は大きくなり得る。

これは地面を観測した部分と補間した部分を分けて考えた話にもつながる。標高の絵だけでなく、どこで接続の判断が揺れるかも見たい。

もう一つ、標高0m以上を選ぶだけでは陸地を選んだことにならない。海面下の陸地を捨て、0mの海セルを残す可能性がある。初期実装のCLIにはこの標高フィルタがあるが、陸地限定の議論には独立した陸海マスクが必要になる。

同じ層に入ることと、移動できること

地図ができると、その色をそのまま移動可能な範囲として読みたくなる。でも、例えば15m、19m、16mと並ぶ道は、すべて同じ5m帯に入りながら、途中に3mの下りを含んでいる。

「一度も下らずに到達できるか」を調べるなら、標高帯の連結とは別の条件が必要になる。経路を $\gamma(t)$ 、経路上の位置を指定するパラメータを $t$ とすると、

h(\gamma(t_2))\ge h(\gamma(t_1)) \qquad (t_2\ge t_1)

が途中のすべての区間で成立しなければならない。これは追加の経路条件を置いた思考実験で、現在の連結成分抽出では検証していない。同じ経路をこの条件で往復できるなら、その経路の標高は一定である必要がある。

標高帯として連続する範囲と、条件付きで移動できる範囲は違う。道路や橋、傾斜、通行可能な幅まで知りたければ、そのためのデータを足す必要がある。

次は、切り方を変えて比べたい

最初の実行では、成分のセル数合計、独立に行別集計した球面面積との一致、接触辺の上向き性、8近傍の成分数が4近傍以下になることを検証した。実験スクリプト

これは実装内部の整合性の確認であって、DEMが現実の地面をどれだけ正確に表すかの検証とは別になる。次は、条件を一つずつ変えてみたい。

同じ範囲で層幅を5m・10m・20mに変え、成分の統合を追う。
水平解像度を変え、1セル成分の割合や大きな成分の形がどこまで残るかを見る。
独立した陸海マスクを入れ、標高だけで選んだ結果と比べる。
DEMの誤差情報をもとに標高を揺らし、消えたり現れたりする接続を調べる。
解析範囲を広げ、境界で切れた成分が領域外でつながるか確かめる。

同じ標高原点で5m帯を10m帯、20m帯へまとめ、対象セルと接続規則を固定するなら、既存の成分は分裂せず、残るか統合される。したがって全層を合計した成分数は増えない。この性質は次の計算を確かめる手掛かりにもなる。

全球に広げるには、日付変更線付近の接続、極、タイル境界の成分統合も必要になる。現在のラベリングはセル数 $N$ と占有層数 $K$ に対して $O(KN)$ の走査を行う地域プロトタイプなので、そのまま全球運用へ持っていく設計ではない。

世界地図を5mごとに切れるだろうか、という疑問には、まず関東周辺の一枚で答えられた。すると次に気になったのは、切れたように見える場所が、地形のどこで切れているのか、それとも地図の粗さによって切れて見えるのか、ということだった。

高さで地図を切ってみたら、今度はその地図の目の粗さが見えてきた。

Sources

MozyOk, 関東周辺を5m標高帯に分ける — 初回実データ実験、実行条件・入力SHA-256・集計値、再現スクリプト。標高帯の図と実データ集計の出典。統計図は同コミットのlayers-4.csv、layers-8.csv、components-4.csvから、再現スクリプトの描画条件に沿って再描画した。参照コミット：2e8f44384b6664bcc38a5a41d937d547e74849d1。確認日：2026-09-20 JST。
MozyOk, 数理モデル、解析コード、テスト。層・連結成分・球面面積・接触グラフの定義。確認日：2026-09-20 JST。
NOAA NCEI, ETOPO Global Relief Model、使用したETOPO 2022 Surface / 60秒角GeoTIFF。製品概要と実験入力。製品説明の確認日：2026-09-20 JST。入力の同定は上記実験のprovenanceに基づく。
SciPy, scipy.ndimage.label。連結成分と接続規則。確認日：2026-09-20 JST。
2×2格子と15→19→16mは説明用の仮想例。最大飛び地面積比率と誤差境界の議論は本文の検討案。成分数の減少率と1セル成分の割合は、公開された実験集計値から算出した。