エクセルで作る地形解析図(3)

目次参考

3.作図の流れ

 閲覧コンバートソフトで表示させるためには、大きく3段階の作業を行う。
  1)10m間隔の標高値が入っているXML文書をダウンロードする。
  2)その標高値をエクセルを使って書き換える。
  3)再びXML文書にする。
 GRASS GISを入手していれば、3)でXML文書にするより、書き換えたエクセルデータをそのままテキストエディタにコピペし、位置座標等の情報を加えてGRASS GISに取り込み、GeoTiff画像にしたほうが簡単だし、そのままGISソフトで表示・加工できるようになる。(2014.8.12追記)


3.1 XML文書のダウンロード

 基盤地図情報には、JPGIS2.x形式とJPGIS2.x(GML)形式の2つある。これが技術的にどう違うのかはさっぱりわからないが、とりあえずファイルサイズはGML形式のほうが小さく、使いやすいので、ここでは『基盤地図情報 数値標高モデル』 JPGIS2.x(GML)形式をダウンロードし、利用する (2014.7.31以降、ダウンロードファイルはJPGIS2014(GML)形式だけになった)
 最後に作成するのもGML形式のほうがよさそうなのだが、試行錯誤してみたがGML形式で作成してもなぜか閲覧コンバートソフトで読み込めなかった。ずいぶん長い間、原因がわからなかったのだが、ある時、ファイル名とファイルサイズに制約があることを知った。zip形式での圧縮率も全然違うので、作成するのはGML形式ではない方が良い。ただし、今後閲覧コンバートソフトがバージョンアップした場合、読み込めなくなる恐れがあるように思う (とりあえず2014.7.31のバージョンアップでは大丈夫のようだ)。

以下、ダウンロード画面も変わっているが、変更画面の説明は省略

 ダウンロード画面は、2次メッシュを1つずつインデックスマップで指定してダウンロードするようになっている。
 例えば、利尻島を含むファイルは、 6741 という1次メッシュ番号がついている。
 前2桁は範囲の南端の緯度の1.5倍の数値となっている。この例である 67 であれば、南端は北緯44.666667度(44度40分)となる。後2桁は範囲の西端の経度から100引いた数値となっている。この例である 41 であれは、東経141度となる。
 1次メッシュ番号をクリックすると2次メッシュマップが出てくるので、ここでファイルを指定する。
 余談だが、数字の順番というのは1、2、3、・・・と1から始まるのが普通だと思う。しかしコンピュータの世界では0から始まるのが常識らしい。2次メッシュ番号でも、8列に対して0〜7の数値が割り振られている。下記のXML文書内の記述でも、行、列の最大値は1124、749となっていて、最初はなぜ1つ小さいのかわからなかった。
【追記】 画像解析の本を見ていたらこんな記述があった。
『 コンピュータプログラムでは、通常、C の場合は0,1,2,・・・とし、Fortranの場合は1,2,3,・・・とすることが多いが、特に決まったものがあるわけではない。多くの場合は本質的な違いはないので、様々な表現があるということを知っておけば実務上困ることはない。』


3.2 標高値の書き換え

