（COPC Viewerで可視化）

みなさんこんにちわ！！COPC触ってますか！！！

COPCというのはクラウド上で大規模点群データを扱うことを想定して作られたCloud OptimizedなLAZ形式のデータです。

HTTP Range Requestと組み合わせることで必要な部分だけをオンデマンドに取得可能なのが特徴です。

LAZ形式のデータや点群データの描画については他の記事で熱く語っているのでぜひご覧ください。

COPCには以下のような特徴があります。

8分木（Octree）の階層構造をメタデータに持つ
データ全体をボクセル状に区切り、それぞれを小さな圧縮チャンクとして保持している。
HTTP Range Requestへの最適化
ファイルの特定のバイト範囲だけを取得すれば、任意のボクセルに含まれる点群のみを読み出すことができる。
LAZ形式として取り扱える
点群そのものはLAZ形式で圧縮されており、部分的に解凍が可能。従来のLAS/LAZライブラリとも親和性がある。

ただのLAZ形式の場合でも、ローカルマシンで操作したり、ダウンロード目的でデータを配信する分には全く問題ないです。しかし、Web上の配信・可視化や、任意の領域・解像度に応じて効率的に利用するためには、適切にランダムアクセス可能なフォーマットが望ましく、そのために便利なのがCOPCです。

今回はこんなCOPCを可視化していこうと思います！

COPCのファイル構成

細かい仕様は全て公式サイトに書いています。

COPCファイルは拡張子は .copc.laz となっている場合が多いですが、必ずしもそうする必要はなく、内部的には「標準的なLAS/LAZのヘッダー + VLR(可変長レコード) + 階層構造 + チャンクデータ」という形になっています。

大まかな構成は以下のとおりです。

LASヘッダー
- LAS 1.4に対応したヘッダ情報が含まれます(ファイルシグネチャ、バージョン情報、ポイント数、スケール・オフセットなど)。
- COPCの特徴として、必ずLAS 1.4で、ポイントフォーマットが6〜8のいづれかと定められています。
VLR(Variable Length Record)群
- 通常のLAS/LAZ同様、拡張用のVLRが並んでいます。
- COPC特有のVLRとして、Copc Info VLRやCopc Hierarchy VRLなどがあります。
  - 特にCopc Info VLRはCOPC内部の階層データのルートを特定するために必要なオフセットやサイズを格納しています。
階層ヒエラルキー（Hierarchy）
- ルートからリーフノードまでの空間分割（Octree）が格納されています。
- 各ノードにおいて「どのバイト範囲に圧縮チャンクがあり、点数がいくつ含まれているか」などの情報がまとまっています。
LAZ圧縮チャンク
- 実際の点群データ（ポイントレコード）は、複数の小さなチャンクに分割され、LAZ圧縮された状態で格納されています。
- 階層ノードごとに参照できるようになっており、対応するバイト範囲を取り出して解凍すれば、そこに含まれる点群が取得できる仕組みになっています。

なぜレンジリクエストが有効なのか

とても巨大なCOPCファイルがあるとして、すべてダウンロードしようとすると、GB〜TB級のサイズになることも珍しくありません。しかし、COPCは階層構造を持ち、チャンクごとに「どこにどの点が含まれるか」が整理されています。

HTTP Range Requestを用いて、必要なヘッダーやVLR、ヒエラルキー、チャンクだけをオンデマンドで取得すれば、全体をダウンロードせずとも「特定レベルのデータだけ」「特定ボクセル(ノード)だけ」などのアクセスが可能になります。これにより、クラウド上の膨大な点群データを扱う際に、通信コストと待ち時間を大幅に削減できます。

特にAWS S3やCloudFrontでは、HTTP Range Requestに対応しているため、byte-range指定をすれば効率よく必要な範囲を取ってくることができます。

ヘッダーやVLRの取得

プロジェクトの立ち上げ

では早速データの解析をしていきましょう！
Python/Jupyter Labを利用して解析を行います。

今回は（あまり）外部モジュールに頼らず自力でパースしていくことを目的としています。このため、基本的にはstructモジュールなどを利用していきますが、部分的にライブラリを利用するためライブラリのインストールなどを行います。

パッケージマネージャーには超高速なuvを利用します。

OSによって差分はありますが、macOSであれば以下のようにするだけで良いです。

curl -LsSf https://astral.sh/uv/install.sh | sh
source $HOME/.cargo/env
mkdir my-project
cd my-project
uv init
uv python pin 3.12
uv add jupyterlab numpy copclib open3d requests
uv run --with jupyter jupyter lab

