※　本アーティクルは，北海道大学理学研究院在職時に行った検討を元にしたものです．

－ その１ 基礎編 －

DEM データとは？

　“DEM (Digital Elevation Model)”とはその名の通り，ある地域の標高をデジタルデータで表現したものです．そのデータ形式にはいろいろありますが，もっとも単純なのは XYZ データ，つまり XY 座標と標高という三つの変数で一つの地点の標高を表しています．通常は XY 座標というのは緯度経度ですが，もちろん他の相対座標とかでも構いません．下は，あるカンマ区切り DEM データの一部で，もっとも単純な経度緯度の XYZ 形式になっています（左参照）．

　これが何万何十万行も延々と続くわけです．ということで，要するに DEM データというのは，単に数値が並んでいるだけの実に無味乾燥なデータに過ぎません．人間がこのデータを見ても何も見えてきません．これをどうやって，地質・地形に関連させていけばよいのでしょうか？　また，そういうデータのメリットはどこにあるんでしょうか？

　それに対する答えとしては；　地形を単純なメッシュ標高データとしてデジタルに表現してしまえば，さまざまなコンピュータプログラム処理によって，その可視化や演算・変換が（スケーラブルに）可能になる，ということになります（右図参照）．

　なお，低精度の DEM は写真測量や既存地形図の等高線からデジタル化されていますが，近年，航空機搭載などのレーザ測量機器による高精度の DEM が一般的になりつつあります（LiDAR data: Light Detection And Ranging data）．

北海道の DEM データ

　DEM データと一口に言っても，地球上のいろいろな地域について様々な形式のものがあります．ここでは，国土地理院が公開している北海道の DEM データについて簡単に紹介したいと思います．

①　国土数値情報（1992年頃）
　・３次メッシュ（1 km 間隔）ごとの数値化情報．データサイズは約 30 MB．
　・標高は３次メッシュを16分割＝ 250 mメッシュ DEM（精度：10 m）．
　・北海道全域というスケールでは，十分な表現力があった．
②　基盤地図情報（2009～2010 年頃からダウンロード提供開始）
　・5 m メッシュデータ：　標高データ精度^（下注）は 1 cm．全域のデータはない（下図）．国土地理院整備の『公共測量成果』と，国土交通省提供の『基本測量成果』がある．後者のカバーエリアはきわめて限定的．データ取得方法は，写真測量とレーザ測量．　　注）：データの表現桁を示すもので，誤差ということではない．

　・10 m メッシュデータ：　標高データ精度は 1 m．日本全域のデータがある．データサイズは，北海道のみで約 11 GB．データ取得方法は，おもに 1/25,000 地形図の等高線から．

データ形式：　国土地理院等で配布されている DEM データは XML (GML) 形式になっています．これは，データ前後に階層化されたメタデータがくっついており，その中でデータ両隅の緯度経度やデータの行と列の数などが与えられているものです．下のデータ例を参照してください．

<gml:lowerCorner>42.525 140.775</gml:lowerCorner>
<gml:upperCorner>42.533333333 140.7875</gml:upperCorner>
　　　　(snip)
<gml:low>0 0</gml:low>
<gml:high>224 149</gml:high>

　データ自身はカンマ区切りで右のようになっており，XYZ 形式のデータにはなっていません．最初にあるのはそのデータの属性で，２番目が標高（この場合，単位は m）です．

　座標がどこにも見当たりませんが，上に示したメタデータ中に記述されているデータの縦横グリッド数や両隅の座標から incremental に各データの座標が得られるようになっています．

DEM データの可視化

　DEM データを可視化する方法として，一般に以下のような三つが考えられます．

　　①　XY 座標点の“色（輝度）”として標高を表わす．単なる平面図．
　　②　陰影処理を施し地形が立体的に見える画像として表わす．
　　③　地形を３Ｄポリゴンとして表わす．視点などを変えて見ることができる．

　おそらくこれ以外にもあるんでしょうけど，私にはこれしか思いつきません．①はまあ私のような素人のプログラム処理能力でもなんとかなります（次項参照）が，②③は到底無理です．
　それではどうするのか？というと．．．市販ソフトの力を借ります．私が使ったのは，この筋では有名どころの，Golden Software という大げさな名前（失礼！）のメーカの Surfer というソフトでした．（⇒付記参照）この使用例については別項で述べます．
　結構高いソフトなので，気軽に買えるようなものではありません．最新のバージョン 16 の single-user license が $849 となっています．国内で買えば，10万円コースですね．．．有名な HULINKS で扱っているようですが，\150,120 （税込み）でした．た，高い．これは無収入者の私にはもう買えません．これに相当するフリーソフトはずっと探しているのですが，まだ見つかっていません．なにか良いのないですかね．．．

