はじめに

前回の記事で、ドローンの画像から点群データ（LAS形式）を作りました。今回は、作った点群データから、Voxel のベクトルタイルを作って、ウェブ地図（MapLibre GL JS）として表示してみることにします。（前回の記事は下のものです）

初心者向けに書こうとおもったのですが、思ったよりも長くなり、また私のやり方がややこしいのでそんなに簡単ではないかもしれません。すみません。

手順

optgeo/kid-cのレポジトリと、YouTube動画（青山学院大学古橋研究室, 2021）を参考にして作業を進めます。kid-cにあるようにRakefileを完璧に作ってボクセルタイルを作成するのはハードルが高いので、今回の作業ではRubyファイルの実行とコマンドラインの実行を組み合わせて、何をしているか理解しながらデータを作成できればと思います。

Step 1. 点群データの確認　

そもそも点群データがどんなものか見てみました。WebODMからCSVファイルを出力してエクセルで見てみます。
X,Yはウェブメルカトルの座標で、Zは標高の値のようですね。また、画像のRGBがRed, Green, Blueとして入っています。RGBはそれぞれ0-255までの8bitのようです。

Step 2. 作業環境の検討と構築（Dockerでやってみます。）

点群データ（LAS）からVoxelを作るために、YouTubeの動画を見ると、以下のツールの実行環境が必要になりそうです。

pdal (Point Data Abstraction Library)
ruby
tippecanoe
unvt/charites

これらをWindowsで作るのは少し難しいので（特にtippecanoeのインストール）、Dockerで実行することにしました。unvt/rgbifyやunvt/nanbanのDockerfileを参考にして、Dockerfileを以下の通り作ってみました。（ファイルはこちら　https://github.com/ubukawa/voxel-br/blob/main/Dockerfile ）
なお、unvt/charitesを入れるためにはnodejsや周辺ツールが必要になるので、まだ入れないことにします。

Dockerfile

FROM osgeo/gdal:ubuntu-small-3.4.0
RUN apt update \
&& apt upgrade -y \
&& apt install -y \
gdal-bin \
git \
libgdal-dev \
libproj-dev \
pdal \
ruby \
sqlite3 \
automake \
cmake \
g++ \
&& git clone https://github.com/mapbox/tippecanoe \
&& cd tippecanoe \
&& make -j3 LDFLAGS="-latomic" \
&& make install \
&& cd .. \

これをunvt/voxelという名前でビルドしてみます。

docker build -t unvt/voxel .

そして実行します。

docker run -it --rm -v ${PWD}:/data unvt/voxel
cd /data

以下の図のように、pdal、ruby、tippecanoeが入っていることを確認しました。（一応、unvt/voxelとしてプッシュしておきました。）

Step 3. 点群データ（LAS）からベクトルタイルの作成

YouTube動画の以下の画面前後の説明を繰り返し聞きました。このYouTubeレクチャーのソースコードは https://github.com/optgeo/kid-c にあります。私はrubyが書けない（≒読めない）ですが、必要最低限を調べながらファイルの内容を理解するように努めました。

(LASデータはウェブメルカトルで準備できてるので、私の作業では投影変換は不要です。また、最後にmbtilesをpbfに展開するtile-joinも使った経験があるので大丈夫です。今回のチャレンジングなところは、lasからGeoJSONを作ってベクトルタイルに変換するというところがメインになります。)

YouTubeでは、Rubyでpdalの処理を記述したあと、パイプ（｜）を使ってstdinとしてpdalに渡しているとういう説明がありました。
作業全体の流れとしては、①投影法変換、②フィルタリング（ボクセルのNearest Neighbourの値を取得）、③GeoJSONへの出力と加工、④ベクトルタイルに変換、という流れになっています。

私はpdalを使うのもこれが初めてです。Rubyでpdalの作業を記述して、コマンドラインのパイプ（｜）でわたす、ということがちょっとユニークで面白いなぁと思いました。pdal pilelineを見てみましたが、例に倣ってrubyを使った方が早そうでした。

