長崎県全域の点群データが自由に使えるらしいからみんなで使っちゃおうぜ！？

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/203944/5a3c706f-8d25-2a59-1f61-a7199a7c7b31.png)

（※[オープンナガサキ](https://opennagasaki.nerc.or.jp/)のトップページ）

くぅ〜〜〜〜痺れるトップページですね〜。

点群データ先進国家「SHIZUOKA」に引き続き長崎県も全域の点群データを公開してくれるそうです！
(一部地域はまだ調整中ですが…)

https://twitter.com/nokonoko_1203/status/1635458469137625088?s=20

（VIRTUAL SHIZUOKAの点群データは[こちら](https://www.geospatial.jp/ckan/organization/shizuokapref?q=&organization=shizuokapref&sort=views_recent+desc))

こうやってデータをバシバシ出してくれている以上は、我々エンジニアとしては有効活用させていただかねばなりません！！

と言うことで早速データをダウンロードしてGoogle Earthみたいな3D地球儀を作成できるJavaScriptライブラリである「[Cesium](https://www.cesium.com/)」を利用して地球上に載せてみました！

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/203944/c3377175-6fad-a590-c794-6b38a0c2505c.png)

（※この記事はPDALやpy3dtiles・Cesiumなど点群操作に慣れている人向けプログラマー向けですが、そうでない方も出来るようになるべく丁寧に記載してはいます。）

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

何はともあれまずはデータをダウンロードしましょう！

[オープンナガサキ](https://opennagasaki.nerc.or.jp/)のページに移動しスクロールしていくと「ダウンロードページ」のリンクが出てきます！

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/203944/6b75b464-d442-dc09-1967-02194fe23b78.png)

クリックすると利用規約が出てきます。
> データの座標参照系は、世界測地系(JGD2011)／平面直角座標系第1系です。

ここ、テストに出るので覚えておいてくださいね。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/203944/fd2f4333-d17f-843d-9813-a65f8dbb45e1.png)

ダウンロードページへ進むと、良い感じに長崎にフォーカスされた地図が出てきますね。
地図から公開データがある箇所とない箇所がわかります。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/203944/3fb40e18-35f3-1ced-f02a-d913d1cf382f.png)

画面左の+ボタンを押してズームしていってみましょう。
すると以下のような8文字からなる英数字がいくつも表示されると思います。
（この英数字は`国土基本図図郭`と呼ばれる図郭です。参考: [国土基本図図郭とは｜図郭コードの読み方](https://club.informatix.co.jp/?p=1293)）

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/203944/6216fbf3-e478-c66e-7055-e80652b0035c.png)

これらのどれか一個をクリックするとサイドバーに表示されるので「ダウンロード」をクリックしましょう。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/203944/f54ad1fa-e822-b3cd-0a51-013a6829093c.png)

重たいのでしばらく待ちます…

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/203944/e532d439-2960-8c3d-657e-6d532b35f454.png)

`01KE9821.zip`のようなzipファイルはダウンロードできましたか？
これを解凍すると`01ke9821_org.las`のようなファイルが出てくると思いますが、これが点群データ本体になります。

## PDALでメタデータの確認・CRSの付与・XY座標の入れ替え

