もう少し広範囲の３次元メッシュ（xml）

　国土地理院の10mメッシュを使った八丈島(はちじょうじま)と八丈小島(はちじょうこじま)です。

　人間の欲望は、きりがありません。１区画の３次元メッシュでは物足りないので、もう少し広い範囲を作ります。

　xmlファイルが複数なので、次のように処理することにします。

　作業用のフォルダに10mメッシュか5mメッシュのどちらか統一したxmlファイルを入れておきます。xmlファイルの先頭がFG-GML-で始まるxmlファイルを処理対象とします。

　今回は、xmlファイルはどうせプログラムの自動出力で中身は間違いないだろうということで、項目チェックを簡素化しています。さっそく、BlenderのPythonスクリプトを見てみましょう。ここに生pyファイルがあります。

import bpy import math from xml.etree.ElementTree import ElementTree import glob # # もう少し広範囲の３次元メッシュ（xml） # demmulti.py # 2016/11/17 Shizuka-Toikanbetsu # XmlFile = "FG-GML-*.xml" # # 既存カーブの削除 for item in bpy.context.scene.objects: if item.type == 'MESH': bpy.context.scene.objects.unlink(item) for item in bpy.data.objects: if item.type == 'MESH': bpy.data.objects.remove(item) for item in bpy.data.meshes: bpy.data.meshes.remove(item) # # namespaceを定義 ns = { 'xml': 'http://fgd.gsi.go.jp/spec/2008/FGD_GMLSchema', 'gml': 'http://www.opengis.net/gml/3.2' } # fact = 36000.0 # # 世界座標の基準緯度と世界座標の目標座標の距離をメートル単位で返す def GetLatOffset(pBaseLat, pLat): return ( pLat - pBaseLat ) * 40007000.0 / 360.0 # # 世界座標の基準経度と世界座標の目標座標の距離をメートル単位で返す def GetLonOffset(pBaseLon, pLat, pLon): return ( pLon - pBaseLon ) * 40075000.0 * math.cos( math.radians( pLat ) ) / 360.0 # class CMinMax: def __init__(self): self.Xmin = 100000000.0 self.Ymin = self.Xmin self.Xmax = - self.Xmin self.Ymax = self.Xmax self.px = 0 self.py = 0 self.Xspan = 0 self.Yspan = 0 # # ２次メッシュ、３次メッシュ専用緯度経度補正関数 def Adjust2(Item): return round( Item * fact ) # def GetMapRange(pXml): print(pXml) NowRange = CMinMax() tree = ElementTree() tree.parse( pXml ) xpath = ".//gml:lowerCorner" lowerCorner = tree.find( xpath, ns ) tempstr = lowerCorner.text item = tempstr.split( ' ' ) NowRange.Ymin = Adjust2( float( item[ 0 ] ) ) NowRange.Xmin = Adjust2( float( item[ 1 ] ) ) xpath = ".//gml:upperCorner" upperCorner = tree.find( xpath, ns ) tempstr = upperCorner.text item = tempstr.split( ' ' ) NowRange.Ymax = Adjust2( float( item[ 0 ] ) ) NowRange.Xmax = Adjust2( float( item[ 1 ] ) ) xpath = ".//gml:high" high = tree.find( xpath, ns ) tempstr = high.text item = tempstr.split( ' ' ) NowRange.px = int( item[ 0 ] ) + 1 NowRange.py = int( item[ 1 ] ) + 1 NowRange.Xspan = ( NowRange.Xmax - NowRange.Xmin ) / NowRange.px NowRange.Yspan = ( NowRange.Ymax - NowRange.Ymin ) / NowRange.py return NowRange # print( "\nStart Script" ) # destminmax = CMinMax() # get xml files files = glob.glob( XmlFile ) if len( files ): print( 'Now get range from all xml.' ) for file in files: NowMapRange = GetMapRange( file ) destminmax.Xmin = min( NowMapRange.Xmin, destminmax.Xmin ) destminmax.Xmax = max( NowMapRange.Xmax, destminmax.Xmax ) destminmax.Ymin = min( NowMapRange.Ymin, destminmax.Ymin ) destminmax.Ymax = max( NowMapRange.Ymax, destminmax.Ymax ) destminmax.px = round( ( destminmax.Xmax - destminmax.Xmin ) / NowMapRange.Xspan ) destminmax.py = round( ( destminmax.Ymax - destminmax.Ymin ) / NowMapRange.Yspan ) wx = ( destminmax.Xmin + destminmax.Xmax ) / 2.0 / fact wy = ( destminmax.Ymin + destminmax.Ymax ) / 2.0 / fact BlockWidth = ( NowMapRange.Xmax - NowMapRange.Xmin ) BlockHeight = ( NowMapRange.Ymax - NowMapRange.Ymin ) BlockXspan = int( BlockWidth / NowMapRange.px ) BlockYspan = int( BlockHeight / NowMapRange.py ) # bpy.ops.mesh.primitive_grid_add( x_subdivisions = destminmax.px, y_subdivisions = destminmax.py, location = ( 0, 0, 0 ) ) Grid = bpy.data.objects[ 'Grid' ] Grid.scale[0] = 0.001 Grid.scale[1] = 0.001 Grid.scale[2] = 0.001 bpy.context.scene.objects.active = Grid for file in files: print( file ) tree = ElementTree() tree.parse( file ) # NowRange = CMinMax() xpath = ".//gml:lowerCorner" lowerCorner = tree.find( xpath, ns ) tempstr = lowerCorner.text item = tempstr.split( ' ' ) NowRange.Ymin = Adjust2( float( item[ 0 ] ) ) NowRange.Xmin = Adjust2( float( item[ 1 ] ) ) xpath = ".//gml:upperCorner" upperCorner = tree.find( xpath, ns ) tempstr = upperCorner.text item = tempstr.split( ' ' ) NowRange.Ymax = Adjust2( float( item[ 0 ] ) ) NowRange.Xmax = Adjust2( float( item[ 1 ] ) ) xpath = ".//gml:tupleList" temp = tree.find( xpath, ns ) tupleListStr = temp.text xpath = ".//gml:startPoint" startPoint = tree.find( xpath, ns ) tempstr = startPoint.text item = tempstr.split( ' ' ) sx = int( item[ 0 ] ) sy = int( item[ 1 ] ) # # 標高データの文字列を分解し、実数型の準備データを作成 AltSrc = [] AltList = tupleListStr.split( '\n' ) for AltLine in AltList: if AltLine: # 空行でなければ以下を実行 LineData = AltLine.split( ',' ) if LineData: AltSrc.append( float( LineData[ 1 ] ) ) # # 準備データを本番データに転記 Alt = [] for destptr in range( 0, NowMapRange.px * NowMapRange.py ): Alt.append( -9999.0 ) srcptr = 0 destptr = sy * NowMapRange.px + sx for ptr in range( 0, len( AltSrc ) ): Alt[ destptr ] = AltSrc[ srcptr ] destptr = destptr + 1 srcptr = srcptr + 1 Xblock = int( ( NowRange.Xmin - destminmax.Xmin ) / BlockWidth ) Yblock = int( ( NowRange.Ymin - destminmax.Ymin ) / BlockHeight ) AltPtr = 0 for Yptr in range( 0, NowMapRange.py ): Ybase = ( NowMapRange.py - 1 ) - Yptr Ybase = Yblock * NowMapRange.py + Ybase for Xptr in range( 0, NowMapRange.px ): Xbase = Xblock * NowMapRange.px + Xptr destindex = Ybase * destminmax.px + Xbase Xpos = ( Xbase * BlockXspan + destminmax.Xmin ) / 36000.0 Ypos = ( Ybase * BlockYspan + destminmax.Ymin ) / 36000.0 vert = bpy.context.object.data.vertices[ destindex ] vert.co[ 0 ] = GetLonOffset( wx, Ypos, Xpos ) vert.co[ 1 ] = GetLatOffset( wy, Ypos ) vert.co[ 2 ] = Alt[ AltPtr ] AltPtr = AltPtr + 1 # # 未定義データー（標高値 = -9999.0）があればをまとめて削除 print( 'Now delete undefine data.' ) bpy.ops.object.mode_set(mode='EDIT') bpy.ops.mesh.select_all(action='DESELECT') bpy.ops.object.mode_set(mode='OBJECT') for v in bpy.context.object.data.vertices: if v.co[2] == -9999.0: v.select = True bpy.ops.object.mode_set(mode='EDIT') bpy.ops.mesh.delete(type='VERT') bpy.ops.object.mode_set(mode='OBJECT') # print( "End Script" )

