自治体が実施する航空レーザ測量の点群データ

自治体等が実施する航空レーザ測量の点群データは、主に

メッシュデータ
グランドデータ
オリジナルデータ

になります。ほかに航空写真や、分析した結果の等高線、図面などもあります。
成果品として納品される航空レーザの点群データはすべてテキストファイルとなっており、このテキストファイルのままでは、QGIS等では非常に使いづらいものになります。
そのため、このテキストファイルをLASやLAZといったフォーマットに変換することで、GISで使いやすいものになるのですが、このようなフォーマットでオープンデータとして公開されているものは稀です。
大抵は、標高DEMラスタ（TIFF）として、地表面のデータやCS立体図、赤色立体地図などの画像ファイルが公開されていることが多いです。
北海道庁が取得した航空レーザデータは、酪農学園大学の協力により、おもにLAZファイルでオープンデータを公開しています。
北海道オープンデータポータル　航空レーザ測量データ
 北海道の航空レーザデータダウンロードマップ

テキスト点群データをLAZに変換する

成果品として納品される航空レーザデータは、次のようなフォルダで保存されています。

1MCSV：メッシュデータ。はじめの「1M」は1mメッシュという意味
GROUND：メッシュ化する前の地表面のデータ
ORIGINAL：建物や樹木なども含めたオリジナルの点群データ
PHOTO：航空写真TIFFファイル

1MCSV、GORUND、ORIGINALの各フォルダには、テキストファイルで点群データが保存されています。

これらのフォルダ構成がただしければ、次のコード使って、PythonのPDALでGROUNDとORIGINALのテキストファイルをLAZファイルに一括変換できます。
LAZファイルにすると、ファイル容量も10分の1程度になります。

PDALの導入については、次の記事を参照してください。
https://qiita.com/koukita/items/f95f868e07c7cdbdf78d
また、スクリプトの実行は、次の記事を参考にしてください。
https://qiita.com/koukita/items/26bb23bcceffd7460769

このコードは、GROUND、ORIGINALをLAZに変換する際に、座標参照系の設定をします。
また、ORIGINALのLAZには航空写真から色データを取得します。
その後、5000分の1図郭ごとにファイルをまとめてZIP圧縮します。

一括でLAZに変換するコード

Lazer_txt2laz.py

import os
import zipfile
import shutil
import subprocess
import json
import pdal

#==============================
# ベースディレクトリ
base_dir = r"d:\aaa\bbbb" #航空レーザデータのフォルダを指定
d_year = "2020"  #西暦年度 ZIPファイルのファイル名に使用
d_name = "北見北部地区流木"  #省略した事業名　ZIPファイルのファイル名に使用

# 検索対象のフォルダ（優先順位あり）
folder_map = {
    "GROUND_LAZ": ["GROUND_LAZ", "GROUND_LAS", "GROUND"],
    "ORIGINAL_LAZ": ["ORIGINAL_LAZ", "ORIGINAL_LAS", "ORIGINAL"],
    "PHOTO_Jpeg": ["PHOTO_Jpeg", "PHOTO"]
}

# ZIP_CSVのフォルダリスト
subfolders = ["1MCSV","1MDTM", "DEM", "PHOTO_Jpeg", "WATER","1MCSMAP"]

#==============================

