PLATEAU の CityGML Building LOD2 を MVT にしてみた

これは ベクトルタイル Advent Calendar 2021 の24日目の記事です。

※ なお、この記事は筆者の所属するインディゴ株式会社の取組を紹介するものですのでご承知おきください。 ※

はじめに

今年の3月に Project PLATEAU が公開され、日本各地の都市モデルが続々とリリースされていきました。インディゴ株式会社では 6月くらいに 東京都23区の CityGML Building LOD0 の MVT を CC-BY4.0 のデータセットとして GitHub に公開しました。Qiita の記事含め、リファーしていただいたり関心を持っていただいたみなさん、この場にて御礼申し上げます。

さて、この記事ではちょっとレベルを上げて 東京都23区の CityGML Building LOD2 を MVT にしてみた事例 をご紹介します。

できあがり

まずは以下をご覧ください。
ブラウザ上で maplibre-gl-js を動かして地理院淡色地図 を表示、その上に今回作成した MVT を fill-extrusion で表示してみたものです。

都庁周辺

https://indigo-lab.github.io/plateau-lod2-mvt/#15.72/35.689623/139.694231/39.6/60

羽田空港

https://indigo-lab.github.io/plateau-lod2-mvt/#15.17/35.55008/139.7861/25.2/60

東京タワー

https://indigo-lab.github.io/plateau-lod2-mvt/#16/35.658754/139.746405/66/60

大手町周辺

https://indigo-lab.github.io/plateau-lod2-mvt/#15.58/35.685395/139.764087/53.9/49

皇居周辺

https://indigo-lab.github.io/plateau-lod2-mvt/#13.18/35.68207/139.76231/53.9/49

データセット

3D都市モデル（Project PLATEAU）東京都23区（CityGML 2020年度） をソースとし、23区全域について MVT を生成しました。
作成したデータセットは以下のレポジトリで公開しています。ライセンスは CC-BY4.0 です。
MVT の詳細についてもこちらに書いていますのでデータを使用したい方はどうぞ。

デモページの入り口はこちらになります。

技術情報

つくりかた

LOD0 の場合はこんな手順で MVT を作っていました。

1. <bldg:measuredHeight> を持つような <bldg:Building> を見つける
2. <bldg:Building> 配下の <bldg:lod0RoofEdge> あるいは <bldg:lod0FootPrint> から緯度経度ポリゴンを抽出
3. <bldg:measuredHeight> より建物の高さ属性を取得
4. 建物高さを属性として持つような緯度経度ポリゴンを GeoJSON として生成
5. Tippecanoe で処理することによって MVT を得る

※ <hoge:puyo> は CityGML(XML) の中から該当する XML Element を探しましょう、ということです。念のため。

今回は LOD2 をターゲットにしていますが、こんな手順にアップデートしています。

1. <bldg:RoofSurface> と <bldg:GroundSurface> を持つような <bldg:Building> を見つける
2. <bldg:GroundSurface> の三次元ポリゴン座標から Z 座標のみを抽出して接地面の Z 座標の代表値を決定する
3. <bldg:RoofSurface> 配下の各三次元ポリゴンについて、Z 座標の代表値を決定し、また、緯度経度ポリゴンを抽出
4. 緯度経度ポリゴン + 高さ属性 (3. の Z 座標代表値から 2. の Z 座標代表値を引いたもの) を GeoJSON として生成
5. Tippecanoe で処理することによって MVT を得る

<bldg:RoofSurface> というのが要するに LOD2 の屋根です。
ひとつの建物が複数の <bldg:RoofSurface> を持つこともあります。
元データの <bldg:RoofSurface> は必ずしも地面に対して水平ではないのですが、
ここでは LOD1 のように「水平である」ものとして、緯度経度ポリゴン＋Z座標代表値、の形で単純化していることにご注意下さい。

代表値については現状 ポリゴンを構成する各点の Z 座標の最大値と最小値の平均値 を採用しています。

スタイリング

MVT には特に色情報を持たせていません。fill-extrusion-color で任意の着色が可能です。デモでは、せっかく LOD2 でこまかなパーツが描かれているのでそれをなるべく目立たせるために塗り分けをしています。

具体的には以下のようなスタイルを設定することで Z値の小数点以下第一位の数字をキーとした塗り分け を実現しています。ちなみに使用している配色はディック・ブルーナのブルーナカラーです。

{
  "id": "fill-extrusion-bldg",
  "type": "fill-extrusion",
  "source": "bldg",
  "source-layer": "bldg",
  "minzoom": 10,
  "maxzoom": 20,
  "filter": ["has", "z"],
  "paint": {
    "fill-extrusion-color": [
      "interpolate-lab",
      ["linear"],
      ["%", ["round", ["*", 10, ["get", "z"]]], 10],
      0, "#d95d2a",
      1, "#d95d2a",
      2, "#f8df4a",
      3, "#f8df4a",
      4, "#1d508e",
      5, "#1d508e",
      6, "#3c782e",
      7, "#3c782e",
      8, "#845c40",
      9, "#726f62"
    ],
    "fill-extrusion-height": ["*", ["get", "z"], 1]
  }
}