※※付記：　元山形大学の川辺孝幸さんからご教示いただいたのですが．．．この Golden Software という会社名，“金ぴかソフトウェア”ではありませんでした．コロラド州の Golden という市で設立された会社だったんですね．ゴールドラッシュでできた町のようですが，Golden はそういう由来ではなくて，ここに入植した金鉱夫の名前なんだとか（本名ですかね？）．で，Golden にはあの有名なビールメーカー Coors があるんだとか．驚きました．自然に囲まれた素晴らしい町です．ここにはそういう経緯でコロラド州立鉱山大学があるんですが，その学生さん二人によって設立されたソフトウェア会社です．なお，最新版の Surfer はバージョン番号表記がないのですがさらに値上がりしており，Single User license が $999 でした．HULINKS だと 20万円コースでしょうか？　旧バージョンを giveaway してくれないかな．．．？　あと Overview のページの一番最後に用途分野例が書かれているのですが，『地質・地形』にもろ相当するものがないですね．．．Education かな？　ちょっと残念．（2021/07/25）

データ抽出・変換・表示プログラム（XML converter）の開発

　基盤地図情報の XML データを単純な XYZ（緯度・経度・標高）データに変換するプログラムを自作しました．なぜそういうことが必要かというと，上述の Surfer は XML (GML) 形式のデータを直接読むことができず，XYZ 形式のデータしか処理できないからです．もう一つ重要な理由として，XYZ 形式のデータは，行ごとにそれ自身完結しているので，データの連結や一部抽出などが面倒な追加処理なしに可能，ということがあげられます．

　処理対象データは GML 形式の 5 m, 10 m メッシュデータ，プログラム開発環境は Microsoft VisualBasic 6 です．古くて，しかも non-professional な開発環境ですが，Windows10 でもまだ動くし，過不足はありません．私が C 系言語にまったく経験を持っていないプログラミング素人である，ということもありますが．
　なお，VisualBasic 6 とそれによって生成された実行ファイルが現在のＯＳ環境で動作しなくなった時の“保険”として，Microsoft Visual Studio のコンポーネントの一つである VisualBasic 2008 以降に対応したバージョンも念のため作成しています．

　XML コンバータの基本機能は，北海道の任意の区域を緯度経度レンジにより指定し，メッシュデータを XYZ 形式の DAT データとして書き出し，同時にそのデータを標高に応じたグレースケール輝度によりグラフィック表示するものです．このプログラムの主要なアルゴリスムは，非常に簡単なもので，以下の通りです．

　　①　GML ファイルを１行ずつ読み込む．
　　②　そこからデータの左上隅と右下隅の緯度経度（上述）を得る．
　　③　同じく，データの行列数（上述）を得る．
　　④　②と③から緯度経度の行と列あたりの増分を得る．
　　⑤　データから標高を読み込み，その緯度経度を②④で increment する．

　この処理の基本ループを示すと下のようになります．もちろんこれはループのイメージを表現したもので，動作可能なコードにはなってません．

Delta_Row = (Longitude_max - Longitude_min) / Row_num
Delta_Line = (Latitude_max - Latitude_min) / Line_num

Do While Not EOF(1)
　　Line Input #1, LineString
　　Altitude = Val(AltitudeString)
　　Longitude = Longitude + Delta_Row
　　I = I + 1
　　If I >= Row_num Then
　　　　Latitude = Latitude - Delta_Line
　　　　Longitude = Longitude_min
　　　　I = 0
　　End If
　　Write #2, Longitude, Latitude, Altitude
Loop

　これによって得られたデータを XYZ 形式のデータファイルに書き込み，同時に標高を適当にグレースケール輝度（0 以上 255 以下）に変換し，グラフィック表示しています．グラフィック表示は，まあおまけですが，データ変換がうまく行っているかどうかを視覚的に確認できるという実務的なメリットがあります．

　こんな簡単なループを基本とするプログラムですが，あの膨大な数値データがちゃんと北海道の地形表現になっていくのを画面で見ていると，やはり感動しますね．蛇足ですが，それがコンピュータを使ったデータ処理の醍醐味かな？と感じます．

GIS じゃダメなの？！

　さて，ここまで読んできた皆さんの中には，『なんでわざわざそんな？　GIS じゃダメなんですか？』と思っている方がいらっしゃるのではないでしょうか？　よく分かります．なんで GIS じゃダメなんでしょう？

　言わずもがなですが，世の中には『GIS (Geographic Information System)』というプロフェッショナルな仕組み・枠組み（≒ソフトウェア）が既に確固としてあり，DEM 情報を含む各種の地理情報をハンドリングできます．下にあげたのはその一例で，QGIS というオープンソースのフリーウェアで，DEM データをハンドリングしたものです．