(1)XML文書より必要なデータの取り出し

 地理院のHPには 『応用スキーマによって定義された構造を持つデータ』 とあって定義ファイルも公開されているが、見ても結局よくわからなかった。でも、全部を理解できなくても、必要な部分が取り出せればここで作業するのにはそれで問題ない。
 まず、DLしたzipファイルからXMLファイルを取り出し、ファイル名の一番最後に .csv と拡張子を書き足す (なお、windowsの設定で拡張子を表示するようにしておかないとできないかもしれないので、要確認)。そしてこのcsvファイルをダブルクリックすれば普通にExcel2007、2010だったら開けるはずである。Excel2003だと行が足りなくて開けない(2022年現在、csvにして開くと文字化けする。2016年にxmlファイルの文字セットがSHIFT-JISからUTF-8に変更されているが、Excelはそれを認識できていないのではないかと思う)
 最下行が843,811行で ”</Dataset>”となっており、48行から843,797行のA列に ”その他”、B列に ”数値”が入っていれば正常に開けている。この843,750個の数字が標高データである。
 ちなみに、GML形式でないほうだと506万行のデータが入っていてエクセルでは開けない。

 このファイルの中でほかにチェックしておいたほうがよいのはデータ範囲を示す緯度経度の2行だけで、後は気にしなくてよい。
   27行目 <gml:lowerCorner>45.333333333 142</gml:lowerCorner>
   28行目 <gml:upperCorner>45.416666667 142.125 </gml:upperCorner>

 気にしなくてよいと書いておきながら何だが、知っておいたほうがよいのは測地系の知識 (この説明は、例えば私が愛用しているフリーウェアソフトであるカシミールのHPを参照)。緯度経度の前の26行目に <gml:Envelope srsName="fguuid:jgd2000.bl"> とある。すなわち、この緯度経度はJGD2000によっていると記述されている。最近は、おそらくたいていはJGD2000かWGS84が使われていて、実用上はほとんど差がない。しかし、ほんのちょっとのずれが問題になるときは、測地系のチェックも必要かもしれない。
  5mDEMの公開範囲がどんどん広がっている。5mDEMでも、ここで書いていることはそのまま使える。ただ、このGML形式ファイルで10mDEMとやや異なるところがある。
 まず、ひとつのファイルに入っているデータの数が違うので、下の表形式に変換する際には何列のデータなのかチェックが必要。
 そして10mDEMでは、私が見る限りではcsvで開いたとき、データ行のA列はすべて ”その他” となっていたが、5mDEMではA列に ”データなし” とか ”地表面” とか、いろいろあるようだ。でもこの文字は気にしなくてよい。気にしなければならないのは、下から5行目の startPoint タグ。例えば、
   <gml:startPoint>15 36</gml:startPoint>
こうなっている場合、データは37行16列から入っている、ということを示す。
 さらに上記文字と関連しているのかもしれないが、河川などのデータが -9999 になっていることがある。計算する際には例えば、IF文を使って計算対象範囲に -9999 があれば -9999 にする、というような対応が必要である。


(2)データを表形式に変換

 843,750個の数値は、750行×1125列のデータを北西から南東に向かって順に並べられたものである。したがって、843,750行×1列のデータを並べ替えればよい。
   1行目から1125行目の1列を、縦横変換して別なシートの1行目に1125列コピー。
   1126行目から2250行目の1列を、同じく別なシートの2行目にコピー。
 この操作を750回繰り返せばよい。  って、やってられないので、関数を使うなり、マクロを組むなり、ご自由にどうぞ。私はマクロを使い、1行目はメモ行とするために空けて2行目からデータを並べるようにした(4章)。
 表がそのまま地形図の座標になるので、表示を一番縮小すると、地形図らしき模様が見える。
 なお、DEMデータで海は -9999 が入っている。これをそのまま計算すると海岸付近はとんでもない値となる。関数を加えればそれでも対応できるだろうが、私は面倒なので全部 空白 に置換するようにした。ただし、後述する計算式の中には空白セルがあるとエラーとなるものもあるので、その場合は該当セルに 0 を入力するようにした。すなわち、標高0mの真っ平ら平原。最終的な図面にするときは、海岸線からはみ出した部分を画像ソフトでカットすればよい。1章の図もそうしている。
 最初は画像ソフトで処理していたが、その後、5mDEMでは No Data(-9999)が普通にあるので、元DEMの計算範囲に-9999があれば-9999とするようにした。
 =if(countif([計算範囲],-9999)>0,-9999,[各種計算])


 この表の別な利用法として、緯度経度を2点指定することで、その点を通る断面図をエクセルのグラフ機能で即座に作図させることもできる。これはもちろん後述の計算結果の方でも使える。つまり、地形断面と同じ測線の地形解析断面を並べて表示させられる。ほかの情報、例えば断層などの構造線位置なども合わせて表示させたりと、いろいろ遊べる。例えば・・・


