LoginSignup
16
4

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 1 year has passed since last update.

@nokonoko_1203(のこのこ)inMIERUNE

超強力な「PDAL」おすすめコマンド3選・「メタデータの表示」「CRSの付与」「XY座標の入れ替え」

Posted at

image.png

みなさん「PDAL」ってご存知ですか？
知らないですよね。

弊社ではちょこちょこ記事を出していたり言及していたりするんですが、マジでマイナーツールだとは理解しています。

ちなみに「プードル🐩」って読むらしいです。（※諸説あり）

全然関係ないですが、弊社代表から初めてもらったUSBは「pdal」の文字が刻まれていました。
E10CDD10-054C-434D-A2A1-791E6FCB1DD3.JPG

使い方がむずい

で、このpdalですが、何をするツールかというと「点群データ」に対してとにかくなんでも出来るツールです。

「点群データ」と一口に言っても、ファイル形式は様々です。

  • ept
  • obj
  • pcd
  • ply
  • las
  • laz

など…

そんなありとあらゆる点群の読み書きを行うことができるほか、以下のような機能も兼ね備えています。

  • 点群データの間引き
  • クラスタリング
  • ラスター化
  • マージ
  • パイプラインによる自動化

など、その他100個くらいの機能…

詳しくはこの辺りを参照してください。
https://pdal.io/en/latest/stages/filters.html#

で、あまりにもリッチすぎるため、コマンドとオプションが多く使い方が複雑で、そもそもコマンドなんて覚えられんという感じではあるのですが、「どんな操作」を行うと「どうなるのか」を知っておくと、実際に利用するときに使いやすいかと思いますので、pdalでとにかくよく利用するおすすめコマンドを3つだけお伝えします！

データのダウンロード

と、その前にデータをダウンロードしてきましょう。

なんのデータでも良いのですが、僕は長崎県全域の点群データが自由に使えるらしいからみんなで使っちゃおうぜ！？にも書いた通り、オープンナガサキのデータをダウンロードします。

詳細は上記の記事を参照してほしいですが、今回は以下の区画をダウンロードしました。
1区画だけでも400MBを超えるのでダウンロードの際はご注意ください。

image.png

01KE9822.zipというファイルがダウンロードされますので、解凍して01ke9822_org.lasを取り出しておきましょう。

pdalのダウンロード

pdal自体のダウンロードは公式のQuickstartを参照してください。

基本的にはAnacondaでの利用をおすすめしております。
その他のインストール方法はハマりがちなので、ご注意ください。

Anacondaのインストールについては以下のような記事を書いていますので、参考にしてみてください。
MacでGISデータ分析を始めるためにサクッとAnacondaとjupyter labをインストールしてみる

メタデータの表示

01ke9822_org.lasを利用していきます。

ごく基本的な操作ですが、点群のメタ情報の確認にはpdal infoコマンドが利用できます。

