国土数値情報 ダウンロードサービスの利用

試しに1地域の地図(xml)

 Blenderのカーブオブジェクトに地図データーを使えることは確認できました。そこで、とにかく1地域の面を作ってみます。まだ実用的ではありません。このページでは、xmlファイルを利用します。

 例にするのは東京都千代田(ちよだ)区です。

 小笠原諸島なども東京都ですから、結構データー量があります。この中から千代田区だけを取り出すことを考えると、2つの方法があります。

 千代田区という漢字の文字列を使う方法と、地方公共団体コードを使う方法があります。

 漢字で探す方法もありますが、東京23区とそれ以外では、名称の文字列を格納している位置(項目)が異なって収録されています。

 そこで名称の文字列ではなく、地方公共団体コードを検索に使うことにします。千代田区の場合は13101です。団体コードはどのレコードでも項目位置は同じです。

 カーブオブジェクトから面にするには、セグメントを作る必要があるので座標点をぐるっと1周ではなく、最後の座標点は出発点の1つ手前で登録をやめます。use_cyclic_uをTrueにすると、セグメントを持つ図形の意味となり、カーブのタイプを2Dにすると線の内側が面になります。

 座標の定義ですが先のシェープファイルでの例を考えると、ぐるっと1周した定義なのか、簡略した定義なのか保証できません。そこで、面倒でも1周しているか調べ、1周していれば最後の座標点は出発位置と同じなので削除するようにします。関数CreateSegVertで行っています。

#!BPY
import bpy
from xml.etree.ElementTree import ElementTree
import glob
#
# 試しに1地域の浮き出し地図(xml)
#
XmlName = "N03-*.xml"   # 収録されているxmlファイルを指定
AreaCode = "13101"      # 地方公共団体コードを文字列型で指定
#
# 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     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 == "CURVE" or item.type == "MESH":
                bpy.context.scene.objects.unlink( item )
for item in bpy.data.objects:
        if item.type == "CURVE" or item.type == "MESH":
                bpy.data.objects.remove( item )
for item in bpy.data.curves:
        bpy.data.curves.remove( item )
for item in bpy.data.meshes:
        bpy.data.meshes.remove( item )
#
files = glob.glob( XmlName )
if files:
        file = files[ 0 ]       # きりがないので処理は1ファイルに限る
        tree = ElementTree()
        tree.parse( file )
        #
        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 )
                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 ]
                        polyline.points[ptr].co = ( posX, posY, 0, 0 )
                polyline.use_cyclic_u = True
        tree = None     # どうせ終わるので無くてもよい
        #
        # ↓ 複数レコードを1つにまとめたいなら実行
        for Code in AreaCode:
                bpy.ops.object.select_all( action="DESELECT" )
                for curve in bpy.data.objects:
                        if curve.type == "CURVE":
                                NowName = curve.name
                                NowCode = NowName[:5]
                                if NowCode == AreaCode:
                                        curve.select = True
                                        bpy.context.scene.objects.active = curve
                bpy.ops.object.join()
                bpy.ops.object.select_all(action='DESELECT')
                objectdata.name = RecName
        bpy.context.scene.objects.active = None
        # ↑ 複数レコードを1つにまとめたいなら実行
#
print( '\nEnd Script' )

 先の例と同じように、まだ緯度経度の数値のままです。

 団体コードに該当する定義だけを取り出し、その中から面の情報を取り出し、さらに面の中からカーブ座標定義を取り出すという流れです。最後に、もし複数レコードで1つにまとめたいときは、該当部分を実行します。xmlでのデーター値を使った検索、アトリビュートの検索が記述されています。

 もう一度同じ画像を掲載しますが、千代田区や周辺にお住まいの方ならすぐにわかることですが、この千代田区は横太りです。緯度経度のままでは適切な形状になりません。

 さらに、千代田区ではうまくいきましたが、大島町や小笠原村などではうまくいきません。中心から座標の位置が単精度ではすべて表しきれず、同一座標とみなして処理の中断などが発生します。単精度の壁を乗り越える必要があります。Blenderで地図を扱うときに気を付けることを読んでいない方は、そちらを確認してください。

 実用的な意味は全くありませんが、横太りのままで簡単に面を作るには、単精度の有効値を増やします。ソースの最初のほうを改造します。大島町のコードは13361です。goffxとgoffyの2つの変数を追加します。


AreaCode = "13361"      # 地方公共団体コードを文字列型で指定
#
goffx = 120.0
goffy = 20.0
#
# namespace

 次に終わりのほうで、複数レコードをまとめる処理の上のほうで、posXとposYに値を代入している部分を書き換えます。引き算を加筆します。


posX = verts[ ptr ][ 0 ] - goffx
posY = verts[ ptr ][ 1 ] - goffy

 結局、座標値のxから120、yから20を引いただけです。これを実行すると図形の位置が移動します。次のようになります。

 値を少し引いただけですが、単精度では有効桁数が少し増えるため面が作れたというわけです。微妙な世界です。地図データーでは大きな図形も普通に登場しますから、常に単精度の壁を考えながらスクリプトを作ることになります。ただ、Blenderは地図用ソフトではないことをお忘れなく。

 これは、あくまでも手間をかけずに面を作るので値を引いてみたという実験です。実際の本番向けではきちんとXY座標軸に綺麗に対応させ、縦横の比率も正しくする必要があります。

 今回のスクリプトは、xmlでの文字列検索、属性値検索、Blenderのカーブオブジェクトから面を作る方法の基礎になります。

<トップページに戻る>


Copyright© Shizuka-Toikanbetsu.
inserted by FC2 system