エクセルで作る地形解析図（4）
４．作図方法4.1 表形式データの作成 843750行×1列の表を750行×1125列にする表変換マクロはこれだけ。この例ではB48〜B843797セルにあるデータを、F2〜AQL751セルに出力させるようにしている。 マクロが入ったエクセルファイルを開いておく。ダウンロードしたGML形式のxmlファイルの拡張子をcsvに変え、ファイルを開いて表示した状態で、『alt＋F8』キーからこのマクロを実行させると、csvファイルのその画面（シート）に数列が変換される。変換された表を、データファイルにコピーすればよい（2022年現在、csvは文字化けしてそのまま開けない）。 Dim i, j Range("F2:AQL751").Clear For i = 1 To 750 For j = 1 To 1125 Cells(i + 1, j + 5) = Range("B48:B843797").Cells(j + (i - 1) * 1125, 1) Next Next MsgBox "変換終了", vbInformation End Sub 上記マクロでは変換に1分半程度要していたが、以下に変更したら数秒で終わるようになった。 Sub 表変換() ' ' GML形式xmlファイルにて直接変換 ' Dim i, Ri, Ci As Integer Ri = 750 Ci = 1125 Range("F2:AQL751").Clear Application.ScreenUpdating = False For i = 1 To Ri Range(Cells(Ci * (i - 1) + 48, 2), Cells(Ci * i + 47, 2)).Copy Cells(i + 1, 6).PasteSpecial Paste:=xlPasteValues, Operation:=xlNone, SkipBlanks _ :=False, Transpose:=True Next Application.ScreenUpdating = True MsgBox "変換終了", vbInformation End Sub 4.2 作図計算 （1）傾斜量図 データは”10mDEM”シートのA列2行目以降にあるものとする。計算シートのB3セルに以下の計算式を入れ、あとはデータの数だけコピーする。 =SQRT(((SUM(10mDEM!A2:A4)-SUM(10mDEM!C2:C4))/(6*ピクセル幅))^2+((SUM(10mDEM!A2:C2)-SUM(10mDEM!A4:C4))/(6*ピクセル高さ))^2) ということでもいいのだが、実は、このコピペは非常にメモリを食うらしくて、少なくとも私のPCでは一度に全部できない。そこで、以下のマクロを使う。ただこれもかなり時間がかかる。データ量に比例するが、10数分かかったりする。画面上では 『（応答なし）』 とメッセージが出てフリーズしたようになるが、内部ではちゃんと動いているので、あせらずに待つ。なお、以下の”開始・終了位置”は次の4.3 （2）の説明を参照。 なお、これだと各種計算ごとにマクロを作成しなければならないので、セルの数式をコピペするだけのマクロも作った。これはほかでもいろいろ重宝している。 Dim Ci, Ri For Ci = 開始列位置 To 終了列位置 For Ri = 開始行位置 To 終了行位置 ・・・ 計算式 ・・・ Next Next MsgBox "変換終了", vbInformation End Sub ・・・計算式部分・・・ Sx = ((Worksheets(1).Cells(Ri - 1, Ci - 1) + Worksheets(1).Cells(Ri, Ci - 1) + Worksheets(1).Cells(Ri + 1, Ci - 1)) - (Worksheets(1).Cells(Ri - 1, Ci + 1) + Worksheets(1).Cells(Ri, Ci + 1) + Worksheets(1).Cells(Ri + 1, Ci + 1))) / (6 * [ピクセル幅]) Sy = ((Worksheets(1).Cells(Ri - 1, Ci - 1) + Worksheets(1).Cells(Ri - 1, Ci) + Worksheets(1).Cells(Ri - 1, Ci + 1)) - (Worksheets(1).Cells(Ri + 1, Ci - 1) + Worksheets(1).Cells(Ri + 1, Ci) + Worksheets(1).Cells(Ri + 1, Ci + 1))) / (6 * [ピクセル高さ]) Cells(Ri, Ci) = (Sx ^ 2 + Sy ^ 2) ^ 0.5 （2）斜面勾配図 データは”10mDEM”シートのA列2行目以降にあるものとする。計算シートのB3セルに以下の計算式を入れ、あとはデータの数だけコピーする。 =SQRT(((SUM(10mDEM!A2:A4)+10mDEM!A3-SUM(10mDEM!C2:C4)-10mDEM!C3)/(8*ピクセル幅))^2+((SUM(10mDEM!A2:C2)+10mDEM!B2-SUM(10mDEM!A4:C4)-10mDEM!B4)/(8*ピクセル高さ))^2) ［斜面勾配マクロ］ 1番目のシートにデータがあるとする。 ・・・計算式部分・・・ DzDx = ((Worksheets(1).Cells(Ri - 1, Ci - 1) + 2 * Worksheets(1).Cells(Ri, Ci - 1) + Worksheets(1).Cells(Ri + 1, Ci - 1)) - (Worksheets(1).Cells(Ri - 1, Ci + 1) + 2 * Worksheets(1).Cells(Ri, Ci + 1) + Worksheets(1).Cells(Ri + 1, Ci + 1))) / (8 * [ピクセル幅]) DzDy = ((Worksheets(1).Cells(Ri - 1, Ci - 1) + 2 * Worksheets(1).Cells(Ri - 1, Ci) + Worksheets(1).Cells(Ri - 1, Ci + 1)) - (Worksheets(1).Cells(Ri + 1, Ci - 1) + 2 * Worksheets(1).Cells(Ri + 1, Ci) + Worksheets(1).Cells(Ri + 1, Ci + 1))) / (8 * [ピクセル高さ]) Cells(Ri, Ci) = WorksheetFunction.Atan2(1, (DzDx ^ 2 + DzDy ^ 2) ^ 0.5) * 180 / Application.Pi() （3）斜面方位図 基本的な計算式は、傾斜量図とほぼ同じ。 ただ、1章に示した計算式だけでは±90°の180°しか出てこないので、Hyの正負により場合分けする。すなわち、Hy＝0の場合はHxの正負により東（90°）、平坦、西（−90°）を判断し、Hy＜0の場合は算出されたθに180を加える。これにより、北を0°とし、時計回りに−90°≦θ≦270°が区分される。 ［斜面方位マクロ］ 1番目のシートにデータがあるとする。 ・・・計算式部分・・・ Hx = ((Worksheets(1).Cells(Ri - 1, Ci - 1) + Worksheets(1).Cells(Ri, Ci - 1) + Worksheets(1).Cells(Ri + 1, Ci - 1)) - (Worksheets(1).Cells(Ri - 1, Ci + 1) + Worksheets(1).Cells(Ri, Ci + 1) + Worksheets(1).Cells(Ri + 1, Ci + 1))) / (3 * Dx) Hy = ((Worksheets(1).Cells(Ri + 1, Ci - 1) + Worksheets(1).Cells(Ri + 1, Ci) + Worksheets(1).Cells(Ri + 1, Ci + 1)) - (Worksheets(1).Cells(Ri - 1, Ci - 1) + Worksheets(1).Cells(Ri - 1, Ci) + Worksheets(1).Cells(Ri - 1, Ci + 1))) / (3 * Dy) If Hy > 0 Then Cells(Ri, Ci) = WorksheetFunction.Atan2(1, Hx / Hy) * 180 / Application.Pi() ElseIf Hy < 0 Then Cells(Ri, Ci) = WorksheetFunction.Atan2(1, Hx / Hy) * 180 / Application.Pi() + 180 Else If Hx > 0 Then Cells(Ri, Ci) = 90 ElseIf Hx = 0 Then Cells(Ri, Ci) = -9999 Else Cells(Ri, Ci) = -90 End If End If なお、この方位角の表現方法は統一されたものではないらしく、人によって（ソフトによって？）北を基準に時計回りだったり東を0°として反時計回りに表現していたりする。だれかの式を参考にする際には、方角との対応の確認が必要。 （4）地上開度図、地下開度図 これはマクロを使わないと厳しい。こちら データ量とPCスペックしだいだが、数時間単位で時間がかかることもある。まずはPCを省電力設定でスリープしないようにしたほうがよい。次に、試しに10行程度を計算させて処理速度を計っておき、実際に行う行数を処理速度で割って所要時間を予測しておいたほうが安心できる。例えば、10行に30秒かかれば、1200行で3600秒（＝1時間）となる。これまでに最長17時間の計算もさせたてみたが、ほぼ予測通りの所要時間で終わっている。 （5）尾根谷度図 ［尾根谷度マクロ］ 地上開度シートに地上開度計算結果が、地下開度シートに地下開度計算結果があるとする。 ・・・計算式部分・・・ Cells(Ri, Ci) = (Worksheets("地上開度").Cells(Ri, Ci) - Worksheets("地下開度").Cells(Ri, Ci)) / 2 （6）接峰面図、接谷面図、起伏量図、ラプラシアン、エッジ * 方眼法： データの範囲に”DEM”という名前を定義しているものとする。ここでは50m格子（5×5ピクセル）で計算する例を示す。 計算シートのA2セルに以下の計算式を入れ、あとはデータの数の1/5だけコピーする。100m格子の計算をする場合は、 5 を 10 に変える。 接峰面図 =MAX(INDEX(DEM,(ROW()-2)*5+1,(COLUMN()-1)*5+1):INDEX(DEM,(ROW()-2)*5+5,(COLUMN()-1)*5+5)) 接谷面図 =MIN(INDEX(DEM,(ROW()-2)*5+1,(COLUMN()-1)*5+1):INDEX(DEM,(ROW()-2)*5+5,(COLUMN()-1)*5+5)) 起伏量図 =［接峰面図]−［接谷面図］ * 谷埋め法： 中山、隈本（2000）をご参照。 ［マクロ］ 1番目のシートにデータがあるとする。 （厳密な絶対量を求めるものではなく、変化程度の違いを知りたいだけなので、ほぼ正方であればセルの縦横比は気にしなくていいのではないかと思う。） ・・・計算式部分・・・ Cells(Ri, Ci) = WorksheetFunction.Sum(Range(Worksheets(1).Cells(Ri - 1, Ci - 1), Worksheets(1).Cells(Ri + 1, Ci + 1))) - 9 * Worksheets(1).Cells(Ri, Ci) Sub edge() ・・・計算式部分・・・ Cells(Ri, Ci) = Abs(Worksheets("laplacian").Cells(Ri, Ci)) （7）高度分散量図、高度分散異常図 ・ 高度分散量図 データの範囲に”DEM”という名前を定義しているものとする。ここには計算エリアを50m格子（5×5ピクセル）で、計算間隔も50mとして計算する式を示す。 計算シートのA2セルに以下の計算式を入れ、あとはデータの数の1/5 （計算間隔50m ＝ 5ピクセルの場合） だけコピーする。計算エリアを100m格子とする場合は、後半の +5 を +10 に変える。計算間隔を20mとする場合は、 *5 を *2 とする。 高度分散量 =STDEV.P(INDEX(DEM,(ROW()-2)*5+1,(COLUMN()-1)*5+1):INDEX(DEM,(ROW()-2)*5+5,(COLUMN()-1)*5+5)) ”計算エリア” とは、求める点を中心としてどのくらいの範囲の起伏量を算出するか、という意味であり、”計算間隔” とは、求める点の間隔という意味のようだ。すなわち計算エリアより計算間隔が狭ければ、隣接する点同士の計算エリアは重複する。 ちなみに、サンプルエリアとして約20x15kmの範囲の10mDEMにおいて、計算エリアを100m格子（10x10ピクセル）として、計算間隔が50mと100mの場合を比べてみたが、最高標高部以外はほとんど差がなかった。最高標高部はデータ数が少ないので、平均値がその標高の値を代表しないと思われる。つまり、標準曲線を求めるのに、計算間隔は狭くしてもあまり意味がない。 ・ 高度分散量異常図 ： ちゃんとした説明は三箇、藤原（2000）を参照のこと。 【※ 本項を最初に書いた後、間違えた計算をしていることに気付いた。それについて5章に書き加えた。[2012.09.30] 以下は、その後書き直したもの。】 高度分散量を例えば50m格子で求めたとすれば、その50m格子の平均標高を別途求める。 平均標高 =AVERAGE(INDEX(DEM, (ROW()-2)*5+1, (COLUMN()-1)*5+1):INDEX(DEM, (ROW()-2)*5+5, (COLUMN()-1)*5+5)) 50m格子ごとの平均標高と高度分散量を標高25mごとに仕分けする。そしてその標高差25mごとの平均標高のさらに平均、高度分散量の平均、高度分散量の標準偏差を求める。平均標高と高度分散量について標高差25mごとの平均値を求めるには、countifs関数とsumifs関数を使う。すなわち、countifs関数で例えば標高50m以上75m未満のデータ数を数え、sumifs関数でそれに該当する平均標高および高度分散量の合計を求め、平均値を算出する。 =COUNTIFS([平均標高データ], ">=" & 50, [平均標高データ], "<" & 75) =SUMIFS([高度分散量データ], [平均標高データ], ">=" & 50, [平均標高データ],"<" & 75) 標高差25mごとの高度分散量の標準偏差も簡単に求めたいところだが、マクロを使わないなら、少し手間がかかるが、標高差25mごとに高度分散量を取り出しておく方法しか今のところ思いつかない。すなわち、高度分散量を標高順にソーティングする。 標高差25mごとに仕分けるマクロはこちら。 これをエクセルの散布図でプロットする。そしてそのプロットを右クリックすると、下図のようなメニューが出るので、「近似曲線の追加」を選択すると近似式のオプションが表示される。   下部の 「グラフに・・・」 の2つをチェックすると、グラフ内に近似式と決定係数が表示される。近似の種類を変えると決定係数も変わるのでそれを見ながらプロットに合う式を選択する。 私がいくつかやった中では、どうやら4次多項式が一番合うことが多く、次が対数近似のように思う。 この近似式で、ある標高において期待される高度分散量が求められる。したがって、セルごとに標高から計算される高度分散量とその地点の高度分散量との差を取れば高度分散量異常（絶対値）が求まる。 同様に標準偏差の近似式を求めておけば、それで除することにより高度分散量異常（相対値）も求まる。 なお、近似式によっては低標高で期待される高度分散量がマイナスになることがあるので、その場合は 0 にしちゃって構わないと思う。 それと、4次多項式が最善と判断された場合、グラフに表示される高次の係数の有効数字が1桁なので、実際とかなり異なる。このため、別途係数を求めたほうがよい。求めるには、INDEX関数とLINEST関数を使う。 （8）デコボコ図 こんな名前の図があるのかどうか知らない。私が勝手に名づけた。 斜面の凹凸を示す指標に曲率があるが、偏微分方程式を解かねばならないようで、簡単には計算できそうもない。しかし、凹凸をみるだけなら隣接セルとの標高差をとるラプラシアンがある。 ArcGISの曲率ツールのヘルプをみると、隣接セルとの標高差をセル間隔の二乗で除する係数があった。これを参考に以下の式で作図してみたら、緩斜面の微地形がよくわかるものができた。 4近傍ラプラシアンとは、セル間隔で除するかどうかが異なっている。ということは、DEMデータが完全な正方グリッドであれば同じものになる。しかし、国土地理院の10mDEMは0.4秒グリッドであり、北海道の北部になると8m×12mグリッドとなるので、セル間隔を考慮しないラプラシアンだと東西成分がやや強調されることになる。この差が有意なのかどうかが問題なのかもしれないが。 ［4近傍デコボコ図マクロ］ 1番目のシートにデータがあるとする。 ・・・計算式部分・・・ Dk_x = ((Worksheets(1).Cells(Ri, Ci - 1) + Worksheets(1).Cells(Ri, Ci + 1)) / 2 - Worksheets(1).Cells(Ri, Ci)) / Dx Dk_y = ((Worksheets(1).Cells(Ri - 1, Ci) + Worksheets(1).Cells(Ri + 1, Ci)) / 2 - Worksheets(1).Cells(Ri, Ci)) / Dy Cells(Ri, Ci) = - (Dk_x + Dk_y) / 2 ［8近傍デコボコ図マクロ］ 1番目のシートにデータがあるとする。 ・・・計算式部分・・・ Dk_x = ((Worksheets(1).Cells(Ri, Ci - 1) + Worksheets(1).Cells(Ri, Ci + 1)) / 2 - Worksheets(1).Cells(Ri, Ci)) / Dx Dk_y = ((Worksheets(1).Cells(Ri - 1, Ci) + Worksheets(1).Cells(Ri + 1, Ci)) / 2 - Worksheets(1).Cells(Ri, Ci)) / Dy Dk_xy = ((Worksheets(1).Cells(Ri - 1, Ci - 1) + Worksheets(1).Cells(Ri - 1, Ci + 1) + Worksheets(1).Cells(Ri + 1, Ci - 1) + Worksheets(1).Cells(Ri + 1, Ci + 1)) / 4 - Worksheets(1).Cells(Ri, Ci)) / (Dx^2 + Dy^2)^0.5 Cells(Ri, Ci) = - (Dk_x + Dk_y + Dk_xy * 2) / 4 （9）平滑化 この項を最初に作成した頃は SUMPRODUCT関数 を知らなかったので、どこかのサイトを見て以下の前半の計算式を書いた。あとで、後半の計算式にあるWorksheetFunctionなるワークシート関数があることを知った。 ［マクロ］ 1番目のシートにデータがあるとする。 ・・・計算式部分・・・ Cells(Ri, Ci) = (Worksheets(1).Cells(Ri - 1, Ci - 1) + Worksheets(1).Cells(Ri - 1, Ci) + Worksheets(1).Cells(Ri - 1, Ci + 1) + Worksheets(1).Cells(Ri, Ci - 1) + Worksheets(1).Cells(Ri, Ci + 1) + Worksheets(1).Cells(Ri + 1, Ci - 1) + Worksheets(1).Cells(Ri + 1, Ci) + Worksheets(1).Cells(Ri + 1, Ci + 1) + Worksheets(1).Cells(Ri, Ci)) / 9 = WorksheetFunction.Average (Range(Worksheets(1).Cells(Ri - 1, Ci - 1), Worksheets(1).Cells(Ri + 1, Ci + 1))) その後、セルの数式をコピペするだけのマクロを使うようになったので、この9セルの平均を取るだけなら AVERAGE関数 でできちゃう。また、SUMPRODUCT関数を使えば例えば5×5の加重平均化フィルタによる平滑化も簡単にできる。
4.3 XML文書作成 ようは、標高値を各種計算結果の値に置き換える、ということをくどくど書いていた。（2022.03.24記） 4.4 閲覧コンバートソフトの使用法 データを表示するための基盤地図情報閲覧コンバートソフトは、国土地理院の基盤地図情報のダウンロードサービスよりダウンロードする。zipファイルを解凍し、Start.cmdファイルをダブルクリックするとインストールされる。 使い方は人それぞれなので細かいことは書かない。私がよく使うもののみ。 （1）ファイルの読み込み ファイルを読み込む方法は４種類ある。
既存ファイルのデータの一部を入れ替えるのに、バージョン3.10のときはプロジェクト設定で読み込まれているファイルを解除し、新たのを追加することができた。例えば、海岸線や水涯線データを流用して傾斜量図を開度図に入れ替えることができた。ところが、バージョン3.20になったらプロジェクト設定からファイルを解除できなくなり、新規プロジェクト作成で全部をいちいち読み込み直さなくてはならなくなった。 同じ地域のDEMデータを読み込み直すことは想定されていないからなのかな、、、 （2）表示設定 DEM表示の設定は、 設定 ＞ 表示設定 ＞ DEM表示設定 で行う。DEM表示設定には連続段彩とステップ段彩とがあるが、これは目的によってお好きなほうをどうぞ。ここでは連続段彩で出てくる左側の標高グラデーションについて、裏技的な使い方をひとつを書く。 この使い方がわかりにくくて覚えるのに時間がかかったのだが、面白いことがわかった。まず、表示されている数値を書き換えるのは、『選択された「ペグ」標高値(m)』の欄でしかできない。最高・最低標高値を書き換えるには、その数値の左の ■ をクリックして『選択された・・・』欄に表示させる。そしてその欄の数値を書き換えたらクリックしないでそのまま待つ。パソコンのスペックによって表示速度が変るのか変らないのかわからないが、だいたい2秒待ては最高もしくは最低標高値が置き換わっている。この方法は、▼ペグを選択した場合も同じ。 このソフトは基本的には地形図の表示ソフトなので、0.1m単位で表示されるが、『選択された「ペグ」』欄での入力は1m単位でしかできない。『「ペグ」を追加』欄だと0.1m単位で追加できるのだが、それを修正しようと▼ペグをクリックしたとたん、少数は切り捨てられてしまうことに注意。 普通に地形を表示させる場合にはよいのだが、小数点以下2桁以上の設定をしたい場合はそのまま入力できない。例えば最低値が0、最高値が0.5だと、ペグは最大4つしか入れられない。10段階に変化させたいときはどうするか？ その答えは、まず最高値を10とする。そして『「ペグ」を追加』で9つ追加すれば10段階に変化させられる。その後、最高値をもとの0.5に戻す。もっと細かい変化を設定したければ、実際の値は別として、例えば0〜1000mを仮に指定すると千分率で変化点を指定できる。 ここで、確認のために例えば1.5のペグをクリックしてみると、『選択された・・・』欄に 1.5 という数値がでてくるはずである。しかし、すぐに自動的に 1 の表示になり、ペグも1のところに移動してしまう。つまり、確認しようとするとやり直しになってしまうので、慎重にどうぞ。 ステップ段彩を選択すると、1m単位でしか入力できない。しかし、連続段彩で細かく設定した後、もういちどステップ段彩で表示設定画面を開くとその設定でステップ段彩にできる。 （3）ファイル書き出し ”ファイル”メニューには、「表示画面を画像データとして保存」という項目がある。もちろんこれを使ってもよいが、小さな解像度のものしかできない。 同じく、「印刷」の項目では文字通りプリンタ出力されるのだが、イメージ出力もできる。この「イメージを作成して・・・」を使うと、高解像度の画像データが得られる。 | ||||||||||||||||||||||||||||||||||