このような GIS のメリット・デメリットとしては，下のようになると思います．

メリット：　業界標準・高機能・簡単・ready-to-use．．．
デメリット：　高価^（下注）（QGIS のようなフリーのアプリケーションもあるが）・融通きかない．．．全然．

　　注）：おそらく業界標準と思われる ESRI 社の ArcGIS の基本版 ArcGIS Desktop Basic 単独使用ライセンス が 390,000 円！　桁が違う．

　“融通がきかない”と言うのは，もちろん私の主観的な判断です．上に示したように地形傾斜表示のような，結構複雑な計算処理を必要とするものもメニューから簡単に（計算式など知らなくても）結果をゲットできるわけですから，地形図なんかもオーバーレイできるし，しかも無料のソフトもあるわけなので，文句言ったらバチが当たりますよね．
　しかし．．．私は GIS を使いません．理由は，そうですね．．．
　　①　ここぞという痒いところに手の届くようなピンポイントデータ処理ができない．
　　②　上の図のキャプションにあるような問題点が発生した場合でも，サプライヤー任せでユーザはどうにも手が出せない．

　　　　　　注）もちろん QGIS のようなオープンソースソフトウェアの場合は，『自分でソースいじれるだろ？』ということにもなるわけですが，理屈と現実は厳然と違いますから．．．

　他にもいろいろありますけど，結局は『自分で処理してハマった時の快感』ということになるでしょうか．

　なお，②のもう一つの例を下にあげておきます．QGIS の最新版 3.6 で陰影処理を作成表示した例です．平坦部の陰影処理が全然おかしいのですが，なにか私が作成パラメータをミスっているのかもしれませんけど，現時点ではどうにもならない状態です．

Surfer によるDEM データの３Ｄ表示

　Surfer (Golden Software) は，ソフトの名前（surface にかけたんでしょうけど）もメーカの名前もなんだかいまいちだけど，ソフト自体は素晴らしいものです．私はその機能の全部は使っていませんでしたが，それでも『これに代わるソフトは存在しない』と思わせるに十分なものでした．

　Surfer は要するに，XYZ 形式の DEM データを視覚化するソフトです．ただ，XYZ データを直接扱うことはできず，一度メッシュ補間を行いグリッドデータを作成する必要があります．もともと国土地理院の DEM データはグリッドデータなんだから，なんであらためてメッシュ補間を行う必要があるのか？なんて思うのですが，要するに “四角四面の” グリッドデータが必要ということなんでしょうね．確証はないのですが．
　で，グリッドデータになってしまうと，あとは陰影図（Shaded Relief Map）や３Ｄ俯瞰図（3D Surface Map）が簡単に出てきます．下はその一例です．このほかにコンターマップや３Ｄワイアフレーム図を作る機能もありますが，私はほとんど使っていませんでした．

　Surfer の視覚化機能以外に重要な機能として，グリッド演算があります．これによって例えば，測定時期の違う二つの DEM データの差分（＝標高・地形の変化）を求めることが簡単にできます．これについては“その２　応用編”で紹介したいと思います．
　あと Surfer の重要な特徴として，“データ量が莫大になっても処理速度がそれほど落ちない”という点があります．使ってきた経験内では，データ量の増大でメモリ不足になって異常終了とかいうことはまったくありませんでした．そのへんの安定性も大きな長所ではないかと思います．

国土地理院 5 m メッシュデータのパワー

　先に，国土地理院の基盤地図情報として，5 m と 10 m メッシュデータがあるということを述べました．この二つのデータ，メッシュ間隔としてはたったの２倍（＝データ数は４倍）に過ぎないのですが，得られる結果としては非常に大きな差異があります．
　下にあげたのは，千歳市東部のデータですが，特に平野地域の表現力に大きな違いがあります．5 m メッシュデータでは，道路や農地境界，川筋や堤防などがはっきりと確認できます．一方，10 m メッシュデータではそういった部分は解像度がほとんど無く，処理上の artifact と思われる不自然な陰影も認められます．
　この違いは，メッシュ間隔が違うことに加えて，10 m メッシュデータは 1/25,000 地形図のコンターから補間的に求められているのに対して，5 m メッシュデータはレーザ測量による LiDAR データであるということが大きく効いているのだと思われます．

　上の 5 m メッシュデータによる陰影の右上隅の部分を拡大したのが下図の左です．データ欠損領域が大きいのが玉に瑕ですが，中央部付近に川（嶮淵川）を挟んで NNW-SSE 方向に鋭く見えている直線状地形（リニアメント）が，昨年９月の胆振東部地震で有名になった“石狩低地東縁断層帯”の構成要素の一つである『泉郷断層』の断層（撓曲）崖です．これを例えば Google Earth で見ても，リニアメントはまったく分かりません（下図右）．もちろんこのリニアメントは 10 m メッシュデータでも見えてはいるわけですが，拡大してもここまでリアルには見えません．これが 5 m メッシュデータのパワーだと思います．

