ヘッダーが壊れているLASファイルを修正して、XY座標を入れ替えた上でpy3dtilesで3DTilesを作成する

## 初めに

みなさん3DTilesをご存知でしょうか？
3DTilesは大規模な3D地理空間データをWebブラウザ上で高速に、なおかつ地図上で閲覧するための仕組みで、Cesiumという組織が先導して開発したフォーマットになります。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/203944/c0d9e9fa-8db8-9bf7-b14a-50c93846589f.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/203944/81f5f953-3a0e-60ff-ca19-86053a20b01f.png)

（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からインストールすることができます。

```python
pip install py3dtiles
```

このほかにも、LAStoolsというツールやpdalというツールを使いますので、必要に応じで以下を参考にインストールしてください。

- [LASデータ操作のスイスアーミーナイフ「LAStools」をMacOSで利用したい](https://qiita.com/nokonoko_1203/items/183317f35ca1f07d311c)
- [iPhoneでも使えるようになったLiDARの標準ファイル形式「LAS」ってどんなデータなの？を、PDAL/Laspyを使って調べてみる。](https://qiita.com/nokonoko_1203/items/8b8568545148acfe11d3)

## データのダウンロード

今回は[My City Construction](https://mycityconstruction.jp/)という、工事データなどをオンライン上で登録することが可能なサイトから、静岡県さんが公開している点群データを利用していきましょう。
以下からダウンロードすることができます！

https://mycityconstruction.jp/products/36/folders/view/ICON

ダウンロードできたら、とりあえずpy3dtilesで変換してみましょう！
オプションで元データのEPSGコードが`EPSG:6676`であることを明示的に指定し、出力先は`EPSG:4978`であることを示しています。
元データのEPSGコードは`<都道府県名> EPSG`などGoogleで調べると出てきますし、出力先は3DTiles標準ですので`4978`固定でOKです。

コマンドを実行するとなにやらエラーが出て終了してしまいました。

```python
❯ 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などでは読み込み可能です。
なので、メタ情報を見てみましょう。

```bash
❯ 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に関する記事よりも多いので、ググればすぐ見つけられると思います。）

```bash
❯ ll /Applications/CloudCompare.app/Contents/MacOS/ | grep CloudCompare

-rwxr-xr-x@ 1 xxxx  admin   4.4M  5 29 05:00 CloudCompare*
```

確認できたら、お使いのシェルの設定にaliasを設定します。

```bash
echo "alias CloudCompare=/Applications/CloudCompare.app/Contents/MacOS/CloudCompare" >> ~/.zshrc
source ~/.zshrc
```

これでヘッダー修正の準備完了です！

## CloudCompareのCLIを使ってヘッダーを修正

ヘッダーを修正というよりも、CloudCompareで読み込み後、普通にファイルを出力してあげることで、ヘッダーをいい感じに調整してくれるというだけの話なのですが、以下のコマンドで行います。

```bash
CloudCompare \
-SILENT \
-O "遊水地データ.las" \
-C_EXPORT_FMT LAS \
-SAVE_CLOUDS FILE "yusuichi_fix_header.las"
```

出力されたファイルを確認してみましょう。

```bash
❯ 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が出ているので、潰してしまいましょう。

```bash
las2las -i yusuichi_fix_header.las -rescale 0.0001 0.0001 0.0001 -o yusuichi_rescale.las 

lasinfo -i yusuichi_rescale.las -repair
```

warningがすべて無くなりましたね。

```bash
❯ 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に変換しても座標がおかしいです。

```bash
py3dtiles convert --srs_in 6676 --srs_out 4978 --out yusuichi_rescale yusuichi_rescale.las 
```

このように、岐阜県の山奥に表示されてしまいました。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/203944/103b5cae-ca80-c4bc-0cf8-e5b9c88422bc.png)

なのでpdalを利用してXとY座標を入れ替えていきます。
（余談ですが、このようにlasを操作するときはいろんなツールを使いこなさないといけないので結構辛いです。）

まずは、`xy_switch_pipeline.json`のような名前でjsonを作成します。

```bash
[
  {
    "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に記載のファイル名を上書きして、動的に処理することができます。

```bash
pdal pipeline xy_switch_pipeline.json \
--readers.las.filename=yusuichi_rescale.las \
--writers.las.filename=yusuichi_switching_rows.las
```

## 3DTilesに変換してCesiumで表示

これで最初の障害であったヘッダーの問題とXY座標の問題が解決されたので、いよいよ3DTilesに変換することができます！
以下のようにして変換していきましょう（上の方で一度利用したコマンドと同じです。）

```bash
py3dtiles convert --srs_in 6676 --srs_out 4978 --out yusuichi yusuichi_switching_rows.las 
```

しばらくするとフォルダが出来上がります。

その後、Cesiumで閲覧するために、こんなHTMLファイルを用意し、カレントディレクトリに配置しましょう。

```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の核となるファイルを指定することで地図上に可視化することができます。

```js
const tileset = viewer.scene.primitives.add(
    new Cesium.Cesium3DTileset({
        url: './yusuichi/tileset.json'
    })
)
```

CORSの関係でローカルのファイルをそのまま読むことができませんので、別のターミナルを開いてカレントディレクトリでpythonの簡易サーバーを起動しましょう。

```bash
❯ python -m http.server 8000
Serving HTTP on :: port 8000 (http://[::]:8000/) ...
```

ブラウザで`localhost:8000`を開きましょう。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/203944/c4e11968-28e1-a2de-5f4a-13de40faf251.png)

道路沿いにピッタリあった形で正しい位置に表示されました！
（ちなみに、遊水地とは河川が氾濫しそうな時に水を貯めとく池みたいなもんです。なので河川沿いに存在します。）

## 終わりに

ということで諸々細かい調整をして3DTilesを作成してきました。

丁寧に解説された情報はあまり見つかりませんが、現場では割と遭遇する事案で、多分いろんな方がハマっているんだと思いますが、この辺りの作業を行うことで、大体解決するとおもいます。

みなさんもどんどん3DTilesを活用していきましょう！