※上の画像は掛川城の点群データをCesiumというGoogleEarthのようなライブラリを用いてJavaScriptで可視化した「掛川城の点群データってこんな感じだよー」っていうのを示したイメージ画像になります。（この記事ではJSの話までは取り扱いません）
※点群データがどんなものか確認したい方はデモサイトを用意していますので、こちらからどうぞ！
https://demo-3d.mierune.io/
las/lazデータとは
みなさん、lasデータってご存知でしょうか？
主に三次元のポイントデータの集合の交換形式として利用されているバイナリフォーマットで、主にiPhone12などから導入され、手軽に利用できるようになったLiDERデータの標準として利用されており、lazはそのlasデータの圧縮形式になります。
lasデータの仕様: https://www.asprs.org/wp-content/uploads/2010/12/LAS_1_4_r13.pdf
Wikipedia: https://en.wikipedia.org/wiki/LAS_file_format
と…めんどくさい言い方をしましたが、いわゆる点群・ポイントクラウドデータのことであり、土木・建築の分野のみならず、デジタルシティやメタバースなどの領域でも近年一般的に利用されるようになっています。
工事の分野ではドローンを飛ばして点群データを作成し日々の進捗を確認したり、現実にある物体や人物をスキャンしてAR/VR空間に持って行ってみたり、身近なところだと引越し前に新しい住居をスキャンし、適切に家具を配置したりなど…実はその辺で利用されています。
そんな点群データですが、国内では静岡県がめちゃくちゃ整備しておりまして、なんとオープンデータとして誰でも利用できるようになっています。
点群データ公開プラットフォームであるShizuoka Point Cloud DBを整備したり、県で発注した工事のデータはMy City Constructionというサイトを利用してオンラインで保存するようにしていたり…
Shizuoka Point Cloud DB
My City Construction
この中から今回は掛川城の点群データを利用させてもらって、Pythonで可視化していきましょう！
データのダウンロード
Shizuoka Point Cloud DBのデータは一部手軽に利用できるようにG空間情報センターというところにアップロードされているので、こちらからダウンロードしましょう。
掛川城オープンデータ化プロジェクト
詳細ページに飛び…
ダウンロードボタンよりダウンロードしましょう！
ファイルは圧縮後で2GB、解凍すると5GBを超えますので、注意してダウンロードしましょう！
※macで開こうとすると、文字コードの問題で普通には解凍できないかもしれません。
こちら（The Unarchiver）のアプリなどで解凍すると開くと思います。
ダウンロード、解凍が終わったら早速点群を弄ってみましょう。
jupyterの環境構築
今回はPythonのみで可視化していきますので、jupyterなどを利用すると便利です。
こちらを参考に環境構築しましょう。
- MacでGISデータ分析を始めるためにサクッとAnacondaとjupyter labをインストールしてみる
- jupyter labで最低限これだけは入れとけっていう拡張機能の紹介！
- Pythonで地理空間情報（GIS）やるために必要なパッケージ全部入りの「geospatial」が便利すぎた
Leafmapとは
こちらに詳しい利用方法などまとめていますのでご覧ください！
Pythonのコードを1行書くだけで誰でも手軽にインタラクティブな地図アプリを作れるLeafmapを使ってみよう
pythonでlasデータを見てみる
まずはパッケージをインポートします
import leafmap
import numpy as np
import open3d as o3d
次にlasデータをleafmapで読み込んでいきましょう。
Pythonならたったこれだけのコードで点群データを読み込めますし、5GBのデータも5秒以内に読み取れました！
%%time
filename = "20190308掛川城.las"
las = leafmap.read_lidar(filename)
# CPU times: user 1.13 s, sys: 2.48 s, total: 3.61 s
# Wall time: 4.03 s
ポイントの数を確認してみるとなんと1.9億ポイントあります。
lasの内部はメタデータ付きの表形式のような構造なので、2億レコード弱のデータと考えると、5GBも納得ですね。
las.header.point_count
# 192366079
Leafmapでlasデータを読み取るとlaspyというパッケージのインスタンスになるようです。
type(las)
# laspy.lasdata.LasData
lasデータはバージョンにもよりますが、このようなプロパティを持っています。
# lasデータが持っている属性を表示
list(las.point_format.dimension_names)
# ['X',
# 'Y',
# 'Z',
# 'intensity',
# 'return_number',
# 'number_of_returns',
# 'scan_direction_flag',
# 'edge_of_flight_line',
# 'classification',
# 'synthetic',
# 'key_point',
# 'withheld',
# 'scan_angle_rank',
# 'user_data',
# 'point_source_id',
# 'red',
# 'green',
# 'blue']
各列はnumpyのndarryaになっているようで、1次元の配列です。
# 各列はndarrayで保持
type(las.X)
# numpy.ndarray
las.X.shape
# (192366079,)
Leafmapで可視化
ではこのデータを可視化していきましょう！可視化にはLeafmapとopen3dを利用します！
Leafmapを利用すればたった1行で可視化できます。
（裏ではOpen3Dというライブラリが動いています）
%%time
leafmap.view_lidar("20190308掛川城.las", backend="open3d")
# CPU times: user 2min 32s, sys: 15.4 s, total: 2min 48s
# Wall time: 5min 35s
この通り！
と、言いたいところなのですが、あまりにもデータが大きすぎて読み込みに5分かかっている上、重たくてまともにプレビューが動かせません…
その上、高度に応じて勝手に色が決められてしまっていますね。
なので、点群を間引いたり、別の表示方法も試していきましょう！
点群を間引く
現状、読み込んだ点群のXYZ座標は以下のように1次元の配列、つまり通常のリストのような形式になっています。
[0, 1, 2, 3] # X座標
[4, 5, 6, 7] # Y座標
[8, 9, 10, 11] # Z座標
まずはこれをまとめて2次元配列にしましょう。
numpyのvstackメソッドを使えば簡単にまとめられます。
こちらの点群は色情報（RGB）を持っているので、そちらもまとめちゃいましょう。
points = np.vstack((las.X, las.Y, las.Z))
colors = np.vstack((las.red, las.green, las.blue))
ただし、このままだとカラムが1.9億列で、3行の2次元配列になってしまいます。
# こんなイメージ
array([[0, 1, 2, 3...x],
[0, 1, 2, 3...y],
[0, 1, 2, 3...z]])
XとYとZの3カラムを持つポイントが、1.9億レコードあるのが望ましい姿なので、行と列を入れ替えましょう。
今度はnumpyのtransposeメソッドを利用していきます。
points = points.transpose()
colors = colors.transpose()
そうすると2次元配列はこのようになります
# こんなイメージ
array([[0, 0, 0],
[1, 1, 1],
[2, 2, 2],
[3, 3, 3],
...
[x, y, z]])
こうすることでレコード（ポイント）がXYZの座標をもつ構造になりましたね。
実行するときは下のように配列の合成と行列の置換をまとめてやってしまいましょう。
%%time
points = np.vstack((las.X, las.Y, las.Z)).transpose()
colors = np.vstack((las.red, las.green, las.blue)).transpose()
# CPU times: user 2.3 s, sys: 832 ms, total: 3.13 s
# Wall time: 3.18 s
さて、これで二次元配列の中には要素が1.9億個ある状態になりました。
これを1/100に間引いていきましょう。
ただの配列なので、以下のようにブラケットでサクッと間引けます。
%%time
step = 100
decimated_points = points[::step]
decimated_colors = colors[::step]
# CPU times: user 7 µs, sys: 0 ns, total: 7 µs
# Wall time: 10 µs
配列の長さを確認すると、ちゃんと間引けて190万ポイントになっていることがわかります。
%%time
len(decimated_points)
# CPU times: user 3 µs, sys: 1e+03 ns, total: 4 µs
# Wall time: 5.96 µs
# 1923661
ポイントを間引いたら、open3dという3Dデータ分析用のパッケージを利用して点群のジオメトリを作成していきます。
ジオメトリ作成後、ポイントとカラーの情報を投入していきましょう。
%%time
# open3dで点群用のジオメトリを作成する
point_cloud = o3d.geometry.PointCloud()
# ポイント属性とカラー属性に配列を格納する
point_cloud.points = o3d.utility.Vector3dVector(decimated_points)
point_cloud.colors = o3d.utility.Vector3dVector(decimated_colors / 65535)
# CPU times: user 3.46 s, sys: 49.4 ms, total: 3.51 s
# Wall time: 3.5 s
Open3Dで可視化する
いよいよ可視化していきましょう！
今度はカラーの情報をちゃんとopen3dのインスタンスに組み込んだので点群が本来持っている色で可視化されます。
また、ファイルサイズも1/100に軽量化されているので、プレビューもサクサク動きます！
Leafmapのインスタンスでは無くなってしまったので、Open3Dを直接利用して可視化していきましょう。こちらも1行で可視化できます！
%%time
o3d.visualization.draw_geometries([point_cloud])
# CPU times: user 3.18 s, sys: 3.92 s, total: 7.1 s
# Wall time: 7.1 s
出てきましたね！
さて、この段階で、一度plyファイルに吐き出しておきましょう。
ここで一つ注意点があります。
open3dはlas形式で出力することができません。
http://www.open3d.org/docs/0.9.0/tutorial/Basic/file_io.html
というより、その他多くのパッケージがlasでの出力をサポートしていないように見えます。
おそらく、lasがheaderに複雑なメタ情報を持っているためかと思われますが、詳しくはわかっていません。
ただ、逆に言えば多くのパッケージが、plyの読み込み、書き出しに対応していますので、plyに変換しておいて困ることはほとんどありません。またplyファイルはblenderでそのまま読み込むことができるので、何かと便利です。
%%time
# 多分、lasをスムーズに読み取れるのはlaspyくらいしかないと思う
o3d.io.write_point_cloud("decimated_kakegawa.ply", point_cloud)
# CPU times: user 170 ms, sys: 144 ms, total: 314 ms
# Wall time: 314 ms
# True
これでplyファイルが出力されましたね。お疲れ様でした！
終わりに
ということで、Leafmap/Open3Dを使って点群データ（las）を可視化する方法を紹介していきました！
が、加工したデータはlasで出力することができませんでした。
plyファイルも汎用的に利用できる便利なフォーマットなので、このままでも問題ないといえばないのですが、lasで出力したい時もあると思いますので、そちらに関しては別の記事にまとめたいと思います！