　ちなみに，この泉郷断層のリニアメント＝断層・撓曲地形は，実際に現地に行って見てみると下の写真のようなことになっています．畑の右側に見えている低い段差状部がそれです．このような微妙な地形でも 5 m メッシュデータのパワーは明瞭に描き出すことができるんだということになるでしょう．

（2019/04/20）

－ その２ 応用編 －

地質構造と地形

　DEM の視覚化で見えてくるもっともスぺクタクルなものはやはり，地質体の構造が形作るマクロ地形でしょう．
　下に示したのは，産総研の地質図 Navi による日高山脈南部の西麓部の地質概略図です．

　10 m メッシュ DEMデータから作成したこの部分の３Ｄ俯瞰図（下図）を見ると，まさに（私の口癖ですが）『地形と地質は 1：1 に対応している』ことが如実に示されています．
　東側の標高が高く起伏の険しい部分が，空知－エゾ帯の付加変成体であるイドンナップ・神居古潭帯に相当しています．両者の間に白亜紀前弧海盆堆積体の蝦夷層群が分布する蝦夷東帯がありますが，幅が非常に狭いため地形特徴はあまり明瞭ではありません．
　西側の低標高・低起伏の部分は新第三系の分布地域です．その南部，海岸線とほぼ平行な構造が明瞭に見える部分は，新第三系中の顕著な褶曲・断層帯で，三石断層帯（青矢印）には下位地質体である蛇紋岩・高圧変成岩体が狭長に露出しています．その特徴的な地形が見えています（下図・下）．

高精細 DEM から見えるもの　－変動地形

　石狩低地東縁断層帯（活断層研究会，1991）は，石狩平野東縁部を南北に走る活断層帯で，延長 50 km 以上に達するものです．

　高精細（5 m メッシュ）DEM から作成した３Ｄ陰影図（下図）には，これらが明瞭に表現されています．特に，北縁部の岩見沢断層の作る撓曲地形と思われる部分（赤矢印）が目を引きます．この部分とそのすぐ西側には国土地理院の活断層図で二（～三）本の活断層線が引かれていますが，地形図や空中写真ではほとんど認識不可能なものです．

　いまさらながら，右のような地形陰影図を見ると，高精細 DEM データの“パワー”をひしひしと感じてしまいます．私が専ら野外地質調査に従事していた 1970-1990 年代には，少なくともこういう情報は身の回りに存在していませんでした．紙に印刷された 1/25,000 地形図（調査地域の書店で購入していた！）がもっとも高精度な『地形情報』だったわけで．．．隔世の感があります．

高精細 DEM から見えるもの　－正体不明地形

　以下の話は私自身真面目な確信が持てないもので，いわば“冗談半分（７分３分？）”です．そんなこと書いていいのか？と言われそうですが，可能性がゼロでないうちは書きたくなるのが地質屋の悲しい性（さが）と言うか，先の短い年寄りの特権と言うか．．．

　『泉郷_{（いずみさと）}断層』は，上の項で述べた石狩低地東縁断層帯の南縁部における西限断層です．しかし，泉郷断層からさらに西に約 5 km 離れた千歳市祝梅_{（しゅくばい）}付近の低地上に，その断層帯にほぼ平行な NNW-SSE 方向の不明瞭リニアメントが認められます（下図左）．この起伏の高さは最大で 10 m 程度です．その正体は不明ですが，開析された断層崖あるいは撓曲地形の可能性があるのではないでしょうか？　平川ほか（2008）はこの近傍に活断層の存在を推定しており，国土地理院都市圏活断層図でも，それが図示されています（下図右）．

　この祝梅周辺の現地確認を（野次馬根性で？！）行ってみました（2013/11）．リニアメント南部の最大起伏点付近では比高 10 m 以下の緩やかな丘状地形が確認されましたが，地形改変もあり，（半）素人が地上で見た限りでは，その全体的な形状や由来はまったく分かりませんでした．
　なおリニアメント北部では，リニアメント方向に対応した直線状の小起伏も視認されました．その正体はまったく不明で，あまりに小規模なので人工地形の可能性も高く，単に畑を拡幅した時の名残なのかもしれませんが，周辺に推定されている同方向の活断層・活撓曲・活褶曲（上図参照）との関係が興味深いところです．

高精細 DEM から見えるもの　－短期地形変動

