LoginSignup
52
43

More than 1 year has passed since last update.

@nokonoko_1203(のこのこ)inMIERUNE

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

Last updated at Posted at 2023-03-15

image.png

（※オープンナガサキのトップページ）

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

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

（VIRTUAL SHIZUOKAの点群データはこちら)

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

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

image.png

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

データのダウンロード

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

オープンナガサキのページに移動しスクロールしていくと「ダウンロードページ」のリンクが出てきます！

image.png

クリックすると利用規約が出てきます。

データの座標参照系は、世界測地系(JGD2011)／平面直角座標系第1系です。

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

image.png

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

image.png

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

image.png

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

image.png

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

image.png

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

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

ダウンロードしてきた点群データをPDALというツールを利用して確認・調整していきます。
インストールの方法はこちらを参照してください。

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

% 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"
}

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

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

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

（参考: EPSGコード一覧表/日本でよく利用される空間座標系（座標参照系））

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

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

%  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\"のような記述も見つかり、ちゃんと座標系が付与されていそうです。

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

% 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であることがわかります。

国土地理院が緯度、経度への換算というサービスを用意しているので、これを利用して座標を経緯度に変換し、どんな位置にあるデータなのか確認してみましょう。

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

image.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を作成します。

[
  {
    "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を利用して軸を入れ替えます。

pdal pipeline xy_switch_pipeline.json

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

% 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という形式に変換するのがスムーズです。

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

※とはいえ、地形情報などを表示させる際にはCesium ionのトークンを利用しますので無料アカウントを作成し、トークンを発行しておくことをお勧めします。
Cesium ion Access Tokens

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

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

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

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

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

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

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
<!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
body {
    padding: 0;
    margin: 0;
    overflow: hidden;
    height: 100%;
}

#cesiumContainer {
    display: flex;
    width: 100%;
    height: 100%;
}
  • main.js
    • YOUR_TOKENのところはCesium ionのトークンを記述してください。
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を起動しましょう！

% python -m http.server 8000

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

image.png

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

52
43
2

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
52
43