縦横比が正しい浮き出し地図を作る（シェープファイル）

　今度はシェープファイルで同じものを作ってみます。

　シェープファイルではdbfの扱いが面倒でソースが長くなります。そこで、どさくさまぎれに別ファイルに分離しました。ソースを分けるだけが目的で、本格的なモジュールにはしていません。少し苦労しますが、処理が速くなります。

　作業用のフォルダには、２つのスクリプトを置くことになります。

　まずは、シェープファイルアクセス担当のスクリプト CShpR.py。ここに生pyファイルもあります。

#!BPY import struct # class CShpR: # class CDbfFieldRec: def __init__(self): self.name = "" self.Top = 0 self.Btm = 0 # def __init__(self): self.fshp = None self.fshx = None self.fdbf = None # def GetBinDword(self, Item): return struct.unpack( "L", Item )[ 0 ] # def GetBinDouble(self, Item): return struct.unpack( "d", Item )[ 0 ] # def GetBinDwordRev(self, Item): return struct.unpack( "L", Item[::-1] )[ 0 ] # def close(self): if self.fdbf: self.fdbf.close() self.fdbf = None if self.fshx: self.fshx.close() self.fshx = None if self.fshp: self.fshp.close() self.fshp = None # def open(self, pFile): self.close() self.fshp = open( pFile, "rb" ) # self.fshp.seek( 32, 0 ) tempbin = self.fshp.read( 4 ) self.VolType = self.GetBinDword( tempbin ) tempbin = self.fshp.read( 32 ) self.VolXmin = self.GetBinDouble( tempbin[ 0 : 8 ] ) self.VolYmin = self.GetBinDouble( tempbin[ 8 : 16 ] ) self.VolXmax = self.GetBinDouble( tempbin[ 16 : 24 ] ) self.VolYmax = self.GetBinDouble( tempbin[ 24 : 32 ] ) self.VolWx = ( self.VolXmin + self.VolXmax ) / 2.0 self.VolWy = ( self.VolYmin + self.VolYmax ) / 2.0 # temp = pFile.split( "." ) self.fshx = open( temp[ 0 ] + ".shx", "rb" ) self.fdbf = open( temp[ 0 ] + ".dbf", "rb" ) self.fdbf.read( 4 ) # skip tempbin = self.fdbf.read( 4 ) self.DbfRecCount = self.GetBinDword( tempbin ) tempbin = self.fdbf.read( 2 ) self.DbfHeaderSize = struct.unpack( "H", tempbin )[ 0 ] tempbin = self.fdbf.read( 2 ) self.DbfRecSize = struct.unpack( "H", tempbin )[ 0 ] self.DbfFieldCount = int( ( self.DbfHeaderSize - 33 ) / 32 ) # フィールドのレコード内位置を割り出す self.fdbf.seek( 32, 0 ) self.DbfField = [] CodePosE = 1 for fldptr in range( 0, self.DbfFieldCount ): self.DbfFieldRec = self.CDbfFieldRec() tempbin = self.fdbf.read( 32 ) temp = tempbin[:10] temp = temp.decode('shift_jis') self.DbfFieldRec.name = temp.rstrip( ' \0' ) temp = tempbin[16:17] temp = struct.unpack( "B", temp ) temp = temp[ 0 ] CodePosE = CodePosE + temp CodePosS = CodePosE - temp self.DbfFieldRec.Top = CodePosS self.DbfFieldRec.Btm = CodePosE self.DbfField.append( self.DbfFieldRec ) # def readdbf(self, pRecNo): self.fdbf.seek( self.DbfRecSize * pRecNo + self.DbfHeaderSize ) return self.fdbf.read( self.DbfRecSize ) # def getdbfvalue(self, pDbfRec, pIndex): temp = pDbfRec[ self.DbfField[ pIndex ].Top : self.DbfField[ pIndex ].Btm ] temp = temp.rstrip( b' ' ) return temp.decode('shift_jis') # def readrec(self, pRecNo): self.fshx.seek( pRecNo * 8 + 100, 0 ) tempbin = self.fshx.read( 8 ) shppos = self.GetBinDwordRev( tempbin[:4] ) * 2 shplen = self.GetBinDwordRev( tempbin[4:8] ) * 2 # self.fshp.seek( shppos + 8, 0 ) return self.fshp.read( shplen ) # def CreateSegVert(self, pText): self.NumParts = self.GetBinDword( pText[36:40] ) self.NumPoints = self.GetBinDword( pText[40:44] ) tparts = [] ptr = 44 for temp in range( 0, self.NumParts ): NowPart = self.GetBinDword( pText[ptr:ptr+4] ) tparts.append( NowPart ) ptr = ptr + 4 tparts.append( self.NumPoints ) # self.XY = [] for temp in range( 0, self.NumPoints ): buff = [] dx = self.GetBinDouble( pText[ptr:ptr+8] ) buff.append( dx ) ptr = ptr + 8 dy = self.GetBinDouble( pText[ptr:ptr+8] ) buff.append( dy ) ptr = ptr + 8 self.XY.append( buff ) # self.parts = [] for temp in range( 0, self.NumParts ): btm = temp + 1 buff = [] buff.append( tparts[ temp ] ) if self.XY[ tparts[ temp ] ] == self.XY[ tparts[ btm ] - 1 ]: buff.append( tparts[ btm ] - 2 ) else: buff.append( tparts[ btm ] - 1 ) self.parts.append( buff )

