サイトマップ
お知らせ、メモ
案内板
うちのヘッドライン
 




トップ  >  基盤地図情報をBlenderに取り込んで3D表示

基盤地図情報をBlenderに取り込んで3D表示

この記事は FOSS4G Advent Calendar 2012 の12月11日の記事です。 今年は論文が投稿されたりして昨年にも増してレベルが高いので、 筆者がレベルを下げて調整しようとしましたが、さらに追い討ちをかけるかのように複数の方向からネタつぶしの妨害が、、、。 というわけで、高いハードルの下をくぐるかのような、いまさらな内容になってしまいました。

単なるデータ変換とそれを利用した絵の作り方という内容で、FOSS4GどころかGISの話ですらないかもしれません。。。 まあ、、中だるみということでご勘弁いただければと思います。

中禅寺湖から男体山を望む絵のつもり
中禅寺湖から男体山を望む絵のつもり

基盤地図情報がGDAL/OGRで取り扱えるようになったわけですが、それによってホイホイと取り扱いやすいフォーマットに変換して色々と加工することが容易になりました。 これを利用してBlenderに取り込んでCG的に表示してみました。 この手法はVisualising QGIS data with Blenderとよく似ていますが、 上記のリンクがDEM画像をモディファイアに適用してオブジェクトに高さを付与するのに対して、こちらはBlender Pythonスクリプトで直接取り込んでいます。

もっとも、AutoCAD Map 3Dを使えばこちらのようにもっとしっかりした絵が作れるとは思いますが、、、 さらに言えばカシミール 3Dでいいのではと思われますが、 まあ一応CGソフトに取り込む利点としては任意のオブジェクトを追加することができる点が挙げられるかと思います。

以下では、この絵ができるまでの工程を説明します。

DEMを取り込む

まずは取り込むDEMを作成します。 ここではDEMをBlender Pythonスクリプトで取り扱いやすい(めんどくさかった)のでArcInfo Ascii Grid形式に変換します。

# とりあえず必要なDEMデータをまとめてGeoTIFFに変換する
> for %i in (*.xml) do gdal_translate -of GTiff -ot Float32 %i %i.tif 

# マージとクリップ
> gdal_merge -o merge.tif -of GTiff -ot Float32 -ul_lr 139.4136 36.7889 139.5322 36.7047 "ファイル1" "ファイル2"

# 座標変換
# Blender等CGソフトに取り込む場合は経緯度では都合が悪い
# ここでは平面直角座標9系にする
# また、gdalwarpコマンドではArcInfo Ascii Grid形式には出力できない
> gdalwarp -r bilinear -s_srs EPSG:4612 -t_srs EPSG:2451 merge.tif merge_w.tif

# 最後にArcInfo Ascii Grid形式に変換
> gdal_translate -of "AAIGrid" -ot Float32 merge_w.tif merge.txt

今回使用したデータでは元となるデータにNodataが無いのでこれで大丈夫ですが、Nodataがある場合は対策が必要です。 Nodataがある場合、基盤地図情報ドライバではデータ無しは-9999.0が割り当てられていますが、 この値をそのままCGソフトに取り込んでしまうわけにはいかないので、 この後の取り込みスクリプトか上記のコマンドのどこかでNodataを0.0ぐらいに値を変えておく必要があります。 上記のコマンドで値を変えるとすればgdalwarpあたりで-srcnodata -dstnodata オプションを指定すればいいのではないでしょうか。

ArcInfo Ascii Gridは簡単なテキストフォーマットなので、取り扱いも簡単です。 Blenderに取り込むスクリプトは以下のとおりです。 CGソフトはGISソフトではないので(当たり前や)、編集画面の視界の中心座標値は0中心で動かせません。 また、Blenderは編集画面の視界範囲が狭い座標範囲で固定されています。 そこで、DEMの座標をスケール/オフセットして取り込む必要があります。 以下のスクリプトでは、Xの範囲が-5.0~5.0になるように調節しています。

Blender Pytho APIは2.5から大幅に変更しました。なので筆者の豆知識は全てパーです。 以下のスクリプトは以前公開していたメッシュ取り込みスクリプトの新API版のものです。

