關閉      標題：全臺地形坡度圖20m
內容：

資料來源：https://www.tgos.tw/MDE/MetaData/CRUD/ViewMetaData?GUID=74c6a7b7-1816-499c-bb60-89246065390f&VIEW_TYPE=only&SOURCE=1
不分幅_全台及澎湖DEM


<?php

function cleanSlopePng($inFile, $outFile)
{
    $img = imagecreatefrompng($inFile);
    imagesavealpha($img, true);
    imagealphablending($img, false);

    $w = imagesx($img);
    $h = imagesy($img);

    // 先準備顏色
    $mark1 = imagecolorallocate($img, 251, 7, 6);   // 紅
    $mark2 = imagecolorallocate($img, 26, 150, 78);  // 綠
    $black = imagecolorallocate($img, 0, 0, 0); // 黑
    $white = imagecolorallocate($img, 255, 255, 255); // 白
    $trans = imagecolorallocatealpha($img, 0, 0, 0, 127);

    // 四角都處理
    $m = [[0,0],[$w-1,0],[$w-1,$h-1],[0,$h-1]];
    foreach($m as $v){

      // 1️⃣ 從左上角 flood fill 紅色
      imagefill($img, $v[0], $v[1], $mark1);

      // 2️⃣ 再從左上 flood fill 綠色（避免孤島）
      imagefill($img, $v[0], $v[1], $mark2);
      
      // 3️⃣ 從左上角 flood fill 紅色
      imagefill($img, $v[0], $v[1], $mark1);

      // 4️⃣ 再從左上 flood fill 綠色（避免孤島）
      imagefill($img, $v[0], $v[1], $mark2);
      
      // 5️⃣ 再從左上 flood fill 白色（避免孤島）
      imagefill($img, $v[0], $v[1], $white);

      // 6️⃣ 再 flood 成黑色，作為最後的「要刪除標記」
      imagefill($img, $v[0], $v[1], $black);
    }

    // 7️⃣ 掃描整張圖，把黑色改透明
    for ($y = 0; $y < $h; $y++) {
        for ($x = 0; $x < $w; $x++) {
            $rgb = imagecolorat($img, $x, $y);

            $r = ($rgb >> 16) & 0xFF;
            $g = ($rgb >> 8) & 0xFF;
            $b = $rgb & 0xFF;

            if ($r === 0 && $g === 0 && $b === 0) {
                imagesetpixel($img, $x, $y, $trans);
            }
        }
    }

    imagepng($img, $outFile);
    imagedestroy($img);
}


$SP = DIRECTORY_SEPARATOR;
$PD = dirname(__FILE__);
$tif_file            = "{$PD}{$SP}不分幅_全台及澎湖DEM{$SP}dem_20m.tif";
$tif_percent_file    = "{$PD}{$SP}不分幅_全台及澎湖DEM{$SP}dem_20m_percent.tif";
$slope_colors_file   = "{$PD}{$SP}slope_colors.txt";
$tif_finial_result = "{$PD}{$SP}不分幅_全台及澎湖DEM{$SP}dem_20m_finial.tif";
// 色表
//-9999 255 255 255
$slope_colors_data = "
0    26  150  78
5   140  202  93
15  215  235 140
30  251  252 194
40  248  222 138
55  211   49  34
100 251    7   6  
";
file_put_contents($slope_colors_file, trim($slope_colors_data));

// 1️⃣ 生成坡度百分比
shell_exec("gdaldem slope \"{$tif_file}\" \"{$tif_percent_file}\" -p -compute_edges");

// 2️⃣ 彩色 PNG，-alpha 會自動把 NoData 透明
shell_exec("gdaldem color-relief \"{$tif_percent_file}\" \"{$slope_colors_file}\" \"{$PD}{$SP}slope_20m_color.png\" -of PNG -alpha");

// 5. 刪色
cleanSlopePng("{$PD}{$SP}slope_20m_color.png", "{$PD}{$SP}slope_20m_color_clean.png");

// 6. 取得原本的 坐標
$meta = json_decode(shell_exec("gdalinfo -json " . escapeshellarg($tif_file)), true);

$cc = $meta['cornerCoordinates'];

$corners3826 = [
    'upper_left'  => $cc['upperLeft'],
    'upper_right' => $cc['upperRight'],
    'lower_right' => $cc['lowerRight'],
    'lower_left'  => $cc['lowerLeft'],
];

// 7. 取得 PNG 尺寸
[$w, $h] = getimagesize("{$PD}{$SP}slope_20m_color_clean.png");


// 8. 將 png 轉回 tif
$ulx = $corners3826['upper_left'][0];
$uly = $corners3826['upper_left'][1];
$lrx = $corners3826['lower_right'][0];
$lry = $corners3826['lower_right'][1];
$cmd = sprintf(
    'gdal_translate -of GTiff -a_srs EPSG:3826 -a_ullr %f %f %f %f -co TILED=YES -co COMPRESS=DEFLATE "%s" "%s"',
    $ulx, $uly, $lrx, $lry,
    "{$PD}{$SP}slope_20m_color_clean.png",
    $tif_finial_result
);

shell_exec($cmd);

echo "Done...{$tif_finial_result}\n";