縦横比が正しい浮き出し地図を作る（xml）

　緯度経度の値で作っただけの地図は横太りです。Blenderには単精度の壁があって、正しい図形をきちんと保持することは難しいのですが、それでも一応、読み込んだ直後の地図くらいはできるだけ正確にしておきたいものです。Blender上でのマスターとなります。

　緯度経度の値から一般的な座標に変換するのは結構難しいです。まじめにやろうとすればするほど難しくなります。国土地理院では対象地域に応じたパラメーターを使い分けながら計算する方法を使っています。しかし、変換用のデーターを用意するのも大変です。

　そこで、とりころーるきゃっとさんが掲載された、フライトシミュレーターではこうやって計算しているはずという式をPythonに移植して使うことにします。

import math # # 世界座標の基準緯度と世界座標の目標座標の距離をメートル単位で返す 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

　定数項はあらかじめ計算しておき、まとめて掛け算もできますが、このリストでは意味が分かるように掛けたり割ったりしています。40007000は地球を縦に１周した距離、40075000は赤道を１周した距離で単位はメートルです。うまく考えられた計算式です。Blenderの図形オブジェクトアクセスは大変遅いので、この部分で定数を使って処理速度を上げようとしても、全く体感できる効果はありません。

　メートル単位で作りますから、このままだとオブジェクトが大きくなります。そこで、オブジェクトのスケール値を1000分の1にすることにします。

　オブジェクトの中心を原点ゼロに配置するように作ります。

　なお、Blenderで地図を扱うときに気を付けることを読んでいない方は、そちらを読んでおいてください。

　スクリプトを文字でも掲載しますが、ここに生pyファイルがありますので、ダウンロードできます。