　広範囲の場合は、開始前に一度全部のxmlを読み、範囲を取り出す必要があるので処理時間も長くなります。座標点は300万以上あるので待つしかありません。すべてのxmlをまるごとメモリに読むのではなく、１つずつ処理するようにしています。

　処理中は数ギガバイトの領域を確保しているようです。Windows32ビット版では、やめたほうが良いと思います。

　座標点が多いため、10mメッシュでも４区画程度にしないと色々と処理が重くなります。レンダリングはかなり時間がかかります。

　少しでも精度をよくするため、座標点の緯度と経度は原則的に整数型で保持し、実際に必要になった時点で本物の実数型にしています。36000という値がところどころにあるのはそのためです。

　海に見立てた濃い青色の大きな平面オブジェクトを配置し、ワールドの背景色を空の色にすると次のようになります。水平線が本当に水平なので、少しウソくさいです。

　サンプルではオブジェクトが原点ゼロになるよう移動しています。オブジェクトを移動する方法については、前節を参照してください。

　立体模型として眺めるだけでなく、ちょっとした建物を自分で作ると箱庭として遊べます。

　5mメッシュの三宅島の三宅高校がある平地部分です。白い箱を置くと、高校ができそうです。5mメッシュでは、水の部分に穴が開いているので、必要に応じて埋めるなり、水を張る必要があります。

　立体模型として眺めるだけでなく、ちょっとした建物を自分で作ると箱庭として遊べます。

　一部の地域では5mメッシュも配信されています。表現力は増しますが、データー量が多くなります。東京都の利島(としま)で、港の防波堤から振り返ると次のようになります。

10mメッシュ	5mメッシュ

　立体模型として眺めるだけでなく、ちょっとした建物を自分で作ると箱庭として遊べます。にしても、5mメッシュの威力すごいこと。レンダリングが重いです。三宅島の頂上付近です。色を塗っていないので雪山に見えます。

　行ったことのない土地でも作れるのも魅力です。

　世界遺産の母島で、左の画像は南から、右の画像は乳房山の頂上から南を向いています。左の画像で手前にある３つの島は左からスリムな姉島、太った妹島、姪島です。

　最後は、雑誌か何かのアンケートで日本人が行きたくない温泉でダントツ１位だったという大分県の湯布院。韓国語と中国語しか聞こえてこないそうです。

　本格的な地図の処理だとtifファイルを使いますが、Blenderに取り込む程度であればxmlをそのまま読めば余計な手間もかからず、簡単ですね。

<トップページに戻る>