|
国土地理院 基盤地図情報の利用 |
１区画の３次元メッシュ（xml）
今度はとりあえず、３次元メッシュの地図をBlenderのPythonで作ってみます。国土地理院の数値標高データーを利用します。上記の画像は東京都の
他サイトやブログでは、GDALが出てきてxmlからtif形式に変換して、何かのソフトを通して、また変換してと、随分と回りくどいことをやっています。
BlenderではPythonでそのままxmlのデーターをメッシュにすればよいだけです。
また、XY座標の配置ルールも自分で好きに定義できます。
なお、ここで紹介する方法は、単純に１枚のメッシュオブジェクトに高さを持たせるもので、簡単ゆえにオブジェクトサイズが小さい利点があります。ただし、3Dプリンターで使われるような台が付いている立体をお望みのときはメッシュに高さを加えるか、ボックスを用意し上面だけをdivideしながら作る必要があります。当然オブジェクトサイズが大きくなります。
この画像は北海道の
つまり、３次元メッシュ上に建物などのオブジェクトを作りこむと、箱庭を作ることもできます。
国土地理院の基盤地図情報はxmlで配信されています。
この例では東京都の鳥島をつくってみます。国土地理院のサイトで数値標高モデルから10mメッシュを選び、地図の中から選択するよう指定した後、伊豆諸島を南下し鳥島が収録されている454052というエリアを探し、ダウンロードします。なお、初めての方は登録（無料）が必要になります。
zipファイルを解凍し、ファイル名がFGで始まり、途中にdemという３文字が入っているxmlにデーターが収録されています。なお、10mメッシュでは、火山のデーターもあわせて提供されている地域があり、場合によっては同一地域で年度違いとかタイプ違いで３枚、４枚と収録されていることがあります。demの b というタイプが通常提供のものになります。a タイプは火山です。
Pythonのスクリプトは次のようになります。ここに生pyファイルがあります。
import bpy import math import xml.etree.ElementTree as ET # # 国土地理院 数値標高モデルのxmlから１区画の３次元メッシュ（xml） # demsingle.py # 2016/11/14 Shizuka-Toikanbetsu # demtree = ET.parse('FG-GML-4540-52-dem10b-20161001.xml') # ここに xml ファイルを指定 # # 世界座標の基準緯度と世界座標の目標座標の距離をメートル単位で返す 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 # # ２次メッシュ、３次メッシュ専用緯度経度補正関数 def Adjust2(Item): return round( Item * 36000.0 ) / 36000.0 # print( '\nStart Script.' ) # # 既存メッシュの削除 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) # demroot = demtree.getroot() # # namespaceを定義 ns = { 'xml': 'http://fgd.gsi.go.jp/spec/2008/FGD_GMLSchema', 'gml': 'http://www.opengis.net/gml/3.2' } demdem = demroot.find( 'xml:DEM', ns ) # メッシュコードを得る。ここから本格的 for demmesh in demdem.findall( 'xml:mesh', ns ): meshstr = demmesh.text # オブジェクト名に使うので文字列型のままで mesh = int( meshstr ) # こちらは整数型（実は未使用なので不要） for demcoverage in demdem.findall( 'xml:coverage', ns ): for demboundedBy in demcoverage.findall( 'gml:boundedBy', ns ): for demEnvelope in demboundedBy.findall( 'gml:Envelope', ns ): for demlowerCorner in demEnvelope.findall( 'gml:lowerCorner', ns ): lowerCornerStr = demlowerCorner.text temp = lowerCornerStr.split( ' ' ) Ymin = Adjust2( float( temp[ 0 ] ) ) # 南端の緯度 Xmin = Adjust2( float( temp[ 1 ] ) ) # 西端の経度 for demupperCorner in demEnvelope.findall( 'gml:upperCorner', ns ): upperCornerStr = demupperCorner.text temp = upperCornerStr.split( ' ' ) Ymax = Adjust2( float( temp[ 0 ] ) ) # 北端の緯度 Xmax = Adjust2( float( temp[ 1 ] ) ) # 東端の経度 for demgridDomain in demcoverage.findall( 'gml:gridDomain', ns ): for demGrid in demgridDomain.findall( 'gml:Grid', ns ): for demlimits in demGrid.findall( 'gml:limits', ns ): for demGridEnvelope in demlimits.findall( 'gml:GridEnvelope', ns ): for demhigh in demGridEnvelope.findall( 'gml:high', ns ): demhighStr = demhigh.text temp = demhighStr.split( ' ' ) px = int( temp[ 0 ] ) + 1 # 横の座標数 py = int( temp[ 1 ] ) + 1 # 縦の座標数 for demrangeSet in demcoverage.findall( 'gml:rangeSet', ns ): for demDataBlock in demrangeSet.findall( 'gml:DataBlock', ns ): for demtupleList in demDataBlock.findall( 'gml:tupleList', ns ): tupleListStr = demtupleList.text for demcoverageFunction in demcoverage.findall( 'gml:coverageFunction', ns ): for demGridFunction in demcoverageFunction.findall( 'gml:GridFunction', ns ): for demstartPoint in demGridFunction.findall( 'gml:startPoint', ns ): demstartPointStr = demstartPoint.text temp = demstartPointStr.split( ' ' ) sx = int( temp[ 0 ] ) # 開始点の横位置 sy = int( temp[ 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, px * py ): Alt.append( -9999.0 ) srcptr = 0 destptr = sy * px + sx for ptr in range( 0, len( AltSrc ) ): Alt[ destptr ] = AltSrc[ srcptr ] destptr = destptr + 1 srcptr = srcptr + 1 # # メッシュオブジェクトを作って頂点の配置と標高の書き込み wx = ( Xmin + Xmax ) / 2.0 wy = ( Ymin + Ymax ) / 2.0 span = ( Xmax - Xmin ) / px bpy.ops.mesh.primitive_grid_add(x_subdivisions=px, y_subdivisions=py, radius=1, location=(0, 0, 0)) Grid = bpy.data.objects['Grid'] Grid.name = meshstr Grid.scale[0] = 0.001 Grid.scale[1] = 0.001 Grid.scale[2] = 0.001 bpy.context.scene.objects.active = Grid dpx = 0 dpy = py - 1 for PointAlt in Alt: xpos = dpx * span + Xmin + span / 2.0 ypos = dpy * span + Ymin - span / 2.0 dx = GetLonOffset( wx, ypos, xpos ) dy = GetLatOffset( wy, ypos ) bpy.context.object.data.vertices[ dpy * px + dpx ].co[ 0 ] = dx bpy.context.object.data.vertices[ dpy * px + dpx ].co[ 1 ] = dy bpy.context.object.data.vertices[ dpy * px + dpx ].co[ 2 ] = PointAlt dpx = dpx + 1 if dpx >= px: dpx = 0 dpy = dpy - 1 # # 未定義データー（標高値 = -9999.0）があればをまとめて削除 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( 'Exit' )
３次元の鳥島が現れます。
Pythonの処理内容ですが、値の読み込み位置を少しまじめにチェックしているので項目が長くなっています。
標高データの実体の文字列の扱いですが、値だけの配列を準備配列として作成し、次に本番の配列を-9999.0で埋め、その後準備配列から該当位置に値を転記します。値の文字列値は配列の途中区間だけという図面もあるのでこのような仕様になります。
本番用の配列ができたらメッシュオブジェクトを作り、座標の再配置と標高の設定をしますが、メッシュオブジェクトの配列の順序は一般的な画像データーと同じで左下が原点、国土地理院の配列の順序は左上が原点と異なっているので注意が必要です。横方向は左から右と同じですが、縦方向は上下を入れ替えないとメッシュのポリゴンが反転します。
最後に未定義データー、標高値が-9999.0のものは不要なのでまとめて削除します。bmeshを使うともっと簡潔に記述できますが、余計な取り込みをしたくないのでbpyの範囲だけで記述しています。
ところで、鳥島が中央に表示されていません。これは、図面全体の中心を原点にしているためです。音威子府村のように全域が地面だと問題ありませんが、鳥島は１枚の図面の左上のほうに存在するのでこのようになります。
この例では鳥島以外に陸地はないので、鳥島を中央に移動することにします。オブジェクトを選択した状態であれば次の処理で移動できます。
print( '\nStart Script.' ) # Xmin = 100000000.0 Ymin = Xmin Xmax = -Xmin Ymax = -Ymin for v in bpy.context.object.data.vertices: Xmin = min( Xmin, v.co[ 0 ] ) Ymin = min( Ymin, v.co[ 1 ] ) Xmax = max( Xmax, v.co[ 0 ] ) Ymax = max( Ymax, v.co[ 1 ] ) Xoff = ( Xmin + Xmax ) / 2.0 Yoff = ( Ymin + Ymax ) / 2.0 for v in bpy.context.object.data.vertices: v.co[ 0 ] = v.co[ 0 ] - Xoff v.co[ 1 ] = v.co[ 1 ] - Yoff # print( 'Exit' )
選択に関係なく、オブジェクト名で指定するときは、次の処理で移動できます。
print( '\nStart Script.' ) # print( '\nStart Script.' ) # Grid = bpy.data.objects['454052'] # ここにオブジェクト名 # Xmin = 100000000.0 Ymin = Xmin Xmax = -Xmin Ymax = -Ymin for v in Grid.data.vertices: Xmin = min( Xmin, v.co[ 0 ] ) Ymin = min( Ymin, v.co[ 1 ] ) Xmax = max( Xmax, v.co[ 0 ] ) Ymax = max( Ymax, v.co[ 1 ] ) Xoff = ( Xmin + Xmax ) / 2.0 Yoff = ( Ymin + Ymax ) / 2.0 for v in Grid.data.vertices: v.co[ 0 ] = v.co[ 0 ] - Xoff v.co[ 1 ] = v.co[ 1 ] - Yoff # print( 'Exit' )
これで中央に配置されてすっきりです。
標高データを作るときは座標点の位置に注意が必要です。例えばこの図面だと左上は経度が140.25、緯度が30.5ですが、その位置は標高点の座標ではありません。その位置を左上とする 小さなマス目の中心の位置が実際の図形の座標点 になります。
10mメッシュでは１区画で1125×750で84万の座標点がありますから、データーは大きくなります。5mメッシュだとなおさらです。
このスクリプトは１枚の図面専用ですから、複数図面には使えません。例えば音威子府村の隣の図面も用意しようと別々に作ると、接続部分のメッシュは無いので、その部分だけ空間になります。手作業で面を張るには数が多すぎます。
BlenderのPythonで複数の図面に対応するためには、全体の大きさになる大きなメッシュを１つ作り、例えば左と右の２枚だったら１回目の作業で左半分の値を設定し、２回目で右半分を設定するという分割作業になります。分割作業でば、佐渡島のような大きなものも一応作れます。大きいと何をやるにも処理が重くて大変です。
複数地域については次節で行います。
現役で火山の鳥島でお湯を沸かすこともできます。Blenderでは物理演算もありますから、色々といたずらできます。
東京都で１区画だけで済む場所だと、
Copyright© (C) Shizuka-Toikanbetsu.