はじめに

この記事は FOSS4G Advent Calendar 2023 3日目の記事です。

GIS界隈ではCloud Optimizedが重要なキーワードになっています。Cloud Optimizedとは、特別なサーバ実装を必要とせず、大きな位置情報データの一部分を配信することを可能とするファイル形式の総称です。本記事は、Cloud Optimized形式のうち、Cloud Optimized Point Cloud（以下、COPC）の生成方法に焦点を当てた記事になります。具体には、東京都（東京都デジタルツイン実現プロジェクト）がG空間情報センターを通じてオープンデータとして公開している、点群データ（オリジナルデータ（DSM）及びグラウンドデータ（DEM））からCOPCを生成する方法について説明します。

アウトプットイメージ

デモサイト

Ctrlを押しながらドラッグで回転できます。
色はデフォルトでElevationになっているのでパネルからRGBに変更すると実際の色で表示できます。
※PCでの閲覧専用です

Cloud Optimized Point Cloud（COPC）の特徴

COPC（Cloud Optimized Point Cloud）はクラウド最適化点群仕様に基づくデータ形式です。
ポイントクラウドデータがOctress構造で保存されることによって、クライアントの必要に応じて一部のデータを抜き出すことができます。

PC環境

OS：Windows 10 Pro
CPU：Intel Core i7-9700
メモリ：32GB

前提条件

OSGeo4W（PDAL）がインストール済みであること。
PDALのバージョンは、pdal 2.5.2です。

G空間情報センターから点群データ（LAS形式）を取得済みであること。

・本記事では、東京都奥多摩町の小河内ダム（おごうちだむ）周辺の点群データ（LAS）を使用します。

LASのマージ

LASのマージには、OSGeo4WのPDALのパイプラインを使用します。
LASは、下記のとおり、lasフォルダ内に格納されているとします。

├─las
│      09KC7262.las
│      09KC7263.las
│      09KC7272.las
│      09KC7273.las

まず、パイプラインのjsonを下記のとおり作成します。

merge-pipeline.json

{
  "pipeline": [
    {
      "type": "readers.las",
      "filename": "las/*.las"
    },
    {
      "type": "writers.las",
      "filename": "ogochi-dam.las"
    }
  ]
}

OSGeo4W Shellを起動して、下記のコマンドを実行して、lasをマージします。

pdal pipeline merge-pipeline.json

座標参照系の付与

G空間情報センターから取得した、点群データ（LAS形式）には座標参照系が設定されていないため、PDALのtranslateを使用して、座標参照系を付与します。
なお、付与する座標参照系は、日本測地系2011／平面直角座標系第9系（EPSG：6677）になります。
OSGeo4W Shellを起動して、下記のコマンドを実行して座標参照系を付与します。

pdal translate -i ogochi-dam.las -o ogochi-dam-translated.las --writers.las.a_srs="EPSG:6677"

また、下記のコマンドを実行して、座標参照系を付与されているかGeoTIFFのメタ情報を確認します。

pdal info --metadata ogochi-dam-translated.las

メタ情報の出力結果は下記のとおりになります。

