※引き続きColaboratory上で実行していきます。以下のコードは、それぞれセルに入力して実行する前提です。

準備　三角形化ライブラリの導入

ポリゴンを三角形に分割するためのライブラリを導入します。

こちらのEarcutアルゴリズムのライブラリを使います。

このソースコードを丸ごとセルに貼り付けて実行するとearcutが使えるようになります。

準備　座標変換

WGS84の緯度経度から、分割したタイルの中心を原点として、高さ方向がちゃんと上を向くような座標系に変換します。
以下のプログラムをセルに入力して実行しておきます。

def BLH2XYZ(b,l,h):
  # WGS84
  a = 6378137.0
  f = 1.0 / 298.257223563

  b = math.radians(b)
  l = math.radians(l)

  e2 = f * (2.0 - f)
  N = a / math.sqrt(1.0 - e2 * math.pow(math.sin(b), 2.0))

  X = (N + h) * math.cos(b) * math.cos(l);
  Y = (N + h) * math.cos(b) * math.sin(l);
  Z = (N * (1 - e2) + h) * math.sin(b);

  return (X,Y,Z)

# spcx,spcy,zoom = actualy target tile value
def originTrans(spcx : int,spcy : int,zoom : int,poslist):
  # calc center point
  (cx,cy) = tile2deg(2*spcx+1,2*spcy+1,zoom+1)
  ox,oy,oz = BLH2XYZ(cy,cx,0)

  retposes = []
  for pos in poslist:
    positions = pos.split(' ')

    retpos = []
    for i in range(int(len(positions)/3)):
      _y = float(positions[3*i])
      _x = float(positions[3*i+1])
      _z = float(positions[3*i+2])

      x,y,z = BLH2XYZ(_y,_x,_z+getGeoidValue(_y,_x))
      x = x - ox
      y = y - oy
      z = z - oz

      s = math.radians(-cx)
      rx = x * math.cos(s) - y * math.sin(s)
      ry = x * math.sin(s) + y * math.cos(s)

      s = math.radians(-cy)
      rxx = rx * math.cos(s) - z * math.sin(s)
      rz  = rx * math.sin(s) + z * math.cos(s)

      retpos.append(ry)
      retpos.append(rz)
      retpos.append(rxx)
    retposes.append(retpos)
  return retposes

OBJ変換実行

以下のプログラムを実行するとタイルごとに分割された座標データファイルを統合してタイルごとのOBJファイルを作ってくれます。

x0s = os.listdir(path='data')

print(x0s)
for x0 in x0s:
  y0s = os.listdir(path=f"data/{x0}/")
  for y0 in y0s:
    x1s = os.listdir(path=f"data/{x0}/{y0}/")
    for x1 in x1s:
      y1s = os.listdir(path=f"data/{x0}/{y0}/{x1}/")
      for y1 in y1s:
        xy2s = os.listdir(path=f"data/{x0}/{y0}/{x1}/{y1}/")
        for xy2 in xy2s:
          spcx = int(x0)*1000+int(x1)*10+int(int(xy2)/10)
          spcy = int(y0)*1000+int(y1)*10+int(int(xy2)%10)
          print(spcx,spcy)

          objfile = open(f"obj/{spcx}_{spcy}.obj",'w')
          index = 1
          vs = []
          fs = []

          filelist = os.listdir(f"data/{x0}/{y0}/{x1}/{y1}/{xy2}/")
          for file in filelist:
            f=open(f"data/{x0}/{y0}/{x1}/{y1}/{xy2}/{file}")

            gmlid = f.readline()
            bldg_id = f.readline()
            measured_height = float(f.readline())
            type = f.readline().strip()
            if type!='Roof':
              print(type)
            poslist = f.readlines()

            positions = originTrans(spcx,spcy,16,poslist)

            for poses in positions:
              vids = earcut(poses,dim=3)
              numv = int(len(poses)/3)
              for i in range(numv):
                x = poses[3*i]
                z = poses[3*i+1]
                y = poses[3*i+2]
                vs.append(f"v {x} {y} {z}\n")
                if type == 'Roof':
                  vs.append(f"v {x} {y-measured_height} {z}\n")
                else:
                  vs.append(f"v {x} {y+measured_height} {z}\n")
              numf = int(len(vids)/3)
              
              for i in range(numf):
                v1 = 2*vids[3*i]+index + (0 if type=='Roof' else 1)
                v2 = 2*vids[3*i+1]+index + (0 if type=='Roof' else 1)
                v3 = 2*vids[3*i+2]+index + (0 if type=='Roof' else 1)
                fs.append(f"f {v1} {v2} {v3} {v1}\n")
              
              for j in range(numv):
                p1 = 2*j
                p2 = 2*j+1
                p3 = 2*j+2
                p4 = 2*j+3

                if p3 >= numv*2:
                  p3 = p3 - numv*2
                if p4 >= numv*2:
                  p4 = p4 - numv*2
                
                fs.append(f"f {p1+index} {p2+index} {p3+index} {p1+index}\n")
                fs.append(f"f {p4+index} {p3+index} {p2+index} {p4+index}\n")
              index = index+numv*2
            f.close()

          for vline in vs:
            objfile.write(vline)
          for fline in fs:
            objfile.write(fline)
          objfile.close()

作成されたOBJが以下です。