みなさん！！！点群データ使ってますかーーーー！！！

というのは、毎回点群記事を書くたびに聞いている気がして、前回「LAZ形式の点群データ」について記事を出したときも同じようなことを言っていました。

前回記事
- LAZ形式の点群データがLASやCSVと比較してなぜ軽量なのか、をまとめてみる

前回は神奈川県が点群データを出してくれていたんですが、今回は本命！
東京都が点群データが公開されました！

https://www.geospatial.jp/ckan/dataset/tokyopc-23ku-2024

今まで見てきたオープン点群データの中で最高品質な気がしていて

強度や分類などの属性も埋め込まれている
座標系の設定がされている
脅威の16point/m2の高密度点群データ

となっています！
（ただし、やはりv1.2…なぜ…）

❯ pdal info --metadata ~/Downloads/pointcloud/data/09LD1729.las
{
  "file_size": 216896854,
  "filename": "~/Downloads/pointcloud/data/09LD1729.las",
  "metadata":
  {
    "comp_spatialreference": "PROJCS[\"JGD2011 / Japan Plane Rectangular CS IX\",GEOGCS[\"JGD2011\",DATUM[\"Japanese_Geodetic_Datum_2011\",SPHEROID[\"GRS 1980\",6378137,298.257222101004,AUTHORITY[\"EPSG\",\"7019\"]],AUTHORITY[\"EPSG\",\"1128\"]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"6668\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",36],PARAMETER[\"central_meridian\",139.833333333333],PARAMETER[\"scale_factor\",0.9999],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH],AUTHORITY[\"EPSG\",\"6677\"]]",
    "compressed": false,
    "copc": false,
    "count": 6379285,
    "creation_doy": 303,
    "creation_year": 2023,
    "dataformat_id": 3,
    "dataoffset": 1164,
    "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      GTRasterTypeGeoKey (Short,1): RasterPixelIsArea\n      GTModelTypeGeoKey (Short,1): ModelTypeProjected\n      ProjectedCSTypeGeoKey (Short,1): Code-6677 (JGD2011 / Japan Plane Rectangular CS IX)\n      ProjLinearUnitsGeoKey (Short,1): Linear_Meter\n      GeogAngularUnitsGeoKey (Short,1): Angular_Degree\n      GeogAngularUnitSizeGeoKey (Double,1): 0.0174532925199433\n      GeographicTypeGeoKey (Short,1): Code-6668 (JGD2011)\n      GeogGeodeticDatumGeoKey (Short,1): Code-1128 (Japanese Geodetic Datum 2011)\n      GeogEllipsoidGeoKey (Short,1): Ellipse_GRS_1980\n      GeogPrimeMeridianGeoKey (Short,1): PM_Greenwich\n      GeogSemiMajorAxisGeoKey (Double,1): 6378137          \n      GeogInvFlatteningGeoKey (Double,1): 298.257222101    \n      GeogPrimeMeridianLongGeoKey (Double,1): 0                \n      GTCitationGeoKey (Ascii,34): \"PCS Name = JGD_2011_Japan_Zone_9\"\n      GeogCitationGeoKey (Ascii,85): \"GCS Name = GCS_JGD_2011|Datum = D_JGD_2011|Ellipsoid = GRS_1980|Primem = Greenwich|\"\n      PCSCitationGeoKey (Ascii,428): \"ESRI PE String = PROJCS[\"JGD_2011_Japan_Zone_9\",GEOGCS[\"GCS_JGD_2011\",DATUM[\"D_JGD_2011\",SPHEROID[\"GRS_1980\",6378137.0,298.257222101]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"False_Easting\",0.0],PARAMETER[\"False_Northing\",0.0],PARAMETER[\"Central_Meridian\",139.8333333333333],PARAMETER[\"Scale_Factor\",0.9999],PARAMETER[\"Latitude_Of_Origin\",36.0],UNIT[\"Meter\",1.0]]\"\n      End_Of_Keys.\n   End_Of_Geotiff.\n",
    "header_size": 227,
    "major_version": 1,
    "maxx": -8000.001,
    "maxy": -33600.001,
    "maxz": 151.226,
    "minor_version": 2,
    "minx": -8400,
    "miny": -33900,
    "minz": 3.059,
    "offset_x": 0,
    "offset_y": 0,
    "offset_z": 0,
    "point_length": 34,
    "project_id": "00000000-0000-0000-0000-000000000000",
    "scale_x": 0.001,
    "scale_y": 0.001,
    "scale_z": 0.001,
    "software_id": "ArcGIS",
    "spatialreference": "PROJCS[\"JGD2011 / Japan Plane Rectangular CS IX\",GEOGCS[\"JGD2011\",DATUM[\"Japanese_Geodetic_Datum_2011\",SPHEROID[\"GRS 1980\",6378137,298.257222101004,AUTHORITY[\"EPSG\",\"7019\"]],AUTHORITY[\"EPSG\",\"1128\"]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"6668\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",36],PARAMETER[\"central_meridian\",139.833333333333],PARAMETER[\"scale_factor\",0.9999],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH],AUTHORITY[\"EPSG\",\"6677\"]]",
    "srs":
    {
      "compoundwkt": "PROJCS[\"JGD2011 / Japan Plane Rectangular CS IX\",GEOGCS[\"JGD2011\",DATUM[\"Japanese_Geodetic_Datum_2011\",SPHEROID[\"GRS 1980\",6378137,298.257222101004,AUTHORITY[\"EPSG\",\"7019\"]],AUTHORITY[\"EPSG\",\"1128\"]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"6668\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",36],PARAMETER[\"central_meridian\",139.833333333333],PARAMETER[\"scale_factor\",0.9999],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH],AUTHORITY[\"EPSG\",\"6677\"]]",
      "horizontal": "PROJCS[\"JGD2011 / Japan Plane Rectangular CS IX\",GEOGCS[\"JGD2011\",DATUM[\"Japanese_Geodetic_Datum_2011\",SPHEROID[\"GRS 1980\",6378137,298.257222101004,AUTHORITY[\"EPSG\",\"7019\"]],AUTHORITY[\"EPSG\",\"1128\"]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"6668\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",36],PARAMETER[\"central_meridian\",139.833333333333],PARAMETER[\"scale_factor\",0.9999],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH],AUTHORITY[\"EPSG\",\"6677\"]]",
      "isgeocentric": false,
      "isgeographic": false,
      "json": {
  "type": "ProjectedCRS",
  "name": "JGD2011 / Japan Plane Rectangular CS IX",
  "base_crs": {
    "name": "JGD2011",
    "datum": {
      "type": "GeodeticReferenceFrame",
      "name": "Japanese Geodetic Datum 2011",
      "ellipsoid": {
        "name": "GRS 1980",
        "semi_major_axis": 6378137,
        "inverse_flattening": 298.257222101004
      }
    },
    "coordinate_system": {
      "subtype": "ellipsoidal",
      "axis": [
        {
          "name": "Geodetic latitude",
          "abbreviation": "Lat",
          "direction": "north",
          "unit": "degree"
        },
        {
          "name": "Geodetic longitude",
          "abbreviation": "Lon",
          "direction": "east",
          "unit": "degree"
        }
      ]
    },
    "id": {
      "authority": "EPSG",
      "code": 6668
    }
  },
  "conversion": {
    "name": "unnamed",
    "method": {
      "name": "Transverse Mercator",
      "id": {
        "authority": "EPSG",
        "code": 9807
      }
    },
    "parameters": [
      {
        "name": "Latitude of natural origin",
        "value": 36,
        "unit": "degree",
        "id": {
          "authority": "EPSG",
          "code": 8801
        }
      },
      {
        "name": "Longitude of natural origin",
        "value": 139.833333333333,
        "unit": "degree",
        "id": {
          "authority": "EPSG",
          "code": 8802
        }
      },
      {
        "name": "Scale factor at natural origin",
        "value": 0.9999,
        "unit": "unity",
        "id": {
          "authority": "EPSG",
          "code": 8805
        }
      },
      {
        "name": "False easting",
        "value": 0,
        "unit": "metre",
        "id": {
          "authority": "EPSG",
          "code": 8806
        }
      },
      {
        "name": "False northing",
        "value": 0,
        "unit": "metre",
        "id": {
          "authority": "EPSG",
          "code": 8807
        }
      }
    ]
  },
  "coordinate_system": {
    "subtype": "Cartesian",
    "axis": [
      {
        "name": "Easting",
        "abbreviation": "",
        "direction": "east",
        "unit": "metre"
      },
      {
        "name": "Northing",
        "abbreviation": "",
        "direction": "north",
        "unit": "metre"
      }
    ]
  },
  "id": {
    "authority": "EPSG",
    "code": 6677
  }
},
      "prettycompoundwkt": "PROJCS[\"JGD2011 / Japan Plane Rectangular CS IX\",\n    GEOGCS[\"JGD2011\",\n        DATUM[\"Japanese_Geodetic_Datum_2011\",\n            SPHEROID[\"GRS 1980\",6378137,298.257222101004,\n                AUTHORITY[\"EPSG\",\"7019\"]],\n            AUTHORITY[\"EPSG\",\"1128\"]],\n        PRIMEM[\"Greenwich\",0],\n        UNIT[\"degree\",0.0174532925199433,\n            AUTHORITY[\"EPSG\",\"9122\"]],\n        AUTHORITY[\"EPSG\",\"6668\"]],\n    PROJECTION[\"Transverse_Mercator\"],\n    PARAMETER[\"latitude_of_origin\",36],\n    PARAMETER[\"central_meridian\",139.833333333333],\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[\"Easting\",EAST],\n    AXIS[\"Northing\",NORTH],\n    AUTHORITY[\"EPSG\",\"6677\"]]",
      "prettywkt": "PROJCS[\"JGD2011 / Japan Plane Rectangular CS IX\",\n    GEOGCS[\"JGD2011\",\n        DATUM[\"Japanese_Geodetic_Datum_2011\",\n            SPHEROID[\"GRS 1980\",6378137,298.257222101004,\n                AUTHORITY[\"EPSG\",\"7019\"]],\n            AUTHORITY[\"EPSG\",\"1128\"]],\n        PRIMEM[\"Greenwich\",0],\n        UNIT[\"degree\",0.0174532925199433,\n            AUTHORITY[\"EPSG\",\"9122\"]],\n        AUTHORITY[\"EPSG\",\"6668\"]],\n    PROJECTION[\"Transverse_Mercator\"],\n    PARAMETER[\"latitude_of_origin\",36],\n    PARAMETER[\"central_meridian\",139.833333333333],\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[\"Easting\",EAST],\n    AXIS[\"Northing\",NORTH],\n    AUTHORITY[\"EPSG\",\"6677\"]]",
      "proj4": "+proj=tmerc +lat_0=36 +lon_0=139.833333333333 +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 IX\",GEOGCS[\"JGD2011\",DATUM[\"Japanese_Geodetic_Datum_2011\",SPHEROID[\"GRS 1980\",6378137,298.257222101004,AUTHORITY[\"EPSG\",\"7019\"]],AUTHORITY[\"EPSG\",\"1128\"]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"6668\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",36],PARAMETER[\"central_meridian\",139.833333333333],PARAMETER[\"scale_factor\",0.9999],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH],AUTHORITY[\"EPSG\",\"6677\"]]"
    },
    "system_id": "CONVERSION",
    "vlr_0":
    {
      "data": "HvM+ZQAAAAA=",
      "description": "Local Time: 2023 10 30 09 04 46",
      "record_id": 8,
      "user_id": "ESRI_Export"
    },
    "vlr_1":
    {
      "data": "AQABAAAAEAABBAAAAQABAAAEAAABAAEAAAwAAAEAFRoEDAAAAQApIwYIAAABAI4jBwiwhwEAAAAACAAAAQAMGgIIAAABAGgECAgAAAEAaxsDCAAAAQDFIgkIsIcBAAEACwiwhwEAAgANCLCHAQADAAIEsYchAAAAAQixh1QAIQABDLGHrAF1AA==",
      "description": "",
      "record_id": 34735,
      "user_id": "LASF_Projection"
    },
    "vlr_2":
    {
      "data": "UENTIE5hbWUgPSBKR0RfMjAxMV9KYXBhbl9ab25lXzkAR0NTIE5hbWUgPSBHQ1NfSkdEXzIwMTF8RGF0dW0gPSBEX0pHRF8yMDExfEVsbGlwc29pZCA9IEdSU18xOTgwfFByaW1lbSA9IEdyZWVud2ljaHwARVNSSSBQRSBTdHJpbmcgPSBQUk9KQ1NbIkpHRF8yMDExX0phcGFuX1pvbmVfOSIsR0VPR0NTWyJHQ1NfSkdEXzIwMTEiLERBVFVNWyJEX0pHRF8yMDExIixTUEhFUk9JRFsiR1JTXzE5ODAiLDYzNzgxMzcuMCwyOTguMjU3MjIyMTAxXV0sUFJJTUVNWyJHcmVlbndpY2giLDAuMF0sVU5JVFsiRGVncmVlIiwwLjAxNzQ1MzI5MjUxOTk0MzNdXSxQUk9KRUNUSU9OWyJUcmFuc3ZlcnNlX01lcmNhdG9yIl0sUEFSQU1FVEVSWyJGYWxzZV9FYXN0aW5nIiwwLjBdLFBBUkFNRVRFUlsiRmFsc2VfTm9ydGhpbmciLDAuMF0sUEFSQU1FVEVSWyJDZW50cmFsX01lcmlkaWFuIiwxMzkuODMzMzMzMzMzMzMzM10sUEFSQU1FVEVSWyJTY2FsZV9GYWN0b3IiLDAuOTk5OV0sUEFSQU1FVEVSWyJMYXRpdHVkZV9PZl9PcmlnaW4iLDM2LjBdLFVOSVRbIk1ldGVyIiwxLjBdXQA=",
      "description": "",
      "record_id": 34737,
      "user_id": "LASF_Projection"
    },
    "vlr_3":
    {
      "data": "Op1SokbfkT8AAABAplRYQaj565QdpHJAAAAAAAAAAAA=",
      "description": "",
      "record_id": 34736,
      "user_id": "LASF_Projection"
    }
  },
  "now": "2024-12-03T22:28:09+0900",
  "pdal_version": "2.8.0 (git-version: Release)",
  "reader": "readers.las"
}

