PLATEAUのデータ全部処理する その2

CityGMLから必要なデータを抽出する

引き続き、Colaboratoryでやっていきます。以下のコードは、セルに入力する前提です。

import

セルに以下を入力し実行します。

import xml.etree.ElementTree as ET
import math
import os

とりあえず使うものをimportします。

タイル計算

def deg2tile(lon,lat,zoom):
  return (math.floor(((lon + 180) / 360) * 2**zoom), math.floor(((1 - math.log(math.tan((lat * math.pi) / 180) + 1 / math.cos((lat * math.pi) / 180)) / math.pi) / 2) * 2**zoom), zoom)

def tile2deg(x,y,zoom):
  n = math.pi-2*math.pi*y/2**zoom
  return ((x/2**zoom)*360-180, 180/math.pi*math.atan(0.5*(math.exp(n)-math.exp(-n))))

タイルの座標計算のための関数を用意します。

LOD0の座標値だけをタイルごとに仕分けしてファイルに保存

def writeoutcoord(filename):
  context = ET.iterparse(filename,events=("start","end"))
  event,root = next(context)

  com_count=0

  for event,elem in context:
    if event == 'start' and elem.tag.endswith('cityObjectMember'):
      poslist = []
      postype = ''
    if event == 'start' and elem.tag.endswith('Building'):
      gmlid=elem.get('{http://www.opengis.net/gml}id')
    if event == 'start' and elem.tag.endswith('stringAttribute'):
      if elem.get('name') == '建物ID':
        next(context)
        _, _elem = next(context)
        bldgid = _elem.text
    if event == 'end' and elem.tag.endswith('measuredHeight'):
      measured_height = elem.text
    if event == 'start' and elem.tag.endswith('lod0RoofEdge'):
      if postype=='':
        postype='Roof'
      elif postype!='Roof':
        print('Type mixed with Roof and ',postype)
      for _event, _elem in context:
        if _event=='end' and _elem.tag.endswith("posList"):
          poslist.append(_elem.text)
        if _event=='end' and _elem.tag.endswith('lod0RoofEdge'):
          break
    if event == 'start' and elem.tag.endswith('lod0FootPrint'):
      if postype=='':
        postype='Foot'
      elif postype!='Foot':
        print('Type mixed with Foot and ',postype)
      for _event, _elem in context:
        if _event=='end' and _elem.tag.endswith("posList"):
          poslist.append(_elem.text)
        if _event=='end' and _elem.tag.endswith('lod0FootPrint'):
          break
    if event == 'end' and elem.tag.endswith('cityObjectMember'):
      lat = 0
      lon = 0
      count = 0
      for str in poslist:
        vallist = str.split(' ')
        for i in range(int(len(vallist)/3)):
          lat += float(vallist[3*i])
          lon += float(vallist[3*i+1])
          count = count+1
      #print(lat/count)
      #print(lon/count)
      (x,y,_) = deg2tile(lon/count,lat/count,16)
      #print(x,y)

      x0 = int(x/1000)
      y0 = int(y/1000)
      x1 = int((x%1000)/10)
      y1 = int((y%1000)/10)
      xy = int(x%10 * 10 + y%10)
      #print(x0,y0,x1,y1,xy)

      dir = f"data/{x0}/{y0}/{x1}/{y1}/{xy}"
      os.makedirs(dir,exist_ok=True)

      f = open(f"{dir}/{bldgid}.dat","x")
      f.write(gmlid)
      f.write("\n")
      f.write(bldgid)
      f.write("\n")
      f.write(measured_height)
      f.write("\n")
      f.write(postype)
      f.write("\n")

      for line in poslist:
        f.write(line)
        f.write("\n")
      f.close()
      com_count=com_count+1

      root.clear()
  print(com_count)

ElementTreeでiterparseを使って逐次的に読み込んで処理していきます。通常のparseだとすべてメモリに読み込んで処理しますが、CityGMLは結構サイズが大きいので、逐次読み込みでやってみました。

また、一つのCityObjectMemberごとに全頂点の平均座標を計算し、それをもとに、どのタイルに分類するかを決めています。
さらに、タイルに分類＝フォルダ分けする処理をしていますが、大量のファイルがフォルダごとにあると処理が重いので、タイル番号を2桁ずつ区切った番号で階層的にフォルダを分類しています。

以上でタイルごとに分類したLOD0の座標値が書かれたファイルが作成されました。

実際のファイルの中身は以下のような感じになっています。

BLD_d5d3f900-7995-43e1-9bcb-005af22cf053
13117-bldg-63990
16.6
Roof
35.787588338578324 139.72745887135878 15.1315 35.787594839140176 139.7274711420614 15.1315 35.787604671206104 139.7274794258648 15.1315 35.78761639153886 139.72748250781189 15.1315 35.78762810680783 139.72747983728817 15.1315 35.78763819496447 139.72747185894042 15.1315 35.787644944408505 139.72745968128325 15.1315 35.78764736549978 139.72744540749738 15.1315 35.787645189487534 139.72743091855986 15.1315 35.78763868892162 139.7274186478528 15.1315 35.78762885685157 139.72741036405338 15.1315 35.787617136517014 139.72740728211642 15.1315 35.78760542134707 139.72741006327485 15.1315 35.787595423329925 139.7274180415075 15.1315 35.787588583657204 139.7274301086558 15.1315 35.78758616266511 139.72744449305597 15.1315 35.787588338578324 139.72745887135878 15.1315