關閉
標題: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";