　国土地理院で公開されている高精細 DEM データは，一般に“ある時取得された”もので，それが一定期間を置いて更新されるということはめったにないと思います．．．というか，そういう例を私は知りません．ある意味，固定データというか．しかし，なにか特定の目的（研究・防災．．．etc.）をもって取得された DEM データはその例外でしょう．

　その良い例が，『2000年有珠山噴火に伴って取得された DEM データ』（仲野ほか，2001：砂防学会誌）だと思います．この時有珠山は 3/27 頃から噴火の可能性が予測され，3/31 13:07 に噴火が始まりました（Wikipedia による）．
　下の画像は，右が 3/31 の噴火直前に取得された DEM データから作成したもので，左が約１か月後（4/26）に取得されたデータです（土木研究所提供）．3/31 のフライトが何時ころだったのかは分かりませんが，まさに“紙一重”のデータで，このような例はなかなかないと思います．ちなみにデータ（グリッド）数は約 580 万点です．

　この噴火後の DEM 画像には，当たり前ですが，噴火で形成されたいくつかの火口群が見えています．よく見ると西山火口群付近の国道 230 号線周辺には，地形変動によって形成された WNW-ESE 方向の小断層群が発達しているのも分かります．しかしそれ以上のことは画像をじっくり見比べてみてもあまりよくわかりません（見えません）．
　そこで，DEM データがデジタルデータであることを利用して，両者の差分を求めてみました．方法としては，Surfer のグリッド演算機能を使用しました．得られた差分グリッドデータを可視化して見ると下のようになります．

　このような“変動量図”は，仲野ほか（2001）で詳細に示されており，全然私のオリジナルな発想・処理の結果ではありませんが，まあ『自分でも（自前で）やってみたかった』という自己満足ということになるでしょうか．

　で，この図から明らかなことは，西山火口群周辺の顕著な隆起域の存在です．その最大隆起量は 71.6 m です．上昇がこの期間中ずっとリニアに起こっていたと仮定すると，ほぼ 2.7 m / day の日上昇量となります．マグマの供給は 2000/08 頃まで続いていた（Wikipedia）らしいので，この上昇量がそこまで継続したとすれば隆起量はさらに 300 m を越えてしまいますが，もちろんそんなことは実際には起きていません．
　ちなみにデータ中の最大沈降量は -26.7 m ですが，これは金毘羅火口（群）の形成によるもので，“沈降”と言うのは適当ではありませんね．“損失”というべきなのかも．西山火口でも当然同じような“損失”があるはずですが，こちらは隆起によって相殺されています．この隆起の原因は，言うまでもありませんが，マグマ上昇によって潜在円頂丘（cryptodome）が形成されたということでしょう．
　この図を見ていてもう一つ気づくのは，“NNW-SSE 方向に走る線状パターン”（矢印）です．これについては特に誰かによって言及されているわけでもなく，その正体はよく分かりません．データ処理の過程で入った artifact である可能性もゼロではありません．単純に『変動量がある線を境界として不連続になっている』と考えると当然，断層の可能性が考えられますが，実際にこのような断層地形が地表で認められたというわけでもないので，謎です．なお，岡田（2001）はこのパターンに一部相当する部分を，“洞爺湖温泉西端のカルデラ壁に対応したヒンジ型の食い違い変動”と表現していますが，その詳細は不明です．

　なお，西山北麓部には沢地形に対応した数 m 程度の標高減少部分がありますが，これは 3/31 時点で沢型中に残存していた積雪が４月の融雪により消滅したものと考えられます．また，同じ程度の標高減少部分が温泉街西方に広く認められますが，この意味はまったく不明です．

　これらの変動量を頻度分布図にしてみました（下図）．600 万件にも達するデータを Excel で一度に集計するのは普通無理なので，自前の処理プログラムを使って集計し頻度を求めています．当然ですが，変動量ゼロ付近のデータがもっとも多いのですが，隆起を示すプラス変動値もかなりの割合を占めています．ちなみに隆起を示すデータ部分では，隣り合う階級の頻度比がほぼ 0.8 - 0.9 と一定しており，べき乗分布を示すと判断されます．この意味は残念ながら私には分かりません．特定の変動量にピークを持たないということにはなるんでしょうけど．

　ここまでくると当然，次のような興味が浮かんできます．『その後どうなったのかな？』　これについては，おそらく火山・地形屋さんがよってたかって GPS データかなんかでやっているはずなので，私みたいな素人の出る幕ではないのですが，たいして意味なくても思いつくとやらないでは済まない私の性格が．．．そう，国土地理院基盤地図情報と比較してみたら？という発想です．
　そこで，国土地理院 5 m DEM（2009/10/28測量）と上の有珠山噴火直後（2000/04/26）の標高データとの差分を求めてみました．下図は二つの DEM データから作成したそれぞれの陰影図ですが，さすがにと言うべきなのか，両者の間には“精細さ”にちょっと違いがあります．これが実際のデータ精度の違いなのか，作成時のパラメータ設定とかの問題なのかは，よく分かりません．

