LoginSignup
25
23

More than 1 year has passed since last update.

@nokonoko_1203(のこのこ)inMIERUNE

iPhoneでも使えるようになったLiDARの標準ファイル形式「LAS」ってどんなデータなの？を、PDAL/Laspyを使って調べてみる。

Last updated at Posted at 2022-06-15

lasデータってなんや

image.png
※掛川城の点群データをopen3Dで可視化したやつ

以前こんな記事を書きました。
ざっくり要約すると「lasデータをOpen3Dで可視化し、plyで出力する記事」なのですが、この記事で元データとなった「LAS」という形式についてちょろっとだけ触れています。
Leafmap/Open3Dを使って掛川城の大規模点群データ（5GB）をPythonで可視化してみよう！

ただ、「点群データだよ」と「LiDARの標準形式だよ」ということ以外にはほとんど触れなかったので今回はもう少し踏み込んだ上で、Pythonから色々触ってみて、さらに点群データの超軽量化までやってしまいましょう！

どんなデータなん？

仕様はこちらになりますが、概要としてはLiDAR（やその他の）点群データを保持することを目的としているオープンなバイナリのファイル形式になります。
https://www.ogc.org/standards/LAS

LASファイルは主に以下の3つセクションで構成されています。

  • Public header block
    • ポイント総数や点群の範囲などを格納
  • Variable length records(VLR)
    • 空間参照システム（CRS）など、任意のオプションが格納
  • Point data records
    • ポイントやRGBなど実際のスキャンデータを格納

そのほかExtended variable length records(EVLR)というVLRの拡張セクションも存在します。

要はXYZやRGBのようなただの数値データにヘッダーとVLRに付与されている様々なメタ情報を組み合わせて利用される形式だということですね！

とはいえ、これだけ言われてもよくわからないと思うので、laspyを使って実際のlasデータの中身を見ていきましょう！

laspyのインストール

何はともあれ、jupyter/Pythonを利用してLASに触ってみましょう！

まずはlaspyというライブラリを利用していきます。
laspy: https://github.com/laspy/laspy

以前書いたこちらの記事に沿って環境を構築してくれた方はすでにlaspyが入っていると思います。

そうではない方は以下のコマンドでlaspyだけcondaでインストールすることができます！

conda install -c conda-forge laspy

また、PDALもインストールしていきましょう！
（こちらは上の記事で紹介していないので、皆さんダウンロードお願いします。）
割と時間がかかるので、気長にお待ちください！
pdal: https://github.com/PDAL/PDAL

conda create -c conda-forge pdal python-pdal gdal

データのダウンロード

適当なlasをお持ちの方はそちらを利用しても良いのですが、Leafmap/Open3Dを使って掛川城の大規模点群データ（5GB）をPythonで可視化してみよう！ の記事から来ていただいた方はそのまま「掛川城」の点群データをお使いください。

「掛川城の点群データがあるの！？」ってなった方は上の記事内にリンクが貼ってあるのでそちらを参照してください！

laspyでlasのメタ情報を読みとる

インストールが終わったら早速lasデータを触っていきましょう！

まずはインポートして…

import laspy

readしてみましょう！
5GBもあるんですが、3秒ほどでサクッと読み込めます！

# https://www.geospatial.jp/ckan/dataset/kakegawacastle/resource/61d02b61-3a44-4a5b-b263-814c6aa23551
las = laspy.read("./20190308掛川城.las")
# CPU times: user 1.04 s, sys: 1.69 s, total: 2.73 s
# Wall time: 2.72 s

早速ヘッダーの情報を覗いてみましょう。

header = las.header
header
# <LasHeader(1.2, <PointFormat(2, 0 bytes of extra dims)>)>