(3)各種計算

 本題の計算内容については4章へ。
 基本的に、1章に例示したように、着目するセルと隣接するセルとの間で足したり引いたり掛けたりするだけである。


3.3 XML文書の作成

(1)表形式のデータを変換

 計算が終わったら、今度は表形式のデータを1列のデータに並べ替える。そして元のXML文書のデータ行と入れ替えればよい。
 のだが、これが単純ではない。データを並べ替えるのは上記の逆で単純でよいのだが、問題は計算が終わった表は必ずしも750行×1125列ではないということ。例えば、傾斜量など、隣接セルとの計算を行うと、縁のセルは隣接セルがないのでその部分のデータはなくなる。開度図であれば、検索距離を例えば250mとすると、隣接28セルを計算対象とするので、28列のデータがなくなる(【参考】参照、開度計算ではなくならないように改良)。
 XML文書では、緯度経度で範囲を指定し、そこに何行何列のデータがあるか、ということで示している。したがって、左側が28列なくなれば、データの範囲は 28列×0.4秒 だけ東に寄ることになる。元のXML文書の緯度経度および何行何列あるかを指定している箇所を書き換えることになる。これが理解できていれば、50mDEMでも、その他自分で作成した任意のピクセルサイズのデータであっても、同様に閲覧コンバートソフトで表示できる (その後、なくなったセルには-9999を入れて、緯度経度や行列数を変えないようにする方が簡単かつ確実と考えを変えた)


(2)XML文書として保存

 元のGML形式の文書で、データの位置と、緯度経度、行・列数を入れる位置はすでにわかっているので、ここに計算された数値を入れればよいはず、、、 なのだが、どこかに私の気づかないミスがあるのか、閲覧コンバートソフトは読み込んでくれない。それで、とりあえずGML形式で作成するのは頓挫した。いずれまた、試してみようとは思ってる。
【時々気になってあれこれやっていたのだが、ようやく分かった。以下追記。】
【GML形式での保存】 -- 2013.12.21追記 --

 これまで読み込めなかった理由がわかった。
 まずは、コピーミス。オリジナルのXMLファイルから最初の48行ほどをコピーしてエクセルに貼り付け、座標等の数値を書き換えてからテキストエディタにコピーしていた。その時、2行ほどになぜかタブと引用符が余計に加わっていた。それでGML形式ファイルとして不完全なものとなっていた。
 それと、XML文書ファイルの名前とサイズに制約があった。ファイル名はこうじゃないとGML形式のファイルとして読み込んでくれないらしい。

   FG-GML-文字列-文字列-DEM-数字.xml  もしくは
   FG-GML-数字-文字列-DM-文字列.xml

アルファベットは固定文字で、ハイフンでつなぐ。数字部分がコンバートソフトではレイヤーリストの名前として表示される。実はデータファイル名の定義も今年になって地理院HPに載っていたことに、今頃になって気付いた。数字は基本的にメッシュ番号など位置を示すものとなっている。
 さらに、ファイルサイズが21MB程度より大きいと読み込まないこともわかった。

 最近、2475行×4950列(約25km×50km)のデータで作図していたので、ファイル名とサイズを変えてGML形式でも作成してみた。
 GMLではない形式だと、テキストエディタの都合でXML文書ファイルは3分割することになり、ファイルサイズは全部で1200MB程度となる。これをzipファイルにすると圧縮率がだいたい97%で35MB程度になった。
 GML形式だと、コンバートソフトの読み込みファイル容量の都合で8分割しなければならないが、ファイルサイズは全部で150MB程度と小さい。しかしzipファイルにすると圧縮率がだいたい84%で24MB程度にしかならない。
 作成手間を考えると、必ずしもGML形式の方が便利とは言えない。

 また、閲覧コンバートソフトがかなり制限をもっていることがわかった。現在、わりといい加減なファイル名を付けているのだが、今後バージョンアップがあった場合、ファイル名が合わない等により読み込めなくなる恐れがありそうだ。



目次参考