MIERUNE Advent Calendar 2025

数値地図（国土基本情報）メッシュ標高情報データファイル (lem & csv)を高速にGeotiffに変換する

Last updated at 2025-12-20Posted at 2025-12-19

これは MIERUNE Advent Calendar 2025 の19日目の記事です。
昨日は @chizutodesign さんによる CLIが怖いデザイナーが、AIとの対話だけで「東京メトロの指令所っぽいモニタ」を作った話でした。

数値地図（国土基本情報）メッシュ標高情報データファイル (lem & csv)は日本独自の規格のようで国土地理院にその仕様が書いてある。

csvファイル
LEMファイルのメタデータ (座標系・グリッド情報)

測量年,2018
修正年,
東西方向の点数,4000
南北方向の点数,3000
東西方向のデータ間隔,0.5
南北方向のデータ間隔,0.5
区画左下の緯度,425332.207
区画左下の経度,1415425.842
区画右下の緯度,425332.461
区画右下の経度,1415553.995
区画右上の緯度,425421.074
区画右上の経度,1415553.745
区画左上の緯度,425420.820
区画左上の経度,1415425.573
図名,12od033
記録レコード数,3000
平面直角座標系番号,12
区画左下X座標,-12300000
区画左下Y座標,-2800000
区画右上X座標,-12150000
区画右上Y座標,-2600000
コメント,
...
...

lemファイル
標高データが固定幅で格納されている

2799- 2069 2068 2067 2066 2065 2063 2062 2061 ... .. 
2800- 2101 2098 2096 2094 2093 2091 2089 2087 ... ..   
...
...

Geotiffに変換できるWEBでしらべてみると、ググったところ以下のサイトが見つかった

最後はWindowsのネイティブ実装らしい、他はPythonでプログラムが作成されていた
実際に使ってみたがやはり大量のデータやサイズが大きいデータでは処理が遅かった
そこで、ネイティブ実装すればもっと早くできそうなのでやってみた。

Rustでの実装も検討してみたが、結局Geotiff変換にGDALのラッパーを使うので
それであれば全部C/C++でよい感じがするので、結局C++で実装、CUDAも使える

変換工程は以下

Lem、csvが入っているフォルダーを指定（サブフォルダー内にあってもオッケーとする）
csvメタデータを解析　（lemデータを数値化するのに必要）
lemデータを解析
Geotiffファイルを作成
複数のGeotiffファイルをまとめて一個のGeotiffにする

高速化手法

SIMDによる高速化
Intel TBB による並列化
NVIDIAのGPUがあればCUDAによるさらなる高速化

高速化ポイント

固定行データなので改行位置の検出をSIMDで行う、SIMDが使えない場合は単純ループ、その他にも利用

参考コード

 // AVX2 SIMDを使用して行オフセット（バイト位置）を検出
inline std::vector<size_t> find_line_offsets(const char* data, size_t size, size_t max_lines = 0) {
    std::vector<size_t> offsets;
    offsets.reserve(max_lines > 0 ? max_lines + 1 : size / 50);

    const char* end = data + size;
    const char* pos = data;

    // 最初の行は常にオフセット0から開始
    offsets.push_back(0);

#if HAS_AVX2
    const __m256i newline_mask = _mm256_set1_epi8('\n');

    while (pos + 32 <= end) {
        __m256i chunk = _mm256_loadu_si256(reinterpret_cast<const __m256i*>(pos));
        __m256i cmp = _mm256_cmpeq_epi8(chunk, newline_mask);
        uint32_t mask = static_cast<uint32_t>(_mm256_movemask_epi8(cmp));

        while (mask != 0) {
            int bit_pos = __builtin_ctz(mask);
            size_t offset = static_cast<size_t>(pos - data) + bit_pos + 1;

            if (offset < size) {
                offsets.push_back(offset);

                if (max_lines > 0 && offsets.size() > max_lines) {
                    return offsets;
                }
            }

            mask &= mask - 1;
        }

        pos += 32;
    }

    // 残りのバイトを処理
    while (pos < end) {
        if (*pos == '\n') {
            size_t offset = static_cast<size_t>(pos - data) + 1;
            if (offset < size) {
                offsets.push_back(offset);

                if (max_lines > 0 && offsets.size() > max_lines) {
                    return offsets;
                }
            }
        }
        ++pos;
    }

#else
    // フォールバック
    while (pos < end) {
        if (*pos == '\n') {
            size_t offset = static_cast<size_t>(pos - data) + 1;
            if (offset < size) {
                offsets.push_back(offset);

                if (max_lines > 0 && offsets.size() > max_lines) {
                    return offsets;
                }
            }
        }
        ++pos;
    }
#endif

    return offsets;
}

Intel TBB による並列化、パイプライン化は遅くなったので使わない

参考コード

        // 全ての.lemファイルを並列で.tifに変換（TBB enumerable_thread_specific使用）
        tbb::enumerable_thread_specific<std::vector<fs::path>> local_converted;

        tbb::parallel_for_each(dir_info.lem_files.begin(), dir_info.lem_files.end(),
                               [&](const fs::path& lem_file) {
                                   if (process_lem_file(lem_file, out_dir, nodata_value)) {
                                       fs::path tif_path = out_dir / (lem_file.stem().string() + ".tif");
                                       local_converted.local().push_back(tif_path);
                                   }
                               });

        // スレッドローカル結果をマージ
        std::vector<fs::path> converted_files;
        for (const auto& local : local_converted) {
            converted_files.insert(converted_files.end(), local.begin(), local.end());
        }

NVIDIAのGPUがあればCUDAによるさらなる高速化

参考カーネルコード

// LEMデータパース用CUDAカーネル
// 各スレッドが1つのセルを処理
__global__ void parse_lem_kernel(const char* file_data, const size_t* line_offsets, int nx, int ny,
                                 float nodata_value, float* output) {
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int col = blockIdx.x * blockDim.x + threadIdx.x;

    if (row >= ny || col >= nx)
        return;

    // 行開始位置を取得
    size_t line_start = line_offsets[row];
    size_t line_end = (row + 1 < ny) ? line_offsets[row + 1] : line_offsets[row] + nx * 5 + 10 + 2;

    // この列のデータ位置を計算
    // 各行は10文字のプレフィックスで始まり、その後値ごとに5文字
    int data_offset = 10 + col * 5;
    size_t char_pos = line_start + data_offset;

    // 位置が有効かチェック
    if (char_pos + 5 > line_end) {
        output[row * nx + col] = nodata_value;
        return;
    }

    // 5文字の整数値をパース
    int int_value = device_fast_parse_int(&file_data[char_pos], 5);

    // 適切なスケーリングでfloatに変換
    if (int_value == -1111 || int_value == -9999) {
        output[row * nx + col] = nodata_value;
    } else {
        output[row * nx + col] = static_cast<float>(int_value) * 0.1f;
    }
}