關閉
標題：dem 20 m RTree gdr 轉資料庫 可作 dem 地形高查詢
內容：
把 *.gdr / *.grd 轉成 csv，再合併成 taiwan_gdr.db (EPSG:3826)
// 並建立快速最近點查詢用的 RTree 索引
<?php
// build_taiwan_gdr_db.php
// 需求：把 *.gdr / *.grd 轉成 csv，再合併成 taiwan_gdr.db (EPSG:3826)
// 並建立快速最近點查詢用的 RTree 索引
ini_set('memory_limit', '-1');
set_time_limit(0);
function my_pathinfo($path)
{
$dirname = '';
$basename = '';
$extension = '';
$filename = '';
$pos = strrpos($path, '/');
if ($pos !== false) {
$dirname = substr($path, 0, strrpos($path, '/'));
$basename = substr($path, strrpos($path, '/') + 1);
} else {
$basename = $path;
}
$ext = strrchr($path, '.');
if ($ext !== false) {
$extension = substr($ext, 1);
}
$filename = $basename;
$pos = strrpos($basename, '.');
if ($pos !== false) {
$filename = substr($basename, 0, $pos);
}
return array(
'dirname' => $dirname,
'basename' => $basename,
'extension' => $extension,
'filename' => $filename
);
}
function mainname($fname)
{
$pathinfo = my_pathinfo($fname);
return $pathinfo['filename'];
}
function rrmdirFilesByPattern($pattern)
{
foreach (glob($pattern) as $f) {
if (is_file($f)) {
@unlink($f);
}
}
}
function runCmd($cmd)
{
echo "\n[CMD] {$cmd}\n";
exec($cmd . " 2>&1", $output, $ret);
if (!empty($output)) {
echo implode("\n", $output) . "\n";
}
if ($ret !== 0) {
throw new Exception("指令執行失敗，exit code={$ret}\n{$cmd}");
}
}
function grdToCsv($srcFile, $csvFile, $sourceName)
{
$in = fopen($srcFile, 'r');
if (!$in) {
throw new Exception("無法開啟檔案: {$srcFile}");
}
$out = fopen($csvFile, 'w');
if (!$out) {
fclose($in);
throw new Exception("無法寫入 CSV: {$csvFile}");
}
// 加欄位，保留來源檔名
fputcsv($out, array('X', 'Y', 'Z', 'SOURCE_FILE'));
$rowCount = 0;
while (($line = fgets($in)) !== false) {
$line = trim(str_replace(array("\r", "\n"), '', $line));
if ($line === '') {
continue;
}
// 允許空白 / tab / 多個空白
$parts = preg_split('/\s+/', $line);
if (count($parts) < 3) {
continue;
}
$x = trim($parts[0]);
$y = trim($parts[1]);
$z = trim($parts[2]);
if (!is_numeric($x) || !is_numeric($y) || !is_numeric($z)) {
continue;
}
fputcsv($out, array($x, $y, $z, $sourceName));
$rowCount++;
}
fclose($in);
fclose($out);
return $rowCount;
}
function buildCombinedCsv($csvFiles, $outputCsv)
{
$out = fopen($outputCsv, 'w');
if (!$out) {
throw new Exception("無法建立 {$outputCsv}");
}
fputcsv($out, array('X', 'Y', 'Z', 'SOURCE_FILE'));
$totalRows = 0;
foreach ($csvFiles as $csv) {
$in = fopen($csv, 'r');
if (!$in) {
continue;
}
$isFirst = true;
while (($row = fgetcsv($in)) !== false) {
if ($isFirst) {
$isFirst = false; // 跳過 header
continue;
}
if (count($row) < 4) {
continue;
}
fputcsv($out, $row);
$totalRows++;
}
fclose($in);
}
fclose($out);
return $totalRows;
}
// =========================================================
// 主流程
// =========================================================
try {
$dbFile = 'taiwan_gdr.db';
$allCsv = 'all_points.csv';
$vrtFile = 'all_points.vrt';
$tableName = 'points';
// 清理舊檔
@unlink($dbFile);
@unlink($allCsv);
@unlink($vrtFile);
rrmdirFilesByPattern('*.csv');
// 同時支援 .gdr / .grd
$files = array_merge(glob('*.gdr'), glob('*.grd'));
$files = array_values(array_unique($files));
if (count($files) === 0) {
throw new Exception("目前目錄找不到 *.gdr 或 *.grd");
}
echo "找到 " . count($files) . " 個來源檔\n";
$csvFiles = array();
$grandTotal = 0;
$maxK = count($files);
foreach ($files as $k => $srcFile) {
$mn = mainname($srcFile);
$csvFile = $mn . '.csv';
if(is_file($csvFile)) {
continue;
}
//echo $csvFile;
//exit();
//echo "WTF";
$rows = grdToCsv($srcFile, $csvFile, basename($srcFile));
$csvFiles[] = $csvFile;
$grandTotal += $rows;
echo sprintf("[%d/%d] %s -> %s, rows=%d\n", $k + 1, $maxK, $srcFile, $csvFile, $rows);
}
echo "全部點數: {$grandTotal}\n";
// 合併成單一 CSV
$mergedRows = buildCombinedCsv($csvFiles, $allCsv);
echo "合併完成: {$allCsv}, rows={$mergedRows}\n";
// 建立 VRT，讓 ogr2ogr 可把 X/Y/Z 變成 3D Point
$csvLayerName = mainname($allCsv); // all_points
$vrt = <<<EOF
<OGRVRTDataSource>
<OGRVRTLayer name="{$tableName}">
<SrcDataSource>{$allCsv}</SrcDataSource>
<SrcLayer>{$csvLayerName}</SrcLayer>
<GeometryType>wkbPoint25D</GeometryType>
<LayerSRS>EPSG:3826</LayerSRS>
<GeometryField encoding="PointFromColumns" x="X" y="Y" z="Z"/>
</OGRVRTLayer>
</OGRVRTDataSource>
EOF;
file_put_contents($vrtFile, $vrt);
// 匯入 SQLite
// 保留 X/Y/Z/SOURCE_FILE 欄位 + Geometry
$cmd = 'ogr2ogr -f SQLite '
. escapeshellarg($dbFile) . ' '
. escapeshellarg($vrtFile)
. ' -nln ' . escapeshellarg($tableName)
. ' -a_srs EPSG:3826'
. ' -lco SPATIAL_INDEX=NO';
runCmd($cmd);
// 建索引：一般索引 + RTree
$pdo = new PDO('sqlite:' . $dbFile);
$pdo->setAttribute(PDO::ATTR_ERRMODE, PDO::ERRMODE_EXCEPTION);
// 一般索引
$pdo->exec("CREATE INDEX IF NOT EXISTS idx_points_x ON {$tableName}(X)");
$pdo->exec("CREATE INDEX IF NOT EXISTS idx_points_y ON {$tableName}(Y)");
$pdo->exec("CREATE INDEX IF NOT EXISTS idx_points_xy ON {$tableName}(X, Y)");
$pdo->exec("CREATE INDEX IF NOT EXISTS idx_points_source_file ON {$tableName}(SOURCE_FILE)");
// RTree for nearest search
$pdo->exec("DROP TABLE IF EXISTS points_rtree");
$pdo->exec("
CREATE VIRTUAL TABLE points_rtree USING rtree(
id,
minX, maxX,
minY, maxY
)
");
$pdo->exec("
INSERT INTO points_rtree(id, minX, maxX, minY, maxY)
SELECT
OGC_FID,
CAST(X AS REAL), CAST(X AS REAL),
CAST(Y AS REAL), CAST(Y AS REAL)
FROM {$tableName}
");
// 統計資訊
$pdo->exec("ANALYZE");
$count = $pdo->query("SELECT COUNT(*) FROM {$tableName}")->fetchColumn();
echo "DB 建立完成: {$dbFile}\n";
echo "總筆數: {$count}\n";
echo "座標系統: EPSG:3826\n";
echo "資料表: {$tableName}\n";
echo "已建立 RTree: points_rtree\n";
} catch (Exception $e) {
echo "\n[ERROR] " . $e->getMessage() . "\n";
exit(1);
}
$POSTS_STRING = "x,y"; // x6 y7
$POSTS = getGET_POST($POSTS_STRING,"POST");
$x = (float) $POSTS['x'];
$y = (float) $POSTS['y'];
$pdoGDR = new PDO('sqlite:dem/taoyuan/taiwan_gdr.db');
$pdoGDR->setAttribute(PDO::ATTR_ERRMODE, PDO::ERRMODE_EXCEPTION);
//$result = findNearestPoint($pdoGDR, 250000.0, 2650000.0);
$result = findNearestPoint($pdoGDR, $x, $y, 1, 5000, 16);
$OUTPUT=[];
$OUTPUT['status'] = 'OK';
$OUTPUT['data'] = $result;
// 計算高度
$OUTPUT['x'] = $x;
$OUTPUT['y'] = $y;
if(is_array($result)){
$OUTPUT['z'] = getGridInterpolatedZ($result,$x,$y); // 內插高度
}
echo json_encode($OUTPUT);
exit();