#LasからLazを作成する関数　GROUND、ORIGINAL共用
def las_to_laz(d_folder,las_folder,laz_folder,photo_path):
    las_files = [f for f in os.listdir(las_folder) if f.lower().endswith('.las')]

    for las_file in las_files:
        input_path = os.path.join(las_folder, las_file)
        base_name, _ = os.path.splitext(las_file)
        output_path = os.path.join(laz_folder, f"{base_name}.laz")

        if os.path.exists(output_path):
            print(f"スキップ: {output_path} は既に存在します。")
            continue
        
        #CRSの設定　北海道の場合でJGD2000平面直角座標系を設定
        if las_file.startswith("11"):
            epsg = 2453
        elif las_file.startswith("12"):
            epsg = 2454
        elif las_file.startswith("13"):
            epsg = 2455
        else:
            print(f"警告: {las_file} の EPSG を判定できませんでした。スキップします。")
            continue

        # PDAL パイプラインの構築
        match d_folder:
            case "GROUND": 
                pipeline = {
                    "pipeline": [
                        {"filename": input_path, "type": "readers.las", "spatialreference":f"EPSG:{epsg}"},
                        {"type": "writers.las", "filename": output_path}
                    ]
                }
            case "ORIGINAL": 
                pipeline = {
                    "pipeline": [
                        {"filename": input_path, "type": "readers.las", "spatialreference":f"EPSG:{epsg}"},
                        {"type": "filters.decimation", "step": 2}, #点数を1/2にする
                        {"type": "writers.las", "filename": output_path}
                    ]
                }

        json_file = os.path.join(base_dir, "pipeline.json")
        with open(json_file, "w") as f:
            json.dump(pipeline, f, indent=4)

        print(f"処理中: {las_file} -> {output_path}")
        pdal_command = ["pdal", "pipeline", json_file]
        result = subprocess.run(pdal_command, stdout=subprocess.PIPE, stderr=subprocess.PIPE)

        if result.returncode == 0:
            print(f"成功: {las_file} -> {output_path}")
        else:
            print(f"エラー: {las_file}\n{result.stderr.decode()}")

        #ORIGINALの場合は、航空写真から色を取得
        if d_folder == "ORIGINAL":
            tif_file_name = os.path.splitext(las_file)[0].replace("_org", "") + ".tif"
            tif_file_path = os.path.join(photo_path, tif_file_name)
                
            if os.path.exists(tif_file_path):
                print(f"Processing {las_file} with {tif_file_path}")
                pipeline_json = create_pipeline(output_path, tif_file_path, output_path)
                pipeline = pdal.Pipeline(pipeline_json)
                pipeline.execute()
            else:
                print(f"Warning: No matching TIF file for {las_file}")

#txtからLazを作成する関数　GROUND、ORIGINAL共用
def txt_to_laz(d_folder,txt_folder,laz_folder,photo_path):

    txt_files = [f for f in os.listdir(txt_folder) if f.lower().endswith('.txt')]

    for txt_file in txt_files:
        input_path = os.path.join(txt_folder, txt_file)
        base_name, _ = os.path.splitext(txt_file)
        output_path = os.path.join(laz_folder, f"{base_name}.laz")

        if os.path.exists(output_path):
            print(f"スキップ: {output_path} は既に存在します。")
            continue

        add_header_to_txt(d_folder,input_path)  # ヘッダー追加処理

        if txt_file.startswith("11"):
            epsg = 2453
        elif txt_file.startswith("12"):
            epsg = 2454
        elif txt_file.startswith("13"):
            epsg = 2455
        else:
            print(f"警告: {txt_file} の EPSG を判定できませんでした。スキップします。")
            continue

        # PDAL パイプラインの構築
        match d_folder:
            case "GROUND": 
                pipeline = {
                    "pipeline": [
                        {"type": "readers.text", "filename": input_path, "separator": ",", "spatialreference": f"EPSG:{epsg}"},
                        {"type": "writers.las", "filename": output_path, "compression": "true"}
                    ]
                }       
            case "ORIGINAL": 
                pipeline = {
                    "pipeline": [
                        {"type": "readers.text", "filename": input_path, "separator": ",", "spatialreference": f"EPSG:{epsg}"},
                        {"type": "filters.decimation", "step": 2},
                        {"type": "writers.las", "filename": output_path, "compression": "true"}
                    ]
                }       

        json_file = os.path.join(base_dir, "pipeline.json")
        with open(json_file, "w") as f:
            json.dump(pipeline, f, indent=4)

        print(f"処理中: {txt_file} -> {output_path}")
        pdal_command = ["pdal", "pipeline", json_file]
        result = subprocess.run(pdal_command, stdout=subprocess.PIPE, stderr=subprocess.PIPE)

        if result.returncode == 0:
            print(f"成功: {txt_file} -> {output_path}")
        else:
            print(f"エラー: {txt_file}\n{result.stderr.decode()}")

        if d_folder == "ORIGINAL":
            tif_file_name = os.path.splitext(txt_file)[0].replace("_org", "") + ".tif"
            tif_file_path = os.path.join(photo_path, tif_file_name)
                
            if os.path.exists(tif_file_path):
                print(f"Processing {txt_file} with {tif_file_path}")
                pipeline_json = create_pipeline(output_path, tif_file_path, output_path)
                pipeline = pdal.Pipeline(pipeline_json)
                pipeline.execute()
            else:
                print(f"Warning: No matching TIF file for {txt_file}")