　実はお分かりのように，この二つの DEM データは共通の XY 座標系を持っていません．位置が合っていないと差分も何もないわけで，そこは非常に苦労した点です．単純な計算だけでは不可能で，手作業含めて多大な（部分的にアバウトな）変換作業が必要でした．正直言って，検討に値するほどの正確さが実現できたかどうかは自信がないのですが，結果はそれなりに reasonable なものになりましたので，そこそこの座標変換ができたのではないかと．
　座標変換さえできれば，あとは上の差分を求める手法をそのまま使って，約９年半の間の地表変動が見えてくる．．．はずです（下図）．

　人為的改変^下注）を含めて状況は複雑ですが，西山火口周辺に弱上昇域が確認されました．その上昇量は 4 ～ 8 m で，噴火時の隆起域と形態的にほぼオーバーラップしています．なお右下のブランク部は，データ欠損により比較の不可能な領域です．

　この結果の学問的な妥当性・意義は（素人作業なので）明瞭ではありませんが；
　　可能性１：　潜在円頂丘の隆起が 2000/04/26 直後にこのレベルで停止した．
　　可能性２：　潜在円頂丘の速い隆起は 2000/04/26 以前に停止したが，その後も非常に緩やかな隆起がずっと継続していた．
　　可能性３：　潜在円頂丘の隆起が 2000/04/26 以降もしばらく進行しピークに達した後に沈降したが 2000/04/26 のレベルよりやや高いところで停止した．

．．．といったいくつかの可能性が想定されるのではないかと思われます．

注）標高減少部は，おもに褐色（-10 ～ -15 m）～濃オレンジ色（- 15 ～ -20 m 以上）で示されている．その想定要因は図中にキャプションで示した．温泉街周辺の減少域は，切土造成地と撤去された建造物である．西山火口西方の二つの矩形の減少部は道道（？）の整備による切土部，その間の増加部は谷埋めの盛り土部である．画像左下隅の減少部は噴火後に稼行した土採り場であろう．金毘羅火口とその対面にある火口は標高減少が著しく，金毘羅火口は内部に見かけ上の上昇部を持つが，その詳細は不明である．

（2019/04/27：2019/05/01 加筆）

－ その３ 解析編 －

はじめに：

“解析編”と言うのは実に大げさなタイトルで，以下の内容はそのようなたいしたものではありません．“応用編”が DEM データを地形陰影図に変換したものを主体にしていたので，それ以外のものをまとめてみたということになります．

DEM 派生データ　－標高分布

　標高データのセットが手の中にあってそれを可視化した後，なんとなく私の頭に浮かんだのは『標高の分布ってどうなってる？』でした．右にあげた図は，実に古い話ですが，川村（1993）（山岸ほか『北海道地すべり地形データベース』）に掲載した北海道の標高とその頻度図です．元データは国土地理院国土数値情報の３次（= 250 m）メッシュデータで，データ数は 1,340,246 = 約 134 万件です．
　ちなみにこの数値情報データですが，1992 年ころ（？）国土地理院から巨大な磁気テープのセットとして供給されたもので，私のＰＣ環境にはもちろんそれを読めるディバイスなど無く，当時同じ講座の助教授だったＦ先生に大学の計算機センターで変換してもらった記憶があります．たしか全部で 200 MB くらいのデータだったと思うので，そんな大げさなと思うんですが，そういう時代だったんでしょうね．

　まあこれでも必要十分な情報が得られているわけで，標高のべき乗的分布や，その中にあるいくつかのスパイク・凸凹などがすぐに見て取れます．

　現在は 10 m メッシュデータとして北海道の精細な標高データが配布されているわけなので，あらためてそこから標高分布を求めてみるとどうなるでしょうか？　そのデータ数は，単純に考えて 1340246 x (250 / 10) ^ 2 = 837653750 = 約８億件となります．まさに“ビッグデータ”と言うにふさわしいでしょう．このようなビッグデータの処理は EXCEL のような on-memory の表計算ソフトでは不可能で，大量のデータファイルを逐次行処理するプログラムを作成して初めて可能になります．しかしその処理はというと，実に簡単な一つのループ処理で “小学生でも書ける” というか．．．

Do While Not EOF(1)
　　Input #1, AltitudeString
　　Altitude = Val(AltitudeString)
　　If Altitude >= 0 Then
　　　　Elev_Freq(Int(Altitude / 5)) = Elev_Freq(Int(Altitude / 5)) + 1
　　End If
Loop