詳しい利用方法は公式サイトを参考にしてみてください。

COPCの作成

今回はこちらの記事で紹介した東京都の点群データを利用します。

任意の箇所のデータをダウンロードしてください。

その後、以下のようなコマンドでCOPCに変換します。

pdal translate -i 09LD1895.las -o 09LD1895.copc.laz --writers.copc.forward=all

変換後はS3などhttpリクエスト可能な任意のストレージに保存してください。

もしストレージを用意できない場合は、公式サイトにサンプルがいくつか配置されているのでそちらを利用してください。
（https://hobu-lidar.s3.amazonaws.com/sofi.copc.laz などが良いと思いますが、クリックした途端に2.3GBのCOPCのダウンロードが始まるのでご注意ください。）

（copc.io）

COPCの取得

まずは通常のLAS/LAZ共通のヘッダーを取得します。urlには任意のCOPCのURLを格納してください。
COPCではLAS 1.4と同様のヘッダーが利用されるため、375バイト固定長のヘッダーを持ちます。
requestモジュールを利用して先頭375バイトのみ取得してみましょう。

その後、必要な情報を収集していきます。

import numpy as np
import struct
import copclib as copc
from copclib import VectorChar
import open3d as o3d
import requests
import array


url = "https://xxxxxxxxxxx.cloudfront.net/09LD1895.copc.laz"

# 375バイト (LAS 1.4ヘッダー)を取得する
request_headers = {"Range": "bytes=0-374"}
response = requests.get(url, headers=request_headers)
las_header = response.content

# LASヘッダーを読み取る
# 最初の4バイトを読み、"LASF"か確認
magic = las_header[:4].decode("ascii")
if magic != "LASF":
    raise ValueError("Not a valid LAS file")
# バージョン情報を取得（位置24から2バイト）
version_major, version_minor = struct.unpack("<2B", las_header[24:26])
# ヘッダーサイズの取得（位置94）
header_size = struct.unpack("<H", las_header[94:96])[0]
print("header_size:", header_size)

# ポイントデータの開始オフセット（位置96）
offset_to_point_data = struct.unpack("<I", las_header[96:100])[0]
print("offset_to_point_data:", offset_to_point_data)

# VLRの個数
num_vlrs = struct.unpack("<I", las_header[100:104])[0]
print("num_vlrs:", num_vlrs)

# 点データ形式と点レコード長の取得（LAS 1.4ではタイプ6,7,8のみ）
point_format = struct.unpack("<B", las_header[104:105])[0]
point_record_length = struct.unpack("<H", las_header[105:107])[0]
print("point_format:", point_format)
print("point_record_length:", point_record_length)

# スケールとオフセットの取得
x_scale, y_scale, z_scale = struct.unpack("<3d", las_header[131:155])
x_offset, y_offset, z_offset = struct.unpack("<3d", las_header[155:179])
print(f"X scale: {x_scale}, Y scale: {y_scale}, Z scale: {z_scale}")
print(f"X offset: {x_offset}, Y offset: {y_offset}, Z offset: {z_offset}")

大体、以下のような情報が得られていると思います。
点群データをよく触っている方にとっては馴染みのある情報ではないかと思います。

header_size: 375
offset_to_point_data: 1368
num_vlrs: 3
point_format: 7
point_record_length: 36
X scale: 0.001, Y scale: 0.001, Z scale: 0.001
X offset: 0.0, Y offset: 0.0, Z offset: 0.0

今回特に重要な情報として

VLRが3つ
1つの点のバイト数が36バイト
ポイントフォーマットが7

という情報が得られました！

VLRの取得

次に、VLR（可変長レコード）を取得しますが、VLRはその名の通り可変長のため、取得が多少手間です。
まずは先ほど取得したヘッダー部分以降で、かつポイントレコード部直前までのデータを取得しましょう。

# VLRを取得する
request_headers = {"Range": f"bytes={header_size}-{offset_to_point_data - 1}"}
response = requests.get(url, headers=request_headers)
vlr_data = response.content

VLRはヘッダー部とデータ部の2項目に別れていますが、固定されているのはヘッダー部のみで、データ部はある程度自由に格納することができてしまいます。

幸いヘッダーの情報によりVLRは3つあること、VLRのヘッダー部は54バイトで固定されていることなどを参考に、以下のように解析できます。