エラー処理などは省略していますので注意!

   1: import bpy
   2: import bmesh
   3: import bpy_extras
   4: 
   5: list_h = []
   6: res = 0.0
   7: rows = 0
   8: cols = 0
   9: min_h = 10000.0
  10: scale = 1.0
  11: 
  12: def readMesh(filename):
  13:     
  14:     i = 0
  15:     j = 0
  16:     c = 0
  17:     nodata = 0
  18:     
  19:     global list_h, res, rows, cols, scale, min_h
  20:     
  21:     for lines in open(filename, 'r'):
  22:         list_c = lines[:-1].split()
  23: 
  24:         if list_c[0] == "ncols":
  25:             cols = int(list_c[1])
  26:         elif list_c[0] == "nrows":
  27:             rows = int(list_c[1])            
  28:         elif list_c[0] == "xllcorner":
  29:             west = float(list_c[1])
  30:         elif list_c[0] == "yllcorner":
  31:             south = float(list_c[1])
  32:         elif list_c[0] == "cellsize":
  33:             res = float(list_c[1])
  34:         elif list_c[0] == "NODATA_value":
  35:             nodata = int(list_c[1])
  36:         else:
  37:             for c in list_c:
  38:                 list_h.append(float(c))
  39:     
  40:     wx = res*cols
  41:     wy = res*rows
  42:     scale = 10.0 / wx
  43:     res = res * scale
  44:     print(rows, cols, west, south, res)
  45: 
  46:     for h in list_h:
  47:         if h < min_h:
  48:             min_h = h
  49: 
  50: class ImportOperator(bpy.types.Operator):
  51:     bl_idname = "object.import_aaigrid"
  52:     bl_label = "Import ArcInfo ASCII Grid"
  53:     
  54:     global res, list_h
  55:     
  56:     filepath = bpy.props.StringProperty( subtype="FILE_PATH" )
  57: 
  58:     @classmethod
  59:     def poll(cls, context):
  60:         pass # 間違えを指摘して頂きました
  60          return True # ここは本来コンテキストを調べてプラグインが起動可能かどうかを判断するべき
  61: 
  62:     def execute(self, context):
  63:         readMesh(self.filepath)
  64:         mesh = bpy.data.meshes.new("Terrain")
  65: 
  66:         u = 0
  67:         v = 0
  68:         verts = []
  69:         faces = []
  70:         
  71:         for h in list_h:
  72:             vert = [-5.0+res*u, 5.0-res*v, (h - min_h)*scale]
  73:             verts.append(vert)
  74:             
  75:             u = u+1
  76:             if u == cols:
  77:                 v = v+1
  78:                 u = 0
  79:         
  80:         for y in range(rows-1):
  81:             for x in range(cols-1):
  82:                 face = [y*cols+x, y*cols+x+1, (y+1)*cols+x+1, (y+1)*cols+x]
  83:                 faces.append(face)
  84:         
  85:         mesh.from_pydata(verts, [], faces)
  86:         ob = bpy.data.objects.new('Terrain', mesh)
  87:         scn = bpy.context.scene
  88:         scn.objects.link(ob)
  89:         
  90:         return {'FINISHED'}
  91: 
  92:     def invoke(self, context, event):
  93:         wm = context.window_manager
  94:         wm.fileselect_add(self)
  95:         
  96:         return {'RUNNING_MODAL'}
  97: 
  98: def register():
  99:     bpy.utils.register_class(ImportOperator)
 100: 
 101: 
 102: def unregister():
 103:     bpy.utils.unregister_class(ImportOperator)
 104: 
 105: 
 106: if __name__ == "__main__":
 107:     register()
 108:     bpy.ops.object.import_aaigrid('INVOKE_DEFAULT')

実行すると取り込むファイルを選択する画面が表示され、ファイルを選択すると以下のように取り込まれます。

実行結果
実行結果

点数が多過ぎて少し間引きしたい場合は、上記のgdal_translateコマンドで-outsizeオプションを指定するとよいでしょう。

Blender上で色々と加工

ここまでくれば後は好きなようにオブジェクトを加工するといいでしょう。 まずDEMのテクスチャを以下のように設定しました。

DEMテクスチャ設定
ここでは簡単に遠景の木を表現するために、プロシージャルテクスチャのボロノイを適用
サイズはNoiseのサイズかMappingのサイズで調整

DEMマテリアル設定
マテリアル設定
少しだけ鏡面反射成分を残しておいたほうがよさそう

ひとつ注意点ですが、今回使用した中禅寺湖のように内水面がある場合、基盤地図情報にも内水面の部分に点があります。 水面を設定する上ではこの点は邪魔なのですが、消すのもめんどくさい場合は水面に使用する平面をひとつ追加して、 基盤地図情報の水面よりもわずかに高いところに配置するといいでしょう。

以下は水面のテクスチャ/マテリアル設定です。ちゃんと波を表現するにはWaveモディファイアを適用しますが、 ここではプロシージャルテクスチャで法線ベクトルを設定して波を表現しています。

水面テクスチャ設定1
Cloudテクスチャを作成、サイズで波の大きさを設定

水面テクスチャ設定2
Geometryの"Normal"にチェックを入れて、Cloudテクスチャを法線ベクトルの値に設定

水面マテリアル設定1
マテリアル設定
Diffuseを0にして鏡面反射のみに設定

水面マテリアル設定2
"Mirror"にチェックを入れ、反射強度を調節して水らしさを表したつもり


降雪バージョン
こちらは降雪バージョン
実際の木を植えていないため近景が残念なことに、、、

