Code for History

"Code for History"はIT技術を歴史学上の問題の解決に使うコミュニティです。強調したいのは、我々にとってIT技術は「手段」であって「目的」ではありません。「目的」は歴史学上の問題を解決する事であって、必要であればITでない手段も活用します。常に最優先なのは、問題を解決することです。

国土地理院DEMから段彩陰影図作成の手順

国土地理院基盤地図情報の5m DEM提供範囲段彩陰影図作成*1が盛り上がってたりしたのですが、元がジオメディア系出身の私は正直あまり価値が判ってませんでした。
が、歴史国土*2の企画でGeoアクティビティフェスタの優秀賞をいただいたり、それを使って作った自分向けの居住地周辺旧版地形図を持って息子と近所探検等をするうち、かなり地形への興味が湧いてきています。

それで自分でもDEMから地形図を作ってやろうと思ったんですが、
東京と違って地方なので地元のデータ作って話題振ってもあまり盛り上がりませんし、サービス系出身の事もありどうしても限られた地域のデータだけ作るよりは全国作りたいよねーとか思ったり、5m DEMと10m DEMの精度以外の違いがあまり判ってなかった事もあって、
まあ5m DEMはすごいけど4分の1のデータ密度でも全国あればそれなりに面白いもんできるんじゃないのかなと思い、4日ほどかかって全国の10mDEM落としてきました。
並行して段彩陰影図の作り方をあちこちのサイト調べたり、ジオな方々に聞いて教えてもらったりして、なんとかスクリプト等でTMSタイルを作成自動化できる形での手順を身につけました。
ズバリでなければあちこちに手順あるのですが、パラメータやオプションも含め手順出せばそれなりに意味もあるのではと思い、共有します。*3

基盤地図情報の取得

基盤地図情報はの元データは、先にも述べた国土地理院のサイトから落としてください。
アカウントを登録し、基盤地図情報ダウンロードサービス > ログイン > 数値標高モデル(JPGIS形式、GMLじゃない方)選択 > 欲しいデータタイプ(5mか10mか等)と地域選択 > その地域に含まれるファイル一覧が出るのでダウンロード、です。
ファイル一覧からさらに欲しい地域を絞り込みたい場合は、地域メッシュコードというものをベースにファイル名が名付けられているので、それを元に欲しい地域を絞り込んでみてください。
地域メッシュコードは、こちら等で調べられます。
この記事では、陸域と海域がいい感じで含まれている三重県志摩市の地域メッシュ5136-36付近で試してみたいと思います。*4

基盤地図情報からのDEMの作成

基盤地図情報からのDEMの作成は、OSGeo財団のプロジェクト*5基盤地図情報を読み込めるように改造したGDAL/OGRが公開されているので、それを利用させてもらいます。
コンパイルされたものはWindows版しかないですが、Linuxでは普通にコンパイル通るようですので、普通に./configure 、make、 make installでOKです。
ただバージョンが1.9.0ベースなので、既に1.9.1のGDALを入れてしまっていたりした場合は、libディレクトリ下でのシンボリックリンク張り直し等が適宜必要です。

mext版GDALをインストールすると、gdal_translateというデータ形式変換プログラムが、JPGIS形式を読み取れるようになります。
これを使って基盤地図情報を、GeoTiffというGISの標準データフォーマットを用いて、各画素に標高を表すデータが入っているDEMというデータに変換します。
走らせるコマンドは

gdal_translate FG-JPS-5136-36-dem10b-20090201.xml wgs_dem.tif

これだけでOKです。*6
ファイルの拡張子や、なかのフォーマットから、入力データ形式と出力データ形式はGDALが推論してくれます。*7

これを行うと、DEM形式のGeoTiffファイルができます。

が、単なる白黒ファイルに見えますね…。

でも大丈夫です。
QGISという、GISファイルを読めるソフトで、
レイヤ > ラスタレイヤを追加 > 作ったtiffを選択 > レイヤのプロパティでカラーマップを原色グレースケールから原色に変更、をすると、

(c)国土地理院
こんな感じでちゃんと地形データになってます。

基本、ちゃんとできるので大丈夫なんですが、たまにオプション等を変えたりしてると本当に失敗する時があり、それを確かめるためにも、この作業をするならQGISは入れておいた方がいいと思います。

球面メルカトル座標に変換

ここでできたDEMファイルは、WGS84世界測地系の経緯度そのまま)で記述されています。
これをTMSにしたいので、Googleメルカトルへの変換をここで行ってやります。

