4.作図方法4.1 表形式データの作成 843750行×1列の表を750行×1125列にする表変換マクロはこれだけ。この例ではB48〜B843797セルにあるデータを、F2〜AQL751セルに出力させるようにしている。 マクロが入ったエクセルファイルを開いておく。ダウンロードしたGML形式のxmlファイルの拡張子をcsvに変え、ファイルを開いて表示した状態で、『alt+F8』キーからこのマクロを実行させると、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) = 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) - 8 * 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)) ”計算エリア” とは、求める点を中心としてどのくらいの範囲の起伏量を算出するか、という意味であり、”計算間隔” とは、求める点の間隔という意味のようだ。すなわち計算間隔が狭ければ、隣接する点同士の計算エリアは重複する。 ・ 高度分散量異常図 : ちゃんとした説明は三箇、藤原(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ごとの平均高度のさらに平均、高度分散量の平均、高度分散量の標準偏差を求める。これをエクセルの散布図でプロットする。そしてそのプロットを右クリックすると、下図のようなメニューが出るので、「近似曲線の追加」を選択すると近似式のオプションが表示される。 ![]() ![]() 下部の 「グラフに・・・」 の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関数 を知らなかったので、どこかのサイトを見て以下の計算式を書いた。 [マクロ] 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 その後、セルの数式をコピペするだけのマクロを使うようになったので、この9セルの平均を取るだけなら AVERAGE関数 でできちゃう。また、SUMPRODUCT関数を使えば例えば5×5の加重平均化フィルタによる平滑化も簡単にできる。
4.3 XML文書作成 XML文書は3つの部分で構成されている。 (1)前段部分 行数が86行と多いのでテキストファイルをリンクする。 この中で、62〜65行目のデータ範囲を示す緯度・経度と、84行目のデータの行数・列数を書き換え、テキストエディタにコピーする。なお、データ数がエクセルの行数に収まるのであれば3.3(2)の記述のようにそのままエクセルを使ってもかまわない。 緯度・経度は、ダウンロードしたデータファイルのデータをそのまま使うのであれば、それと変らない。しかし、必要な範囲だけ取り出したり、あるいは開度図のように縁のデータがなくなってピクセル数が減った場合は、その分を計算しなければならない。10mDEMであれば1ピクセルは0.4秒であるのでピクセル数をカウントしてもよいし、50m格子の場合でもやはりその分のピクセル数をカウントしてもよいが、行・列数と端の緯度・経度があるのだから計算したほうが楽だ。 次の(2)に示す開始・終了位置により、緯度経度の範囲を計算する。以下の緯度・経度範囲、データ行・列数は各種計算前のDEMデータ部分のもの。すなわち、傾斜量や開度図でデータがなくなったピクセルも含んでいる。かつ、DEMデータはA2セル、すなわち2行目以下にあるものとしている。 <jps:westBoundLongitude>行 =西端経度+(東端経度-西端経度)/データ列数*(開始列位置-1) <jps:eastBoundLongitude>行 =東端経度-(東端経度-西端経度)/データ列数*(データ列数-終了列位置) <jps:southBoundLatitude>行 =南端緯度+(北端緯度-南端緯度)/データ行数*(データ行数-終了行位置+1) <jps:northBoundLatitude>行 =北端緯度-(北端緯度-南端緯度)/データ行数*(開始行位置-2) 行数・列数のカウントは、0 から始まっていることと、列 行 の順になっていることに注意。すなわち、10行×20列だとすると、ここに記入するのは <jps:coordValues>19 9</jps:coordValues> となる。 (2)データ部分 ![]() データを複数列に変換した場合は、当然ながらこの数式でデータ行を複数列作成してから、順にテキストエディタにコピーして下に付け加えていく。 以下のマクロは1セルずつ指定しているので、例えば2000行×2000列のデータの場合10分程度を要していた。これを、先にまず30万行×13列+1列に変換(行単位で縦横変換)し、その後文字列を追加するようにしたら、数10秒で終わるようになった。 また当然ながら、この部分を例えば ポイント番号、緯度、経度、標高 の順で並べてテキスト保存すればPNEZ形式データとなるので、CADやGISソフトで読み込めることになる。 ”開始・終了位置”とは、 「データの入っているWorksheetのうち、20行50列目のセルから80行70列目のセルまでの範囲を変換する」 という意味の行・列のセル位置である。PCのスペックによると思うのだが、データ数が多いとメモリ不足で処理できなくなるので、分割してファイルを作成する必要がある。この行・列は、表データ部分のものではなく、エクセルシートの行・列番号としている。すなわち、データが B5セル から始まっていれば、5行2列目が開始行・列となる。 データは1番目のシートに入っているものとする。 Dim xb, yb, xe, ye As Integer Dim i, j, gyo, retu, K, Ra As Integer Ra = [集約列数] '1列にまとめる列の数 xb = [開始行位置] yb = [開始列位置] xe = [終了行位置] ye = [終了列位置] gyo = xe - xb + 1 retu = ye - yb + 1 ' For i = 1 To gyo For j = 1 To retu K = (i - 1) Mod Ra Cells(j + 4 + retu * K, Int((i - 1) / Ra) + 11) = "<jps:values><jps:memberValue><type>その他</type><alti>" & Range(Worksheets(1).Cells(xb, yb), Worksheets(1).Cells(xe, ye)).Cells(i, j) & "</alti></jps:memberValue></jps:values>" Next Next MsgBox "変換終了", vbInformation End Sub (3)後段部分 ここは変更する箇所がないので、以下をそのままデータ部分に付け加える。 <jps:type>linear</jps:type> <jps:scanDirection>+x -y</jps:scanDirection> </jps:sequencingRule> <jps:startSequence> <jps:coordValues>0 0</jps:coordValues> </jps:startSequence> </jps:CV_GridValuesMatrix> </jps:valueAssignment> </coverage> </DEM> </dataset> </GI> 4.4 閲覧コンバートソフトの使用法 データを表示するための基盤地図情報閲覧コンバートソフトは、国土地理院の基盤地図情報のダウンロードサービスよりダウンロードする。zipファイルを解凍し、Start.cmdファイルをダブルクリックするとインストールされる。 使い方は人それぞれなので細かいことは書かない。私がよく使うもののみ。 (1)ファイルの読み込み ファイルを読み込む方法は4種類ある。
既存ファイルのデータの一部を入れ替えるのに、バージョン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)ファイル書き出し ”ファイル”メニューには、「表示画面を画像データとして保存」という項目がある。もちろんこれを使ってもよいが、小さな解像度のものしかできない。 同じく、「印刷」の項目では文字通りプリンタ出力されるのだが、イメージ出力もできる。この「イメージを作成して・・・」を使うと、高解像度の画像データが得られる。 |