{
  "file_size": 886509750,
  "filename": "ogochi-dam-translated.las",
  "metadata":
  {
    "comp_spatialreference": "PROJCS[\"JGD2011 / Japan Plane Rectangular CS IX\",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\",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[\"Northing\",NORTH],AXIS[\"Easting\",EAST],AUTHORITY[\"EPSG\",\"6677\"]]",
    "compressed": false,
    "copc": false,
    "count": 26073803,
    "creation_doy": 0,
    "creation_year": 2023,
    "dataformat_id": 3,
    "dataoffset": 448,
    "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,40): \"JGD2011 / Japan Plane Rectangular CS IX\"\n      GeogCitationGeoKey (Ascii,8): \"JGD2011\"\n      GeogAngularUnitsGeoKey (Short,1): Angular_Degree\n      ProjectedCSTypeGeoKey (Short,1): Code-6677 (JGD2011 / Japan Plane Rectangular CS IX)\n      ProjLinearUnitsGeoKey (Short,1): Linear_Meter\n      End_Of_Keys.\n   End_Of_Geotiff.\n",
    "header_size": 227,
    "major_version": 1,
    "maxx": -70400,
    "maxy": -22800,
    "maxz": 625.95,
    "minor_version": 2,
    "minx": -71200,
    "miny": -23400,
    "minz": 382.82,
    "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.2 (57c4e7)",
    "spatialreference": "PROJCS[\"JGD2011 / Japan Plane Rectangular CS IX\",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\",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[\"Northing\",NORTH],AXIS[\"Easting\",EAST],AUTHORITY[\"EPSG\",\"6677\"]]",
    "srs":
    {
      "compoundwkt": "PROJCS[\"JGD2011 / Japan Plane Rectangular CS IX\",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\",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[\"Northing\",NORTH],AXIS[\"Easting\",EAST],AUTHORITY[\"EPSG\",\"6677\"]]",
      "horizontal": "PROJCS[\"JGD2011 / Japan Plane Rectangular CS IX\",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\",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[\"Northing\",NORTH],AXIS[\"Easting\",EAST],AUTHORITY[\"EPSG\",\"6677\"]]",
      "isgeocentric": false,
      "isgeographic": false,
      "json": "{\n  \"type\": \"ProjectedCRS\",\n  \"name\": \"JGD2011 / Japan Plane Rectangular CS IX\",\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\": 36,\n        \"unit\": \"degree\",\n        \"id\": {\n          \"authority\": \"EPSG\",\n          \"code\": 8801\n        }\n      },\n      {\n        \"name\": \"Longitude of natural origin\",\n        \"value\": 139.833333333333,\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\": 6677\n  }\n}",
      "prettycompoundwkt": "PROJCS[\"JGD2011 / Japan Plane Rectangular CS IX\",\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\",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[\"Northing\",NORTH],\n    AXIS[\"Easting\",EAST],\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.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\",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[\"Northing\",NORTH],\n    AXIS[\"Easting\",EAST],\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.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\",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[\"Northing\",NORTH],AXIS[\"Easting\",EAST],AUTHORITY[\"EPSG\",\"6677\"]]"
    },
    "system_id": "PDAL",
    "vlr_0":
    {
      "data": "AQABAAAABwAABAAAAQABAAEEAAABAAEAAgSxhygAAAABCLGHCAAoAAYIAAABAI4jAAwAAAEAFRoEDAAAAQApIw==",
      "description": "GeoTiff GeoKeyDirectoryTag",
      "record_id": 34735,
      "user_id": "LASF_Projection"
    },
    "vlr_1":
    {
      "data": "SkdEMjAxMSAvIEphcGFuIFBsYW5lIFJlY3Rhbmd1bGFyIENTIElYfEpHRDIwMTF8AA==",
      "description": "GeoTiff GeoAsciiParamsTag",
      "record_id": 34737,
      "user_id": "LASF_Projection"
    }
  },
  "now": "2023-10-06T20:36:05+0900",
  "pdal_version": "2.5.2 (git-version: 57c4e7)",
  "reader": "readers.las"
}

COPCの生成

LASからCOPCへの変換にはPDALのtranslateを使用します。
OSGeo4W Shellを起動して、下記のコマンドを実行してCOPCを生成します。

pdal translate -i ogochi-dam-translated.las -o ogochi-dam-translated.copc.laz --writers.copc.forward=all

COPCの表示

上記の方法で生成したCOPCをCOPC Viewerで表示します。

以下のように直接URLを指定してアクセスすることができます。
https://viewer.copc.io/?copc=https://xs489works.xsrv.jp/copc-data/tokyo-digitaltwin/ogochi-dam-translated.copc.laz
Ctrlを押しながらドラッグで回転できます。
色はデフォルトでElevationになっているのでパネルからRGBに変更すると実際の色で表示できます。