3-1. 投影変換　reprojection　

私のデータはウェブメルカトルになっているので、本当は必要ないのですが、Rubyとpdalのテストもかねて、ファイルを読み込んで（そのまま）書き出すということをテストしました。optgeo/kid-cのreproject_pipeline.rbを見ながら、微調整をしつつ、まず 1.rb というファイルを準備します。

1.rb

# for reprojection. (actually it is not needed)
require './constants'
require 'yaml'
require 'json'

basename = ENV['BASENAME']
src_path = "#{TMP_DIR}/#{basename}.las"
dst_path = "#{TMP_DIR}/#{basename}-3857.las"

pipeline = <<-EOS
pipeline: 
  - 
    filename: #{src_path}
    type: readers.las
    spatialreference: "EPSG:3857" 
  -
    type: writers.las
    filename: #{dst_path}
EOS

# Use the following if reprojection is needed.
=begin
pipeline = <<-EOS
pipeline: 
  - 
    filename: #{src_path}
    type: readers.las
    spatialreference: "EPSG:3857" #modify it if it is not in webmercator.
  -
    type: filters.reprojection
    out_srs: "EPSG:3857"
  -
    type: writers.las
    filename: #{dst_path}
EOS
=end

print JSON.dump(YAML.load(pipeline))

このときのconstants.rbは以下のものだけです。BASENAMEはコマンドラインで変数を与えます。

そして以下のように実行します。

BASENAME=sample ruby 1.rb | pdal pipeline --stdin

作業の結果、新しいデータが出てきました。中身はそのままですが、Rubyとpdalが動いていることがわかりました。

3-2. リサンプル

次にリサンプルを準備します。optgeo/kid-cのresample_text_pipeline.rbを見ながら、2.rbというファイルを作ります。pdalの filters.voxelcenternearestneighbor というものを使ってリサンプルしているんですね。このときにCellのspacingを与えてあげるということのようです。

2.rb

require './constants'
require 'yaml'
require 'json'

z = ENV['Z'].to_i
spacing = (BASE ** (Z_ONE_METER - z)).to_f
basename = ENV['BASENAME']
src_path = "#{TMP_DIR}/#{basename}-3857.las" #unless File.exist?(src_path)
dst_path = "#{TMP_DIR}/#{basename}-#{z}.las"
$stderr.print "#{spacing}m for #{dst_path} from #{src_path}\n"

pipeline = <<-EOS
pipeline: 
  - 
    filename: #{src_path}
    type: readers.las
  -
    type: filters.voxelcenternearestneighbor
    cell: #{spacing}
  -
    type: writers.text
    format: csv
    filename: STDOUT
EOS

print JSON.dump(YAML.load(pipeline))

このときのconstants.rbは、以下の通りです。Z_ONE_METERとBASEを足しています。zは環境変数としてあたえるのでconstanrsには入れません。Z_ONE_METERとBASEはSpacingを与えるために入れているようですが、環境変数で与える z が Z_ONE_METER と一致するときは、spacing が 1 になることがわかります。

そして、以下のように実行します。

Z=19 BASENAME=sample ruby 2.rb | pdal pipeline --stdin

STDOUTで出てくるので、以下のような感じです。CSV形式で出ていますね。これを次のステップの作業でgeojsonにしつつ、少し加工することになります。

3-3. 加工、そしてGeoJSONに

上で出てきたものをGeoJSONにすることをイメージしながら3.rbを作ります。optgeo/kid-cのtogeojson.rbをよく見ながら作業します。