# VLRヘッダーの構造（固定長54バイト）
vlt_header_struct = struct.Struct("<H 16s 2H 32s")

copc_info_vlr_data = None
laszip_vlr_data = None

# 順にVLRを取得していく
offset = 0
for _ in range(num_vlrs):
    vlr_header_bytes = vlr_data[offset : offset + vlt_header_struct.size]
    reserved, user_id_bytes, record_id, rec_length, description_bytes = (
        vlt_header_struct.unpack(vlr_header_bytes)
    )
    user_id = user_id_bytes.decode("ascii", errors="ignore").rstrip("\x00")
    description = description_bytes.decode("ascii", errors="ignore").rstrip("\x00")
    print(
        f"VLR {i}: User ID: {user_id}, Record ID: {record_id}, Length: {rec_length}, Description: {description}"
    )

    # それぞれ、仕様に応じたVLRヘッダーを持つ
    # VLR本体はVLRヘッダーの後に続き、rec_lengthバイト分のデータを持つ
    if user_id.lower() == "copc" and record_id == 1:
        print(f"Found COPC info VLR (record ID 1) with length {rec_length} bytes")
        copc_info_vlr_data = vlr_data[
            offset + vlt_header_struct.size : offset
            + vlt_header_struct.size
            + rec_length
        ]
        offset += vlt_header_struct.size + rec_length
    elif record_id == 22204:
        print(f"Found laszip VLR with length {rec_length} bytes")
        laszip_vlr_data = vlr_data[
            offset + vlt_header_struct.size : offset
            + vlt_header_struct.size
            + rec_length
        ]
        offset += vlt_header_struct.size + rec_length
    else:
        offset += vlt_header_struct.size + rec_length

# COPC仕様のVLR
if copc_info_vlr_data:
    info_struct = struct.Struct("<5d 2Q 2d 11Q")
    (
        center_x,
        center_y,
        center_z,
        halfsize,
        spacing,
        root_hier_offset,
        root_hier_size,
        gps_time_min,
        gps_time_max,
        *reserved_vals,
    ) = info_struct.unpack(copc_info_vlr_data)
    print("COPC Info VLR:")
    print("  Octree center XYZ:", (center_x, center_y, center_z))
    print("  Octree half-size:", halfsize)
    print("  Root spacing:", spacing)
    print("  Root hierarchy page offset:", root_hier_offset)
    print("  Root hierarchy page size:", root_hier_size)

# LAZデータに必須のVLR
if laszip_vlr_data:
    special_vlr_contents_fmt = "<2H 2B H I I 2q H"
    contents_size_special = struct.calcsize(special_vlr_contents_fmt)
    special_contents = struct.unpack(
        special_vlr_contents_fmt, laszip_vlr_data[:contents_size_special]
    )
    (
        compressor,
        coder,
        ver_major,
        ver_minor,
        ver_revision,
        options,
        chunk_size,
        num_special_evlrs,
        offset_special_evlrs,
        num_items,
    ) = special_contents

    print("LAZ Special VLR Contents:")
    print("  Compressor:", compressor)
    print("  Coder:", coder)
    print("  Version: {}.{}.{}".format(ver_major, ver_minor, ver_revision))
    print("  Options:", options)
    print("  Chunk Size (points per chunk):", chunk_size)
    print("  Number of special EVLRs:", num_special_evlrs)
    print("  Offset of special EVLRs:", offset_special_evlrs)
    print("  Number of Item records:", num_items)

    item_fmt = "<3H"
    item_size = struct.calcsize(item_fmt)
    items = []
    offset = contents_size_special
    for i in range(num_items):
        item_data = laszip_vlr_data[offset : offset + item_size]
        if len(item_data) < item_size:
            break
        item_type, item_size_val, item_version = struct.unpack(item_fmt, item_data)
        items.append((item_type, item_size_val, item_version))
        offset += item_size
    print("LAZ Special VLR Item Records:")
    for idx, (it_type, it_size, it_version) in enumerate(items):
        print(
            f"  Item {idx}: Type={it_type}, Size={it_size} bytes, Version={it_version}"
        )

実行すると、以下のような情報が出力されます。

VLR 6545844: User ID: copc, Record ID: 1, Length: 160, Description: COPC info VLR
Found COPC info VLR (record ID 1) with length 160 bytes
VLR 6545844: User ID: laszip encoded, Record ID: 22204, Length: 46, Description: lazperf variant
Found laszip VLR with length 46 bytes
VLR 6545844: User ID: LASF_Projection, Record ID: 2112, Length: 625, Description: 
COPC Info VLR:
  Octree center XYZ: (-5800.0005, -35800.0005, 187.60449999999992)
  Octree half-size: 199.9994999999999
  Root spacing: 2.72108163265306
  Root hierarchy page offset: 72685703
  Root hierarchy page size: 10240