※例によってこの VB コードは単なるイメージで動作可能なコードではありません．集計される頻度階級の幅は 5 m．

　要するにこのループが８億回繰り返されるだけなわけです．いつも思うのですが，ここがコンピュータ処理の良いところと言うか，つまらないところと言うか．．．要するにコンピュータは人間と違って，いくらループ回してもスピードも正確さもぜんぜん変わらないわけですね．疲れないし飽きないし．その結果をグラフにすると下のようになります．

　上図を見るとお分かりのように，標高分布は，ほぼ『べき乗分布（power-law distribution）』に近い形態を示しますが；
　　　　①　400 m 以下の部分に何本かのスパイクがある．
　　　　②　200 m 以下の部分で“段付き”が認められる．
　①はいくつかの湖の湖面標高によるスパイクと考えられます．参考としていくつかの湖の湖面標高は；　洞爺湖：84 m，屈斜路湖：121 m，支笏湖：247 m，朱鞠内湖：282 m，摩周湖：355 m，阿寒湖：420 m となっています（by Wikipedia）．これらは上図のスパイクの位置とよく対応しています．
　②は，詳細にみると標高 40 m, 130 m の部分に存在し，段丘面（・崖）の発達による影響の可能性が考えられますが，詳細は不明です．

　この標高頻度分布を頻度軸を対数にした片対数グラフで見ると，標高 2,000 m 以下の部分ではほぼ直線となり，べき乗分布をしているものと考えられます．
　隣の階級との頻度比を求めてみると，やはり標高 2,000 m 以下の部分でほぼ一定（≒ 0.8 ～ 0.98）となっています（下図）．これはこの分布がべき乗分布であることの一つの根拠だと思われます．

　じゃあ，標高分布が基本的にべき乗分布をしているとして，それはどういう意味を持っているのでしょうか？　私にはスケール則と地形形成の営力との関係とかよく分からないのですが，標高（＝地形）の形成は『スケールに依存しない何らかの法則』に支配されているということになるんでしょうか？　それ以上はまったくの no idea です．いずれにせよ私如きが首を突っ込めるような領域ではないと思いますので，事実だけ書いて，ここまでにしておきます．

DEM 派生データ　－傾斜量

　次に誰しもの頭に浮かぶのはやはり『地形の傾斜はどうなってる？』でしょう．左にあげた図は，前項で述べましたが 1992 年ころに入手した国土地理院国土数値情報の３次メッシュデータから作成した北海道の傾斜図です．

　このデータにはメッシュノードの傾斜値が最初から入っているため，なんの計算も必要なく，この図はそれを単にメッシュ上に表示しているだけです．言うまでもなく 250 m 間隔の地点の傾斜が色分けされて表示されているだけですから，分かるのは北海道全体のマクロ的な傾斜の特徴だけで，それ以上のことはなにも分かりません．まあそれでも，日高山脈がとびぬけて傾斜の急な地域であることは良く分かりますが．
　なおこの３次メッシュ国土数値情報のドキュメントには，この傾斜がどのように計算（・取得）されたものかは，記述されていなかったと記憶しています．

　そこで，北海道全体で８億点もある 10 m メッシュデータから傾斜を求めることを考えました．その方法ですが，もちろんある一点の標高から傾斜が分かるわけはないので．．．とりあえず独自に次のような方法を考えました．これを私は勝手に『 Max - Min 法』と呼んでいます．
　まず，ある標高点を中心とした 3 x 3 = 9 点の標高を取ります．その中から最大値と最小値を求め，その差を標高点の水平距離で割った値の逆 tan を傾斜とします．数式風に書くと，SlopeAngle = Atn ((Altitude_max - Altitude_min) / Separation) となります．処理の詳細は省略しますが，この方法では傾斜方向も（離散値になりますが）得ることができます．

　と・こ・ろ・が．．．こんなことをしなくても，傾斜とその方向を得る方法が既に世の中にあることが，その後分かりました．それは，傾斜量図と地形・地質との関係をおそらく日本では初めて指摘した神谷ほか（1999, 2000：情報地質）に紹介されていますが，右図のような方法です．

　ちゃんとした方法がとっくに世の中に示されていた（≒私は無駄骨）ということで，非常にがっかりしました．悔しいので，私の独自方法の結果との比較をやってみました（下図）．

　やはりと言いますか，私の Max - Min 法は，“切れ味”が鈍いようで，遠目に見るとちょっとコントラストが低いかな？という程度ですが，等倍で見ると明らかにデータの分離度が低いです．まあでも，それなりに妥当な結果を得ていたんだな，と少し満足しました．こんなものでしょう．