x, y, h を属性として作っています。点の(x, y, h)ではなくて、それぞれをspacingで割った剰余を引いていて、ボクセルのスタート地点を調整しているようですね（to_fはストリングにするファンクション、%は剰余の計算みたいです。）。
私がODMから作ったデータはoptgeo/kid-cとは別のソースデータでしたが、RGBはそれぞれ13,14,15になっていたのでそのままで大丈夫そうです。だけど色の深度が違うかもしれません。私のデータは最初からRGBがそれぞれ 8 bitみたいなので、4096で割らないようにします。
minzoom = (z == MINZOOM) ? MINCOPYZOOM : z　はどうしてそうやっているのかがよくわかりません。constraints.rbで与えたMINZOOMがzと同じならMINCOPYZOOM（＝MINZOOM）にして、そうでないならZにしています。
- maxzoomはzなので、minzoomがzよりも大きくならないようにしようということなのかもしれません。
- あるいはzを変えながら何度か実行することで、ズームレベルごとに違うスペーシングでデータをつくるということなのかもしれません。

3.rb

require './constants'
require 'yaml'
require 'json'

z = ENV['Z'].to_i
minzoom = (z == MINZOOM) ? MINCOPYZOOM : z
maxzoom = z
spacing = (BASE ** (Z_ONE_METER - z)).to_f
n = 0

first = true
start_time = Time.now
while gets
  if first
    first = false
    next
  else
    n += 1
  end
  r = $_.strip.split(',')
  x = r[0].to_f - r[0].to_f % spacing
  y = r[1].to_f - r[1].to_f % spacing
  h = r[2].to_f - r[2].to_f % spacing
  h = h.to_i
  color = '#' + r[13..15].map{|v| sprintf('%01x', v.to_i )}.join #check your rgb location
#  color = '#' + r[13..15].map{|v| sprintf('%01x', v.to_i / 4096 )}.join #if your RGB is more
  g = <<-EOS
type: Polygon
coordinates: 
  -
    -
      - #{x}
      - #{y}
    -
      - #{x + spacing}
      - #{y}
    -
      - #{x + spacing}
      - #{y + spacing}
    -
      - #{x}
      - #{y + spacing}
    -
      - #{x}
      - #{y}
  EOS
  g = YAML.load(g)
  f = <<-EOS
type: Feature
properties: 
  color: '#{color}'
  h: #{h}
  spacing: #{spacing}
tippecanoe:
  minzoom: #{minzoom}
  maxzoom: #{maxzoom}
  layer: #{LAYER}
  EOS
  f = YAML.load(f)
  f[:geometry] = g
  print JSON.dump(f), "\n"
end

def hms(s)
  min, sec = s.to_i.divmod(60)
  hour, min = min.divmod(60)
  "%02d:%02d:%02d" % [hour, min, sec]
end

$stderr.print "\n[#{z}] #{hms(Time.now - start_time)}s\n"

このときのconstants.rbは、以下の通りです。Tippecanoeでの変換に備えるために、ズームレベルやレイヤー名も入れるようになってきています。

まえのステップのSTDOUTをパイプするので、以下のようにコマンドを追加してテストしました。もう一つコマンドが追加されましたね。

Z=19 BASENAME=sample ruby 2.rb | pdal pipeline --stdin | Z=19 ruby 3.rb

出力の結果が以下の図です。属性として、色、h、スペーシングが入っていますし、tippecanoe用のパラメータも適切に入っていそうな感じです。このままTippecanoeでの変換にすすんでよさそうです。

3-4. Tippecanoeに流し込む

これまででGeoJSONが吐き出されることが分かったので、このままTippecanoeにしてみます。

Z=19 BASENAME=sample ruby 2.rb | pdal pipeline --stdin | Z=19 ruby 3.rb | tippecanoe --maximum-zoom=19 --minimum-zoom=10 --projection=EPSG:3857 --force --output output.mbtiles --no-tile-size-limit --no-feature-limit

Z=19では、252KBのMBtilesファイルができました。

3-5. ファイルシステムに展開

mbtiles形式のベクトルタイルが得られたので、以下のようにファイルシステムに展開してみます。（展開してみるとZ=19しかなかったので、GeoJSONの中のMAXとMINが優先されたのだと思います。）

