台つき １区画の３次元メッシュ（xml）

　台がついた10mメッシュ音威子府村の一部です。3Dプリンターで出力して、子供たちに色を塗ってねーとかありますね。ただ、座標点が多くなるので個人的には全く使っていません。

　Blenderのメッシュオブジェクトの場合、厚み付けは下方向に伸びます。座標点のインデックス番号ですが、上面の座標は平面のときと全く同じで下面のインデックス番号が同じ並びで続きます。つまり、全部の座標点のうち、前半は上面、後半は下面になります。

　台が付いていると未定義データー部分がそのまま欠損では形が悪いので、ありがちな標高ゼロに置き換えています。

　Pythonのスクリプトは次のようになります。ここに生pyファイルがあります。

import bpy import math import xml.etree.ElementTree as ET # # 台つき１区画の３次元メッシュ（xml） # table_1.py # 2016/11/16 Shizuka-Toikanbetsu # demtree = ET.parse('FG-GML-6742-02-dem10b-20161001.xml') # ここに xml ファイルを指定 tblheight = 2 # 台の部分の高さ # # 世界座標の基準緯度と世界座標の目標座標の距離をメートル単位で返す 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.' ) # # 既存メッシュの削除 bpy.ops.object.mode_set( mode = "OBJECT" ) 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 # bpy.ops.object.modifier_add(type='SOLIDIFY') bpy.context.object.modifiers["Solidify"].thickness = tblheight * 1000 bpy.context.object.modifiers["Solidify"].thickness_clamp = 0 bpy.ops.object.modifier_apply(apply_as='DATA', modifier="Solidify") # 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 ) ndx = dpy * px + dpx bpy.context.object.data.vertices[ ndx ].co[ 0 ] = dx bpy.context.object.data.vertices[ ndx ].co[ 1 ] = dy if PointAlt == -9999.0: PointAlt = 0 bpy.context.object.data.vertices[ ndx ].co[ 2 ] = PointAlt ndx = px * py + ndx bpy.context.object.data.vertices[ ndx ].co[ 0 ] = dx bpy.context.object.data.vertices[ ndx ].co[ 1 ] = dy dpx = dpx + 1 if dpx >= px: dpx = 0 dpy = dpy - 1 # Grid.location[2] = tblheight # print( 'Exit' )

　ただ、鳥島のような小さな陸地は切り出して使ったほうがよさそうです。

<トップページに戻る>