　ところで，こうやって傾斜量を可視化して見ると，左の図のように “脳状模様” というかそんな風な独特なパターンに見えるのが特徴的です．この表現は神谷ほか（1999, 2000）で使われており，尾根線の発達する山地地域における傾斜量図の特徴とされています．

　さて次に，傾斜量図と地形・地質との切っても切れない？関係について簡単に触れておきます．

　右の図は，北見地域のある区域の 10 m メッシュ標高データを Surfer で地形陰影図にしたものです．地形特徴が見事に見えていますが，地表面の傾斜分布を直感的に把握することはこの図からは（不可能ではありませんが）案外難しいです．

　この地域の地形図を国土地理院の電子国土 Web から取得してみると，下の図のようになっています．単なる等高線図ですが，実は我々『地図判読スキル保有者』は，この地図から傾斜量を直感的に読み取ることができます．

　傾斜量というのは要するに標高の微分ですが，等高線図における傾斜 = tan^-1 (等高線の垂直間隔／水平間隔) なので，言うまでもありませんが垂直間隔は一定ですから，水平間隔が小さいほど傾斜は大きくなります．そして我々の眼は，等高線の水平間隔を等高線を遠目に見た時の“色の濃さ”として感じ取っているというわけです．

　つまり，左の図で中央部から下部にかけての部分が傾斜量が大きいことは，我々には『一目で見て分かる』わけです．等高線図というのは偉大な地形表現システムだなとつくづく思います．しかし，だからと言って傾斜量図が不要というわけではありません．なぜかと言うと，“計算して求めればデジタルな定量的データとなる” からです．

　上の 10 m メッシュデータから作成した傾斜量図を下に示します．ちなみに現在では，このような傾斜量図も国土地理院電子国土Web から取得できるようになっています．そうなんですか．．．また無駄骨？　まあいいか．

　この傾斜量図で分かることは，等高線図でもわかりますけど，中央部から南部にかけての傾斜の急な領域があり，それは北部と西部に明瞭な境界をもって緩傾斜部と接しているということです．この理由は，地質図を見ると一目瞭然です（下図）．

　急傾斜部は，新第三系達媚_{（たつこぶ）}層の分布（オレンジ・紫・濃青色部）とほぼ 1:1 に対応しています．達媚層は硬質頁岩を主体とした地層で非常に硬質な岩相を持っています．急傾斜部の北縁はゆるくカーブした境界を持っていますが，これは WSW-ENE 方向の栄森断層に一致しています．栄森断層の北側には砂岩礫岩主体の古第三系栄森層（青色部）が分布しています．堆積年代は古いのですが，比較的軟らかい地層のようです．また急傾斜部の西側には白亜紀付加体仁頃層群（濃い茶色部）が分布していますが，新第三系とは NS 方向の二又断層で接しています．この地質境界も地形急傾斜部の境界として明瞭に表れています．
　新第三系中には NNW-SSE 方向の軸を持つ褶曲構造が発達しています．中央にタツコブ背斜，その西にタツコブ向斜があります．傾斜量図をよく見ると，この褶曲構造に対応するパターンが見えています．

　なおここでは詳細は示しませんが，地形急傾斜部の中には，緩傾斜を示す色の明るい不規則な楕円形の部分が，最大で数百 m の規模でいくつか散在しています．これらは，地すべり地形の分布・形状に対応しています．実はここで図示した部分を含む北見南方山地は，北海道でも有数の地すべり地密集地域となっている区域です（川村，1993）．

DEM 派生データ　－起伏量

　これは正直言って蛇足です．まあ一応やってみたことなので．．．それは地形の『起伏量（ruggedness）』です．これがどのようなパラメータなのか具体的にはどうもイメージできないのですが，私は勝手に“ ≒ 傾斜量の微分”つまり“標高の２回微分”と考えています．一般には『地形の粗さ』（凸凹度？）を示すパラメータとされているようです．その計算式（Riley et al., 1999）は下の通りです．

　なんと言うか．．．こうやって可視化してみても，傾斜量図とほとんど同じパターンが表示されるだけで，私にはその意義を理解することはいまいちできませんでした．処理がどこかで間違っているのかな？
　そのほかに，ラプラシアン（Laplacian）というパラメータもあります（岩橋・神谷，1995：情報地質）．ある点の標高と，その周辺４点の標高の和から計算され（４近傍法 4-neighbors Method），標高の２回微分に相当し，地形の凹凸度を表現するとされるものですが．．．これも一応計算して可視化してみたのですが，起伏量以上になんだかわけのわからないものでした．ここでは省略します．

（2019/05/01）（2022/12/26ページ統合）

mak-kawa ウェブの目次へ戻る