合併內政部20M DEM台灣本島與澎湖圖幅，並用30M SRTM補上被挖除的軍事區域

DEM_merge.py

import os
import urllib.request
import zipfile
import math
import numpy as np
import rasterio
from rasterio.merge import merge
from rasterio.warp import transform_bounds, reproject, Resampling, calculate_default_transform
from rasterio.transform import Affine
from rasterio.io import MemoryFile

def merge_and_fill_bbox(moi_taiwan_path, moi_penghu_path, output_path):
    # --- 步驟 1：下載範圍內所需的所有 SRTM ---
    srtm_hgt_files = []
    print("[1/6] 正在檢查並下載範圍內的 SRTM 瓦片 (這包含海洋與陸地)...")
    for lat in range(21, 26):
        for lon in range(119, 123):
            tile_name = f"N{lat:02d}E{lon:03d}"
            hgt_filename = f"{tile_name}.hgt"
            zip_path = f"{tile_name}.zip"
            url = f"https://step.esa.int/auxdata/dem/SRTMGL1/{tile_name}.SRTMGL1.hgt.zip"
            
            if os.path.exists(hgt_filename):
                srtm_hgt_files.append(hgt_filename)
                continue
            
            try:
                req = urllib.request.Request(url, method='HEAD')
                urllib.request.urlopen(req, timeout=5)
                print(f"  -> 下載 {tile_name}...")
                urllib.request.urlretrieve(url, zip_path)
                with zipfile.ZipFile(zip_path, 'r') as zip_ref:
                    zip_ref.extractall(".")
                os.remove(zip_path)
                srtm_hgt_files.append(hgt_filename)
            except Exception:
                print(f"  -> {tile_name} 無檔案 (可能全為海洋)，已跳過。")

    # --- 步驟 2：合併 SRTM ---
    print("\n[2/6] 正在合併 SRTM 瓦片...")
    srtm_srcs = [rasterio.open(f) for f in srtm_hgt_files]
    srtm_merged, srtm_transform = merge(srtm_srcs)
    srtm_crs = srtm_srcs[0].crs
    for src in srtm_srcs:
        src.close()

    # --- 步驟 3：處理座標系 (CRS) 衝突並合併內政部 DTM ---
    print("\n[3/6] 正在合併內政部 DTM (自動處理 EPSG:3826 與 3825 衝突)...")
    taiwan_src = rasterio.open(moi_taiwan_path)
    penghu_src = rasterio.open(moi_penghu_path)

    moi_srcs_to_merge = [taiwan_src]
    memfile = None
    penghu_reprojected = None

    if taiwan_src.crs != penghu_src.crs:
        print(f"  -> 發現座標系不同！本島: {taiwan_src.crs}, 澎湖: {penghu_src.crs}")
        print("  -> 正在將澎湖動態投影至本島座標系...")
        
        # 計算澎湖轉移到本島座標系後的變形參數
        transform, width, height = calculate_default_transform(
            penghu_src.crs, taiwan_src.crs, penghu_src.width, penghu_src.height, *penghu_src.bounds
        )
        
        kwargs = penghu_src.meta.copy()
        kwargs.update({
            'crs': taiwan_src.crs,
            'transform': transform,
            'width': width,
            'height': height
        })
        
        # 在記憶體中建立一個虛擬的 TIFF 來裝載轉換後的澎湖
        memfile = MemoryFile()
        penghu_reprojected = memfile.open(**kwargs)
        
        reproject(
            source=rasterio.band(penghu_src, 1),
            destination=rasterio.band(penghu_reprojected, 1),
            src_transform=penghu_src.transform,
            src_crs=penghu_src.crs,
            dst_transform=transform,
            dst_crs=taiwan_src.crs,
            resampling=Resampling.nearest
        )
        moi_srcs_to_merge.append(penghu_reprojected)
    else:
        moi_srcs_to_merge.append(penghu_src)

    print("  -> 執行無縫拼接...")
    moi_merged, moi_transform = merge(moi_srcs_to_merge)
    moi_crs = taiwan_src.crs
    moi_nodata = taiwan_src.nodata
    moi_meta = taiwan_src.meta.copy()

    # 清理檔案釋放記憶體
    if penghu_reprojected: penghu_reprojected.close()
    if memfile: memfile.close()
    taiwan_src.close()
    penghu_src.close()

    # --- 步驟 4：計算目標大網格的範圍 ---
    print("\n[4/6] 正在建立目標範圍 (119~122.6E, 21~26N) 的標準大網格...")
    left_wgs, bottom_wgs, right_wgs, top_wgs = 119.0, 21.0, 122.6, 26.0
    left_moi, bottom_moi, right_moi, top_moi = transform_bounds(
        'EPSG:4326', moi_crs, left_wgs, bottom_wgs, right_wgs, top_wgs
    )

    col_min = math.floor((left_moi - moi_transform.c) / moi_transform.a)
    col_max = math.ceil((right_moi - moi_transform.c) / moi_transform.a)
    row_min = math.floor((top_moi - moi_transform.f) / moi_transform.e)
    row_max = math.ceil((bottom_moi - moi_transform.f) / moi_transform.e)

    out_width = col_max - col_min
    out_height = row_max - row_min
    out_transform = moi_transform * Affine.translation(col_min, row_min)

    print(f"  -> 生成矩陣大小: {out_width} x {out_height} 像素")

    # --- 步驟 5：將資料投影至標準畫布 ---
    print("\n[5/6] 正在將資料繪製到大網格上 (預計佔用較多記憶體)...")
    canvas_moi = np.full((1, out_height, out_width), fill_value=(moi_nodata if moi_nodata is not None else np.nan), dtype=np.float32)
    canvas_srtm = np.full((1, out_height, out_width), fill_value=-32768, dtype=np.float32)

    reproject(
        source=moi_merged, destination=canvas_moi,
        src_transform=moi_transform, src_crs=moi_crs,
        dst_transform=out_transform, dst_crs=moi_crs,
        resampling=Resampling.nearest
    )

    reproject(
        source=srtm_merged, destination=canvas_srtm,
        src_transform=srtm_transform, src_crs=srtm_crs,
        dst_transform=out_transform, dst_crs=moi_crs,
        resampling=Resampling.bilinear
    )

    # --- 步驟 6：填補缺值與輸出 ---
    print("\n[6/6] 正在將 SRTM 填入缺口並輸出...")
    if moi_nodata is None or np.isnan(moi_nodata):
        missing_mask = np.isnan(canvas_moi[0])
    else:
        missing_mask = (canvas_moi[0] == moi_nodata)

    valid_srtm_mask = (canvas_srtm[0] > -100) & (canvas_srtm[0] < 4500)
    fill_mask = missing_mask & valid_srtm_mask
    
    final_data = np.copy(canvas_moi[0])
    final_data[fill_mask] = canvas_srtm[0][fill_mask]

    # 將海洋/無資料區域統一補為 0
    final_data[np.isnan(final_data)] = 0
    if moi_nodata is not None and not np.isnan(moi_nodata):
        final_data[final_data == moi_nodata] = 0

    moi_meta.update({
        "height": out_height,
        "width": out_width,
        "transform": out_transform,
        "compress": "lzw",
        "dtype": "float32",
        "nodata": 0
    })

    with rasterio.open(output_path, "w", **moi_meta) as dest:
        dest.write(final_data, 1)

    print("\n🎉 大功告成！檔案儲存完畢。")

if __name__ == "__main__":
    INPUT_MOI_TAIWAN = "DEM_tawiwan_V2025.tif"
    INPUT_MOI_PENGHU = "DEM_Penghu_V2025.tif"
    OUTPUT_FINAL_DTM = "DEM_merged_V2025.tif"
    
    merge_and_fill_bbox(INPUT_MOI_TAIWAN, INPUT_MOI_PENGHU, OUTPUT_FINAL_DTM)