なお、Googleメルカトルへの変換を必ずしもここでやる必要はありません。
好みによって、段彩陰影図のGeoTiffまでWGS84で作ってやって、TMSにする際にgdal2tilesにメルカトル化を任せるという方向性もありだとは思います。
が、私としてはいろいろ試してみた結果、ある程度画像化してから座標系変換するとブロックノイズみたいなのがどうしても載ってしまうので、早いうちに変換した方がノイズは少ないかと思い、こうしています。

変換コマンドはこんな感じです。

gdalwarp -s_srs EPSG:4326 -t_srs EPSG:3785 -r bilinear wgs_dem.tif merc_dem0.tif
gdalwarp -s_srs EPSG:4326 -t_srs EPSG:3785 -r bilinear -srcnodata "" -dstnodata -9999 wgs_dem.tif merc_dem9.tif

2つコマンドを走らせていますが、データのない(=海面)領域のデータを0とするか、-9999でnodataフラグをつけるかの2パターンを作るために行っています。
なぜ2バージョン作っているかというと、後工程で海面領域が0のデータ(前者)と-9999のデータ(後者)の2通りが必要になるからです。
それぞれ別の名前をつけて、保存しておきます。

カラーレリーフ(段彩図)を作成する

段彩陰影図というのは、段彩図(高度によって高さを変えた地形図)と陰影図(地形に光を当てた時の影を描いた図)の2種を混ぜた図な訳ですが、これの前者の図を作成します。

変換コマンドは、前行程での-9999nodataの画像側を使って、

gdaldem color-relief -alpha -co ALPHA=YES merc_dem9.tif ramp.txt colorrelief.tif

というふうに行います。
が、ちょっと待ってください。
このコマンド、入力ファイルと出力ファイルの指定はいいのですが、合間に「ramp.txt」とかいう謎のファイルが挟まっています。
これはなんでしょうか?

実はこれ、高度と色分けの関係を定義するファイルです。
私は、こちらで公開されていたデータを参考に、

-9999 0 0 0 0
-50 0 0 205
0 0 191 191
0.1 57 151 105
100 117 194 93
200 230 230 128
500 202 158 75
1000 214 187 98
2000 185 154 100
3000 220 220 220
3800 255 255 255

という値を今のところ使っていますが、当然ながらこれの出来によって後の段彩陰影図の出来が変わってきてしまいます。
いろいろ変更して、納得のいくものを作るとよいでしょう。
また、-9999だけ4つの値が指定されていますが、この4つ目の値0が大切です。
これは透過フラグで、ここが0になっているピクセルは、タイル化の際に透過となるので、海領域を透過にしたい場合は、この指定が必須になります。

出来た画像のサンプルはこんな感じ。

(c)国土地理院

ヒルシェード(陰影図)を作成する

段彩図ができれば、今度は陰影図です。
こちらは前工程での海域0のファイルを使い、次のようにします。

gdaldem hillshade -compute_edges merc_dem0.tif hillshade.tif

海域が0のファイルを使う理由は、どうもこのgdaldem hillshadeというコマンドは透過情報を読み取ってくれないらしく、-9999nodataのデータを使うと海岸に断崖絶壁を作ってくれるため、それを防ぐために0にしています。
透過は、後でカラーレリーフとマージする事で、そっちから透過データを写し取ります。

なお、compute_edgesというオプションは、これを入れないと四辺隅に1ピクセル幅の枠が出来てしまうため、加えています。

また、上記コマンドは光源の高さや角度、方角等を制御していないため、デフォルトの設定で陰影図ができてしまいますが、詳細は省略しますがオプションで変更できます。
いろいろ試して、好みの陰影を作るとよいでしょう。*8

出来た画像のサンプルはこんな感じ。

(c)国土地理院

段彩図と陰影図をマージして、段彩陰影図を作る

マージするには、hsv_merge.pyというpythonスクリプトを使います。
これはgdalに標準では入ってないので、こちら等からダウンロードしてやります。
また、pythonはGDALモジュール拡張のインストール等が事前に必要なので、その準備はしておきます。

走らせるコマンドは次の通り。

./hsv_merge.py colorrelief.tif hillshade.tif colorshade.tif

これで、段彩図と陰影図がマージされ、段彩陰影図が作られます。
Google Maps等で使えるTMSにするならばまだ一工程必要ですが、QGIS等で用いるのであれば、ここまでで完了です。
お疲れ様でした!

(c)国土地理院

タイル化してGoogle Maps等で使える形にする

タイル化するには後もう一工程必要です。

gdal2tiles.pyというgdal同梱のスクリプトを使って、タイル化するのですが、このスクリプト、そのまま普通に実行すると、図郭全体に薄いグレーの1ピクセル枠ができてしまいます。
こんな枠があると、隣のメッシュと並べて表示した際に、邪魔になってしまいます。

