初めに
みなさん3DTilesをご存知でしょうか？
3DTilesは大規模な3D地理空間データをWebブラウザ上で高速に、なおかつ地図上で閲覧するための仕組みで、Cesiumという組織が先導して開発したフォーマットになります。
（https://cesium.com/blog/2015/08/10/introducing-3d-tiles/ より）
3DTilesは主にCesium.jsというライブラリを利用して、Google Earthのような地球儀の上に建物など地理空間情報を載せることのできるライブラリで、建築土木部門などで日本でも急速に利用者が増えています。
（というか、地図上に大規模な3Dのデータを載せようとすると、ほぼ一強です。）
そんな3DTilesは、py3dtilesというツールを利用するとコマンド一発でLASデータから変換可能ですが、データの一部が壊れていたりすると変換できません。
また、日本の測量時に利用される平面直角座標系と呼ばれる座標は縦軸がXで横軸がYです。
が、GIS（Cesiumでも同様だと思います）は一般的に数学座標を利用し、僕らが見慣れているようなXが横軸でYが縦軸として取り扱います。
座標系のＸ軸は、座標系原点において子午線に一致する軸とし、真北に向う値を正とし、座標系のＹ軸は、座標系原点において座標系のＸ軸に直交する軸とし、真東に向う値を正とする。
（https://www.gsi.go.jp/LAW/heimencho.html より）
一般的に日本で公開されている点群データ（特に、構造物など詳細な建物のデータなどではなく、航空機から取るようなデータの場合）はXYが逆なことが多いです。
このように、3DTiles作成作業というのはデータが壊れていて変換できなかったり、変換しても変な位置に表示されたりとハマりポイントが多いタスクなんですね。
今回はこれらの問題に一挙に対処していきます。
py3dtilesのダウンロード
まずは肝心のツールをダウンロードしていきましょう。
pypiからインストールすることができます。
pip install py3dtiles
このほかにも、LAStoolsというツールやpdalというツールを使いますので、必要に応じで以下を参考にインストールしてください。
- LASデータ操作のスイスアーミーナイフ「LAStools」をMacOSで利用したい
- iPhoneでも使えるようになったLiDARの標準ファイル形式「LAS」ってどんなデータなの？を、PDAL/Laspyを使って調べてみる。
データのダウンロード
今回はMy City Constructionという、工事データなどをオンライン上で登録することが可能なサイトから、静岡県さんが公開している点群データを利用していきましょう。
以下からダウンロードすることができます！
ダウンロードできたら、とりあえずpy3dtilesで変換してみましょう！
オプションで元データのEPSGコードがEPSG:6676であることを明示的に指定し、出力先はEPSG:4978であることを示しています。
元データのEPSGコードは<都道府県名> EPSGなどGoogleで調べると出てきますし、出力先は3DTiles標準ですので4978固定でOKです。
コマンドを実行するとなにやらエラーが出て終了してしまいました。
❯ py3dtiles convert --srs_in 6676 --srs_out 4978 --out 遊水地データ 遊水地データ.las
Error opening 遊水地データ.las. Skipping.
Incoherent header size
Traceback (most recent call last):
File "/Users/xxxx/opt/anaconda3/envs/jupyter/lib/python3.9/site-packages/py3dtiles/command_line.py", line 39, in main
convert.main(args)
File "/Users/xxxx/opt/anaconda3/envs/jupyter/lib/python3.9/site-packages/py3dtiles/convert.py", line 303, in main
return convert(args.files,
File "/Users/xxxx/opt/anaconda3/envs/jupyter/lib/python3.9/site-packages/py3dtiles/convert.py", line 397, in convert
infos['aabb'][0][0], infos['aabb'][0][1], infos['aabb'][0][2])))
TypeError: 'NoneType' object is not subscriptable
usage: py3dtiles [-h] [--verbose VERBOSE] {convert,info,merge,export} ...
Read/write 3dtiles files
positional arguments:
{convert,info,merge,export}
convert Convert .las files to a 3dtiles tileset.
info Extract informations from a 3DTiles file
merge Merge several pointcloud tilesets in 1 tileset
export Generate a tileset from a set of geometries
optional arguments:
-h, --help show this help message and exit
--verbose VERBOSE Print logs (-1: no logs, 0: progress indicator, 1+: increased verbosity) (default: 0)
Incoherent header sizeという文字が表示されていますが、これに関する情報がまっっっっっっっっっっっっったく存在せず、長い長い地獄のheader調査沼に落ちていくことになるのですが、この話はまた別の記事で…
ともあれ、lasデータのheaderに異常があり、データを変換できないことがわかりました。
ヘッダーの情報をほぼ維持したままデータを修正する
py3dtilesでは読み取れない（laspyなどでも読み取れない）ことはわかりましたが、pdalやLAStoolsなどでは読み込み可能です。
なので、メタ情報を見てみましょう。
❯ lasinfo 遊水地データ.las
WARNING: for LAS 1.3 header_size should at least be 235 but it is only 227
lasinfo (220613) report for '遊水地データ.las'
reporting all LAS header entries:
file signature: 'LASF'
file source ID: 0
global_encoding: 0
project ID GUID data 1-4: 00000000-0000-0000-0000-000000000000
version major.minor: 1.3
system identifier: 'FC'
generating software: 'TREND-POINT'
file creation day/year: 0/0
header size: 227
offset to point data: 227
number var. length records: 0
point data format: 2
point data record length: 26
number of point records: 2301317
number of points by return: 0 0 0 0 0
scale factor x y z: 0.000000044428938 0.000000040671672 0.000000002527102
offset x y z: -8955.918810275501528 -109804.676442435113131 3.5446009174646
min x y z: -8955.918810275501528 -109804.676442435113131 3.5446009174646
max x y z: -8860.508391625911827 -109717.334692435106263 8.971512104833167
start of waveform data packet record: 0
reporting minimum and maximum for all LAS point record entries ...
X 0 2147483646
Y 0 2147483647
Z 0 2147483647
intensity 0 0
return_number 0 0
number_of_returns 0 0
edge_of_flight_line 0 0
scan_direction_flag 0 0
classification 1 1
scan_angle_rank 0 0
user_data 0 0
point_source_ID 29 29
Color R 0 64005
G 0 64005
B 0 64005
number of first returns: 2301317
number of intermediate returns: 0
number of last returns: 2301317
number of single returns: 2301317
WARNING: there are 2301317 points with return number 0
WARNING: there are 2301317 points with a number of returns of given pulse of 0
histogram of classification of points:
2301317 unclassified (1)
上記の通り、以下の情報が読み取れます。
- lasのバージョンが1.3でpoint formatが2のデータ
-
WARNING: for LAS 1.3 header_size should at least be 235 but it is only 227というwarningが出ている- ヘッダーが8バイト足りてないらしい
ほかにも情報はたくさんありますが、現時点で必要な情報はこれだけです。
こちらを修正するためにものすごい労力をかけて調査したのですが、別記事にまとめるとして、結果的にはCloudCompareというツールをCLIから利用してlasをI/Oしてあげることで、ヘッダーが修正されました。
CloudCompareのダウンロード
こちらからダウンロード可能なので、ご自分のOSに合わせてインストールしてみましょう。
http://www.danielgm.net/cc/release/
インストールが終わったら以下の場所にCLIツール本体（実行ファイル）があることを確認しましょう。
（Windowsだと全然違う場所にあるかと思いますが、Windowsについてまとまった日本語記事の方がMacに関する記事よりも多いので、ググればすぐ見つけられると思います。）
❯ ll /Applications/CloudCompare.app/Contents/MacOS/ | grep CloudCompare
-rwxr-xr-x@ 1 xxxx admin 4.4M 5 29 05:00 CloudCompare*
確認できたら、お使いのシェルの設定にaliasを設定します。
echo "alias CloudCompare=/Applications/CloudCompare.app/Contents/MacOS/CloudCompare" >> ~/.zshrc
source ~/.zshrc
これでヘッダー修正の準備完了です！
CloudCompareのCLIを使ってヘッダーを修正
ヘッダーを修正というよりも、CloudCompareで読み込み後、普通にファイルを出力してあげることで、ヘッダーをいい感じに調整してくれるというだけの話なのですが、以下のコマンドで行います。
CloudCompare \
-SILENT \
-O "遊水地データ.las" \
-C_EXPORT_FMT LAS \
-SAVE_CLOUDS FILE "yusuichi_fix_header.las"
出力されたファイルを確認してみましょう。
❯ lasinfo yusuichi_fix_header.las
lasinfo (220613) report for 'yusuichi_fix_header.las'
reporting all LAS header entries:
file signature: 'LASF'
file source ID: 0
global_encoding: 0
project ID GUID data 1-4: 00000000-0000-0000-0000-000000000000
version major.minor: 1.3
system identifier: 'PDAL'
generating software: 'PDAL 2.3.0 (ff11ae)'
file creation day/year: 179/2022
header size: 235
offset to point data: 235
number var. length records: 0
point data format: 2
point data record length: 26
number of point records: 2301317
number of points by return: 2301317 0 0 0 0
scale factor x y z: 0.0000001 0.0000001 0.0000001
offset x y z: -8955.918810275499709 -109804.676442434996716 3.5446009174646
min x y z: -8955.9189453 -109804.6796875 3.5446010
max x y z: -8860.5087891 -109717.3359375 8.9715118
WARNING: stored resolution of min_x not compatible with x_offset and x_scale_factor: -8955.9189453125
WARNING: stored resolution of min_y not compatible with y_offset and y_scale_factor: -109804.6796875
WARNING: stored resolution of min_z not compatible with z_offset and z_scale_factor: 3.544600963592529
WARNING: stored resolution of max_x not compatible with x_offset and x_scale_factor: -8860.5087890625
WARNING: stored resolution of max_y not compatible with y_offset and y_scale_factor: -109717.3359375
WARNING: stored resolution of max_z not compatible with z_offset and z_scale_factor: 8.971511840820312
start of waveform data packet record: 0
reporting minimum and maximum for all LAS point record entries ...
X -1350 954100212
Y -32451 873405049
Z 0 54269109
intensity 0 0
return_number 1 1
number_of_returns 1 1
edge_of_flight_line 0 0
scan_direction_flag 0 0
classification 1 1
scan_angle_rank 0 0
user_data 0 0
point_source_ID 29 29
Color R 0 64000
G 0 64000
B 0 64000
WARNING: 2 points outside of header bounding box
number of first returns: 2301317
number of intermediate returns: 0
number of last returns: 2301317
number of single returns: 2301317
overview over number of returns of given pulse: 2301317 0 0 0 0 0 0
histogram of classification of points:
2301317 unclassified (1)
WARNING: real min y smaller than header min y by 0.000000
WARNING: real min z smaller than header min z by 0.000000
「ヘッダーが不整合やで」ってエラーが消えてますね！
scaleなどを微調整する
その他、ちょこちょこwarningが出ているので、潰してしまいましょう。
las2las -i yusuichi_fix_header.las -rescale 0.0001 0.0001 0.0001 -o yusuichi_rescale.las
lasinfo -i yusuichi_rescale.las -repair
warningがすべて無くなりましたね。
❯ lasinfo yusuichi_rescale.las
lasinfo (220613) report for 'yusuichi_rescale.las'
reporting all LAS header entries:
file signature: 'LASF'
file source ID: 0
global_encoding: 0
project ID GUID data 1-4: 00000000-0000-0000-0000-000000000000
version major.minor: 1.3
system identifier: 'LAStools (c) by rapidlasso GmbH'
generating software: 'las2las (version 220613)'
file creation day/year: 179/2022
header size: 235
offset to point data: 235
number var. length records: 0
point data format: 2
point data record length: 26
number of point records: 2301317
number of points by return: 2301317 0 0 0 0
scale factor x y z: 0.0001 0.0001 0.0001
offset x y z: -8955.918810275499709 -109804.676442434996716 3.5446009174646
min x y z: -8955.9189 -109804.6796 3.5446
max x y z: -8860.5088 -109717.3359 8.9715
start of waveform data packet record: 0
reporting minimum and maximum for all LAS point record entries ...
X -1 954100
Y -32 873405
Z 0 54269
intensity 0 0
return_number 1 1
number_of_returns 1 1
edge_of_flight_line 0 0
scan_direction_flag 0 0
classification 1 1
scan_angle_rank 0 0
user_data 0 0
point_source_ID 29 29
Color R 0 64000
G 0 64000
B 0 64000
number of first returns: 2301317
number of intermediate returns: 0
number of last returns: 2301317
number of single returns: 2301317
overview over number of returns of given pulse: 2301317 0 0 0 0 0 0
histogram of classification of points:
2301317 unclassified (1)
XYを入れ替える
ヘッダーが修正できましたが、このまま3dtilesに変換しても座標がおかしいです。
py3dtiles convert --srs_in 6676 --srs_out 4978 --out yusuichi_rescale yusuichi_rescale.las
このように、岐阜県の山奥に表示されてしまいました。
なのでpdalを利用してXとY座標を入れ替えていきます。
（余談ですが、このようにlasを操作するときはいろんなツールを使いこなさないといけないので結構辛いです。）
まずは、xy_switch_pipeline.jsonのような名前でjsonを作成します。
[
{
"type": "readers.las",
"filename": "rescale.las",
"spatialreference": "EPSG:6676"
},
{
"type": "filters.reprojection",
"in_srs": "EPSG:6676",
"out_srs": "EPSG:6676",
"in_axis_ordering": "2, 1"
},
{
"type": "writers.las",
"filename": "translated.las",
"forward": "header,scale,vlr",
"offset_x": "auto",
"offset_y": "auto",
"offset_z": "auto"
}
]
これはpdalのpipeline(https://pdal.io/apps/pipeline.html )という機能で処理をするためのフローを記載したjsonになります。
このJSONを利用すると、以下のような処理を行うことができます。
- lasの読み込み・投影法の設定
- XとY軸の入れ替え
- ヘッダーとスケールと可変長レコードを維持したままファイルの吐き出し（offsetは軸を入れ替えるので自動で調整）
以下のように利用することで、jsonに記載のファイル名を上書きして、動的に処理することができます。
pdal pipeline xy_switch_pipeline.json \
--readers.las.filename=yusuichi_rescale.las \
--writers.las.filename=yusuichi_switching_rows.las
3DTilesに変換してCesiumで表示
これで最初の障害であったヘッダーの問題とXY座標の問題が解決されたので、いよいよ3DTilesに変換することができます！
以下のようにして変換していきましょう（上の方で一度利用したコマンドと同じです。）
py3dtiles convert --srs_in 6676 --srs_out 4978 --out yusuichi yusuichi_switching_rows.las
しばらくするとフォルダが出来上がります。
その後、Cesiumで閲覧するために、こんなHTMLファイルを用意し、カレントディレクトリに配置しましょう。
<!DOCTYPE html>
<html>
<head>
<meta charset="UTF-8">
<title>Cesium</title>
<script src="https://cesium.com/downloads/cesiumjs/releases/1.72/Build/Cesium/Cesium.js"></script>
<link href="https://cesium.com/downloads/cesiumjs/releases/1.72/Build/Cesium/Widgets/widgets.css" rel="stylesheet">
<style>
#cesiumContainer {
position: absolute;
top: 0;
left: 0;
height: 100%;
width: 100%;
margin: 0;
overflow: hidden;
padding: 0;
font-family: sans-serif;
}
html {
height: 100%;
}
body {
padding: 0;
margin: 0;
overflow: hidden;
heght: 100%;
}
</style>
</head>
<body>
<div id="cesiumContainer"></div>
<script>
const raster_tile = new Cesium.OpenStreetMapImageryProvider({
url: 'https://cyberjapandata.gsi.go.jp/xyz/std/',
credit: new Cesium.Credit('地理院タイル', '', 'https://maps.gsi.go.jp/development/ichiran.html'),
maximumLevel: 18,
})
const viewer = new Cesium.Viewer('cesiumContainer', {
imageryProvider: raster_tile,
baseLayerPicker: false,
geocoder: false,
homeButton: false
})
viewer.camera.setView({
destination : Cesium.Cartesian3.fromDegrees(140.00, 36.14, 20000000.0)
})
const tileset = viewer.scene.primitives.add(
new Cesium.Cesium3DTileset({
url: './yusuichi/tileset.json'
})
)
tileset.style = new Cesium.Cesium3DTileStyle({
pointSize: 5
})
viewer.zoomTo(tileset)
</script>
</body>
</html>
この中に以下のようなコードがありますが、ここにtileset.jsonと呼ばれる3DTilesの核となるファイルを指定することで地図上に可視化することができます。
const tileset = viewer.scene.primitives.add(
new Cesium.Cesium3DTileset({
url: './yusuichi/tileset.json'
})
)
CORSの関係でローカルのファイルをそのまま読むことができませんので、別のターミナルを開いてカレントディレクトリでpythonの簡易サーバーを起動しましょう。
❯ python -m http.server 8000
Serving HTTP on :: port 8000 (http://[::]:8000/) ...
ブラウザでlocalhost:8000を開きましょう。
道路沿いにピッタリあった形で正しい位置に表示されました！
（ちなみに、遊水地とは河川が氾濫しそうな時に水を貯めとく池みたいなもんです。なので河川沿いに存在します。）
終わりに
ということで諸々細かい調整をして3DTilesを作成してきました。
丁寧に解説された情報はあまり見つかりませんが、現場では割と遭遇する事案で、多分いろんな方がハマっているんだと思いますが、この辺りの作業を行うことで、大体解決するとおもいます。
みなさんもどんどん3DTilesを活用していきましょう！