基盤地図情報DEMをOCADにインポートする方法と雑学
はじめに
（公開後に問題点を一つ追記と、本文も修正しました）
実はOCAD最新版では基盤地図情報のDEMが直接インポートできるんですよ。いつの間にかそんなアップデートがされています。いつの間にかというか実は、2月くらいにそのアップデートが入って、ヨルクさんから試してみて！ってメールが来てたのですが、すっかり忘れてGWになってしまいました。
で、早速試してみたところ、すんなりインポートできました。OCAD最新版をお持ちの方は試してみてください。
↑ファイル形式でGMLファイルを選びます
で、で、これまで通りQGISで作ったデータと比較しようとしたのですが、PCの買い替えなどもあり、QGISの使い方もすっかり忘れていたので、頑張って思い出しました。
最近は基盤地図情報もたまにしか使う機会がなくなり、久しぶりに使おうと思うとなんだかいろいろ手間取るものです。よってもう誰かが作っていたかもしれないのですが、備忘のために手順を書き留め、せっかくなのでそれをシェアしておこうというのが本記事の趣旨その1です。
趣旨その2は、その変換方法によって生成される等高線に微妙な差が生じること。この辺、実用的には大きな問題にならないのですが、せっかく得られた知見をシェアしたくなった、というというところです。
記事中で、現状のOCADの基盤地図情報直接インポート機能の問題点を指摘しますが、これはOCAD側にフィードバックするのできっとそのうち改善されるはず。また、実用上はぶっちゃけそんなに大きな問題ではないです。とすると、QGISを噛ませた七面倒臭い手順など完全に不要になってしまい、趣旨その１は無意味になる可能性もあるのですが、まぁ話の種の雑学として残しておこう、という思いです。GWで暇だしね。
（追記：実用上問題ある場合もあり、後ろに一節追記しました。）
また、最後のところにリンクを張りましたが基盤地図情報の変換に関しては別の方がわかりやすく記事にしていたりします。本記事で紹介する手順は、それとは微妙に異なるかもしれません。いわば「別解」みたいなものです。そういう意味でも、雑学的な意味の方が大きい記事だと思います。
生成方法によるコンタの違い
上に３つ基盤地図情報で生成したコンタの比較を示します。コンタ自体はすべてOCADの機能で生成したもので、違いはOCADに読み込ませるデータの下準備だけです。
①はなんか遠目に見るとブロックっぽく？カクカクしてるような感じがします。②は①より少しカクカク感が減ってるような気もします。③はだいぶなめらかに見えます。（ちなみに、基盤地図情報をFGDVでXYZに変換してOCADにインポートするのも、同じような感じになると思います。試してないけど）
また、①と③を比べると、コンタの位置が少しずれています。図には出していませんが、②と③の位置はほぼ同じです。VectorMapMakerでコンタを生成することもできますが、そのコンタも②、③とほぼ同一の位置に来ます。要はOCADで直接基盤地図情報をインポートすると、なんか位置がずれるんです。これはなんででしょうか？
基盤地図情報の変換って何をやってるの？
まぁファイル形式を変換してるってのもあるんですが、もっと肝心な部分として座標系の変換をしています。基盤地図情報は緯度経度座標です。まぁ日本全土を統一でかカバーしようと思えば普通そうなります。一方、O-MAPなどの地図では平面直角座標系を使いますで、この座標変換をしなければいけません。
これは↑のようなイメージになります。基盤地図情報では緯度経度座標できりのいい間隔（大体５m）で標高データが並んでいます。これを平面直角座標に直すと標高データの位置がずれることになります。これは上手いことやれば同じ点を使うように回避できそうなもんですが、緯度経度は極座標でこれを平面座標に変換したら厳密に言えばどうやったってズレが生じるという本質的な問題があります。
で、このようなことをするときに、どうやって新しい座標系での値を決定するか、というところですが、いくつかの方法があります。
・最近傍
最も近い点の値をそのまま持ってくる。
・バイリニア
隣接する４つの値の加重距離平均値を取る
・キュービック
隣接する4×4セルから曲線で補完して値を取る
他にもありますが、省略します。おそらくOCADでは最近傍を使っているのではないかと思います。このアルゴリズムの問題点は、セルサイズの最大1/2分だけ水平位置がずれてしまうことです。先に示したQGISで変換した場合との位置のズレの原因はこれです。ちなみに、バイリニアとキュービックは自分が試した限りほとんど差が無いように思います。
ちなみに、なんでOCADが最近傍を使ってるのか、ですが、OCADでは基盤地図情報のDEMデータをグリッドに区切られたデータではなく、生データ扱いしているからではないかと思います。生データというのは、標高データの位置がランダムに散らばってるようなデータです。こういうデータからだと、最近傍が使いやすいのかな？と思います。あんま詳しくないので、もしかしたらそういうデータからでも補間できるのかもしれないけど。
追記：海岸エリアのコンタのズレ方を見るに、「最近傍」を使っているとも言い切れないかもしれない…分からなくなってきた。
QGISを使用した基盤地図情報DEMの変換
では、ここから具体的にQGISを使用したデータ変換の手順を示します。
①QGISの用意
面倒なことに、使用するプラグインが最新版のQGISに対応していないので、敢えて古いバージョンを使用します。まぁ最新版を使用したい場合は別のソフトを使う手もあるんですけどね。
下記URLからQGIS2系列のものをダウンロードし、インストールします。なんでもいいですがとりあえず最後の2系列である2.18.28にしておきましょう。
qgis.org/downloads/
②基盤地図情報変換プラグインを入れる
https://github.com/minorua/fgddemImporter/releases
からzipをダウンロードして解凍します。それをC:\Users\(ユーザー名)\.qgis2\python\pluginsにフォルダごとコピーします。
③QGIS desktopを起動
ツールバーのプラグインに「基盤地図情報DEMインポーター」というのがあるはず。なかったら、「プラグインの管理とインストール」でインストール済みのところにfgddemimporterがあって、チェックが付いているかを確認。なかったら、コピーするフォルダか何かが間違っていると思います。
④インポートするファイルに基盤地図情報からダウンロードしたzipファイルを追加。
出力先フォルダに全角文字があるとだめなので、適当なフォルダを指定。うまくいけばtiffファイルが出てくるはずです。
⑤座標変換
ラスター→投影法→ワープ
保存ファイル名、変換先の座標系を入力。あと、下のコマンドのところの”-of ERS”を消さないと何故かうまくいきません。
⑥ascに変換
ラスター→変換→形式変換
で、画像の通り指定して.ascに変換します。このとき出力サイズを200%程度にすると、ほんのちょっとなめらかなコンタが作れます（その代わりデータが4倍になる）。最初の画像の③はこれをしています。
⑦OCADのDEMインポートウィザートで.ascファイルを読み込む
あとは、通常の手順通りに…
もう一つの問題点：欠損値の扱い
もう一つ、問題なのが海や川の部分、その他データが欠損しているエリアの扱いです。下図の青で塗りつぶされていないエリアがデータがない川の部分なのですが、こういうところのコンタが汚くなってしまっています。これはデータがない部分を無理やり補間しようとしているためです。
QGIS経由の場合はこうはなりません。（ちゃんとデータ無しをデータ無しをとして扱い無理に補完したりしない）
もっと問題なのは、データ無しの補間に釣られて、データのある部分のコンタまでおかしくなっていることです。さすがにこれは実用上も問題があります。下図では堤防のコンタや鉄道の下のコンタの位置がズレて道路との位置関係がおかしくなっています。これは困ったことです。
そんなこんなで
久々にやると、人間いろいろ忘れるもんですね。なかなかQGISが上手く使えなくて苦戦しました。プラグインがQGIS3で使えないというのも罠でした。一方で、ほんのちょっとコンタを滑らかにする裏技もみつかりました。まぁ実用上は誤差みたいなもんですが。でもちょっとだけ見栄えが良くなるから、無修正でそのまま地図にしちゃうときにはいいですよね。
OCADも積極的に基盤地図情報とか対応してくれようとしているようです。ぶっちゃけそこより作図する上での微妙なストレスを軽減してほしいなぁとか思ってるんですが、それもちゃんと具体的に要望を出せば改善してくれるかもしれません。いわゆる「誰も消防車を呼んでいない」状態なのかもしれません。
補足
・DEM以外も含めた基盤地図情報の活用について、↓の記事がわかりやすいと思います。（DEMに関しては工事中みたい）。
今回はQGIS2でしか使えないプラグインを使用しましたが、別のソフトでtiffに変換する手もあります。（↑記事で紹介されてます）
・OCADを使わずにコンタを作るにはこの辺をどうぞ。OCADでコンタを作るメリットとしては、計曲線、補助曲線などのかき分けが簡単なこと、OCAD側で標高値が確認できて、コース設定時に登距離が簡易的に計算できたりすることかな。
・区域によってはDEMデータのA（航空レーザー測量）とB（写真測量）が混在している場所があって、それらをマージする方法も今回確認しました。QGISのラスター→その他→結合で合体できます。QGIS便利ですね。