高さで色分けする

地図帳なんかでよく見かける、低いところが緑で高いところが茶色のような絵にしたい場合は、 DEM画像を利用してマテリアルを設定するといいでしょう。 DEM画像は上記の変換の過程で副産物としてできているので、それを利用します。

# DEM画像を8bit Bitmapに変換する
# -scaleオプションで指定する高さの範囲は色を割り当てる範囲に応じて設定する
# -scaleオプションの最大値をDEMの最大値にすると山頂付近だけが茶色の絵になってしまうので、
# 少し低いところを最大値にするとよさそう
> gdal_translate -of BMP -ot Byte -scale 1000 2000 merge_w.tif merge_8bit.tif

できた8bitグレースケール画像をGIMP等で編集してテクスチャ画像を作成してもいいですが、ここではBlenderらしく以下のように設定してみます。

ちなみに、GIMPで緑→茶画像を作るとしたら、まず背景を緑で塗りつぶした画像を作成し、 上記で作成した8bitグレースケールDEM画像をマスクレイヤとして割り当てて、前景色を茶色にして塗りつぶし、とかですかね? そちらの方法であれば、茶色の割り当て具合はマスクレイヤーの濃度で調節できそうです。

ノードエディタ
ノードエディタでマテリアルを設定
高所、低所用のマテリアルをそれぞれ用意して、上記で作成した8bit高さ画像を混合比率として使用する


地図帳バージョン
雰囲気に合わせて水面なども調節するといいでしょう
テキスト要素に日本語を使えるといいんですけどね、、、
あと細かいことですが、この絵ではテキスト要素は影無しに設定していますが、あってもまあいいかもしれません。

ベクタデータをオーバーレイする

Visualising QGIS data with Blenderにもありますが、 同じ区域のオルソ画像があればそれをテクスチャ画像として貼り付けることは簡単にできます。 ですが筆者の手元には無いので、代わりに基盤地図情報のベクタデータをテクスチャ画像に変換して貼り付けてみます。

本来はDEM画像と貼り付ける画像の座標が完全に一致していなくても、スクリプトで計算してUVマッピングすればいいのですが、 残念ながら筆者はまだそのスクリプトを作れないので、以下の例では完全に一致した画像を作成してはりつけます。

基盤地図情報対応GDAL/OGRであれば、入力ファイルにxmlファイルを指定することができます。

# まずはベクタデータをDEM画像の範囲でクリップ
# 複数の項目がある場合はそれぞれ以下のコマンドを実行する
> ogr2ogr -f "ESRI Shapefile" -clipsrc 139.4136 36.7889 139.5322 36.7047 "出力ファイル".shp "入力ファイル"

# ラスタライズ、複数ファイルを1つのファイルに書き込むことができる。
# 最初の1ファイルだけ-tr -init -teオプションを付ける
# -initオプションは省略すると背景が黒になる
# -teオプションは省略すると画像範囲はベクタデータを包含する最小矩形になってしまうため、今回の例では必ず指定する必要がある
# 以下の例では背景を薄緑色に、ベクタデータを橙色に設定している
> gdal_rasterize -of "GTiff" -ot Byte -init 200 255 200 -burn 255 -burn 192 -burn 0 
  -tr 0.000111 0.000111 -te 139.4136 36.7889 139.5322 36.7047 -l clip_rd1 clip_rd1.shp tex.tif

# 2枚目以降は上記のファイルに書き込んでいく。-bオプションで書き込むバンドを指定する
# 以下の例では書き込む色を青に指定している
> gdal_rasterize -b 1 -b 2 -b 3 -burn 0 -burn 0 -burn 255 -l clip_wl clip_wl.shp tex.tif

# 最後にワープ
> gdalwarp -s_srs EPSG:4612 -t_srs EPSG:2451 tex.tif tex_w.tif

あとは、できたテクスチャファイルをBlenderに取り込んだDEMオブジェクトに貼り付けます。 テクスチャマッピングはデフォルトでは"Generated"ですが、範囲が完全に一致していればこのままでぴったり重なります。

ベクタオーバーレイバージョン
視点を東側のいろは坂、華厳の滝方面からに設定
これもいまいち残念なできばえ、、、
一応位置はあっているのがわかると思います

終わりに

筆者のCGソフトの腕前がいまひとつなので、できた絵もいまひとつですが、上手に使えればきれいな絵ができるでしょう。 今回は間に合いませんでしたが、CGソフトに取り込んでしまえば、木を植えたり建物を建てたり、トンネルを掘ったり橋を架けたり、隕石をぶち込んだりなど、いろいろなことができるようになります。 Blenderは想像以上に恐ろしいソフトです。極めれば大概のことはできるでしょう。


クリエイティブ・コモンズ・ライセンス
This documents by Yamate,N is licensed under a Creative Commons 表示 - 継承 3.0 非移植 License.
login