LAZ Special VLR Contents:
  Compressor: 3
  Coder: 0
  Version: 3.4.3
  Options: 0
  Chunk Size (points per chunk): 4294967295
  Number of special EVLRs: -1
  Offset of special EVLRs: -1
  Number of Item records: 2
LAZ Special VLR Item Records:
  Item 0: Type=10, Size=30 bytes, Version=3
  Item 1: Type=11, Size=6 bytes, Version=3

VLRのヘッダーを解析すると、User IDやRecord ID、Record Length After Headerなどが取得できます。
COPCの場合は仕様として、必ず先頭のVLRのUser IDがcopcになっており、Record IDが1になっています。
また、Record Length After Headerが160バイトで固定されているのも特徴です。

LAZデータの場合は、Record IDが22204になっているようでした。

VLR以降のバイトで、なおかつRecord Length After Headerの範囲がVLRのデータ部になるため、その特徴を利用して、VLRのヘッダー部とデータ部を解析することができました。

ヒエラルキーVLRの取得

COPCはその特徴として、「階層情報」を持っているため、cloudに最適化されていると言えます。
先ほどは触れませんでしたが、COPC Info VRLのデータ部にはroot_hier_offsetとroot_hier_sizeの情報が含まれています。

このオフセット・サイズのバイナリを取得することでヒエラルキーVLRを得ることができます。

# ルートの階層ページを取得する
start_hier = root_hier_offset
end_hier = root_hier_offset + root_hier_size - 1

request_headers = {"Range": f"bytes={start_hier}-{end_hier}"}
response = requests.get(url, headers=request_headers)

hier_data = response.content
# 各ノードは32バイトのため、割ると数を取得できる
num_entries = len(hier_data) // 32

# 全てのノードを取得する
entries = []
for i in range(num_entries):
    rec = hier_data[i * 32 : (i + 1) * 32]
    level, ix, iy, iz, off, bsize, pcount = struct.unpack("<4i Q 2i", rec)
    entries.append(
        {
            "level": level,
            "x": ix,
            "y": iy,
            "z": iz,
            "offset": off,
            "byte_size": bsize,
            "point_count": pcount,
        }
    )

ポイントレコードの取得

LAS/LAZ（およびCOPC）での点群データは、「ポイントフォーマットID」や「1点あたりのバイト数（point_record_length）」などの情報によって構造が決まっています。今回のコードでは、まずは「どのノードのチャンクを取得するか」をヒエラルキーから調べたうえで、そのノードが持つチャンク（圧縮されたバイト範囲）を取得します。

ヒエラルキーの取得

COPCの場合、CopcInfoVlrから得られるroot_hier_offsetとroot_hier_sizeに基づき、階層データのルート部分を読み取ります。その中に32バイトごとの「ノード情報」が複数並んでいて、各ノードに対して以下のようなデータが含まれます。

level, x, y, z: ノードの空間的インデックス
offset, byte_size: このノードが格納するLAZ圧縮チャンクのファイル上のオフセットとバイトサイズ
point_count: このノードが含む点の数

# ヒエラルキーを取得
start_hier = root_hier_offset
end_hier = root_hier_offset + root_hier_size - 1
range_header_hier = {"Range": f"bytes={start_hier}-{end_hier}"}
res_hier = requests.get(url, headers=range_header_hier)
hier_data = res_hier.content

# 32バイト毎にノード情報をパース
num_entries = len(hier_data) // 32
entries = []
for i in range(num_entries):
    rec = hier_data[i * 32 : (i + 1) * 32]
    level, ix, iy, iz, off, bsize, pcount = struct.unpack("<4iQ2i", rec)
    entries.append(
        {
            "level": level,
            "x": ix,
            "y": iy,
            "z": iz,
            "offset": off,
            "byte_size": bsize,
            "point_count": pcount,
        }
    )

指定ノードのLAZ圧縮チャンクを取得

データは8分木になっているため、特定のチャンクを指定してデータを取得します。

（PLATEAU Learning TOPIC 24）

今回の例では(level=0, x=0, y=0, z=0)というルートのノードを探し、そのチャンク（LAZ圧縮部分）を取り出しています。ノードには圧縮チャンクへのオフセットとサイズが格納されているため、それを利用します。