LasHeader(1.2, <PointFormat(2, 0 bytes of extra dims)とありますね。
これはlasのバージョンが1.2、ポイントフォーマットが2であることを示しています。

さっぱり意味がわかりませんね。

よくわからないので仕様を確認してみましょう！
https://pythonhosted.org/laspy/tut_background.html

これ見ているとLASにはバージョン1.0〜1.4まで存在し、ポイントフォーマットは0~10まで存在することがわかります。
Field Nameに()がついているものは必須ではない属性かと思います。

image.png

つまり、バージョンによってヘッダーの情報が異なるし、ポイントフォーマットによって属性（カラム）の項目が異なるめんどくさいデータということがわかりますね！！！

ちなみに今回のデータはバージョン1.2のポイントフォーマット2なので、ヘッダーはこちら

image.png

属性はこちら

image.png

ということになります。

実際に情報が存在するか見てみましょう。
プロパティを見てみると、laspyのheaderインスタンスが持つメソッドなども一緒に表示され見づらいですが、おおむねバージョン1.2のヘッダー情報が存在していそうなことがわかるのではないでしょうか。

dir(header)
# [...
#  'xtra_dim',
#  'add_extra_dims',
#  'are_points_compressed',
#  'creation_date',
#  'extra_header_bytes',
#  'extra_vlr_bytes',
#  'file_source_id',
#  'generating_software',
#  'global_encoding',
#  'grow',
#  'major_version',
#  'maxs',
#  'minor_version',
#  'mins',
#  'number_of_evlrs',
#  'number_of_points_by_return',
#  'offset_to_point_data',
#  'offsets',
#  'partial_reset',
#  'point_count',
#  'point_format',
#  'read_from',
#  'remove_extra_dim',
#  'remove_extra_dims',
#  'scales',
#  'set_compressed',
#  'set_version_and_point_format',
#  'start_of_first_evlr',
#  'start_of_waveform_data_packet_record',
#  'system_identifier',
#  'update',
#  'uuid',
#  'version',
#  'vlrs',
#  'write_to',
#  'x_max',
#  'x_min',
#  'x_offset',
#  'x_scale',
#  'y_max',
#  'y_min',
#  'y_offset',
#  'y_scale',
#  'z_max',
#  'z_min',
#  'z_offset',
#  'z_scale']

次に、ポイントの属性を見てみましょう。
必須の属性がちゃんとあることが確認できますね。

point_format = las.point_format
list(point_format.dimension_names)
# ['X',
#  'Y',
#  'Z',
#  'intensity',
#  'return_number',
#  'number_of_returns',
#  'scan_direction_flag',
#  'edge_of_flight_line',
#  'classification',
#  'synthetic',
#  'key_point',
#  'withheld',
#  'scan_angle_rank',
#  'user_data',
#  'point_source_id',
#  'red',
#  'green',
#  'blue']

このように、一口にlasと言っても、バージョンやポイントフォーマットによって所有している情報が異なるので、気をつけましょう！

laspyでポイントの情報を読み取る

今度はヘッダーではなく、属性情報を読み取っていきましょう！
ポイントはこのように読み取ります。

las.points.array
# array([( 531578298,  855844897, 1295921591,   513, 0, 1, 0, 0, 29, 35445, 31365, 32640),
#        ( 528493437,  914345567, 1295115473,  2309, 0, 1, 0, 0, 29,  4080,  1530,  3825),
#        ( 539360560,  857806553, 1501347396, 64149, 0, 1, 0, 0, 29, 11475,  4080,  9180),
#        ...,
#        (2133125023, 2079379287,  744315865, 30840, 0, 1, 0, 0, 29, 31110, 30855, 29835),
#        (2132550117, 2078187852,  783591739, 40863, 0, 1, 0, 0, 29, 42585, 40290, 37485),
#        (2129016549, 2086335343,  756228502, 31868, 0, 1, 0, 0, 29, 33915, 31110, 29580)],
#       dtype=[('X', '<i4'), ('Y', '<i4'), ('Z', '<i4'), ('intensity', '<u2'), ('bit_fields', 'u1'), ('raw_classification', 'u1'), ('scan_angle_rank', 'i1'), ('user_data', 'u1'), ('point_source_id', '<u2'), ('red', '<u2'), ('green', '<u2'), ('blue', '<u2')])

ここで「あれ！？」と思った方は日頃から座標のデータに慣れ親しんでいる方かもしれません…
なんか気がつきましたか？
そうです！！！！！！！なんとXYZが整数値（int32）なんですね！！！！！！

las.X
# array([ 531578298,  528493437,  539360560, ..., 2133125023, 2132550117,
# 2129016549], dtype=int32)

これはlaspyのバグでも何でもなくて、lasデータは座標（というかポイント）をint32でしか持てないのが理由になります。
つまり仕様上の制限になります。

一見しただけではなかなか気づかず、ちょっと込み入った作業をしようとして、生データをそのまま利用してしまうと望まない結果になってしまうので、注意してくださいね。

で、どうやって小数点まで含めた座標を算出するのかというと、ここでヘッダーの情報が役に立ちます。

ヘッダーはスケールとオフセットというパラメーターをXYZごとに持っていて、これを処理することで正確なポイント座標を得ることができます。

# スケール・オフセットの意味
# 座標は整数: 531578298
x_first = las.X[0]

# xのスケール: 7.131602618438667e-08 -> 0.0000007
x_scale = header.x_scale

# xのoffset: -44528.753
x_offset = header.x_offset

# なので実際の座標は 531578298 * 0.0000007 + (-44528.753) = -44490.842948180776
x_coordinate = (x_first * x_scale) + x_offset

print(f"{x_first=}")
print(f"{x_scale=}")
print(f"{x_offset=}")
print(f"{x_coordinate=}")

# x_first=531578298
# x_scale=7.131602618438667e-08
# x_offset=-44528.753
# x_coordinate=-44490.842948180776

このようにして、座標が計算されます。
ではPDALを用いてこのデータが実際に正しいのか確認してみましょう。

PDALはコマンドラインから利用するツールです。（Pythonのラッパーが存在します。）
今回は利用しやすいコマンドラインツールを使っていきましょう。
（弊社代表のPDAL記事がわかりやすいのでこちらもどうぞ: PDALで点群データを処理）

以下は最初のポイントを取得するコマンドです。

% pdal info ./20190308掛川城.las -p 0      
{
  "file_size": 5001518281,
  "filename": "./20190308掛川城.las",
  "now": "2022-06-15T17:24:24+0900",
  "pdal_version": "2.4.1 (git-version: Release)",
  "points":
  {
    "point":
    {
      "Blue": 32640,
      "Classification": 1,
      "EdgeOfFlightLine": 0,
      "Green": 31365,
      "Intensity": 513,
      "NumberOfReturns": 0,
      "PointId": 0,
      "PointSourceId": 29,
      "Red": 35445,
      "ReturnNumber": 0,
      "ScanAngleRank": 0,
      "ScanDirectionFlag": 0,
      "UserData": 0,
      "X": -44490.84295,
      "Y": -135781.1752,
      "Z": 54.58493098
    }
  },
  "reader": "readers.las"
}

laspyから算出: -44490.842948180776
pdalで表示: -44490.84295

ということで、ほぼ一致しましたね！

ちなみにメタデータもがさっとpdalで見れます！
スケールやオフセットのほか、ファイルサイズ・ポイントの数・最大最小値など色々見れますね！便利！

❯ pdal info ./20190308掛川城.las --metadata
{
  "file_size": 5001518281,
  "filename": "./20190308掛川城.las",
  "metadata":
  {
    "comp_spatialreference": "",
    "compressed": false,
    "copc": false,
    "count": 192366079,
    "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": -44375.603,
    "maxy": -135673.849,
    "maxz": 73.5996558,
    "minor_version": 2,
    "minx": -44528.753,
    "miny": -135852.29,
    "minz": 25.648,
    "offset_x": -44528.753,
    "offset_y": -135852.29,
    "offset_z": 25.648,
    "point_length": 26,
    "project_id": "00000000-0000-0000-0000-650078006500",
    "scale_x": 7.13160261843867e-08,
    "scale_y": 8.30930658071832e-08,
    "scale_z": 2.23292297800368e-08,
    "software_id": "TREND-POINT\u0000E��\u0001\u001A��1;�\u0001\u0014\u0000f\u00005",
    "spatialreference": "",
    "srs":
    {
      "compoundwkt": "",
      "horizontal": "",
      "isgeocentric": false,
      "isgeographic": false,
      "prettycompoundwkt": "",
      "prettywkt": "",
      "proj4": "",
      "units":
      {
        "horizontal": "unknown",
        "vertical": ""
      },
      "vertical": "",
      "wkt": ""
    },
    "system_id": "FC\u0000B'�\u001E��\0012\u0002G�\fKdE�\u001B&\u0000\u0001\u0000&\u0000�\u0011"
  },
  "now": "2022-06-15T16:32:04+0900",
  "pdal_version": "2.4.1 (git-version: Release)",
  "reader": "readers.las"
}

PDALで点群を軽量化する

さて、色々触ってきましたが、元データが5GBもあると何をするにも挙動が重たいのでサクッと軽量化してしまいましょう！

Leafmap/Open3Dを使って掛川城の大規模点群データ（5GB）をPythonで可視化してみよう！ の記事ではplyに変換しましたが、点群処理専用ツールであるPDALなら当然las形式でも出力できます！

まずはこのようにしてpipeline.jsonを作成します。
元データのlasの点群を1/100000に間引いて、ヘッダーやスケールの情報を維持したまま、decimated_kakegawa.lasを作成する、というコマンドになります。
詳細はこちらからどうぞ。

[
  "20190308掛川城.las",
  {
    "type": "filters.decimation",
    "step": 100000
  },
  {
    "type": "writers.las",
    "filename": "decimated_kakegawa.las",
    "forward": "header,scale,offset,vlr"
  }
]

作成後、以下のコマンドを実行してください。

pdal pipeline pipeline.json

そうすると、decimated_kakegawa.lasが作成されるはずですので、PDALで情報を見てみましょう！

% pdal info decimated_kakegawa.las -p 0 
{
  "file_size": 50251,
  "filename": "decimated_kakegawa.las",
  "now": "2022-06-15T17:38:59+0900",
  "pdal_version": "2.4.0 (git-version: Release)",
  "points":
  {
    "point":
    {
      "Blue": 32640,
      "Classification": 1,
      "EdgeOfFlightLine": 0,
      "Green": 31365,
      "Intensity": 513,
      "NumberOfReturns": 0,
      "PointId": 0,
      "PointSourceId": 29,
      "Red": 35445,
      "ReturnNumber": 0,
      "ScanAngleRank": 0,
      "ScanDirectionFlag": 0,
      "UserData": 0,
      "X": -44490.84295,
      "Y": -135781.1752,
      "Z": 54.58493098
    }
  },
  "reader": "readers.las"
}

% pdal info decimated_kakegawa.las --metadata
{
  "file_size": 50251,
  "filename": "decimated_kakegawa.las",
  "metadata":
  {
    "comp_spatialreference": "",
    "compressed": false,
    "copc": false,
    "count": 1924,
    "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": -44377.447,
    "maxy": -135676.781,
    "maxz": 72.42308311,
    "minor_version": 2,
    "minx": -44526.901,
    "miny": -135848.501,
    "minz": 26.168,
    "offset_x": -44528.753,
    "offset_y": -135852.29,
    "offset_z": 25.648,
    "point_length": 26,
    "project_id": "00000000-0000-0000-0000-650078006500",
    "scale_x": 7.13160261843867e-08,
    "scale_y": 8.30930658071832e-08,
    "scale_z": 2.23292297800368e-08,
    "software_id": "TREND-POINT\u0000E��\u0001\u001A��1;�\u0001\u0014\u0000f\u00005",
    "spatialreference": "",
    "srs":
    {
      "compoundwkt": "",
      "horizontal": "",
      "isgeocentric": false,
      "isgeographic": false,
      "prettycompoundwkt": "",
      "prettywkt": "",
      "proj4": "",
      "units":
      {
        "horizontal": "unknown",
        "vertical": ""
      },
      "vertical": "",
      "wkt": ""
    },
    "system_id": "FC\u0000B'�\u001E��\0012\u0002G�\fKdE�\u001B&\u0000\u0001\u0000&\u0000�\u0011"
  },
  "now": "2022-06-15T17:39:22+0900",
  "pdal_version": "2.4.0 (git-version: Release)",
  "reader": "readers.las"
}

メタデータが維持された状態で軽量化できました！！！

ファイルサイズ: 5GB -> 50KB
ポイント数: 192366079個 -> 1924個

軽量化した点群を見てみる

最後に、leafmapで前回同様に可視化していきましょう!
可視化はこれだけでできるんでしたね。

import leafmap
leafmap.view_lidar("decimated_kakegawa.las", backend="open3d")

image.png

なにがなんだかわかんねぇ。
（かろうじて城があるのがわかるか…？くらい）

ということで間引きすぎには皆さん気をつけましょう！！！

お疲れ様でした！！

25
23
0

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
25
23