關閉      標題:grd 20M 地形圖轉成 sqlite 可發圖資 與 tif 可轉 cesium 地形資料
內容:

<?php
// grd_to_merged_tif.php

$fps = glob("*.grd");
$max_k = count($fps);

foreach ($fps as $k => $v) {
    $mn = pathinfo($v, PATHINFO_FILENAME);

    $tif_file = "{$mn}.tif";
    $fix_file = "{$mn}_fix.tif";

    if (!is_file($tif_file)) {
        $CMD = 'gdal_translate -of GTiff -a_srs EPSG:3826 "' . $v . '" "' . $tif_file . '"';
        exec($CMD, $output, $return_var);

        if ($return_var !== 0) {
            echo "⚠️ Failed to translate {$v}. Attempting vector-rescue...\n";

            $csv_file = "{$mn}.csv";
            $vrt_file = "{$mn}.vrt";

            $data = "X,Y,Z\n";
            $grd_data = preg_replace("/\s+/", ",", trim(file_get_contents($v)));
            file_put_contents($csv_file, $data . $grd_data);

            $vrt_template = <<<EOF
<OGRVRTDataSource>
  <OGRVRTLayer name="{$mn}">
    <SrcDataSource>{$csv_file}</SrcDataSource>
    <GeometryType>wkbPoint</GeometryType>
    <LayerSRS>EPSG:3826</LayerSRS>
    <GeometryField encoding="PointFromColumns" x="X" y="Y" z="Z"/>
  </OGRVRTLayer>
</OGRVRTDataSource>
EOF;
            file_put_contents($vrt_file, $vrt_template);

            $rescue_cmd = 'gdal_rasterize -a Z -tr 20 20 -a_nodata -32768 -ot Float32 -a_srs EPSG:3826 "' . $vrt_file . '" "' . $tif_file . '"';
            exec($rescue_cmd, $out, $rescue_var);

            @unlink($csv_file);
            @unlink($vrt_file);

            if ($rescue_var === 0 && is_file($tif_file)) {
                echo "✅ Successfully rescued {$v}!\n";
            } else {
                echo "❌ Rescue failed for {$v}. Skipping warp step.\n";
                continue;
            }
        }
    }

    echo sprintf("%d / %d ... %s\n", ($k + 1), $max_k, $v);

    if (is_file($tif_file) && !is_file($fix_file)) {
        $CMD = 'gdalwarp -r near -dstnodata -32768 "' . $tif_file . '" "' . $fix_file . '"';
        exec($CMD, $warp_out, $warp_ret);

        if ($warp_ret !== 0 || !is_file($fix_file)) {
            echo "❌ Warp failed for {$tif_file}\n";
            continue;
        }
    }
}

// 清理舊輸出
$cleanup_files = ["all.vrt", "list.txt", "merged_dem.tif"];
foreach ($cleanup_files as $file) {
    if (is_file($file)) {
        unlink($file);
    }
}

// 建立 all.vrt
$fix_tifs = glob("*_fix.tif");
if (empty($fix_tifs)) {
    die("No valid _fix.tif files generated. Exiting.\n");
}

file_put_contents("list.txt", implode("\n", $fix_tifs));

$CMD = 'gdalbuildvrt -input_file_list list.txt all.vrt';
exec($CMD, $vrt_out, $vrt_ret);

if ($vrt_ret !== 0 || !is_file("all.vrt")) {
    echo implode("\n", $vrt_out) . "\n";
    die("Failed to build all.vrt\n");
}

// 輸出合併後 GeoTIFF
echo "Building merged GeoTIFF...\n";

$CMD = 'gdal_translate -of GTiff '
     . '-a_nodata -32768 '
     . '-co TILED=YES '
     . '-co COMPRESS=LZW '
     . '-co BIGTIFF=IF_SAFER '
     . 'all.vrt merged_dem.tif';

exec($CMD, $merge_out, $merge_ret);

if ($merge_ret !== 0 || !is_file("merged_dem.tif")) {
    echo implode("\n", $merge_out) . "\n";
    die("Failed to build merged_dem.tif\n");
}

echo "✅ Process completed. Output: merged_dem.tif\n";