# ルートノードを選択する
target_node = None
for e in entries:
    if e["level"] == 0 and e["x"] == 0 and e["y"] == 0 and e["z"] == 0:
        target_node = e
        break

if not target_node:
    raise RuntimeError("Could not find node.")

print("Selected node:", target_node)
point_count = target_node["point_count"]
if point_count <= 0:
    print("Warning: This node may not have actual points (point_count <= 0).")

chunk_start = target_node["offset"]
chunk_end = chunk_start + target_node["byte_size"] - 1

request_headers = {"Range": f"bytes={chunk_start}-{chunk_end}"}
response = requests.get(url, headers=request_headers)

chunk_data = response.content
print("Chunk size:", len(chunk_data), "(expected:", target_node["byte_size"], ")")

Selected node: {'level': 1, 'x': 0, 'y': 0, 'z': 0, 'offset': 63480802, 'byte_size': 2390957, 'point_count': 181830}
Chunk size: 2390957 (expected: 2390957 )

LAZの解凍 (外部モジュールを利用)

取得したデータはLAZ形式のポイントレコードとして圧縮されているため、そのままでは利用できません。
LAZの解凍も自力でやってみたかったですが、ちょっとすぐに出来るようなものでもないので外部ライブラリを利用します。色々種類はありますが、今回はcopclibを利用しました。

copclibのDecompressBytesを使って、圧縮チャンクを直接解凍しています。ここで注意したいのは、point_format_idやpoint_record_length が正しく設定されていること、加えてLAS仕様上のコアレコード長との差分(Extra Bytes)をどう扱うかなどです。

# copclibのDecompressBytes()を利用する
# LAS 1.4における標準のコアレコード長は、点形式によって異なる
# point format 7の場合はコアレコード長は30バイト、全体は36バイト → eb_byte_size=6
point_format_id = point_format & 0x3F
# eb_byte_size = point_record_length - 30
# だが、エラーが出るので、eb_byte_sizeは0を指定している
eb_byte_size = 0

arr = array.array("b", chunk_data)
v = VectorChar(arr)

print("point_format_id:", point_format_id)
print("eb_byte_size:", eb_byte_size)
print("point_count:", point_count)

uncompressed_data = copc.DecompressBytes(
    v,
    point_format_id,
    eb_byte_size,
    point_count,
)
print("Decompressed data size:", len(uncompressed_data))

# 解凍されたデータを解析して、座標データを取得
point_data = []
for i in range(0, len(uncompressed_data), point_record_length):
    x, y, z = struct.unpack("<3i", uncompressed_data[i : i + 12])
    x_real = x * x_scale + x_offset
    y_real = y * y_scale + y_offset
    z_real = z * z_scale + z_offset
    point_data.append([x_real, y_real, z_real])
# 中心値を求める
center = np.mean(point_data, axis=0)
print("Center of points:", center)

コメントにも書いている通り、動作はしているが想定外の挙動をする可能性がある、といった状況です。
もう少し丁寧に設定していく必要があるため、引き続き調査していきますが、今回はこのような設定でも結果として、解凍後のバイト列を得ることができました。

解凍されたので、X,Y,Zを取り出します。LASのスケールとオフセットを用いて、実際の座標値（メートル換算など）に変換するのが基本的な流れです。

計算された平均値は、概ねCOPC全体の中心に一致するため、少なくともXYZ座標は取得できていそうです。

point_format_id: 7
eb_byte_size: 0
point_count: 90284
Decompressed data size: 3250224
Center of points: [ -5790.53926006 -35843.95578703     37.93877992]

可視化

最後に、Open3Dなどを用いて、取得した座標群を点群として表示してみます！

points = np.array(point_data)

pcd = o3d.geometry.PointCloud()
pcd.points = o3d.utility.Vector3dVector(points)

o3d.visualization.draw_geometries([pcd])

これにより、S3(CloudFront経由)にあるCOPCファイルを部分的に読み込み、解凍、そしてローカルPC上で点群として表示するまでを一通り実現できました！

おわりに

Pythonでネットワーク上のCOPCに対してRange Requestで部分的にデータを取得し、可視化することができました！

COPCは比較的新しい技術ですが、クラウドを活用した大規模点群データ処理において今後も注目が高まると考えられます。ぜひご自身で試してみてください！

次はCOG（GeoTIFF）とかでもやってみたいですね！