これを消すには、antialiasというオプションをつけるだけで消えるという情報もありましたが、ちょっと確かめられていません。
私は、それより先にこちらの記事をみつけてしまい、それに従ってソースコード改変してしまっていたので、ソースコード改変無しでantialias指定だけで枠が消えるかは確認できません。
リンク先はQ&Aサイトでリンク内容が消えるかもしれないので、念のため変更内容を引用しておきます。

# Scaling by PIL (Python Imaging Library) - improved Lanczos
array = numpy.zeros((querysize, querysize, tilebands), numpy.uint8)
for i in range(tilebands):
array[:,:,i] = gdalarray.BandReadAsArray(dsquery.GetRasterBand(i+1), 0, 0, querysize, querysize)
im = Image.fromarray(array, 'RGBA') # Always four bands
im1 = im.resize((tilesize,tilesize), Image.BICUBIC)
if os.path.exists(tilefilename):
im0 = Image.open(tilefilename)
im1 = Image.composite(im1, im0, im1)
im1AsArray = numpy.array(im1)
alpha = im1AsArray[:,:,3]
semiTransparentIndices = alpha < 255
alpha[semiTransparentIndices] = 0
im1AsArray[:,:,3] = alpha
im1 = Image.fromarray(im1AsArray, 'RGBA')
im1.save(tilefilename,self.tiledriver)

この変更を、gdal2tiles.pyの中に同じコメント文で始まる1ブロックがありますので、置き換えてやります。

その後、下記のコマンドを実行してやります。
上記変更を施しても、antialiasオプションは必須です。

gdal2tiles.py -k -r antialias colorshade.tif

そうすると、タイル分割された画像とともに、それを表示するGoogle Maps APIOpenLayersを使ったサンプルのhtmlが吐かれます。
このhtmlが、いまだにGoogle Maps API v2だったりして意外にしょぼいのですが、まあタイルの確認程度には使えます。

(c)国土地理院, OpenStreetMap CCBYSA
というか、オープンソースなんだからしょぼいとか言うんなら自分で直せって話ですよね。はい。すみません。考えておきます。
また、-kオプションをつけているので、Google Earth で使えるkml等も一緒に吐かれます。
kmlいらなければ-kオプション不要ですね。

(c)国土地理院Google

そんな感じで

一応、段彩陰影図の作り方手順一通りはまとめられたようです。
が、できてから改めてみましたが、結構しょぼいです…。
東京地形地図さんあたりとは、5m DEMと10m DEMの違いを考慮しても、色使いのダイナミックさや美しさで勝負になりません。
もし本当にサービスにするなら、もっと美麗に仕上げられるよう、ramp.txtや陰影図オプションの試行錯誤が必要そうです。

また、レーザー計測の5m DEMと違い、等高線から起こしている10mDEMは、海以外の内水面、川や湖等はnodata領域扱いではないので、内水面も陸地扱いになってしまいます。*9
この辺をしっかりやろうと思うと、基盤地図情報から水涯線のデータ等を落としてきて、ポリゴン化してマスク判定とか、そんな手順まで必要になってきてしまいます。
そこまでできるか判りませんが、とりあえず今現在、段彩陰影図作成について知っている事をまとめてみました。

*1:例えば東京地形地図さんなど

*2:すみません、諸事情ありその後ほとんど進んでいません…

*3:申し訳ありません、本記事いろいろなところで収集した情報の寄せ集め+自己経験+知人に教えていただいた事で作ってるのですが、収集した情報ソースの整理がまだできておりません。取り急ぎリファー不足で公開しますが、追って確認でき次第リファーを追加していきます。

*4:東京や大阪でやりたかったのですが、いい感じで陸と海が混ざってるメッシュがなかったので…この地域の処理はやってないので、ブログ書きながらのぶっつけ本番並行処理になります。ドキドキです。

*5:すみません、間違ってたら教えてください

*6:当然、入力ファイル名、出力ファイル名は適宜変えてください

*7:もちろん明示もできますが、詳しく触れる事は省略します

*8:また、メルカトル図法は座標の単位がmですが、それは赤道上での一致で、実際には緯度が高くなる程、座標系上の距離は実際の距離より長くなります。なので、DEMでの陰影は、係数を掛けないと緯度に従い実際の高度より低い感じで描画される事になります。これを嫌う場合は、各メッシュ毎の緯度に従った係数を高度スケールに加算してやるべきですが、個人的には、高度情報を得るのは段彩の方であり陰影の方は地形の雰囲気を掴むためのものなので、そこまで正確な高度である必要はないと思い、無係数で作成しています。

*9:レーザー計測では、水面では計測できないため、内水面もデータ無しになるらしい

© Code for History