% pdal info --metadata data/01ke9822_org.las
{
  "file_size": 631697631,
  "filename": "data/01ke9822_org.las",
  "metadata":
  {
    "comp_spatialreference": "",
    "compressed": false,
    "copc": false,
    "count": 24296054,
    "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": 35999.999,
    "maxy": -27000.001,
    "maxz": 147.515,
    "minor_version": 2,
    "minx": 35000,
    "miny": -27750,
    "minz": -3.841,
    "offset_x": 35000,
    "offset_y": -27750,
    "offset_z": -3.841,
    "point_length": 26,
    "project_id": "00000000-0000-0000-0000-000000000000",
    "scale_x": 4.65660821863294e-07,
    "scale_y": 3.49245499982147e-07,
    "scale_z": 7.04806298345703e-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-07-19T23:50:12+0900",
  "pdal_version": "2.5.0 (git-version: Release)",
  "reader": "readers.las"
}

上記より、点群のポイント数が24296054と膨大であることがわかったり、x座標の最大値が35999.999であることがわかったりします。

pdal infoはオプションが多くいろんなことができます。

  • 最初のポイントの情報を取得
% pdal info -p 0 data/01ke9822_org.las
{
  "file_size": 631697631,
  "filename": "data/01ke9822_org.las",
  "now": "2023-07-19T23:53:59+0900",
  "pdal_version": "2.5.0 (git-version: Release)",
  "points":
  {
    "point":
    {
      "Blue": 43690,
      "Classification": 1,
      "EdgeOfFlightLine": 0,
      "Green": 44975,
      "Intensity": 3812,
      "NumberOfReturns": 0,
      "PointId": 0,
      "PointSourceId": 29,
      "Red": 42662,
      "ReturnNumber": 0,
      "ScanAngleRank": 0,
      "ScanDirectionFlag": 0,
      "UserData": 0,
      "X": 35011.977,
      "Y": -27715.101,
      "Z": 34.56399993
    }
  },
  "reader": "readers.las"
}
  • 点群に付与された属性（座標・色など）の統計量を取得
% pdal info --stats data/01ke9822_org.las 
{
  "file_size": 631697631,
  "filename": "data/01ke9822_org.las",
  "now": "2023-07-19T23:55:00+0900",
  "pdal_version": "2.5.0 (git-version: Release)",
  "reader": "readers.las",
  "stats":
  {
    "bbox":
    {
      "native":
      {
        "bbox":
        {
          "maxx": 35999.999,
          "maxy": -27000.001,
          "maxz": 147.515,
          "minx": 35000,
          "miny": -27750,
          "minz": -3.841
        },
        "boundary": { "type": "Polygon", "coordinates": [ [ [ 35000.0, -27750.0, -3.841 ], [ 35000.0, -27000.001, -3.841 ], [ 35999.999000000003434, -27000.001, 147.515 ], [ 35999.999000000003434, -27750.0, 147.515 ], [ 35000.0, -27750.0, -3.841 ] ] ] }
      }
    },
    "statistic":
    [
      {
        "average": 35492.88025,
        "count": 24296054,
        "maximum": 35999.999,
        "minimum": 35000,
        "name": "X",
        "position": 0,
        "stddev": 283.5114865,
        "variance": 80378.76295
      },
      {
        "average": -27379.82475,
        "count": 24296054,
        "maximum": -27000.001,
        "minimum": -27750,
        "name": "Y",
        "position": 1,
        "stddev": 216.0244363,
        "variance": 46666.55709
      },
      {
        "average": 40.42072142,
        "count": 24296054,
        "maximum": 147.515,
        "minimum": -3.841,
        "name": "Z",
        "position": 2,
        "stddev": 30.51080463,
        "variance": 930.909199
      },
      {
        "average": 3415.683034,
        "count": 24296054,
        "maximum": 65535,
        "minimum": 47,
        "name": "Intensity",
        "position": 3,
        "stddev": 2963.964707,
        "variance": 8785086.785
      },
      {
        "average": 0,
        "count": 24296054,
        "maximum": 0,
        "minimum": 0,
        "name": "ReturnNumber",
        "position": 4,
        "stddev": 0,
        "variance": 0
      },
      {
        "average": 0,
        "count": 24296054,
        "maximum": 0,
        "minimum": 0,
        "name": "NumberOfReturns",
        "position": 5,
        "stddev": 0,
        "variance": 0
      },
      {
        "average": 0,
        "count": 24296054,
        "maximum": 0,
        "minimum": 0,
        "name": "ScanDirectionFlag",
        "position": 6,
        "stddev": 0,
        "variance": 0
      },
      {
        "average": 0,
        "count": 24296054,
        "maximum": 0,
        "minimum": 0,
        "name": "EdgeOfFlightLine",
        "position": 7,
        "stddev": 0,
        "variance": 0
      },
      {
        "average": 1,
        "count": 24296054,
        "maximum": 1,
        "minimum": 1,
        "name": "Classification",
        "position": 8,
        "stddev": 0,
        "variance": 0
      },
      {
        "average": 0,
        "count": 24296054,
        "maximum": 0,
        "minimum": 0,
        "name": "ScanAngleRank",
        "position": 9,
        "stddev": 0,
        "variance": 0
      },
      {
        "average": 0,
        "count": 24296054,
        "maximum": 0,
        "minimum": 0,
        "name": "UserData",
        "position": 10,
        "stddev": 0,
        "variance": 0
      },
      {
        "average": 29,
        "count": 24296054,
        "maximum": 29,
        "minimum": 29,
        "name": "PointSourceId",
        "position": 11,
        "stddev": 0,
        "variance": 0
      },
      {
        "average": 32901.45604,
        "count": 24296054,
        "maximum": 65535,
        "minimum": 7453,
        "name": "Red",
        "position": 12,
        "stddev": 12734.92859,
        "variance": 162178406.1
      },
      {
        "average": 36864.89409,
        "count": 24296054,
        "maximum": 65535,
        "minimum": 15420,
        "name": "Green",
        "position": 13,
        "stddev": 10848.22869,
        "variance": 117684065.8
      },
      {
        "average": 36389.1289,
        "count": 24296054,
        "maximum": 65535,
        "minimum": 8224,
        "name": "Blue",
        "position": 14,
        "stddev": 9390.640171,
        "variance": 88184122.81
      }
    ]
  }
}
  • BBOXを取得する
% pdal info --boundary data/01ke9822_org.las
{
  "boundary":
  {
    "area": 2276732.212,
    "avg_pt_per_sq_unit": 0.3041666281,
    "avg_pt_spacing": 0.306117462,
    "boundary": "POLYGON ((35745.951 -28477.869,36333.13 -27969.357,36333.13 -26952.333,35011.977 -26698.077,34424.798 -27715.101,35745.951 -28477.869))",
    "boundary_json": { "type": "Polygon", "coordinates": [ [ [ 35745.951123560000269, -28477.869284540000081 ], [ 36333.130422540001746, -27969.357094920000236 ], [ 36333.130422540001746, -26952.332715669999743 ], [ 35011.976999830003479, -26698.076620859999821 ], [ 34424.797700850002002, -27715.10100009999951 ], [ 35745.951123560000269, -28477.869284540000081 ] ] ] },
    "density": 10.6714588,
    "edge_length": 0,
    "estimated_edge": 508.5121896,
    "hex_offsets": "MULTIPOINT (0 0, -146.795 254.256, 0 508.512, 293.59 508.512, 440.384 254.256, 293.59 0)",
    "sample_size": 5000,
    "threshold": 15
  },
  "file_size": 631697631,
  "filename": "data/01ke9822_org.las",
  "now": "2023-07-19T23:56:47+0900",
  "pdal_version": "2.5.0 (git-version: Release)",
  "reader": "readers.las"
}

詳しいことはhelpを確認してみてください。

% pdal info --help                          
usage: pdal info [options] [input]
standard options:
  --developer-debug  Enable developer debug (don't trap exceptions)
  --label            A string to label the process with
  --driver           Override reader driver
  --help             Print help and exit
options:
  --input, -i               Input file name
  --all                     Dump statistics, schema and metadata
  --point, -p               Point to dump --point="1-5,10,100-200" (0 indexed)
  --query                   Return points in order of distance from the specified location (2D or 3D) --query Xcoord,Ycoord[,Zcoord][/count]
  --stats                   Dump stats on all points (reads entire dataset)
  --boundary                Compute a hexagonal hull/boundary of dataset
  --dimensions              Dimensions on which to compute statistics
  --enumerate               Dimensions whose values should be enumerated
  --schema                  Dump the schema
  --pipeline-serialization  Output filename for pipeline serialization
  --summary                 Dump summary of the info
  --metadata                Dump file metadata info
  --stdin, -s               Read a pipeline file from standard input

For more information, see the full documentation for PDAL at http://pdal.io/

CRSの付与

GIS屋さんなら馴染み深いと思いますが、いわゆる「座標系」の設定もできます。

「メタデータの表示」で確認した通り、デフォルトの状態ではCRSが付与されていないデータがほとんどです。

% pdal info --metadata data/01ke9822_org.las
{
  "file_size": 631697631,
  "filename": "data/01ke9822_org.las",
  "metadata":
  {
    "comp_spatialreference": "",
    "compressed": false,
    "copc": false,
    "count": 24296054,
    "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": 35999.999,
    "maxy": -27000.001,
    "maxz": 147.515,
    "minor_version": 2,
    "minx": 35000,
    "miny": -27750,
    "minz": -3.841,
    "offset_x": 35000,
    "offset_y": -27750,
    "offset_z": -3.841,
    "point_length": 26,
    "project_id": "00000000-0000-0000-0000-000000000000",
    "scale_x": 4.65660821863294e-07,
    "scale_y": 3.49245499982147e-07,
    "scale_z": 7.04806298345703e-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-07-19T23:50:12+0900",
  "pdal_version": "2.5.0 (git-version: Release)",
  "reader": "readers.las"
}

特に座標系っぽいものは見つかりませんね。

なので付与していきましょう。

今回のデータは長崎県でJGD2011ですので、EPSG:6669になります。

pdal translate -i data/01ke9822_org.las -o data/01ke9822_crs.las --writers.las.a_srs="EPSG:6669"

生成されたdata/01ke9822_crs.lasのメタデータを確認してみましょう。

% pdal info --boundary data/01ke9822_org.las
{
  "boundary":
  {
    "area": 2276732.212,
    "avg_pt_per_sq_unit": 0.3041666281,
    "avg_pt_spacing": 0.306117462,
    "boundary": "POLYGON ((35745.951 -28477.869,36333.13 -27969.357,36333.13 -26952.333,35011.977 -26698.077,34424.798 -27715.101,35745.951 -28477.869))",
    "boundary_json": { "type": "Polygon", "coordinates": [ [ [ 35745.951123560000269, -28477.869284540000081 ], [ 36333.130422540001746, -27969.357094920000236 ], [ 36333.130422540001746, -26952.332715669999743 ], [ 35011.976999830003479, -26698.076620859999821 ], [ 34424.797700850002002, -27715.10100009999951 ], [ 35745.951123560000269, -28477.869284540000081 ] ] ] },
    "density": 10.6714588,
    "edge_length": 0,
    "estimated_edge": 508.5121896,
    "hex_offsets": "MULTIPOINT (0 0, -146.795 254.256, 0 508.512, 293.59 508.512, 440.384 254.256, 293.59 0)",
    "sample_size": 5000,
    "threshold": 15
  },
  "file_size": 631697631,
  "filename": "data/01ke9822_org.las",
  "now": "2023-07-19T23:56:47+0900",
  "pdal_version": "2.5.0 (git-version: Release)",
  "reader": "readers.las"
}
(ml) 
satoru@[23:56:47]:~/go/src/github.com/nokonoko1203/0002_study/sample% pdal info --help                          
usage: pdal info [options] [input]
standard options:
  --developer-debug  Enable developer debug (don't trap exceptions)
  --label            A string to label the process with
  --driver           Override reader driver
  --help             Print help and exit
options:
  --input, -i               Input file name
  --all                     Dump statistics, schema and metadata
  --point, -p               Point to dump --point="1-5,10,100-200" (0 indexed)
  --query                   Return points in order of distance from the specified location (2D or 3D) --query Xcoord,Ycoord[,Zcoord][/count]
  --stats                   Dump stats on all points (reads entire dataset)
  --boundary                Compute a hexagonal hull/boundary of dataset
  --dimensions              Dimensions on which to compute statistics
  --enumerate               Dimensions whose values should be enumerated
  --schema                  Dump the schema
  --pipeline-serialization  Output filename for pipeline serialization
  --summary                 Dump summary of the info
  --metadata                Dump file metadata info
  --stdin, -s               Read a pipeline file from standard input

For more information, see the full documentation for PDAL at http://pdal.io/

(ml) 
satoru@[23:57:38]:~/go/src/github.com/nokonoko1203/0002_study/sample% pdal translate -i data/01ke9822_org.las -o data/01ke9822_crs.las  
^C
(ml) 
satoru@[0:03:05]:~/go/src/github.com/nokonoko1203/0002_study/sample% pdal translate -i data/01ke9822_org.las -o data/01ke9822_crs.las --writers.las.a_srs="EPSG:6669" 
(ml) 
satoru@[0:03:22]:~/go/src/github.com/nokonoko1203/0002_study/sample% pdal info --metadata data/01ke9822_crs.las                                                      
{
  "file_size": 826066283,
  "filename": "data/01ke9822_crs.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": 24296054,
    "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": 35999.999,
    "maxy": -27000.001,
    "maxz": 147.5150002,
    "minor_version": 2,
    "minx": 35000,
    "miny": -27750,
    "minz": -3.841,
    "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-07-20T00:04:04+0900",
  "pdal_version": "2.5.0 (git-version: Release)",
  "reader": "readers.las"
}

ものすごい情報量が増えましたね！

これでCRSの付与が完了です。

これをしておくことによって、例えばCesium ionなんかに適当にアップロードした時に、正しい位置に表示される3DTilesを作成してくれたりします！

XY座標の入れ替え

日本でオープンに公開されている点群データのほとんどは、平面直角座標系で生成されたデータですが、平面直角座標系は「縦軸（つまり南北方向）がX」という直感とは反する座標軸を持っています。

一般的な3次元表示ソフトは数学座標（つまりX軸が横軸）のため、点群データを読み込ませると、XとYが反転して表示されてしまいます。

これを防ぐために、ほとんどのケースでXY軸を反転させる必要が出てきます。

が、pdalはこれを1コマンドで処理できないので、「パイプライン」を作成する必要があります。

まずはxy_switch_pipeline.jsonのような適当なJSONを作成し、以下のように記述します。

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

その後、以下のコマンドを実行します。

pdal pipeline xy_switch_pipeline.json

変換前後の最初のポイントを比較してみましょう。

  • 変換前
% pdal info -p 0 data/01ke9822_crs.las                                                            
{
  "file_size": 826066283,
  "filename": "data/01ke9822_crs.las",
  "now": "2023-07-20T00:12:52+0900",
  "pdal_version": "2.5.0 (git-version: Release)",
  "points":
  {
    "point":
    {
      "Blue": 43690,
      "Classification": 1,
      "EdgeOfFlightLine": 0,
      "GpsTime": 0,
      "Green": 44975,
      "Intensity": 3812,
      "NumberOfReturns": 0,
      "PointId": 0,
      "PointSourceId": 29,
      "Red": 42662,
      "ReturnNumber": 0,
      "ScanAngleRank": 0,
      "ScanDirectionFlag": 0,
      "UserData": 0,
      "X": 35011.98,
      "Y": -27715.1,
      "Z": 34.56
    }
  },
  "reader": "readers.las"
}
  • 変換後
% pdal info -p 0 data/01ke9822_swap.las 
{
  "file_size": 826066283,
  "filename": "data/01ke9822_swap.las",
  "now": "2023-07-20T00:13:17+0900",
  "pdal_version": "2.5.0 (git-version: Release)",
  "points":
  {
    "point":
    {
      "Blue": 43690,
      "Classification": 1,
      "EdgeOfFlightLine": 0,
      "GpsTime": 0,
      "Green": 44975,
      "Intensity": 3812,
      "NumberOfReturns": 0,
      "PointId": 0,
      "PointSourceId": 29,
      "Red": 42662,
      "ReturnNumber": 0,
      "ScanAngleRank": 0,
      "ScanDirectionFlag": 0,
      "UserData": 0,
      "X": -27715.1,
      "Y": 35011.98,
      "Z": 34.56
    }
  },
  "reader": "readers.las"
}

XとYの値が入れ替わっているのが分かるかと思います！

おわりに

他にもたくさん紹介したいコマンドはあるので、今後もちょこちょこpdal記事を書いていきたいと思います。

点群データを使うならめちゃめちゃ便利なツールですので、みなさんも利用してみてください！

お題は不問！Qiita Engineer Festa 2023で記事投稿！
16
4
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
16
4

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?