def create_pipeline(laz_file, tif_file, output_file):
    """PDALのPipeline JSONを作成する"""
    pipeline = [
        {"filename": laz_file, "type": "readers.las"},
        {
            "type": "filters.colorization",
            "raster": tif_file,
            "dimensions": ["Red:1:1", "Green:2:1", "Blue:3:1"]
        },
        {"type": "writers.las", "filename": output_file}
    ]
    return json.dumps(pipeline)

#LAZに変換するためには、CSVのヘッダが必要なため、ヘッダを作成し上書き保存する
def add_header_to_txt(d_folder,file_path):
    """
    テキストファイルの先頭に "ID,X,Y,Z" を追加して上書き保存
    """
    try:
        with open(file_path, "r", encoding="utf-8") as f:
            lines = f.readlines()

        if lines and lines[0].strip() == "ID,X,Y,Z":
            return  # すでにヘッダーがある場合は何もしない

        if lines and lines[0].strip() == "ID,X,Y,Z,NO":
            return  # すでにヘッダーがある場合は何もしない 
                
        with open(file_path, "w", encoding="utf-8") as f:
            match d_folder:
                case "GROUND":
                    f.write("ID,X,Y,Z\n")

                case "ORIGINAL": 
                    f.write("ID,X,Y,Z,NO\n")

            f.writelines(lines)
    except Exception as e:
        print(f"エラー: {file_path} のヘッダー追加中に問題が発生しました: {e}")


#「GROUND_LAS」フォルダが有るか確認し、あれば「GROUND_LAS」からLAZファイルを作成する
#無ければ「GROUND」フォルダのテキストファイルからLAZを作成する。
las_path = os.path.join(base_dir, "GROUND_LAS")
txt_path = os.path.join(base_dir, "GROUND")
d_folder = "GROUND"
laz_folder = os.path.join(base_dir, "GROUND_LAZ")
os.makedirs(laz_folder, exist_ok=True)
if os.path.exists(las_path):
    #GROUND_LASからLAZファイルを作成する
    las_to_laz(d_folder,las_path,laz_folder,"")
else:
    if os.path.exists(txt_path):
        print(f"GROUND_LASフォルダが見つからないため、GROUNDフォルダからLAZを作成します。")
        #GROUNDからLAZファイルを作成する
        txt_to_laz(d_folder,txt_path,laz_folder,"")
    else:
        print(f"GROUNDフォルダが見つかりません: {las_path}")
    
#「ORIGINAL_LAS」フォルダが有るか確認し、あれば「ORIGINAL_LAS」からLAZファイルを作成する
#無ければ「ORIGINAL」フォルダのテキストファイルからLAZを作成する。
las_path = os.path.join(base_dir, "ORIGINAL_LAS")
txt_path = os.path.join(base_dir, "ORIGINAL")
photo_path = os.path.join(base_dir, "PHOTO")
d_folder = "ORIGINAL"
laz_folder = os.path.join(base_dir, "ORIGINAL_LAZ")
os.makedirs(laz_folder, exist_ok=True)
if os.path.exists(las_path):
    #GROUND_LASからLAZファイルを作成する
    las_to_laz(d_folder,las_path,laz_folder,photo_path)
else:
    if os.path.exists(txt_path):
        print(f"ORIGINAL_LASフォルダが見つからないため、ORIGINALフォルダからLAZを作成します。")
        #GROUNDからLAZファイルを作成する
        txt_to_laz(d_folder,txt_path,laz_folder,photo_path)
    else:
        print(f"ORIGINALフォルダが見つかりません: {las_path}")

#==============================
#ZIPファイルの作成
def m_zipfile(zip_filename, subfolders):
    # ZIPファイルの存在を確認
    if os.path.exists(zip_filename):
        print("同名のZIPファイルが存在します")
        return

    print(f"処理中: {zip_filename}")
    with zipfile.ZipFile(zip_filename, 'w', zipfile.ZIP_DEFLATED) as zipf:
        for folder in subfolders:
            target_folder = os.path.join(base_dir, folder)

            # フォルダが存在しない場合の処理（優先フォルダの切り替え）
            if folder in folder_map:
                for alt_folder in folder_map[folder]:
                    alt_path = os.path.join(base_dir, alt_folder)
                    if os.path.exists(alt_path):
                        target_folder = alt_path
                        break
                else:
                    print(f"{folder} および代替フォルダが見つかりません（スキップ）")
                    continue
            
            # フォルダが実際に存在するか確認
            if not os.path.exists(target_folder):
                print(f"{target_folder} が存在しません（スキップ）")
                continue

            # ファイルを検索してZIPに追加
            for file in os.listdir(target_folder):
                #if identifier in file:  # 6文字識別子を含むファイル
                if identifier.upper() in file.upper():  # 大文字・小文字を区別せず検索
                    file_path = os.path.join(target_folder, file)
                    zipf.write(file_path, os.path.relpath(file_path, base_dir))

    print(f"ZIP作成完了: {zip_filename}")