tile-join --force --output-to-directory=docs/voxel --no-tile-size-limit --no-tile-compression output.mbtiles

3-6　少し小さいZLも作成

せっかくなので、ZL14からZL18もタイルを作ってみました。

Z=18 BASENAME=sample ruby 2.rb | pdal pipeline --stdin | Z=18 ruby 3.rb | tippecanoe --maximum-zoom=18 --minimum-zoom=10 --projection=EPSG:3857 --force --output output-18.mbtiles --no-tile-size-limit --no-feature-limit
....
Z=14 BASENAME=sample ruby 2.rb | pdal pipeline --stdin | Z=14 ruby 3.rb | tippecanoe --maximum-zoom=14 --minimum-zoom=10 --projection=EPSG:3857 --force --output output-14.mbtiles --no-tile-size-limit --no-feature-limit

tile-join --force --output-to-directory=docs/voxel/v18 --no-tile-size-limit --no-tile-compression output-18.mbtiles
...
tile-join --force --output-to-directory=docs/voxel/v14 --no-tile-size-limit --no-tile-compression output-14.mbtiles

この複数のZLでのタイル作成方法は、私はよくわかっていないようなので、今後改善の余地があります。多分もともとのRakefileの実行ではもっとうまくやっているはず。

3-7 このステップのまとめ

ここまででボクセルタイルのベクトルタイルデータができました。

Step 4. ウェブでのホスティング

4-1. GitHubページ

レポジトリのdocsページをホスティングします。これで、 https://ubukawa.github.io/voxel-br/voxel/{z}/{x}/{y}.pbf にてタイルがホストされました。

4-2. MapLibreのダウンロード

最新バージョンが2.4.0になっていましたね。Powershellでdocsフォルダに移動して以下のようにMapLibre GL JSをダウンロードしてきます。

mkdir maplibre-gl@2.4.0
cd maplibre-gl@2.4.0
curl.exe -O https://unpkg.com/maplibre-gl@2.4.0/dist/maplibre-gl.css
curl.exe -O https://unpkg.com/maplibre-gl@2.4.0/dist/maplibre-gl.js
curl.exe -O https://unpkg.com/maplibre-gl@2.4.0/dist/maplibre-gl.js.map

4-3. スタイルファイルの準備

スタイルですが、この地域では背景になるような地図がないのでVOXELだけのスタイルにしておきます。

style.json

{
  "version": 8,
  "center": [
    -91.993929,
    46.842568
  ],
  "zoom": 17,
  "layers": [
    {
      "id": "background",
      "type": "background",
      "paint": {
        "background-color": "rgb(40, 40, 40)"
      }
    },
    {
      "id": "voxel",
      "type": "fill-extrusion",
      "source": "voxel",
      "source-layer": "voxel",
      "paint": {
        "fill-extrusion-base": [
          "get",
          "h"
        ],
        "fill-extrusion-color": [
          "get",
          "color"
        ],
        "fill-extrusion-opacity": 1,
        "fill-extrusion-height": [
          "+",
          [
            "get",
            "h"
          ],
          [
            "get",
            "spacing"
          ]
        ]
      }
    }
  ],
  "sources": {
    "voxel": {
      "type": "vector",
      "attribution": "PointCloud from pierotofy/drone_dataset_brighton_beach",
      "minzoom": 14,
      "maxzoom": 19,
      "tiles": [
        "https://ubukawa.github.io/voxel-br/voxel/{z}/{x}/{y}.pbf"
      ]
    }
  }
}

4-4. ウェブ地図作成

index.htmlを以下のように準備します。（スタイルの切り替えボタンは今後必要に応じて拡張できます。）

index.html

<!doctype html>
<html>
<head>
<meta charset="utf-8" />
<meta name="viewport" content="width=device-width" />
<title></title>
<script type="module" src="module.js"></script>
</head>
<body>
<div id="map"/>
<div id="menu">
  <input id="style.json" type="radio" name="r" checked="checked" />
  <label>Voxel</label>
