關閉      標題：dem grd hdr 轉 tiff mapserver 發地形圖- 圓滑版 -- 這版不好
內容：


<?php
// grd_to_tif.php


 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 subname($fname){
    //$pathinfo=pathinfo($fname);
    //$pathinfo['extension'];
    
    $m=explode(".",$fname);
    return end($m);
  } 
  function mainname($fname){
    $pathinfo=my_pathinfo($fname);
    //print_r($pathinfo);
    //exit();
    return $pathinfo['filename'];           
  }
  
  
$EPSG = "EPSG:3826";
$DX = 20; $DY = 20;
$NODATA = -32768;   // 你原本用的 nodata
$RADIUS = 30;       // 20m DEM 建議 20~30
$fp = glob("*.grd");
$max_k = count($fp);

foreach($fp as $k => $v){
  $bn = basename($v);
  $mn = mainname($v);

  // 1) grd -> csv（逐行轉，避免空白/多欄位）
  $csv = "{$mn}.csv";
  if(!is_file($csv)){
    $in = fopen($v, "rb");
    $out = fopen($csv, "wb");
    fwrite($out, "X,Y,Z\n");
    while(!feof($in)){
      $line = fgets($in);
      if($line===false) break;
      $line = trim(str_replace("\r","", $line));
      if($line==="") continue;
      $parts = preg_split('/\s+/', $line);
      if(count($parts) < 3) continue;
      fwrite($out, $parts[0].",".$parts[1].",".$parts[2]."\n");
    }
    fclose($in); fclose($out);
  }

  // 2) csv -> OGRVRT (強制幾何 + 指定 SrcLayer)
  // CSV 的 layer 名通常就是檔名（不含副檔名）
  $vrt = "{$mn}.vrt";
  if(!is_file($vrt)){
    $vrt_xml = <<<EOF
<OGRVRTDataSource>
  <OGRVRTLayer name="dem">
    <SrcDataSource relativeToVRT="1">{$mn}.csv</SrcDataSource>
    <SrcLayer>{$mn}</SrcLayer>
    <GeometryType>wkbPoint</GeometryType>
    <LayerSRS>{$EPSG}</LayerSRS>
    <GeometryField encoding="PointFromColumns" x="X" y="Y" z="Z"/>
  </OGRVRTLayer>
</OGRVRTDataSource>
EOF;
    file_put_contents($vrt, $vrt_xml);
  }

  // 3) 掃 min/max 範圍（用來給 gdal_grid 的 -txe/-tye/-outsize）
  //    這一步可以省略，但你要格網正確，最好做
  $xmin=$ymin=null; $xmax=$ymax=null;
  $fh = fopen($v, "rb");
  while(!feof($fh)){
    $line = fgets($fh);
    if($line===false) break;
    $line = trim(str_replace("\r","",$line));
    if($line==="") continue;
    $p = preg_split('/\s+/', $line);
    if(count($p) < 3) continue;
    $x = floatval($p[0]); $y = floatval($p[1]);
    if($xmin===null){ $xmin=$xmax=$x; $ymin=$ymax=$y; }
    else{
      if($x<$xmin) $xmin=$x;
      if($x>$xmax) $xmax=$x;
      if($y<$ymin) $ymin=$y;
      if($y>$ymax) $ymax=$y;
    }
  }
  fclose($fh);

  // 4) extent 對齊 20m 格點（避免像素變 19.x）
  $xmin = floor($xmin/$DX)*$DX;
  $ymin = floor($ymin/$DY)*$DY;
  $xmax = ceil($xmax/$DX)*$DX;
  $ymax = ceil($ymax/$DY)*$DY;

  $cols = (int)(($xmax-$xmin)/$DX) + 1;
  $rows = (int)(($ymax-$ymin)/$DY) + 1;

  // 5) gdal_grid 產 tif（這才是正確的 XYZ->raster）
  $tif = "{$mn}.tif";
  if(!is_file($tif)){
    // 注意：-tye 要 ymax ymin，才能確保 PixelSize Y 為負
    $CMD = "gdal_grid -a_srs {$EPSG} -of GTiff -zfield Z ".
           "-a nearest:radius1={$RADIUS}:radius2={$RADIUS}:nodata={$NODATA} ".
           "-txe {$xmin} {$xmax} -tye {$ymax} {$ymin} -outsize {$cols} {$rows} ".
           "\"{$vrt}\" \"{$tif}\"";
    echo $CMD, "\n";
    `$CMD`;
  }

  // 6) warp 統一到 20m 格網（tap），順便把 nodata 固定好
  $fix = "{$mn}_fix.tif";
  if(!is_file($fix)){
    $CMD = "gdalwarp -overwrite -r near -tr 20 20 -tap ".
           "-srcnodata {$NODATA} -dstnodata {$NODATA} ".
           "\"{$tif}\" \"{$fix}\"";
    echo $CMD, "\n";
    `$CMD`;
  }

  // 7) 平滑（cubic）一定要帶 src/dst nodata，否則邊界會抹出怪值
  $smooth = "{$mn}_smooth.tif";
  if(!is_file($smooth)){
    $CMD = "gdalwarp -overwrite -r cubic ".
           "-srcnodata {$NODATA} -dstnodata {$NODATA} ".
           "\"{$fix}\" \"{$smooth}\"";
    echo $CMD, "\n";
    `$CMD`;
  }

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



// Windows 才有 dir /b；若在 Linux 要改成 ls
`dir /b *_fix.tif > list.txt`;
`gdalbuildvrt -input_file_list list.txt all.vrt`;

// 先整張統一（避免縫）
`gdalwarp -overwrite -r near -tr 20 20 -tap -srcnodata -32768 -dstnodata -32768 all.vrt all_aligned.tif`;

// 想圓滑：整張做一次 cubic（比每片做更不容易放大接縫）
`gdalwarp -overwrite -r cubic -srcnodata -32768 -dstnodata -32768 all_aligned.tif all_smooth.tif`;

// 等高線（10m）
`gdal_contour -i 10 -a elev all_smooth.tif contour_10m.gpkg`;

// 匯出 SQLite，幾何欄位 GEOMETRY
`ogr2ogr -f SQLITE 8429.db contour_10m.gpkg -s_srs EPSG:3826 -t_srs EPSG:4326 -nln 8429 -lco GEOMETRY_NAME=GEOMETRY`;