#===1MCSVほか=======================================

csv_dir = os.path.join(base_dir, "1MCSV")
#csv_dir = os.path.join(base_dir, "1MDTM")

# 1MCSVフォルダ内のTXTファイルを取得
if not os.path.exists(csv_dir):
    print("1MCSVフォルダが存在しません")
    exit()

txt_files = [f for f in os.listdir(csv_dir) if f.endswith(".txt")]
#txt_files = [f for f in os.listdir(csv_dir) if f.endswith(".csv")]

if not txt_files:
    txt_files = [f for f in os.listdir(csv_dir) if f.endswith(".TXT")]

if not txt_files:
    print("1MCSVフォルダにTXTファイルがありません")
    exit()


# 各TXTファイルの6文字識別子を取得
identifiers = {f[:6] for f in txt_files}

# すべての識別子に対してZIPを作成
for identifier in identifiers:
    zip_dir = os.path.join(base_dir, "ZIP_CSV")

    # ZIPフォルダを作成（存在しない場合）
    os.makedirs(zip_dir, exist_ok=True)

    zip_filename = os.path.join(zip_dir, f"{identifier}_CSVほか【{d_year}{d_name}】.zip")
    m_zipfile(zip_filename, subfolders)

#===GROUND=======================================
# ZIP_GROUNDのフォルダリスト
subfolders = ["GROUND_LAZ"]

csv_dir = os.path.join(base_dir, "GROUND_LAZ")

# GROUND_LAZフォルダ内のLAZファイルを取得
if not os.path.exists(csv_dir):
    print("GROUND_LAZフォルダが存在しません")
    exit()

txt_files = [f for f in os.listdir(csv_dir) if f.endswith(".laz")]

if not txt_files:
    print("LAZファイルがありません")
    exit()

# 各TXTファイルの6文字識別子を取得
identifiers = {f[:6] for f in txt_files}

# すべての識別子に対してZIPを作成
for identifier in identifiers:
    zip_dir = os.path.join(base_dir, "ZIP_GROUND")

    # ZIPフォルダを作成（存在しない場合）
    os.makedirs(zip_dir, exist_ok=True)

    zip_filename = os.path.join(zip_dir, f"{identifier}_GROUND_LAZ【{d_year}{d_name}】.zip")

    m_zipfile(zip_filename, subfolders)

#===ORIGINAL=======================================
# ZIP_ORIGINALのフォルダリスト
subfolders = ["ORIGINAL_LAZ"]

csv_dir = os.path.join(base_dir, "ORIGINAL_LAZ")

# ORIGINAL_LAZフォルダ内のLAZファイルを取得
if not os.path.exists(csv_dir):
    print("ORIGINAL_LAZフォルダが存在しません")
    exit()

txt_files = [f for f in os.listdir(csv_dir) if f.endswith(".laz")]

if not txt_files:
    print(f"LAZファイルがありません {txt_files}")
    exit()

# 各TXTファイルの6文字識別子を取得
identifiers = {f[:6] for f in txt_files}

# すべての識別子に対してZIPを作成
for identifier in identifiers:
    zip_dir = os.path.join(base_dir, "ZIP_ORIGINAL")

    # ZIPフォルダを作成（存在しない場合）
    os.makedirs(zip_dir, exist_ok=True)

    zip_filename = os.path.join(zip_dir, f"{identifier}_ORIGINAL_LAZ【{d_year}{d_name}】.zip")
    m_zipfile(zip_filename, subfolders)



print("すべてのZIP圧縮が完了しました。")

このコードをテキストエディタなどにコピーして、「Lazer_txt2laz.py」というファイルで保存してください。

スクリプトに実行方法

Windowsでのスクリプトの実行には「Anaconda Prompt」を使います。スタートメニューから起動してください。
プロンプトにconda activate myenvと入力するとPDALが利用できるようになります。

Lazer_txt2laz.pyを保存したディレクトリにCDで移動して、プロンプトにLazer_txt2laz.pyと入力して実行すると、変換処理が開始されます。

ファイルの変換にはかなりの時間がかかるので、はじめは３～４個のファイルで試すのがいいと思います。