</div>
</body>
</html>

index.htmlが参照しているmodule.jsは以下のように書きます。

module.js

const style = href => {
  const e = document.createElement('link')
  e.href = href
  e.rel = 'stylesheet'
  document.head.appendChild(e)
}

const script = src => {
  const e = document.createElement('script')
  e.src = src
  document.head.appendChild(e)
}

const init = () => {
  style('style.css')
  style('maplibre-gl@2.4.0/maplibre-gl.css')
  script('maplibre-gl@2.4.0/maplibre-gl.js')

  const inputs = document.getElementById('menu').getElementsByTagName('input');
  for (const input of inputs) {
    input.onclick = (layer) => {
      map.setStyle(layer.target.id)
    }
  }
}
init()

let map
const showMap = async (texts) => {
  map = new maplibregl.Map({
    container: 'map',
    hash: true,
    style: 'style.json',
    maxZoom: 20,
    maxPitch: 85
  })
  map.addControl(new maplibregl.NavigationControl())
  map.addControl(new maplibregl.ScaleControl({
    maxWidth: 200, unit: 'metric'
  }))

  let voice = null
  for(let v of speechSynthesis.getVoices()) {
    console.log(v.name)
    if ([
      'Daniel',
      'Google UK English Male',
      'Microsoft Libby Online (Natural) - English (United Kingdom)'
    ].includes(v.name)) voice = v
  }

  map.on('load', () => {
    map.on('click', 'voxel', e => {
      let u = new SpeechSynthesisUtterance()
      u.lang = 'en-GB'
      u.text = 'a voxel of ' + e.features[0].properties.spacing + 'meters.'
      if (voice) u.voice = voice
      speechSynthesis.cancel()
      speechSynthesis.speak(u)
    })
    map.on('moveend', e => {
      let fs = map.queryRenderedFeatures(
	[window.innerWidth / 2, window.innerHeight / 2], 
        { layers : ['voxel'] }
      )
      if (fs.length == 0) return
      let height = fs[0].properties.h 
      console.log(height)
    })
  })
}

window.onload = () => {
  showMap()
}

style.cssも。

style.css

body { margin:0; padding:0; font-family: 'Roboto', sans-serif; color: #333333}
#map { top:0; height: 100vh; width:100vw; position: fixed; z-index: 0}
#menu { position: absolute; background: rgba(242, 242, 242, 0.5); padding: 10px; z-index: 2; }

これで、Web地図ができたはずです。

Step 5. 成果の確認

https://ubukawa.github.io/voxel-br/index.html にアクセスします。ズームレベルが大きいと、右ドラッグで傾けにくいので、ズームレベルを下げるか、表示サイズの割合を50%などにするとよいと思います。

例えば、ブラウザの表示を33%にして、ZL19まで見たボクセルがこんな感じです。拡大すると一部表示されない場合もあるような気がします。このあたりはまだ研究の余地があります。

まとめ

なんとかかんとか、VOXELタイル作成と、ウェブ地図での表示ができました。ズームレベルを上げると操作しにくくなるので、もう少し研究します。
（そもそもVoxelタイルは、地形概要をつかむのに使うのであまり高解像度を目指さない方がよいかもしれませんが。）

今回は、全体的な流れがわかるようになったので、今後、余力があればRakefileをつかってエレガントに処理できる方法も試してみたいと思います。

謝辞

optge/kid-cを開発された@hfuさんに感謝します。また、@hfuさんのレクチャーをYouTubeから公開してくれた青山学院大学の古橋研究室に感謝します。

参考　

PDAL　https://pdal.io/en/stable/
YouTube動画　点群LASファイルからボクセルVTのつくり方｜ UNVT　（古橋研究室,2021） https://www.youtube.com/watch?v=LrDk0VFodTE
参考にしたレポジトリ　https://github.com/optgeo/kid-c
今回の作業のレポジトリ　https://github.com/ubukawa/voxel-br