整備範囲もどどーんと23区全域！！！

とんでもないですねー。~~いくらかかっとるんやろ…~~

ということでこんな点群データ触らない訳にはいかないので、Cesiumで表示できる3D Tilesに変換するためのツールを作って、実際に描画してみました！

ツールの紹介

そして、完成したものがこちらです。

以下のように簡単にインストールすることができます。

インストール

curl -sSf https://raw.githubusercontent.com/MIERUNE/point-tiler/main/install.sh | bash

LAS/LAZの変換方法

実行（LAS to 3D Tiles）

point_tiler --input app/examples/data/*.las \
    --output app/examples/data/output \
    --epsg 6677 \
    --min 15 \
    --max 18

2024/12時点ではmacOSのみに対応しており、windowsでは動かすことができません。
（Linuxでは動くかもしれません）

ただし、マシンにRustがインストールされている場合は、リポジトリのソースコードをダウンロードしてきて、フォルダ内で以下のように動かすことが可能です。

cargo run --release --package app \
    --input app/examples/data/sample.las \
    --output app/examples/data/output \
    --epsg 6677 \
    --min 15 \
    --max 18

（昔こんな記事を書きましたが、VSCodeのターミナルで実行するとめっちゃ遅くなる可能性があるので要注意です！）

設定値は以下のようになっています。

input: .las/.laz/.csv/.txtファイルを指定します。複数のファイルはスペースで区切って入力できます。またglobのように*を利用してパターンマッチすることもできます
output: 出力フォルダを指定します。tileset.jsonとglbを出力します
epsg: ファイルのepsgコードを入力します。すべての点群は同じ座標系として認識されます。（現在は、日本の平面直角座標系のみサポートしています
min: 出力したい最小ズームレベルを指定します
max: 出力したい最大ズームレベルを指定します

min/maxはズームレベルを指定します。
大体、中央区全域が画面に表示される程度の領域がズームレベル16とか17程度になります。

出力される3D Tilesは点群が間引かれますが、最大ズームレベルを上げていくとどんどん間引かれなくなっていき、22レベル程度でほぼ間引きがなくなります。

指定するファイルは、1度の実行でファイル形式が統一されていればLASファイルでもLAZファイルでも同様に実行することができます。

追記

記事の最後の方に記載していますが、利用メモリ指定オプションをつけたので以下のようなコマンドになります。
4096と指定すると4GBのメモリを利用します。

max-memory-mb: 変換に使用可能なメモリ（MB）を指定します

point_tiler --input app/examples/data/*.las \
    --output app/examples/data/output \
    --epsg 6677 \
    --min 15 \
    --max 18 \
    --max-memory-mb 4096

cargo run --release --package app \
    --input app/examples/data/sample.las \
    --output app/examples/data/output \
    --epsg 6677 \
    --min 15 \
    --max 18 \
    --max-memory-mb 4096

CSV/TXTの利用

CSVファイルはスキーマレスのため、カラム名を固定してあげなければXYZ・RGBの列を判定することはできません。

列名は「x」、「y」、「z」、および「r」、「g」、「b」、または「red」、「green」、「blue」でなければなりません。
（文字の大文字・小文字は区別されません。）
また、XYZRGB以外の列は無視されます。

例えば、以下のデータは有効です。

“X“,”Y“,”Z“,”Intensity“,”ReturnNumber“,”NumberOfReturns“,”ScanDirectionFlag“,”EdgeOfFlightLine“,”Classification“,”Synthetic“,”KeyPoint“,”Withheld“,”Overlap“,”ScanAngleRank“,”UserData“,”PointSourceId“,”GpsTime“,”Red“,”Green“,”Blue”
-5599.971,-35106.097,3.602,680.000,1.000,1.000,1.000,0.000,1.000,0.000,0.00 0,0.000,0.000,-4.000,0.000,95.000,188552.207,19532.000,20046.000,20560.000
-5599.976,-35111.458,3.440,715.000,1.000,1.000,1.000,0.000,1.000,0.000,0.00 0,0.000,0.000,-4.000,0.000,95.000,188552.287,20560.000,21074.000,21588.000
-5599.978,-35115.695,3.543,311.000,1.000,1.000,1.000,0.000,1.000,0.000,0.00 0,0.000,0.000,-4.000,0.000,95.000,188552.347,19018.000,19275.000,20303.000
-5599.992,-35119.757,3.603,722.000,1.000,1.000,1.000,0.000,1.000,0.000,0.00 0,0.000,0.000,-4.000,0.000,95.000,188552.407,18504.000,20046.000,20817.000
-5599.992,-35129.327,3.431,505.000,1.000,1.000,1.000,0.000,1.000,0.000,0.00 0,0.000,0.000,-4.000,0.000,95.000,188552.547,18504.000,19789.000,21074.000

変換とCesiumでの描画

実際にデータを変換して、描画してみましょう！
今回は東京駅周辺の42ファイル、計8.8GBのデータをダウンロードしました。

https://gic-tokyo.s3.ap-northeast-1.amazonaws.com/2024/dig/lp/09LD1854.zip
https://gic-tokyo.s3.ap-northeast-1.amazonaws.com/2024/dig/lp/09LD1855.zip
https://gic-tokyo.s3.ap-northeast-1.amazonaws.com/2024/dig/lp/09LD1857.zip
https://gic-tokyo.s3.ap-northeast-1.amazonaws.com/2024/dig/lp/09LD1858.zip
https://gic-tokyo.s3.ap-northeast-1.amazonaws.com/2024/dig/lp/09LD1856.zip
https://gic-tokyo.s3.ap-northeast-1.amazonaws.com/2024/dig/lp/09LD1859.zip
https://gic-tokyo.s3.ap-northeast-1.amazonaws.com/2024/dig/lp/09LD1869.zip
https://gic-tokyo.s3.ap-northeast-1.amazonaws.com/2024/dig/lp/09LD1867.zip
https://gic-tokyo.s3.ap-northeast-1.amazonaws.com/2024/dig/lp/09LD1865.zip
https://gic-tokyo.s3.ap-northeast-1.amazonaws.com/2024/dig/lp/09LD1864.zip
https://gic-tokyo.s3.ap-northeast-1.amazonaws.com/2024/dig/lp/09LD1866.zip
https://gic-tokyo.s3.ap-northeast-1.amazonaws.com/2024/dig/lp/09LD1868.zip
https://gic-tokyo.s3.ap-northeast-1.amazonaws.com/2024/dig/lp/09LD1874.zip
https://gic-tokyo.s3.ap-northeast-1.amazonaws.com/2024/dig/lp/09LD1884.zip
https://gic-tokyo.s3.ap-northeast-1.amazonaws.com/2024/dig/lp/09LD1894.zip
https://gic-tokyo.s3.ap-northeast-1.amazonaws.com/2024/dig/lp/09LD2804.zip
https://gic-tokyo.s3.ap-northeast-1.amazonaws.com/2024/dig/lp/09LD2814.zip
https://gic-tokyo.s3.ap-northeast-1.amazonaws.com/2024/dig/lp/09LD2815.zip
https://gic-tokyo.s3.ap-northeast-1.amazonaws.com/2024/dig/lp/09LD2816.zip
https://gic-tokyo.s3.ap-northeast-1.amazonaws.com/2024/dig/lp/09LD2817.zip
https://gic-tokyo.s3.ap-northeast-1.amazonaws.com/2024/dig/lp/09LD2818.zip
https://gic-tokyo.s3.ap-northeast-1.amazonaws.com/2024/dig/lp/09LD2819.zip
https://gic-tokyo.s3.ap-northeast-1.amazonaws.com/2024/dig/lp/09LD2809.zip
https://gic-tokyo.s3.ap-northeast-1.amazonaws.com/2024/dig/lp/09LD1899.zip
https://gic-tokyo.s3.ap-northeast-1.amazonaws.com/2024/dig/lp/09LD1889.zip
https://gic-tokyo.s3.ap-northeast-1.amazonaws.com/2024/dig/lp/09LD1879.zip
https://gic-tokyo.s3.ap-northeast-1.amazonaws.com/2024/dig/lp/09LD1885.zip
https://gic-tokyo.s3.ap-northeast-1.amazonaws.com/2024/dig/lp/09LD1895.zip
https://gic-tokyo.s3.ap-northeast-1.amazonaws.com/2024/dig/lp/09LD2805.zip
https://gic-tokyo.s3.ap-northeast-1.amazonaws.com/2024/dig/lp/09LD1875.zip
https://gic-tokyo.s3.ap-northeast-1.amazonaws.com/2024/dig/lp/09LD1876.zip
https://gic-tokyo.s3.ap-northeast-1.amazonaws.com/2024/dig/lp/09LD1877.zip
https://gic-tokyo.s3.ap-northeast-1.amazonaws.com/2024/dig/lp/09LD1878.zip
https://gic-tokyo.s3.ap-northeast-1.amazonaws.com/2024/dig/lp/09LD1888.zip
https://gic-tokyo.s3.ap-northeast-1.amazonaws.com/2024/dig/lp/09LD1887.zip
https://gic-tokyo.s3.ap-northeast-1.amazonaws.com/2024/dig/lp/09LD1886.zip
https://gic-tokyo.s3.ap-northeast-1.amazonaws.com/2024/dig/lp/09LD1896.zip
https://gic-tokyo.s3.ap-northeast-1.amazonaws.com/2024/dig/lp/09LD1897.zip
https://gic-tokyo.s3.ap-northeast-1.amazonaws.com/2024/dig/lp/09LD1898.zip
https://gic-tokyo.s3.ap-northeast-1.amazonaws.com/2024/dig/lp/09LD2808.zip
https://gic-tokyo.s3.ap-northeast-1.amazonaws.com/2024/dig/lp/09LD2807.zip
https://gic-tokyo.s3.ap-northeast-1.amazonaws.com/2024/dig/lp/09LD2806.zip

すべてG空間情報センターからダウンロードし、解凍し、以下のように変換します。

point_tiler --input pointcloud/data/*.las \
    --output output/3dtiles_tokyo_pointcloud \
    --epsg 6677 \
    --min 15 \
    --max 18

2024-12-11 13:51:38 [INFO] - start parsing...
2024-12-11 13:54:12 [INFO] - finish parsing in 153.707916083s
2024-12-11 13:54:12 [INFO] - start transforming...
2024-12-11 13:59:17 [INFO] - Finish transforming in 305.207464542s
2024-12-11 13:59:17 [INFO] - start tiling at max_zoom (18)...
2024-12-11 14:01:04 [INFO] - Finish tiling at max_zoom in 107.377314583s
2024-12-11 14:01:04 [INFO] - start zoom aggregation...
2024-12-11 14:01:04 [INFO] - aggregating zoom level: 17
2024-12-11 14:01:43 [INFO] - aggregating zoom level: 16
2024-12-11 14:02:24 [INFO] - aggregating zoom level: 15
2024-12-11 14:03:02 [INFO] - Finish zoom aggregation in 117.477323458s
2024-12-11 14:03:02 [INFO] - start exporting tiles (GLB)...
2024-12-11 14:05:47 [INFO] - Finish exporting tiles in 165.417663125s
2024-12-11 14:05:47 [INFO] - write tileset.json: "output/3dtiles_tokyo_pointcloud/tileset.json"
2024-12-11 14:05:47 [INFO] - Elapsed: 849.190854791s
2024-12-11 14:05:47 [INFO] - Finish processing

変換にかかった時間は14分程度でした。
結構早いんではないでしょうか！？
変換後のファイルは423ファイルになり、合計570MBでした。
（上の方で説明した通り、ズームレベル18だと点群の間引きが発生するため、総ファイルサイズが小さくなっています。）

描画するとこんな感じになりました！！！

追記

ちなみに、Xで指摘していただいたんですが、デフォルトのファイル読み込み数が制限されている可能性があるので、大量の点群データを処理しようとするときには以下のようなコマンドでファイルディスクリプタ数を拡張してあげる必要があるようでした。

ulimit -n 65535

実装方法の紹介

3D Tilesへのデータ変換というと、先日Rust.tokyo 2024というイベントでも発表してきたPLATEAU GIS Converterなどがあります。
その時の発表資料も置いておきます。

当然、これらの実装の知識を多分に活用しているので、3D Tilesに関してはそちらの説明をいつかできればと思います。
点群データのパースに関しては地味な話が多いので、多くは語らないのですが、色々考えて実装した点があります。

「点群」を表す構造体

「点群データって、なんだろう」と考えた時に以下のような特徴があるのかなと思いました。

XYZ座標を持っている
RGB情報を持っていたりそうでなかったり
強度や分類などの属性情報を持っていたりそうでなかったり
その他、点群ファイル自体に座標系などのメタ情報を持っていたり

さらに、「属性情報」は実質LAS形式がスタンダードになっていること、そのLAS自体にメタ情報を埋め込めること、を考えると以下のような構造になるんだろうなーと思いました。
（細かい話をすれば、RGBは8bitだろうかそれとも16bitだろうか〜とか、座標値は整数にしておこうかそれとも浮動小数点数の生データを保持しようか〜とか設計判断が必要な箇所はいくらでもありそうです。）

#[derive(Debug, Clone)]
pub struct PointAttributes {
    pub intensity: Option<u16>,
    pub return_number: Option<u8>,
    pub classification: Option<String>,
    pub scanner_channel: Option<u8>,
    pub scan_angle: Option<f32>,
    pub user_data: Option<u8>,
    pub point_source_id: Option<u16>,
    pub gps_time: Option<f64>,
}

#[derive(Debug, Clone, Default)]
pub struct Color {
    pub r: u16,
    pub g: u16,
    pub b: u16,
}

// LAS data coordinates are expressed in u32 format
// The actual coordinates are calculated based on a combination of scale and offset, as follows
// x = (x * scale[0]) + offset[0]
#[derive(Debug, Clone)]
pub struct Point {
    pub x: f64,
    pub y: f64,
    pub z: f64,
    pub color: Color,
    pub attributes: PointAttributes,
}

#[derive(Debug, Clone)]
pub struct PointCloud {
    pub points: Vec<Point>,
    pub metadata: Metadata,
}

// This represents the maximum and minimum values of the original coordinate values obtained by combining the scale and offset.
#[derive(Debug, Clone, Default)]
pub struct BoundingVolume {
    pub min: [f64; 3],
    pub max: [f64; 3],
}

#[derive(Debug, Clone, Default)]
pub struct Metadata {
    pub point_count: usize,
    pub bounding_volume: BoundingVolume,
    pub epsg: EpsgCode,
    pub scale: [f64; 3],
    pub offset: [f64; 3],
    pub other: HashMap<String, String>,
}

今回は、XYZ/RGB/その他属性情報を持つPointと、Pointの可変長配列とメタ情報を持つPointCloud構造体を作成することにしました。

点群データの間引き

3D Tilesは「タイル形式」、つまり元は1つだった点群データを複数に分割します。
本ツールでは、MVTやラスタータイルと同じような領域でタイル分割を行ったため、例えば、東京都のLASデータ1つ変換することを考えると、z15では分割なし・z16では4分割・z17では16分割、のように分割されます。
（つまり、4分木です。）

カメラが遠い（ズームレベルが低い）場合にはかなり広い範囲が表示されますが、その時には必ず点群間引きを行う必要があり、さもなければ表示される点の数が多すぎてブラウザがクラッシュする可能性が上がります。

このため、今回は3D TilesのgeometricErrorに合わせて適切に間引きを行うようにしました。

具体的には点群データのBounding Volumeをボクセル化し、1つのボクセルに1点のみ格納されるようにする、また、geometricErrorに応じてボクセルサイズを可変にするような実装を行いました。

以下のようなイメージです。
（ただし、今回は平均位置に新たな点を配置するのではなく、実際に存在する点を残しました。）

（Open3Dで3D点群のダウンサンプリングを行うプログラム（Open3DとPythonによる実装）から引用）

点群の量子化

生の点群の座標は32bitの浮動小数点数で表現される（精度の問題であり、経緯度を格納する場合では64bitで保持した方が良い）ことが多いため、それだけで4バイト消費してしまいます。

このため、すべての点を原点(0,0,0)付近にオフセットしてあげることで、原点を基準とした相対座標で表現でき、より少ない桁数で精度を確保できるようになります。

尚且つ、小数点以下をすべて「スケール」という値に書き出してあげることで、それを保持しておけば座標値は整数で格納することも可能になります。

このような手法がLAS形式の内部で行われている（詳しくはこちら）のですが、さらにファイルサイズを削減できます。

オフセットされた値を「0~1」の値に正規化し、さらに指定ビット数の値の範囲（例えば16bit整数なら0~65535）に分散させてあげることで、32bit浮動小数点数だったものが16bit整数で表現できるようになり、単純にファイルサイズが1/2になります。
これを量子化といいますが、ツールではこれを実装しています。

pub fn quantize_unsigned_norm(value: f32, bits: i32) -> i32 {
    let max_value = (1i32 << bits) - 1i32;
    let scale = max_value as f32;

    let clamped = value.clamp(0.0, 1.0);

    (clamped * scale + 0.5) as i32
}

3D Tiles（というかglTF）ではKHR_mesh_quantizationという拡張機能が存在し、この拡張機能が付与されたglTFは、量子化していても自動的にデコードして描画してくれます。

今後の展望

正直なところ、完成度は5割程度といったところで、まだまだやることは残っています。

Windowsへの対応
octreeによるタイル化
ストリーミング処理
gzip圧縮
EXT_mesh_features を使用した属性の割り当て
meshoptを使用した圧縮
海外のCRSに対応

やることは多いですが、ストリーミング処理は必須かなと思っています。

ここには記載していませんが、点群データを読み込む際にはすべてのデータを一度メモリに載せています。

タイル化の際には外部ソートしているため、メモリに乗り切るファイルサイズであれば変換可能ですが、東京23区全域などを対象にするとファイルサイズが500GB程度になるはずなので、1度に大規模の変換を行うのは難しいです。
なので、この辺りの対応を優先的にやる必要がありそうです。

ただし、3D Tiles v1.1への変換にしっかりと対応したツールは少ないことや日本で作成されたため、日本のサポートを手厚くすることができる点はメリットかなと思います。

今後は、GeoTIFF（DEM）やLandXMLなどの地形データへの対応も検討しているため、そういった点は良さそうです。

ということで、ひとまず点群データを変換して描画するところまでやっていきました！
皆さんもぜひ利用してみて、なんならPRいただけると大変嬉しいです！！！

20241216_追記

メモリに乗り切るファイルサイズであれば変換可能ですが、東京23区全域などを対象にするとファイルサイズが500GB程度になるはずなので、1度に大規模の変換を行うのは難しいです。

と、書いたんですが、それは不便だろう…と思ったので、対応しました！
ストリーミングで処理するようにしたのと、利用メモリ制限（MBで指定）のオプションをつけています！

182ファイル、33GBのLASデータを試しに変換してみたところ、無事に変換できました！

処理時間は53分ほどでした！
悪くないんじゃないでしょうか！

point_tiler --input pointcloud/data/*.las \
    --output output/3dtiles_tokyo_pointcloud \
    --epsg 6677 \
    --min 15 \
    --max 18 \
    --max-memory-mb 8192
    
2024-12-16 00:17:30 [INFO] - start parse and transform and tiling at max_zoom (18)...
2024-12-16 00:17:30 [INFO] - max_memory_mb_bytes: 8589934592
2024-12-16 00:17:30 [INFO] - one_chunk_mem: 960000000
2024-12-16 00:17:30 [INFO] - channel_capacity: 8
2024-12-16 00:25:10 [INFO] - Finish transforming and tiling in 459.994233209s
2024-12-16 00:25:10 [INFO] - start sorting...
2024-12-16 00:44:28 [INFO] - Finish sorting in 1157.284924292s
2024-12-16 00:44:28 [INFO] - start zoom aggregation...
2024-12-16 00:44:28 [INFO] - aggregating zoom level: 17
2024-12-16 00:49:03 [INFO] - aggregating zoom level: 16
2024-12-16 00:53:18 [INFO] - aggregating zoom level: 15
2024-12-16 00:57:46 [INFO] - Finish zoom aggregation in 797.973805166s
2024-12-16 00:57:46 [INFO] - start exporting tiles (GLB)...
2024-12-16 01:10:53 [INFO] - Finish exporting tiles in 787.9703145s
2024-12-16 01:10:53 [INFO] - write tileset.json: "/output/3dtiles_tokyo_pointcloud/tileset.json"
2024-12-16 01:10:53 [INFO] - Elapsed: 3203.23167475s
2024-12-16 01:10:53 [INFO] - Finish processing

これだけ広いと圧巻ですね！！！