　続いてメインのスクリプト。ここに生pyファイルがあります。

#!BPY import bpy import math import glob from CShpR import CShpR # # 縦横比が正しい浮き出し地図を作る（シェープファイル） # main002.py # 2016/11/17 Shizuka-Toikanbetsu # ShpName = "N03-??_13_*.shp" # 収録されているshpファイルを指定 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 # 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 # 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( ShpName ) if files: file = files[ 0 ] # きりがないので処理は１ファイルに限る shape = CShpR() shape.open( file ) # for fldptr in range( 0, shape.DbfFieldCount ): temp = shape.DbfField[ fldptr ].name if temp == "N03_001": DbfPref = fldptr elif temp == "N03_003": DbfCounty = fldptr elif temp == "N03_004": DbfCity = fldptr elif temp == "N03_007": DbfArea = fldptr # 事前にセンターの位置を割り出す Xmin = 100000000.0 Ymin = Xmin Xmax = -Xmin Ymax = Xmax for RecNo in range( 0, shape.DbfRecCount ): # tempbin = shape.readdbf( RecNo ) temp = shape.getdbfvalue( tempbin, DbfArea ) if( temp == AreaCode ): tempbin = shape.readrec( RecNo ) shptype = shape.GetBinDword( tempbin[:4] ) if shptype == 5: # polygon is 5 lXmin = shape.GetBinDouble( tempbin[4:12] ) lYmin = shape.GetBinDouble( tempbin[12:20] ) lXmax = shape.GetBinDouble( tempbin[20:28] ) lYmax = shape.GetBinDouble( tempbin[28:36] ) Xmin = min( lXmin, Xmin ) Xmax = max( lXmax, Xmax ) Ymin = min( lYmin, Ymin ) Ymax = max( lYmax, Ymax ) wx = ( Xmin + Xmax ) / 2.0 wy = ( Ymin + Ymax ) / 2.0 # # ここから本番。レコードごとに図形を作る for RecNo in range( 0, shape.DbfRecCount ): # dbfを１レコード読む tempbin = shape.readdbf( RecNo ) temp = shape.getdbfvalue( tempbin, DbfArea ) if( temp == AreaCode ): # ↓ オブジェクト名を作る RecName = AreaCode temp = shape.getdbfvalue( tempbin, DbfPref ) RecName = RecName + "_" + temp temp = shape.getdbfvalue( tempbin, DbfCounty ) if len( temp ): RecName = RecName + "_" + temp temp = shape.getdbfvalue( tempbin, DbfCity ) if len( temp ): RecName = RecName + "_" + temp # ↑ オブジェクト名を作る # tempbin = shape.readrec( RecNo ) shptype = shape.GetBinDword( tempbin[:4] ) if shptype == 5: # polygon is 5 shape.CreateSegVert( tempbin ) # part = 0 #for part in range( 0, len( shape.parts ) ): temp = shape.parts[ part ] top = temp[ 0 ] btm = temp[ 1 ] Count = btm - top + 1 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" ) # polyline.points.add( Count - 1 ) ptr = top for xyptr in range( 0, Count ): posX = shape.XY[ ptr ][ 0 ] posY = shape.XY[ ptr ][ 1 ] ptr = ptr + 1 posX = GetLonOffset( wx, posY, posX ) posY = GetLatOffset( wy, posY ) polyline.points[xyptr].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 # shape.close() # MergeRecName() # 複数レコードを１つにまとめたいなら実行 # print( '\nEnd Script' )

　スクリプトを分けたので、確実に動かすには次のようにします。

1. 　２つのソーススクリプトを作業用フォルダに入れます。

2. 　Blenderを新規に起動するか、初期化した状態にし、作業用フォルダに保存します。作業用フォルダにはこれで３つのファイルが入ります。

3. 　エクスプローラーから保存したBlenderのファイルをダブルクリックで起動します。

4. 　テキストエディタから CShpR.pyを読み込みます。

5. 　続いてテキストエディタからメインのスクリプトを読み込んでスクリプトを実行します。

　BlenderのPythonはモジュールの読み込みデコードを１回しか行わない特性があるので、もしCShpR.py側を変更した場合はスクリプトを保存し、テキストエディタから一度開放し、改めて読み直す必要があります。解放した直後に、Blenderのファイルそのものを一度保存して終了しないと、うまくいかないこともあります。全部開放して、解放した状態を保存してから、改めて読み込みなおすという意味です。

　シェープファイルはちょっと扱いが面倒な部分がありますが、島と周囲の岩など、レコード数が非常に多い場合はxmlよりも処理が速くなります。

<トップページに戻る>