まとめ

CityGML の LOD2 から MVT を生成してみました。MVT では水平でない面の表現に難があるので、必ずしも元データに対して忠実ではない部分もありますが、（人間から見て）ある程度建物の特徴をつかめるような再現ができているのではないかと思います。

おまけ

CityGML(XML) から GeoJSON への変換、どうやっているの？という疑問をお持ちと思います。
参考情報として LOD0 を処理するためのコード例を紹介します。

このような XSLT ファイルを用意して

style.xsl

<?xml version="1.0" encoding="UTF-8"?>
<xsl:stylesheet version="1.0" xmlns:xsl="http://www.w3.org/1999/XSL/Transform"
  xmlns:core="http://www.opengis.net/citygml/2.0" xmlns:bldg="http://www.opengis.net/citygml/building/2.0"
  xmlns:gml="http://www.opengis.net/gml" xmlns:str="http://exslt.org/strings" extension-element-prefixes="str">
  <xsl:output method="text" indent="no" encoding="UTF-8" />
  <xsl:template match="/">
    <xsl:apply-templates select="/core:CityModel/core:cityObjectMember/bldg:Building" />
  </xsl:template>
  <xsl:template match="bldg:Building" />
  <xsl:template match="bldg:Building[bldg:lod0FootPrint][bldg:measuredHeight]">
    <xsl:apply-templates select="bldg:lod0FootPrint//gml:Polygon">
      <xsl:with-param name="z" select="number(bldg:measuredHeight)" />
    </xsl:apply-templates>
  </xsl:template>
  <xsl:template match="bldg:Building[bldg:lod0RoofEdge][bldg:measuredHeight]">
    <xsl:apply-templates select="bldg:lod0RoofEdge//gml:Polygon">
      <xsl:with-param name="z" select="number(bldg:measuredHeight)" />
    </xsl:apply-templates>
  </xsl:template>
  <xsl:template match="gml:Polygon">
    <xsl:param name="z" />
    <xsl:text>{"type":"Feature","properties":{"z":</xsl:text>
    <xsl:value-of select="$z" />
    <xsl:text>},"geometry":{"type":"Polygon","coordinates":[</xsl:text>
    <xsl:for-each select=".//gml:posList">
      <xsl:if test="position()!=1">,</xsl:if>
      <xsl:text>[</xsl:text>
      <xsl:variable name="t" select="str:split(.)" />
      <xsl:for-each select="$t[position() mod 3 = 1]">
        <xsl:variable name="i" select="position()" />
        <xsl:if test="$i != 1">,</xsl:if>
        <xsl:value-of select="concat('[',$t[$i*3-1],',',string(.),']')" />
      </xsl:for-each>
      <xsl:text>]</xsl:text>
    </xsl:for-each>
    <xsl:text>]}}&#10;</xsl:text>
  </xsl:template>
</xsl:stylesheet>

こんなかんじでひとつの CityGML ファイルから GeoJSON Text Sequence を生成できます。

$ xsltproc style.xsl 53392643_bldg_6697_op2.gml
{"type":"Feature","properties":{"z":10.6},"geometry":{"type":"Polygon","coordinates":[[[139.79319067619434,35.53833815357543],...
{"type":"Feature","properties":{"z":31.7},"geometry":{"type":"Polygon","coordinates":[[[139.79143959719713,35.54082886939892],...
...
$ 

入力がたくさんあるのであればこんなかんじで。

$ find . -name '*.gml' -exec xsltproc style.xsl {} \;

MVT 生成するならこのように。

$ (find . -name '*.gml' -exec xsltproc style.xsl {} \; ) | tippecanoe --no-tile-compression -ad -an -Z10 -z16 -e dist -l bldg -ai

実態は40行程度のスタイルシートとワンライナー、でなんとかなってしまうレベルですので気軽に試行錯誤してもらえればいいんじゃないかと思います。東京都23区の CityGML がすべて手元にある状態で、10分もあれば東京都23区全域の MVT の生成が完了する、といった速度で動きます。

ただ、XSLT/xsltproc の特性として、たとえば CSV のパースやら LOD2 の Z座標代表値計算みたいなちょっとややこしい処理を書こうとすると、途端に処理時間の桁が変わるというのがありがちです。上の例では座標列の整形に XSLT の拡張である EXSLT を使ってしのいだりもしています。

実際のところは XSLT で大雑把に処理したものを Node.js で調整、みたいな組み方が落とし所かと思います。適材適所でどうぞ。

（以上）