#!BPY import bpy import math from xml.etree.ElementTree import ElementTree import glob # # 縦横比が正しい浮き出し地図を作る（xml） # main001.py # 2016/11/17 Shizuka-Toikanbetsu # XmlName = "N03-??_13_*.xml" # 収録されているxmlファイルを指定 AreaCode = "13101" # 地方公共団体コードを文字列型で指定 ObjScale = 0.001 # 原寸では大きいので1000分の１ # # 世界座標の基準緯度と世界座標の目標座標の距離をメートル単位で返す 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 # # namespace ns = { "ksj": "http://nlftp.mlit.go.jp/ksj/schemas/ksj-app", "gml": "http://www.opengis.net/gml/3.2", "xlink": "http://www.w3.org/1999/xlink" } # def CreateRecName(pAreaCode, pAdministrativeBoundary): # オブジェクト名を作る RecName = pAreaCode tempstr = pAdministrativeBoundary.find( "ksj:prefectureName", ns ).text if tempstr: RecName = RecName + "_" + tempstr tempstr = pAdministrativeBoundary.find( "ksj:countyName", ns ).text if tempstr: RecName = RecName + "_" + tempstr tempstr = AdministrativeBoundary.find( "ksj:cityName", ns ).text if tempstr: RecName = RecName + "_" + tempstr return RecName # def MergeRecName(): Buff = [] for item in bpy.data.objects: if item.type == "MESH": temp = item.name Buff.append( temp[ 0 : 5 ] ) Buff.sort() AreaCode = [] OldCode = "" for temp in Buff: if temp != OldCode: AreaCode.append( temp ) OldCode = temp # for Code in AreaCode: bpy.ops.object.select_all(action='DESELECT') flag = True for curve in bpy.data.objects: if curve.type == 'MESH': NowName = curve.name NowCode = NowName[:5] if NowCode == Code: if flag: RecName = NowName flag = False curve.select = True bpy.context.scene.objects.active = curve bpy.ops.object.join() bpy.ops.object.select_all(action='DESELECT') # 小数点の付く名前は書き換える for item in bpy.data.objects: if item.type == "MESH": temp = item.name try: item.name = temp[ 0 : temp.index( "." ) ] except ValueError: pass # bpy.context.scene.objects.active = None # def CreateSegVert(pText): LineData = pText.split( "\n" ) XY = [] for LineStr in LineData: if LineStr: buff = [] temp = LineStr.lstrip( " \t" ) if len( temp ): temp = temp.split( " " ) dy = float( temp[ 0 ] ) dx = float( temp[ 1 ] ) buff.append( dx ) buff.append( dy ) XY.append( buff ) btm = len( XY ) if btm: btm = btm - 1 if XY[ 0 ] == XY[ btm ]: del XY[ btm ] return XY # # お待たせしました。ここがメインです。 # 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 ) # files = glob.glob( XmlName ) if files: file = files[ 0 ] # きりがないので処理は１ファイルに限る tree = ElementTree() tree.parse( file ) # 事前にセンターの位置を割り出す Xmin = 100000000.0 Ymin = Xmin Xmax = -Xmin Ymax = Xmax # xpath = ".//*/..[ksj:administrativeAreaCode='" + AreaCode + "']" for AdministrativeBoundary in tree.findall( xpath, ns ): # 面の定義を探す bounds = AdministrativeBoundary.find( 'ksj:bounds', ns ) href = bounds.get('{http://www.w3.org/1999/xlink}href') NowSf = href.lstrip('#') # 面の定義情報からカーブを取り出す xpath = ".//*/..[@gml:id='" + NowSf + "']" Surface = tree.find( xpath, ns ) exterior = Surface.find( ".//gml:exterior", ns ) temp = exterior.find( ".//gml:curveMember", ns ) href = temp.get( "{http://www.w3.org/1999/xlink}href" ) NowCv = href.lstrip('#') # 座標定義位置を探してカーブオブジェクトを作成 xpath = ".//*/..[@gml:id='" + NowCv + "']" Curve = tree.find( xpath, ns ) posList = Curve.find( './/gml:posList', ns ) LineData = posList.text.split( '\n' ) Buff = [] for temp in LineData: if temp: Buff.append( temp ) Count = len( Buff ) - 1 for ptr in range( 0, Count ): NowText = Buff[ ptr ].lstrip( ' \t' ) xyStr = NowText.split( ' ' ) posY = float( xyStr[ 0 ] ) posX = float( xyStr[ 1 ] ) Xmin = min( posX, Xmin ) Xmax = max( posX, Xmax ) Ymin = min( posY, Ymin ) Ymax = max( posY, Ymax ) # # ↓ こんなにやって結局欲しかったのはこの２行 wx = ( Xmin + Xmax ) / 2.0 wy = ( Ymin + Ymax ) / 2.0 # # ここからが本番 xpath = ".//*/..[ksj:administrativeAreaCode='" + AreaCode + "']" for AdministrativeBoundary in tree.findall( xpath, ns ): RecName = CreateRecName( AreaCode, AdministrativeBoundary ) # 面の定義を探す bounds = AdministrativeBoundary.find( 'ksj:bounds', ns ) href = bounds.get('{http://www.w3.org/1999/xlink}href') NowSf = href.lstrip('#') # 面の定義情報からカーブを取り出す xpath = ".//*/..[@gml:id='" + NowSf + "']" Surface = tree.find( xpath, ns ) exterior = Surface.find( ".//gml:exterior", ns ) temp = exterior.find( ".//gml:curveMember", ns ) href = temp.get( "{http://www.w3.org/1999/xlink}href" ) NowCv = href.lstrip('#') # 座標定義位置を探してカーブオブジェクトを作成 xpath = ".//*/..[@gml:id='" + NowCv + "']" Curve = tree.find( xpath, ns ) posList = Curve.find( './/gml:posList', ns ) verts = CreateSegVert( posList.text ) # セグメント用の座標点リスト # curvedata = bpy.data.curves.new( name = RecName, type = "CURVE" ) curvedata.dimensions = "2D" objectdata = bpy.data.objects.new( RecName, curvedata ) objectdata.location = ( 0, 0, 0 ) objectdata.scale[0] = ObjScale objectdata.scale[1] = ObjScale objectdata.scale[2] = ObjScale bpy.context.scene.objects.link( objectdata ) polyline = curvedata.splines.new( "POLY" ) # Count = len( verts ) polyline.points.add( Count - 1 ) for ptr in range( 0, Count ): posX = verts[ ptr ][ 0 ] posY = verts[ ptr ][ 1 ] posX = GetLonOffset( wx, posY, posX ) posY = GetLatOffset( wy, posY ) Point = polyline.points[ptr] Point.co = ( posX, posY, 0, 0 ) polyline.use_cyclic_u = True # カーブをメッシュオブジェクトに変換 objectdata.select = True bpy.context.scene.objects.active = objectdata bpy.ops.object.convert( target = "MESH" ) objectdata.select = False # tree = None # どうせ終わるので無くてもよい # MergeRecName() # 複数レコードを１つにまとめたいなら実行 # print( '\nEnd Script' )

　実は座標点を２度読みする無駄があるのですが、処理を追いかけやすくするためにこのようにしています。

　横太りが解消されました。また、図形が原点ゼロに配置されました。この時点で一度保存しておき、以降の加工作業はコピーを使うとよいでしょう。

　これまでは、カーブオブジェクトでしたが、これからはメッシュオブジェクトに変換し、メッシュオブジェクトで管理します。この先座標点を減らしたりとか、穴の開いた図形を作るにはメッシュオブジェクトでなければなりません。

　これから浮き出し地図を作りますが、精度の良い地図を変換した直後で、座標点が多すぎます。座標点を減らしてみます。

　座標を減らすには、オブジェクトを選択してからモディファイアーのポリゴン数削減で、比率を変えます。好みの形になるまで少しずつ試してください。この例ではおもいきりデフォルメしたいので、0.05にします。

　続いて、高さを付けます。再度モディファイアーから厚み付けで1000にしてみます。オブジェクトが1000分の1になっているので、希望の値を1000倍で入力します。

　厚み付けは下に伸びたので、オブジェクトのZ軸を上に引っ張り上げれば定位置にオブジェクトが配置されます。

<トップページに戻る>