ダウンロードしてきた点群データを[PDAL](https://pdal.io/en/latest/)というツールを利用して確認・調整していきます。
インストールの方法は[こちら](https://pdal.io/en/latest/quickstart.html)を参照してください。

以下のように点群データの確認ができます。
ファイルサイズや点群の数などいろんな情報が表示されます。
（`"count": 18175075`から分かる通り点群は数万〜数億のポイントを持つため、基本的にはファイルサイズがとても大きいです。）

```bash
% pdal info --metadata 01ke9821_org.las
{
  "file_size": 472552177,
  "filename": "01ke9821_org.las",
  "metadata":
  {
    "comp_spatialreference": "",
    "compressed": false,
    "copc": false,
    "count": 18175075,
    "creation_doy": 0,
    "creation_year": 0,
    "dataformat_id": 2,
    "dataoffset": 227,
    "filesource_id": 0,
    "global_encoding": 0,
    "global_encoding_base64": "AAA=",
    "header_size": 227,
    "major_version": 1,
    "maxx": 34999.994,
    "maxy": -27000.001,
    "maxz": 100.905,
    "minor_version": 2,
    "minx": 34000,
    "miny": -27750,
    "minz": -2.807,
    "offset_x": 34000,
    "offset_y": -27750,
    "offset_z": -2.807,
    "point_length": 26,
    "project_id": "00000000-0000-0000-0000-000000000000",
    "scale_x": 4.65658493556854e-07,
    "scale_y": 3.49245499982147e-07,
    "scale_z": 4.82946634517492e-08,
    "software_id": "TREND-POINT",
    "spatialreference": "",
    "srs":
    {
      "compoundwkt": "",
      "horizontal": "",
      "isgeocentric": false,
      "isgeographic": false,
      "json": "",
      "prettycompoundwkt": "",
      "prettywkt": "",
      "proj4": "",
      "units":
      {
        "horizontal": "unknown",
        "vertical": ""
      },
      "vertical": "",
      "wkt": ""
    },
    "system_id": "FC"
  },
  "now": "2023-03-15T14:17:30+0900",
  "pdal_version": "2.5.0 (git-version: Release)",
  "reader": "readers.las"
}
```

次に、「座標系」を設定していきます。
点群データ取得時（つまり、航空機などから地上に向けてレーザーを放ってデータを取得する時）には特定の「座標系」に従って計測を行い、その座標系に沿った値が格納されるのですが、「座標系の定義」自体は付与されていないこともあります。
（座標系とは: [座標系とは｜座標系の種類・平面直角座標系について解説](https://club.informatix.co.jp/?p=708)）

今回のデータは明らかに付与されていないので、PDALを利用して付与していきましょう！

以下のサイトを参考にすると「世界測地系(JGD2011)／平面直角座標系第1系」のEPSGコードは`EPSG:6669`なので、こちらを付与していきます。

（参考: [EPSGコード一覧表/日本でよく利用される空間座標系（座標参照系）](https://lemulus.me/column/epsg-list-gis)）

```bash
% pdal translate -i 01ke9821_org.las -o 01ke9821_translated.las --writers.las.a_srs="EPSG:6669"
```

付与されているか確認しましょう。

```bash
%  pdal info --metadata 01ke9821_translated.las
{
  "file_size": 617952997,
  "filename": "01ke9821_translated.las",
  "metadata":
  {
    "comp_spatialreference": "PROJCS[\"JGD2011 / Japan Plane Rectangular CS I\",GEOGCS[\"JGD2011\",DATUM[\"Japanese_Geodetic_Datum_2011\",SPHEROID[\"GRS 1980\",6378137,298.257222101,AUTHORITY[\"EPSG\",\"7019\"]],AUTHORITY[\"EPSG\",\"1128\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"6668\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",33],PARAMETER[\"central_meridian\",129.5],PARAMETER[\"scale_factor\",0.9999],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Northing\",NORTH],AXIS[\"Easting\",EAST],AUTHORITY[\"EPSG\",\"6669\"]]",
    "compressed": false,
    "copc": false,
    "count": 18175075,
    "creation_doy": 0,
    "creation_year": 2023,
    "dataformat_id": 3,
    "dataoffset": 447,
    "filesource_id": 0,
    "global_encoding": 0,
    "global_encoding_base64": "AAA=",
    "gtiff": "Geotiff_Information:\n   Version: 1\n   Key_Revision: 1.0\n   Tagged_Information:\n      End_Of_Tags.\n   Keyed_Information:\n      GTModelTypeGeoKey (Short,1): ModelTypeProjected\n      GTRasterTypeGeoKey (Short,1): RasterPixelIsArea\n      GTCitationGeoKey (Ascii,39): \"JGD2011 / Japan Plane Rectangular CS I\"\n      GeogCitationGeoKey (Ascii,8): \"JGD2011\"\n      GeogAngularUnitsGeoKey (Short,1): Angular_Degree\n      ProjectedCSTypeGeoKey (Short,1): Code-6669 (JGD2011 / Japan Plane Rectangular CS I)\n      ProjLinearUnitsGeoKey (Short,1): Linear_Meter\n      End_Of_Keys.\n   End_Of_Geotiff.\n",
    "header_size": 227,
    "major_version": 1,
    "maxx": 34999.994,
    "maxy": -27000.001,
    "maxz": 100.9050001,
    "minor_version": 2,
    "minx": 34000,
    "miny": -27750,
    "minz": -2.807,
    "offset_x": 0,
    "offset_y": 0,
    "offset_z": 0,
    "point_length": 34,
    "project_id": "00000000-0000-0000-0000-000000000000",
    "scale_x": 0.01,
    "scale_y": 0.01,
    "scale_z": 0.01,
    "software_id": "PDAL 2.5.0 (Releas)",
    "spatialreference": "PROJCS[\"JGD2011 / Japan Plane Rectangular CS I\",GEOGCS[\"JGD2011\",DATUM[\"Japanese_Geodetic_Datum_2011\",SPHEROID[\"GRS 1980\",6378137,298.257222101,AUTHORITY[\"EPSG\",\"7019\"]],AUTHORITY[\"EPSG\",\"1128\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"6668\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",33],PARAMETER[\"central_meridian\",129.5],PARAMETER[\"scale_factor\",0.9999],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Northing\",NORTH],AXIS[\"Easting\",EAST],AUTHORITY[\"EPSG\",\"6669\"]]",
    "srs":
    {
      "compoundwkt": "PROJCS[\"JGD2011 / Japan Plane Rectangular CS I\",GEOGCS[\"JGD2011\",DATUM[\"Japanese_Geodetic_Datum_2011\",SPHEROID[\"GRS 1980\",6378137,298.257222101,AUTHORITY[\"EPSG\",\"7019\"]],AUTHORITY[\"EPSG\",\"1128\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"6668\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",33],PARAMETER[\"central_meridian\",129.5],PARAMETER[\"scale_factor\",0.9999],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Northing\",NORTH],AXIS[\"Easting\",EAST],AUTHORITY[\"EPSG\",\"6669\"]]",
      "horizontal": "PROJCS[\"JGD2011 / Japan Plane Rectangular CS I\",GEOGCS[\"JGD2011\",DATUM[\"Japanese_Geodetic_Datum_2011\",SPHEROID[\"GRS 1980\",6378137,298.257222101,AUTHORITY[\"EPSG\",\"7019\"]],AUTHORITY[\"EPSG\",\"1128\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"6668\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",33],PARAMETER[\"central_meridian\",129.5],PARAMETER[\"scale_factor\",0.9999],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Northing\",NORTH],AXIS[\"Easting\",EAST],AUTHORITY[\"EPSG\",\"6669\"]]",
      "isgeocentric": false,
      "isgeographic": false,
      "json": "{\n  \"type\": \"ProjectedCRS\",\n  \"name\": \"JGD2011 / Japan Plane Rectangular CS I\",\n  \"base_crs\": {\n    \"name\": \"JGD2011\",\n    \"datum\": {\n      \"type\": \"GeodeticReferenceFrame\",\n      \"name\": \"Japanese Geodetic Datum 2011\",\n      \"ellipsoid\": {\n        \"name\": \"GRS 1980\",\n        \"semi_major_axis\": 6378137,\n        \"inverse_flattening\": 298.257222101\n      }\n    },\n    \"coordinate_system\": {\n      \"subtype\": \"ellipsoidal\",\n      \"axis\": [\n        {\n          \"name\": \"Geodetic latitude\",\n          \"abbreviation\": \"Lat\",\n          \"direction\": \"north\",\n          \"unit\": \"degree\"\n        },\n        {\n          \"name\": \"Geodetic longitude\",\n          \"abbreviation\": \"Lon\",\n          \"direction\": \"east\",\n          \"unit\": \"degree\"\n        }\n      ]\n    },\n    \"id\": {\n      \"authority\": \"EPSG\",\n      \"code\": 6668\n    }\n  },\n  \"conversion\": {\n    \"name\": \"unnamed\",\n    \"method\": {\n      \"name\": \"Transverse Mercator\",\n      \"id\": {\n        \"authority\": \"EPSG\",\n        \"code\": 9807\n      }\n    },\n    \"parameters\": [\n      {\n        \"name\": \"Latitude of natural origin\",\n        \"value\": 33,\n        \"unit\": \"degree\",\n        \"id\": {\n          \"authority\": \"EPSG\",\n          \"code\": 8801\n        }\n      },\n      {\n        \"name\": \"Longitude of natural origin\",\n        \"value\": 129.5,\n        \"unit\": \"degree\",\n        \"id\": {\n          \"authority\": \"EPSG\",\n          \"code\": 8802\n        }\n      },\n      {\n        \"name\": \"Scale factor at natural origin\",\n        \"value\": 0.9999,\n        \"unit\": \"unity\",\n        \"id\": {\n          \"authority\": \"EPSG\",\n          \"code\": 8805\n        }\n      },\n      {\n        \"name\": \"False easting\",\n        \"value\": 0,\n        \"unit\": \"metre\",\n        \"id\": {\n          \"authority\": \"EPSG\",\n          \"code\": 8806\n        }\n      },\n      {\n        \"name\": \"False northing\",\n        \"value\": 0,\n        \"unit\": \"metre\",\n        \"id\": {\n          \"authority\": \"EPSG\",\n          \"code\": 8807\n        }\n      }\n    ]\n  },\n  \"coordinate_system\": {\n    \"subtype\": \"Cartesian\",\n    \"axis\": [\n      {\n        \"name\": \"Northing\",\n        \"abbreviation\": \"\",\n        \"direction\": \"north\",\n        \"unit\": \"metre\"\n      },\n      {\n        \"name\": \"Easting\",\n        \"abbreviation\": \"\",\n        \"direction\": \"east\",\n        \"unit\": \"metre\"\n      }\n    ]\n  },\n  \"id\": {\n    \"authority\": \"EPSG\",\n    \"code\": 6669\n  }\n}",
      "prettycompoundwkt": "PROJCS[\"JGD2011 / Japan Plane Rectangular CS I\",\n    GEOGCS[\"JGD2011\",\n        DATUM[\"Japanese_Geodetic_Datum_2011\",\n            SPHEROID[\"GRS 1980\",6378137,298.257222101,\n                AUTHORITY[\"EPSG\",\"7019\"]],\n            AUTHORITY[\"EPSG\",\"1128\"]],\n        PRIMEM[\"Greenwich\",0,\n            AUTHORITY[\"EPSG\",\"8901\"]],\n        UNIT[\"degree\",0.0174532925199433,\n            AUTHORITY[\"EPSG\",\"9122\"]],\n        AUTHORITY[\"EPSG\",\"6668\"]],\n    PROJECTION[\"Transverse_Mercator\"],\n    PARAMETER[\"latitude_of_origin\",33],\n    PARAMETER[\"central_meridian\",129.5],\n    PARAMETER[\"scale_factor\",0.9999],\n    PARAMETER[\"false_easting\",0],\n    PARAMETER[\"false_northing\",0],\n    UNIT[\"metre\",1,\n        AUTHORITY[\"EPSG\",\"9001\"]],\n    AXIS[\"Northing\",NORTH],\n    AXIS[\"Easting\",EAST],\n    AUTHORITY[\"EPSG\",\"6669\"]]",
      "prettywkt": "PROJCS[\"JGD2011 / Japan Plane Rectangular CS I\",\n    GEOGCS[\"JGD2011\",\n        DATUM[\"Japanese_Geodetic_Datum_2011\",\n            SPHEROID[\"GRS 1980\",6378137,298.257222101,\n                AUTHORITY[\"EPSG\",\"7019\"]],\n            AUTHORITY[\"EPSG\",\"1128\"]],\n        PRIMEM[\"Greenwich\",0,\n            AUTHORITY[\"EPSG\",\"8901\"]],\n        UNIT[\"degree\",0.0174532925199433,\n            AUTHORITY[\"EPSG\",\"9122\"]],\n        AUTHORITY[\"EPSG\",\"6668\"]],\n    PROJECTION[\"Transverse_Mercator\"],\n    PARAMETER[\"latitude_of_origin\",33],\n    PARAMETER[\"central_meridian\",129.5],\n    PARAMETER[\"scale_factor\",0.9999],\n    PARAMETER[\"false_easting\",0],\n    PARAMETER[\"false_northing\",0],\n    UNIT[\"metre\",1,\n        AUTHORITY[\"EPSG\",\"9001\"]],\n    AXIS[\"Northing\",NORTH],\n    AXIS[\"Easting\",EAST],\n    AUTHORITY[\"EPSG\",\"6669\"]]",
      "proj4": "+proj=tmerc +lat_0=33 +lon_0=129.5 +k=0.9999 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs",
      "units":
      {
        "horizontal": "metre",
        "vertical": ""
      },
      "vertical": "",
      "wkt": "PROJCS[\"JGD2011 / Japan Plane Rectangular CS I\",GEOGCS[\"JGD2011\",DATUM[\"Japanese_Geodetic_Datum_2011\",SPHEROID[\"GRS 1980\",6378137,298.257222101,AUTHORITY[\"EPSG\",\"7019\"]],AUTHORITY[\"EPSG\",\"1128\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"6668\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",33],PARAMETER[\"central_meridian\",129.5],PARAMETER[\"scale_factor\",0.9999],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Northing\",NORTH],AXIS[\"Easting\",EAST],AUTHORITY[\"EPSG\",\"6669\"]]"
    },
    "system_id": "PDAL",
    "vlr_0":
    {
      "data": "AQABAAAABwAABAAAAQABAAEEAAABAAEAAgSxhycAAAABCLGHCAAnAAYIAAABAI4jAAwAAAEADRoEDAAAAQApIw==",
      "description": "GeoTiff GeoKeyDirectoryTag",
      "record_id": 34735,
      "user_id": "LASF_Projection"
    },
    "vlr_1":
    {
      "data": "SkdEMjAxMSAvIEphcGFuIFBsYW5lIFJlY3Rhbmd1bGFyIENTIEl8SkdEMjAxMXwA",
      "description": "GeoTiff GeoAsciiParamsTag",
      "record_id": 34737,
      "user_id": "LASF_Projection"
    }
  },
  "now": "2023-03-15T14:26:06+0900",
  "pdal_version": "2.5.0 (git-version: Release)",
  "reader": "readers.las"
}
```

メタデータがめちゃくちゃ増えましたね。
`\"JGD2011 / Japan Plane Rectangular CS I\"`のような記述も見つかり、ちゃんと座標系が付与されていそうです。

次に、データの先頭に格納されているポイントの位置を確認してみましょう。

```bash
% pdal info -p 0 01ke9821_translated.las
{
  "file_size": 617952997,
  "filename": "01ke9821_translated.las",
  "now": "2023-03-15T14:27:45+0900",
  "pdal_version": "2.5.0 (git-version: Release)",
  "points":
  {
    "point":
    {
      "Blue": 45232,
      "Classification": 1,
      "EdgeOfFlightLine": 0,
      "GpsTime": 0,
      "Green": 48830,
      "Intensity": 620,
      "NumberOfReturns": 0,
      "PointId": 0,
      "PointSourceId": 29,
      "Red": 48316,
      "ReturnNumber": 0,
      "ScanAngleRank": 0,
      "ScanDirectionFlag": 0,
      "UserData": 0,
      "X": 34000.12,
      "Y": -27664.3,
      "Z": 20.46
    }
  },
  "reader": "readers.las"
}
```

X座標が`34000.12m`、Y座標が`-27664.3m`であることがわかります。

国土地理院が[緯度、経度への換算](https://vldb.gsi.go.jp/sokuchi/surveycalc/surveycalc/xy2blf.html)というサービスを用意しているので、これを利用して座標を経緯度に変換し、どんな位置にあるデータなのか確認してみましょう。

すると、この座標は海の上に存在することがわかりました。
そんなはずありませんね…

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/203944/82923eb5-99a9-5111-98f5-86392e675560.png)

これは日本の領域の点群データを扱うときによくある問題なんですが…

- 平面直角座標系は、Xが縦軸で、Yが横軸になる、いわゆる左手系
- lasデータなどの点群データは一般的に1番目の要素をXと「呼ぶ」
- 一般的には、1番目の要素には数学的なX（つまり横軸）の座標が入る
- が、平面直角座標の点群データは横軸がYを示すので1番目の要素はY
    - 2番目の要素がXになる
- 横軸（平面直角座標のY）は経度を表す
    - 経度は横軸なので、数学的にはX軸

という問題があります。
3Dの軸と、地球上の座標の取り方の違いの話ですね。

このため、地球上の正しい座標に配置したいときはXとYを入れ替えてあげる必要（もしくは利用したいアプリケーション上で、1番目の項目をY座標として扱う必要）があります。

これをPDALでやっていきましょう。

まずは`xy_switch_pipeline.json`のような適当なファイル名のJSONを作成します。

```json
[
  {
    "type": "readers.las",
    "filename": "01ke9821_translated.las",
    "spatialreference": "EPSG:6669"
  },
  {
    "type": "filters.reprojection",
    "in_srs": "EPSG:6669",
    "out_srs": "EPSG:6669",
    "in_axis_ordering": "2, 1"
  },
  {
    "type": "writers.las",
    "filename": "01ke9821_swaped.las",
    "forward": "header,scale,vlr",
    "offset_x": "auto",
    "offset_y": "auto",
    "offset_z": "auto"
  }
]
```

PDALとこのJSONを利用して軸を入れ替えます。

```bash
pdal pipeline xy_switch_pipeline.json
```

変換されたデータを確認してみましょう。
XとYが入れ替わっていますね！

```bash
% pdal info -p 0 01ke9821_swaped.las
{
  "file_size": 617952997,
  "filename": "01ke9821_swaped.las",
  "now": "2023-03-15T14:40:32+0900",
  "pdal_version": "2.5.0 (git-version: Release)",
  "points":
  {
    "point":
    {
      "Blue": 45232,
      "Classification": 1,
      "EdgeOfFlightLine": 0,
      "GpsTime": 0,
      "Green": 48830,
      "Intensity": 620,
      "NumberOfReturns": 0,
      "PointId": 0,
      "PointSourceId": 29,
      "Red": 48316,
      "ReturnNumber": 0,
      "ScanAngleRank": 0,
      "ScanDirectionFlag": 0,
      "UserData": 0,
      "X": -27664.3,
      "Y": 34000.12,
      "Z": 20.46
    }
  },
  "reader": "readers.las"
}
```

## py3dtilesで点群データを3DTilesに変換する。

点群データをCesium上に載せるには[3DTiles](https://cesium.com/why-cesium/3d-tiles/)という形式に変換するのがスムーズです。

点群データを3DTilesに変換可能なOSSはあまり多くなく、[py3dtiles](https://gitlab.com/Oslandia/py3dtiles)というツールを利用するのがベターかと思われます。
Cesium自体が運営するCesium ionというサービスがありますが、一定以上の利用を超えると有料になります。ただし、手軽に変換でき、かつホスティングもされるのでとても優秀なサービスではあります。

※とはいえ、地形情報などを表示させる際にはCesium ionのトークンを利用しますので無料アカウントを作成し、トークンを発行しておくことをお勧めします。
[Cesium ion Access Tokens](https://cesium.com/learn/ion/cesium-ion-access-tokens/)

僕はM1 Macを利用していますが、M1だとインストールがうまくいかないので、[Docker](https://www.docker.com/)を利用しました。
インストールの方法は以下を参考にしてください。
https://docs.docker.com/engine/install/

最新はv5.0.0ですが、M1 Macではv3.0.0以降動かなかったので、このようなDockerfileを作成しました。

```dockerfile
FROM python:3.8.16-slim
RUN pip install py3dtiles==3.0.0
```

その後、`amd64`としてイメージをビルドします。

```bash
docker build --platform=linux/amd64 -t 3dtiles .
```

ビルド後にpy3dtilesを実行し、データを変換しましょう。
（時間がかかるので、気長に待ちましょう。）

```bash
docker run --platform=linux/amd64 -v "$(pwd)":/home 3dtiles py3dtiles convert --srs_in 6669 --srs_out 4978 --out home/01ke9821 home/01ke9821_swaped.las
```

すると`01ke9821`のようなフォルダが作成され、その中に大量の`.pnts`というファイルが生成されているかと思います。

## Cesiumで3DTilesを表示

生成されたフォルダと同じフォルダに`index.html`と`style.css`と`main.js`を作成しましょう。
その後以下のように修正してください。

- index.html
```html
<!DOCTYPE html>
<html>
  <head>
    <title>Cesium Earth</title>
    <script src="https://cesium.com/downloads/cesiumjs/releases/1.101/Build/Cesium/Cesium.js"></script>
    <link
      href="https://cesium.com/downloads/cesiumjs/releases/1.101/Build/Cesium/Widgets/widgets.css"
      rel="stylesheet"
    />
  </head>
  <body>
    <div id="cesiumContainer"></div>
    <script src="main.js"></script>
  </body>
</html>
```

- style.css

```css
body {
    padding: 0;
    margin: 0;
    overflow: hidden;
    height: 100%;
}

#cesiumContainer {
    display: flex;
    width: 100%;
    height: 100%;
}
```

- main.js
    - `YOUR_TOKEN`のところはCesium ionのトークンを記述してください。

```js
const token = 'YOUR_TOKEN`;

Cesium.Ion.defaultAccessToken = token;

const viewer = new Cesium.Viewer("cesiumContainer");

const las = viewer.scene.primitives.add(
  new Cesium.Cesium3DTileset({
    url: "./01ke9821/tileset.json",
    pointCloudShading: {
      attenuation: true,
      maximumAttenuation: 2,
    },
  })
);

viewer.zoomTo(las);
```

最後に、サーバーを立ち上げてCesiumを起動しましょう！

```bash
% python -m http.server 8000
```

`http://localhost:8000/`へブラウザでアクセスすると点群データが表示されるかと思います！！

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/203944/79dbfefa-8f16-b895-a3ba-15795952421a.png)

以上で終わりです！
点群データ、見るだけならとても簡単なので、みんな触ってみてください！！